A discussion on “A selective review of statistical methods using calibration information from similar studies” by Qin, Liu and Li

We congratulate Qin, Liu and Li (QLL) on a thoughtful and much needed review of many interesting methods for combining information from similar studies. We appreciate being given the opportunity to make a discussion. QLL cover a variety of different settings and methods. Based on that, we will provide a brief review on some additional relevant literature with a focus on methods that deal with population heterogeneity, since in practice it is most likely that different studies sample from different populations and whether information should be combined depends on how similar those populations are, among many other considerations. To keep the discussion focussed, we will follow the setting in Section 5 of QLL, although most of these methods can be more broadly applied.

We adopt the notation in Section 5.1 of QLL with some variations. Let $(Y_i,{\mathit{X}}_i^{\mathrm{T}},{\mathit{Z}}_i^{\mathrm{T}}{)}^{\mathrm{T}}$, $i=1,\dots ,n$, denote the individual data based on a random sample from the internal study, where Y is the response and $\mathit{X}$ and $\mathit{Z}$ are vectors of covariates. The model of interest is $f(Y|\mathit{X},\mathit{Z};\text{β})$ for $f(Y|\mathit{X},\mathit{Z})$ with parameter $\text{β}$. The external study fitted a (possibly misspecified) model $f(Y|\mathit{X};\mathit{\theta })$ for $f^{\ast }\left(Y|\mathit{X}\right)$ with covariates $\mathit{X}$ alone and parameter $\mathit{\theta }$. Throughout this discussion we use a *-superscript to denote distributions/expectations/quantities associated with the external study population. The external model fitting information is summarized by (1) $E^{\ast }\left\{\mathit{h}\right(Y,\mathit{X};{\mathit{\theta }}^{\ast }\left)\right\}=\mathbf{0},$(1) where $\mathit{h}(Y,\mathit{X};\mathit{\theta })$ is the score function for $f(Y|\mathit{X};\mathit{\theta })$ and ${\mathit{\theta }}^{\ast }$ is the solution to the score equation based on the external study sample. Individual data from the external sample are not available. We assume the external sample size is very large so that the uncertainty in ${\mathit{\theta }}^{\ast }$ is negligible compared to the internal study, i.e. Case I in QLL.

QLL in Section 5 give an excellent review of some methods and their comparisons in terms of asymptotic efficiency when the internal and external study populations are the same, based on both the empirical likelihood (EL) formulation (Owen, 1988; Qin & Lawless, 1994) and the constrained maximum likelihood (CML) formulation (Chatterjee et al., 2016; Qin, 2000). These two formulations are closely connected (Han & Lawless, 2016). For ease of discussion, we provide the CML formulation here, which is already covered in QLL. When $f(Y,\mathit{X},\mathit{Z})=f^{\ast }(Y,\mathit{X},\mathit{Z})$, (1(1) $E^{\ast }\left\{\mathit{h}\right(Y,\mathit{X};{\mathit{\theta }}^{\ast }\left)\right\}=\mathbf{0},$(1) ) can be transformed into (2) $E\left\{\mathit{\psi }\right(\mathit{X},\mathit{Z};\text{β}\left)\right\}=\mathbf{0},$(2) where (3) $\begin{array}{rl}\mathit{\psi }(\mathit{X},\mathit{Z};\text{β})& =E\left\{\mathit{h}\right(Y,\mathit{X};{\mathit{\theta }}^{\ast })|\mathit{X},\mathit{Z}\}\\ & =\int \mathit{h}(Y,\mathit{X};{\mathit{\theta }}^{\ast })f(Y|\mathit{X},\mathit{Z};\text{β})\mathrm{d}Y.\end{array}$(3) The CML estimator ${\widehat{\text{β}}}_{cml}$ is defined through (4) $\begin{array}{rl}\underset{\text{β},p_1,\dots ,p_n}{max}& \prod _{i=1}^{n}f(Y_i\mid {\mathit{X}}_i,{\mathit{Z}}_i;\text{β})p_i\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}& p_i\ge 0,\sum _{i=1}^n{p}_i=1,\\ & \sum _{i=1}^n{p}_i\mathit{\psi }({\mathit{X}}_i,{\mathit{Z}}_i;\text{β})=\mathbf{0},\end{array}$(4) where $p_i\equiv \mathrm{d}F({\mathit{X}}_i,{\mathit{Z}}_i)$, $i=1,\dots ,n$.

When $f(Y|\mathit{X},\mathit{Z})=f^{\ast }(Y|\mathit{X},\mathit{Z})$ but $f(\mathit{X},\mathit{Z})\ne f^{\ast }(\mathit{X},\mathit{Z})$, applying the transformation (3(3) $\begin{array}{rl}\mathit{\psi }(\mathit{X},\mathit{Z};\text{β})& =E\left\{\mathit{h}\right(Y,\mathit{X};{\mathit{\theta }}^{\ast })|\mathit{X},\mathit{Z}\}\\ & =\int \mathit{h}(Y,\mathit{X};{\mathit{\theta }}^{\ast })f(Y|\mathit{X},\mathit{Z};\text{β})\mathrm{d}Y.\end{array}$(3) ) to (1(1) $E^{\ast }\left\{\mathit{h}\right(Y,\mathit{X};{\mathit{\theta }}^{\ast }\left)\right\}=\mathbf{0},$(1) ) will lead to (5) $E^{\ast }\left\{\mathit{\psi }\right(\mathit{X},\mathit{Z};\text{β}\left)\right\}=\mathbf{0},$(5) where the expectation $E^{\ast }(\cdot )$ is under $f^{\ast }(\mathit{X},\mathit{Z})$ in contrast to the $E(\cdot )$ in (2(2) $E\left\{\mathit{\psi }\right(\mathit{X},\mathit{Z};\text{β}\left)\right\}=\mathbf{0},$(2) ) that is under $f(\mathit{X},\mathit{Z})$. In this case, although (2(2) $E\left\{\mathit{\psi }\right(\mathit{X},\mathit{Z};\text{β}\left)\right\}=\mathbf{0},$(2) ) no longer holds, there are ways to make use of the external study information summarized by (5(5) $E^{\ast }\left\{\mathit{\psi }\right(\mathit{X},\mathit{Z};\text{β}\left)\right\}=\mathbf{0},$(5) ). One way is to collect a small supplementary sample $({\mathit{X}}_j^{\ast \mathrm{T}},{\mathit{Z}}_j^{\ast \mathrm{T}}{)}^{\mathrm{T}}$, $j=1,\dots ,n^{\ast }$, from the external study population (Chatterjee et al., 2016; Han & Lawless, 2019), and define the CML estimator through $\begin{array}{rl}\underset{\text{β},p_1^{\ast },\dots ,p_{n^{\ast }}^{\ast }}{max}& \prod _{i=1}^{n}f(Y_i\mid {\mathit{X}}_i,{\mathit{Z}}_i;\text{β})\prod _{j=1}^{n^{\ast }}{p}_j^{\ast }\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}& p_j^{\ast }\ge 0,\sum _{j=1}^{n^{\ast }}{p}_j^{\ast }=1,\\ & \sum _{j=1}^{n^{\ast }}{p}_j^{\ast }\mathit{\psi }({\mathit{X}}_j^{\ast },{\mathit{Z}}_j^{\ast };\text{β})=\mathbf{0},\end{array}$where $p_j^{\ast }\equiv \mathrm{d}{F}^{\ast }({\mathit{X}}_j^{\ast },{\mathit{Z}}_j^{\ast })$, $j=1,\dots ,n^{\ast }$. Another way is to specify a density ratio model of the form $f^{\ast }(\mathit{X},\mathit{Z})=\exp \left({\mathit{W}}^{\mathrm{T}}\mathit{\alpha }\right)f(\mathit{X},\mathit{Z})$, where $\mathit{W}$ is a subset of $({\mathit{X}}^{\mathrm{T}},{\mathit{Z}}^{\mathrm{T}}{)}^{\mathrm{T}}$ and $\mathit{\alpha }$ is newly introduced parameter, to link the external distribution $f^{\ast }(\mathit{X},\mathit{Z})$ to the internal distribution $f(\mathit{X},\mathit{Z})$ (Sheng et al., 2021). This approach models the heterogeneity between $f^{\ast }(\mathit{X},\mathit{Z})$ and $f(\mathit{X},\mathit{Z})$ by $\exp \left({\mathit{W}}^{\mathrm{T}}\mathit{\alpha }\right)$, and is in the same spirit that an exponential tilting model is specified to link the case distribution to the control distribution as in Section 7.1 of QLL. With the density ratio model, the CML estimator can be defined through $\begin{array}{rl}\underset{\text{β},\mathit{\alpha },p_1,\dots ,p_n}{max}& \prod _{i=1}^{n}f(Y_i\mid {\mathit{X}}_i,{\mathit{Z}}_i;\text{β})p_i\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}& p_i\ge 0,\sum _{i=1}^n{p}_i=1,\\ & \sum _{i=1}^n{p}_i\{\exp ({\mathit{W}}_i^{\mathrm{T}}\mathit{\alpha })-1\}=\mathbf{0},\\ & \sum _{i=1}^n{p}_i\mathit{\psi }({\mathit{X}}_i,{\mathit{Z}}_i;\text{β})\exp \left({\mathit{W}}_i^{\mathrm{T}}\mathit{\alpha }\right)=\mathbf{0},\end{array}$where $p_i\equiv \mathrm{d}F({\mathit{X}}_i,{\mathit{Z}}_i)$, $i=1,\dots ,n$. Here the last two constraints are based on the fact that $f^{\ast }(\mathit{X},\mathit{Z})=\exp \left({\mathit{W}}^{\mathrm{T}}\mathit{\alpha }\right)f(\mathit{X},\mathit{Z})$ is a density and (5(5) $E^{\ast }\left\{\mathit{\psi }\right(\mathit{X},\mathit{Z};\text{β}\left)\right\}=\mathbf{0},$(5) ), respectively. Note that the dimension of $\mathit{\alpha }$ needs to be no larger than the dimension of $\mathit{\psi }$ plus one for $\mathit{\alpha }$ to be identifiable.

When $f(Y|\mathit{X},\mathit{Z})\ne f^{\ast }(Y|\mathit{X},\mathit{Z})$, applying transformation (3(3) $\begin{array}{rl}\mathit{\psi }(\mathit{X},\mathit{Z};\text{β})& =E\left\{\mathit{h}\right(Y,\mathit{X};{\mathit{\theta }}^{\ast })|\mathit{X},\mathit{Z}\}\\ & =\int \mathit{h}(Y,\mathit{X};{\mathit{\theta }}^{\ast })f(Y|\mathit{X},\mathit{Z};\text{β})\mathrm{d}Y.\end{array}$(3) ) to (1(1) $E^{\ast }\left\{\mathit{h}\right(Y,\mathit{X};{\mathit{\theta }}^{\ast }\left)\right\}=\mathbf{0},$(1) ) leads to neither (2(2) $E\left\{\mathit{\psi }\right(\mathit{X},\mathit{Z};\text{β}\left)\right\}=\mathbf{0},$(2) ) nor (5(5) $E^{\ast }\left\{\mathit{\psi }\right(\mathit{X},\mathit{Z};\text{β}\left)\right\}=\mathbf{0},$(5) ) in general, regardless of if $f(\mathit{X},\mathit{Z})$ is the same as $f^{\ast }(\mathit{X},\mathit{Z})$. In this case, the aforementioned CML estimators are biased. However, since the external study sample size is large, the reduction in variance by making use of external summary information may still benefit the internal study parameter estimation from a mean squared error perspective. Based on consideration of such a bias-variance trade-off, Estes et al. (2018) proposed an empirical Bayes shrinkage estimator of the form ${\hat{\mathit{V}}}_0({\hat{\mathit{V}}}_{mle}+{\hat{\mathit{V}}}_0{)}^{-1}{\widehat{\text{β}}}_{mle}+{\hat{\mathit{V}}}_{mle}({\hat{\mathit{V}}}_{mle}+{\hat{\mathit{V}}}_0{)}^{-1}{\widehat{\text{β}}}_{cml}$that is a weighted average of ${\widehat{\text{β}}}_{mle}$, the MLE using the internal study data alone, and ${\widehat{\text{β}}}_{cml}$ defined through (4(4) $\begin{array}{rl}\underset{\text{β},p_1,\dots ,p_n}{max}& \prod _{i=1}^{n}f(Y_i\mid {\mathit{X}}_i,{\mathit{Z}}_i;\text{β})p_i\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}& p_i\ge 0,\sum _{i=1}^n{p}_i=1,\\ & \sum _{i=1}^n{p}_i\mathit{\psi }({\mathit{X}}_i,{\mathit{Z}}_i;\text{β})=\mathbf{0},\end{array}$(4) ). Here ${\hat{\mathit{V}}}_{mle}$ and ${\hat{\mathit{V}}}_0$ are estimated variance matrices for ${\widehat{\text{β}}}_{mle}$ and for the prior normal distribution on $\text{β}$, respectively. This method shrinks the final estimate towards ${\widehat{\text{β}}}_{mle}$ in the presence of population heterogeneity and towards ${\widehat{\text{β}}}_{cml}$ otherwise. Gu et al. (2021) extended the idea in Estes et al. (2018) to the case of multiple external studies, and propose an estimator that is a weighted average of the empirical Bayes estimators resulted from using each external study separately.

To deal with arbitrary population heterogeneity when information is available from multiple external studies but without causing estimation bias, Zhai and Han (2022) developed an estimation procedure that simultaneously selects the studies that give (2(2) $E\left\{\mathit{\psi }\right(\mathit{X},\mathit{Z};\text{β}\left)\right\}=\mathbf{0},$(2) ) and incorporates the corresponding information into internal model fitting. Their method also applies under the current setting of only one external study. When (2(2) $E\left\{\mathit{\psi }\right(\mathit{X},\mathit{Z};\text{β}\left)\right\}=\mathbf{0},$(2) ) does not hold because of an unknown form of heterogeneity, some components of $E\left\{\mathit{\psi }\right(\mathit{X},\mathit{Z};\text{β}\left)\right\}$ may still be zero even if $E\left\{\mathit{\psi }\right(\mathit{X},\mathit{Z};\text{β}\left)\right\}$ is not a zero vector, and these components still contain useful information to improve internal model estimation efficiency. This observation makes sense because the association between the same response and certain covariates may not differ much across populations with certain specific heterogeneity. Some general examples on this observation are given in Zhai and Han (2022). Let $\mathit{\gamma }\equiv E\left\{\mathit{\psi }\right(\mathit{X},\mathit{Z};\text{β}\left)\right\}$, then the components of $\mathit{\psi }(\mathit{X},\mathit{Z};\text{β})$ that correspond to zero-components of $\mathit{\gamma }$ should be selected to compute ${\widehat{\text{β}}}_{cml}$. To shrink the estimate of the zero-components of $\mathit{\gamma }$ to exactly zero for information integration, under the current setting, the estimator in Zhai and Han (2022) is defined through $\begin{array}{rl}\underset{\text{β},\mathit{\gamma }}{max}& [\underset{p_1,\dots ,p_n}{max}\log \left\{\prod _{i=1}^{n}f\right(Y_i|{\mathit{X}}_i,{\mathit{Z}}_i;\text{β}\left)p_i\right\}\\ & -n\sum _{k=1}^{\mathrm{d}\mathrm{i}\mathrm{m}\left(\mathit{\gamma }\right)}{\lambda }_n\frac{|{\gamma }_k|}{|{\tilde{\gamma }}_k{|}^w}]\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}& p_i\ge 0,\sum _{i=1}^n{p}_i=1,\\ & \sum _{i=1}^n{p}_i\left\{\mathit{\psi }\right({\mathit{X}}_i,{\mathit{Z}}_i;\text{β})-\mathit{\gamma }\}=\mathbf{0},\end{array}$with an adaptive Lasso penalty (Zou, 2006) on $\mathit{\gamma }$ with tuning parameter ${\lambda }_n>0$, where ${\tilde{\gamma }}_k$ is the kth component of $\tilde{\mathit{\gamma }}=n^{-1}\sum _{i=1}^n\mathit{\psi }({\mathit{X}}_i,{\mathit{Z}}_i;{\widehat{\text{β}}}_{mle})$ and w>0 is some user-specified positive number typically taken to be 1 or 2. Similar idea has also been considered for integrating external summary information into survival analysis with different penalty (Chen et al., 2021).

All of the aforementioned methods make use of the external summary information through transformation (3(3) $\begin{array}{rl}\mathit{\psi }(\mathit{X},\mathit{Z};\text{β})& =E\left\{\mathit{h}\right(Y,\mathit{X};{\mathit{\theta }}^{\ast })|\mathit{X},\mathit{Z}\}\\ & =\int \mathit{h}(Y,\mathit{X};{\mathit{\theta }}^{\ast })f(Y|\mathit{X},\mathit{Z};\text{β})\mathrm{d}Y.\end{array}$(3) ). Taylor et al. (2022) took a different approach when both $f(Y|\mathit{X},\mathit{Z};\text{β})$ and $f(Y|\mathit{X};\mathit{\theta })$ are generalized linear models (GLM), namely (6) $g\left\{\mathbb{E}\right(Y|\mathit{X},\mathit{Z}\left)\right\}={\text{β}}_0+{\mathit{X}}^{\mathrm{T}}{\text{β}}_X+{\mathit{Z}}^{\mathrm{T}}{\text{β}}_Z$(6) and $l\left\{\mathbb{E}\right(Y|\mathit{X}\left)\right\}={\theta }_0+{\mathit{X}}^{\mathrm{T}}{\mathit{\theta }}_X$, with possibly different link functions $g(\cdot )$ and $l(\cdot )$ (note that the second GLM may be misspecified). Here the notation $\mathbb{E}(\cdot )$ for expectation is generic to simply present the form of the model. For ease of discussion, assume covariates $\mathit{Z}$ have been orthogonalized to $\mathit{X}$, which can be done by taking $\mathit{Z}$ to be the vector of residuals of the least squares regression of each covariate in $\mathit{Z}$ on $\mathit{X}$ using the internal data. Taylor et al. (2022) showed that ${\text{β}}_X\approx c{\mathit{\theta }}_X$ for some unknown constant c when both GLMs are fitted to the same infinitely large sample and when ${\text{β}}_X$, ${\text{β}}_Z$ and ${\mathit{\theta }}_X$ are all close to zero (here in this sentence, with some abuse of notation, ${\text{β}}_X$, ${\text{β}}_Z$ and ${\mathit{\theta }}_X$ are the values after fitting both GLMs to the same infinitely large sample). See also Neuhaus and Jewell (1993). Based on this result, when $f\left(Y|\mathit{X}\right)$ and $f^{\ast }\left(Y|\mathit{X}\right)$ are similar in the sense that the relative effects of $\mathit{X}$ on Y (but not necessarily the absolute magnitudes) are the same between the two populations, the ${\mathit{\theta }}_X^{\ast }$ produced by the external study can still be used to improve the internal estimation efficiency. With $\mathit{Z}$ orthogonalized to $\mathit{X}$, instead of fitting (6(6) $g\left\{\mathbb{E}\right(Y|\mathit{X},\mathit{Z}\left)\right\}={\text{β}}_0+{\mathit{X}}^{\mathrm{T}}{\text{β}}_X+{\mathit{Z}}^{\mathrm{T}}{\text{β}}_Z$(6) ), Taylor et al. (2022) proposed to fit $g\left\{\mathbb{E}\right(Y|\mathit{X},\mathit{Z}\left)\right\}={\text{β}}_0+\alpha \left({\mathit{X}}^{\mathrm{T}}{\mathit{\theta }}_X^{\ast }\right)+{\mathit{Z}}^{\mathrm{T}}{\text{β}}_Z$with coefficients $({\text{β}}_0,\alpha ,{\text{β}}_Z^{\mathrm{T}}{)}^{\mathrm{T}}$ to the internal study data, which is equivalent to letting ${\text{β}}_X=\alpha {\mathit{\theta }}_X^{\ast }$.

In the presence of different study populations, a crucial question to ask before combining information is which population is of the primary interest. Most of the methods reviewed in this discussion, if not all, explicitly or implicitly assume that the internal study population is of the primary interest and the external summary information is used for efficiency improvement without causing (too much) bias. This is reasonable, for example, when the internal study has a clear target population and is based on a careful design with a well controlled sampling. In practice, with data usually collected based on convenience sampling, the internal study sample may not be representative of the target population, or there may even be ambiguity about the target population itself. Therefore, some cautions are always needed when applying those methods to combine information in the presence of population heterogeneity, and more methodological developments are definitely needed.

Disclosure statement

No potential conflict of interest was reported by the author.