Analyzing Treatment Effect by Integrating Existing Propensity Score and Outcome Regressions with Heterogeneous Covariate Sets

Abstract

Analyzing treatment or exposure effect is a major research theme in scientific studies. In the current big-data era where multiple sources of data are available, it is of interest to perform a synthesized analysis of treatment effects by integrating information from different data sources or studies. However, studies may contain heterogeneous and incomplete covariate sets, and individual data therein may not be accessible. We apply and extend the generalized meta-analysis method to integrate summary results (e.g., regression coefficients) of outcome and treatment (propensity score, PS) regression analyses across different datasets that may contain heterogeneous covariate sets. The proposed integrated analysis utilizes a reference dataset, which contains data on the complete set of covariates. The asymptotic distribution for the proposed integrated estimator is established. Simulations reveal that the proposed estimator performs well. We apply the proposed method to obtain the causal effect of waist circumference on hypertension by integrating two existing outcomes and PS regression analyses with different sets of covariates.

1. Introduction

The impact of a treatment or exposure on an outcome is a significant focus in many scientific fields. These include clinical trials [1], education [2], economics [3], and medical [4] fields. In observational studies, treatment assignment often correlates with subject characteristics, potentially leading to systematic differences between treated and control groups. This can result in biased conclusions when directly comparing these groups.

To perform causal inferences in observational studies, Rubin [5] proposed the counterfactual or potential outcome framework. Suppose the treatment assignment Z is binary with $Z=1$ denoting the treated group and $Z=0$ the control group, and $Y_1$ and $Y_0$ are the potential outcome values a subject would have if he/she had a treatment assignment $Z=1$ and $Z=0$, respectively. The outcome Y observed for a subject with treatment Z is then given by

$$Y=Z{Y}_1+(1-Z)Y_0.$$

A causal treatment effect for a subject is obtained as the difference of two potential outcomes, $Y_1-Y_0$. Usually, only one potential outcome canc be observed for a subject, and the causal effect is unobserved. However, under suitable assumptions, the average treatment effect (ATE)

$$E\left(Y_1-Y_0\right)=E\left(Y_1\right)-E\left(Y_0\right)$$

is estimable from observational studies and has been a popular target of causal inferences. For a binary outcome, the ATE amounts to the marginal risk difference $P\left(Y_1=1\right)-P\left(Y_0=1\right)$. A key assumption for estimating ATE with data from observational studies is that the treatment assignment is strongly ignorable [6]:

$$\left(Y_1,Y_0\right)\perp Z\mid \mathbf{X},$$

where “⊥” stands for statistical independence; that is, the potential outcomes $\left(Y_1,Y_0\right)$, and the treatment assignment Z are independent conditioning on the covariate set $\mathbf{X}$. In other words, the treatment assignment Z is irrelevant to the values of potential outcomes $Y_1$ and $Y_0$ once the covariates $\mathbf{X}$ are controlled. Another key assumption is the stable unit-treatment value assumption (SUTVA), which ensures that $Y_1\left(Y_0\right)$ is the unique outcome associated with $Z=1(Z=0)$ [6]. Under (1) and SUTVA, consistent estimation of the ATE can then be achieved by the technique of matching, stratification, covariate adjustment (CA), or inverse probability of treatment weighting (IPTW) [6,7,8,9,10].

Rosenbaum and Rubin [6] further proposed using propensity score (PS) for conducting causal inferences. A PS is the probability $P(Z=1\mid \mathbf{X})$ of a subject being assigned to the treated group conditioning on his/her observed covariates and is the coarsest balancing score such that the covariate distributions are the same between treatment groups once the score is matched. Accordingly, a PS can replace the full covariate set $\mathbf{X}$ used for matching, stratification, CA, or IPTW to conduct causal inferences on treatment effect when the assumption (1) holds Rosenbaum and Rubin [6]. Since the PS is a univariate minimal sufficient statistic, the PS-based method can also be viewed as an effective dimension-reduction technique, since matching, stratification, CA, or IPTW can be simply performed on a scalar PS rather than on a multi-dimensional covariate vector [6,7,8,9,10]. Since its invention, the PS method has gained tremendous popularity in observational studies; see Figure 1 of Simoneau et al. [11].

In the big data era, there is growing interest in conducting integrated analyses of treatment effects by integrating information from databases across various observational studies [12,13]. However, combining data from different databases presents challenges due to variations in covariate variables, even if they share common outcome and treatment variables, and discrepancies in baseline data distributions. Moreover, stringent data privacy regulations such as the European Union’s General Data Protection Regulation (GDPR) and the Health Insurance Portability and Accountability Act (HIPAA) pose barriers to directly sharing, exchanging, and synthesizing individual-level data from disparate sources.

Analyzing treatment or exposure effects is a major research theme in scientific studies. In the current big-data era where multiple data sources are available, it is interesting to perform a synthesized analysis of treatment effects by integrating information from different data sources or studies. However, studies may contain heterogeneous and incomplete covariate sets, and individual data may not be accessible. This motivates us to consider a framework for integrated regression analysis that can utilize summary results (e.g., regression coefficients) of outcome and treatment (PS) regression analyses across different datasets that may contain heterogeneous covariate sets.

This study is specifically motivated by two recent investigations into the relationship between waist circumference (WC) and hypertension (HT): Ren et al. [14] and Hu et al. [15]. Each study explored how WC influences the risk of HT while accounting for different sets of covariates (detailed in Section 4). The aim here is to examine the impact of WC on HT while controlling for the comprehensive set of covariates (combining the respective sets from both studies). Specifically, this study integrates the regression analysis results from two existing studies on the relationship between WC and HT, considering their respective sets of covariates. Additionally, it incorporates the findings from propensity score analyses of WC as a treatment variable in both studies. This integrated approach aims to conduct a synthesized analysis of the causal effect of WC on HT, utilizing a comprehensive covariate dataset sourced from the Taiwan BioBank database. Specifically, we apply and extend the generalized meta-analysis method developed by Kundu et al. [16] to integrate the existing PS and outcome regression analyses with heterogeneous covariate variables using a dataset on the complete covariate set, which provides a reference to the covariate distribution and is termed “reference dataset” by Kundu et al. [16]. We obtain consistent estimates for the parameters of the full PS and outcome regression models and for the ATE. Asymptotic theory for the estimators of the regression parameters and the ATE is also established. The proposed integrated analysis allows researchers to integrate summary results from existing outcome and PS regression analyses with heterogeneous covariates sets, and hence to obtain enhanced statistical efficiency and power from existing studies.

In Section 2, we describe the proposed estimation methods for the full PS and outcome regression models and the ATE. In Section 3, we report simulation results to show the satisfactory performance of the proposed method. Section 4 presents the empirical study on the effect of waist circumference on hypertension based on two existing studies and a reference sample of covariate data. Finally, concluding discussions are provided in Section 5.

2. The Proposed Method

In this section, we present an inferential framework that integrates results from existing studies and results from treatment (i.e., propensity score, PS) regression analyses, in which different and incomplete covariates may be present. We use the reference set on the full set of covariates to convert existing outcome and PS regression analyses into consistent estimates of the full outcome and PS regression model based on the full set of covariates, which in turn generate causal inferences about treatment effects.

2.1. The Existing Outcome and PS Regression Analyses

Suppose that there exist K independent studies, where the kth study contains independent data on $\left\{Y,Z,{\mathbf{X}}_{\left(k\right)}\right\}$ of size $n_k,k=1,2,\dots ,K$. Here, Y and Z are, respectively, the common outcome and binary treatment indicator of interest in the K studies, and ${\mathbf{X}}_{\left(k\right)}$ is the covariate set considered in study k and is some subset of $\mathbf{X}$, the complete covariate set satisfying the strong ignorability condition (1). Suppose that $g_k\left(Y\mid Z,{\mathbf{X}}_{\left(k\right)};{\mathit{\gamma }}_k\right)$ and $h_k\left(Z\mid {\mathbf{X}}_{\left(k\right)};{\mathit{\theta }}_k\right)$ are, respectively, the “reduced” outcome and PS regression models used in the kth study, where ${\mathit{\gamma }}_k$ and ${\mathit{\theta }}_k$ are the respective vectors of regression coefficients in the reduced outcome and the PS models based on the reduced covariate set ${\mathbf{X}}_{\left(k\right)}$. Note that since the reduced PS and outcome models are from existing studies based on reduced sets of the covariate variables, to make the proposed method applicable even when models used in such studies are wrongly specified, we consider the general case where these models can be misspecified. That is, $g_k\left(Y\mid Z,{\mathbf{X}}_{\left(k\right)}\right)$ and $h_k\left(Z\mid {\mathbf{X}}_{\left(k\right)}\right)$ may not equal to the true conditional distributions $P\left(Y\mid Z,{\mathbf{X}}_{\left(k\right)}\right)$ and $P\left(Z\mid {\mathbf{X}}_{\left(k\right)}\right)$, respectively. Assume that the estimates of ${\mathit{\gamma }}_k$ and ${\mathit{\theta }}_k$ are available for study $k,k=1,\dots ,K$. Individual data in these studies are not required in the proposed method.

Further, assume that the true distribution of Z given the complete covariate set $\mathbf{X}$ is given by the full PS model $f_Z(Z\mid \mathbf{X};\mathit{\alpha })$, where $\mathit{\alpha }$ is the vector of regression coefficients of $\mathbf{X}$, and hence the full PS is given by $e\left(\mathbf{X}\right)=f_Z(Z=1\mid \mathbf{X};\mathit{\alpha })$. Also, assume that the true distribution of Y given Z and $\mathbf{X}$ is given by the full outcome regression model $f_Y(Y\mid Z,\mathbf{X};\mathit{\beta })$, where $\mathit{\beta }$ is the vector of regression coefficients of $(Z,\mathbf{X})$. To fix ideas, in the following we consider the underlying distributions $f_Z(Z\mid \mathbf{X};\mathit{\alpha })$ and $f_Y(Y\mid Z,\mathbf{X};\mathit{\beta })$ that follow the generalized linear models (GLMs) of Nelder and Wedderburn [17], although the ideas may simply extend to more general models. We assume that the underlying treatment and outcome distributions $f_Z(Z\mid \mathbf{X};\mathit{\alpha })$ and $f_Y(Y\mid Z,\mathbf{X};\mathit{\beta })$ are essentially the same across different studies, except that their baselines, namely their intercept parameters in $\mathit{\alpha }$ and $\mathit{\beta }$ are allowed to vary across different studies to accommodate possible differences in baseline data distributions among existing studies; however, we do not make this explicit in the notation for the regression parameters to keep the notation simple.

2.2. Reference Data

In addition to the existing regression results, the proposed method needs a reference dataset with data on the complete covariate set.

Specifically, let ${\left\{{\mathbf{X}}_i\right\}}_{i=1}^n$ be the reference sample of n independent observations of $\mathbf{X}$. Since the intercept parameters in the treatment and outcome distributions $f_Z(Z\mid \mathbf{X};\mathit{\alpha })$ and $f_Y(Y\mid Z,\mathbf{X};\mathit{\beta })$ are already allowed to be different among the existing studies, the covariate distributions among the existing studies and the reference sample are assumed to be the same up to location shifts, namely these covariate distributions can have different means but become the same after mean removal. In fact, under the above assumptions, it is possible just to assume the covariate distributions among the studies and the reference dataset are the same, and to adjust the intercept parameters of the underlying treatment and outcome distributions $f_Z(Z\mid \mathbf{X};\mathit{\alpha })$ and $f_Y(Y\mid Z,\mathbf{X};\mathit{\beta })$ in each study, such that $f_Y(Y\mid Z,\mathbf{X};\mathit{\beta })\times f_Z(Z\mid \mathbf{X};\mathit{\alpha })\times f_{\mathbf{X}}\left(\mathbf{X}\right)$ equals the joint distribution of $(Y,Z,\mathbf{X})$ in each study population. It is such adjusted full models $f_Z(Z\mid \mathbf{X};\mathit{\alpha })$ and $f_Y(Y\mid Z,\mathbf{X};\mathit{\beta })$ that are our estimation targets.

2.3. Estimation of the Full Propensity and Outcome Regression Models

The generalized meta-analysis method of Kundu et al. [16] is a data integration method for combining information on parameters of various outcome regression models with disparate covariate sets. Assuming a common underlying data distribution, this method integrates the estimating equations for parameter estimates of various outcome regression models using the generalized method of moments approach [18] and a reference dataset on the covariates to yield consistent estimation for the full outcome regression model.

We apply and extend the generalized meta-analysis method to estimate the parameters in both the full PS and outcome models $f_Z(Z\mid \mathit{X};\mathit{\alpha })$ and $f_Y(Y\mid Z,\mathit{X};\mathit{\beta })$ using the available estimates ${\widehat{\mathit{\theta }}}_k$ and ${\widehat{\mathit{\gamma }}}_k,k=1,\dots K$, from the existing K studies, as well as the reference dataset ${\left\{{\mathbf{X}}_i\right\}}_{i=1}^n$. Following the usual practice, for $k=1,\dots K$, we assume the estimates ${\widehat{\mathit{\theta }}}_k$ and ${\widehat{\mathit{\gamma }}}_k$ are obtained by the maximum likelihood estimation (of the GLMs) based on the (reduced) models $h_k\left(Z\mid {\mathbf{X}}_{\left(k\right)};{\mathit{\theta }}_k\right)$ and $g_k\left(Y\mid Z,{\mathbf{X}}_{\left(k\right)};{\mathit{\gamma }}_k\right)$, $k=1,\dots K$, respectively. Let $t_k\left(Z\mid {\mathbf{X}}_{\left(k\right)};{\mathit{\theta }}_k\right)={\displaystyle \frac{\partial }{\partial {\mathit{\theta }}_k}}log\left(h_k\left(Z\mid {\mathbf{X}}_{\left(k\right)};{\mathit{\theta }}_k\right)\right)$ and $s_k\left(Y\mid Z,{\mathbf{X}}_{\left(k\right)};{\mathit{\gamma }}_k\right)={\displaystyle \frac{\partial }{\partial {\mathit{\gamma }}_k}}log\left(g_k\left(Y\mid Z,{\mathbf{X}}_{\left(k\right)};{\mathit{\gamma }}_k\right)\right)$ be the score functions of the kth reduced PS and outcome models, respectively, and consider the expected scores

$$\begin{array}{c}{t}_k\left(\mathbf{X};\mathit{\alpha },{\mathit{\theta }}_k\right)=E_{Z\mid \mathbf{X}}\left\{t_k\left(Z\mid {\mathbf{X}}_{\left(k\right)};{\mathit{\theta }}_k\right)\right\}=\sum _{Z=0,1}{t}_k\left(Z\mid {\mathbf{X}}_{\left(k\right)};{\mathit{\theta }}_k\right)f_Z(Z\mid \mathbf{X};\mathit{\alpha }),\mathrm{and}\hfill \\ s\left(\mathbf{X};\mathit{\alpha },\mathit{\beta },{\mathit{\theta }}_k,{\mathit{\gamma }}_k\right)=E_{(Y,Z)\mid \mathbf{X}}\left\{s_k\left(Y\mid Z,{\mathbf{X}}_{\left(k\right)};{\mathit{\theta }}_k,{\mathit{\gamma }}_k\right)\right\}\hfill \\ =\sum _{Z=0,1}\int s_k\left(Y\mid Z,{\mathbf{X}}_{\left(k\right)};{\mathit{\theta }}_k,{\mathit{\gamma }}_k\right)f_Y(Y\mid Z,\mathbf{X};\mathit{\beta })f_Z(Z\mid \mathbf{X};\mathit{\alpha })dY.\hfill \end{array}$$

Let

$$\begin{array}{c}T(\mathbf{X};\mathit{\alpha })=\left(\begin{array}{c}{t}_1\left(\mathbf{X};\mathit{\alpha },{\widehat{\mathit{\theta }}}_1\right)\\ ⋮\\ t_K\left(\mathbf{X};\mathit{\alpha },{\widehat{\mathit{\theta }}}_K\right)\end{array}\right),S(\mathbf{X};\mathit{\alpha },\mathit{\beta })=\left(\begin{array}{c}{s}_1\left(\mathbf{X};\mathit{\alpha },\mathit{\beta },{\widehat{\mathit{\theta }}}_1,{\widehat{\mathit{\gamma }}}_1\right)\\ ⋮\\ s_K\left(\mathbf{X};\mathit{\alpha },\mathit{\beta },{\widehat{\mathit{\theta }}}_K,{\widehat{\mathit{\gamma }}}_K\right)\end{array}\right),\mathrm{and}\\ U(\mathbf{X};\mathit{\alpha },\mathit{\beta })=\left(\begin{array}{c}T(\mathbf{X};\mathit{\alpha })\\ S(\mathbf{X};\mathit{\alpha },\mathit{\beta })\end{array}\right).\end{array}$$

The estimator $(\widehat{\mathit{\alpha }},\widehat{\mathit{\beta }})$ for $(\mathit{\alpha },\mathit{\beta })$ is obtained by minimizing the objective function $U_n{(\mathit{\alpha },\mathit{\beta })}^{\mathrm{T}}\mathit{C}{U}_n(\mathit{\alpha },\mathit{\beta })$, or equivalently solving the estimating equation ${\dot{U}}_n{(\mathit{\alpha },\mathit{\beta })}^{\mathrm{T}}\mathit{C}{U}_n(\mathit{\alpha },\mathit{\beta })=0$, where $U_n(\mathit{\alpha },\mathit{\beta })={\sum }_{i=1}^{n}U\left({\mathbf{X}}_i;\mathit{\alpha },\mathit{\beta }\right)$, ${\dot{U}}_n(\mathit{\alpha },\mathit{\beta })=$$\partial U_n(\mathit{\alpha },\mathit{\beta })/\partial (\mathit{\alpha },\mathit{\beta })$. The matrix $\mathbf{C}$ is an arbitrary positive semi-definite weighting matrix and the optimal (minimal asymptotic covariance of the resulting estimator) choice of $\mathbf{C}$ is $\mathbf{C}={(\mathbf{\Delta }+\mathbf{\Lambda })}^{-1}$[16], where $\mathbf{\Delta }$ is the covariance matrix of $U(\mathbf{X};\mathit{\alpha },\mathit{\beta })$, and $\mathbf{\Lambda }$ arises from the covariances of $\sqrt{n_k}\left({\widehat{\mathit{\theta }}}_k,{\widehat{\mathit{\gamma }}}_k\right)$, $k=1,\dots ,K$, which may be available from the existing studies or estimated using the reference dataset. The way of estimation of the matrices $\mathbf{\Delta }$ and $\mathbf{\Lambda }$ is given in Appendix A.

Following similar arguments in Kundu et al. [16], $\sqrt{n}\{(\widehat{\mathit{\alpha }},\widehat{\mathit{\beta }})-(\mathit{\alpha },\mathit{\beta })\}$ has the asymptotic normal distribution with zero mean and the covariance matrix given by ${\left\{{\mathbf{\Gamma }}^{\mathrm{T}}{(\mathbf{\Delta }+\mathbf{\Lambda })}^{-1}\mathbf{\Gamma }\right\}}^{-1}$, where $\mathbf{\Gamma }=E\{\partial U(\mathbf{X};\mathit{\alpha },\mathit{\beta })/\partial (\mathit{\alpha },\mathit{\beta }\left)\right\}$ whose estimation can be based on its sample analog in the reference dataset, and the matrices $\mathbf{\Delta }$ and $\mathbf{\Lambda }$ are estimated in the way mentioned in Appendix A using the final estimator $(\widehat{\mathit{\alpha }},\widehat{\mathit{\beta }})$ for $(\mathit{\alpha },\mathit{\beta })$.

We conclude this section by noting that, the covariance matrices of the regression parameter estimates from the existing studies are not necessary for implementing the proposed method. However, when the covariance matrix for the regression parameter estimates is unavailable and the outcome regression model in an existing study contains a dispersion parameter not fixed to 1 (e.g., the normal regression model), the estimate for that dispersion parameter is required for the proposed method to implement estimation of the covariance matrix of the regression parameters (see Appendix A for detail).

2.4. Estimation of the Average Treatment Effect (ATE)

Let $Q\left(\mathbf{X}\right)=E(Y\mid Z=1,\mathbf{X})-E(Y\mid Z=0,\mathbf{X})$. By the SUTVA and the strong ignorability assumption (1), $Q\left(\mathbf{X}\right)=E\left(Y_1\mid Z=1,\mathbf{X}\right)-E\left(Y_0\mid Z=0,\mathbf{X}\right)=E\left(Y_1\mid \mathbf{X}\right)-E\left(Y_0\mid \mathbf{X}\right)$ corresponds to the conditional average treatment effect (CATE) at $\mathbf{X}$, and $E_{\mathbf{X}}\left\{Q\left(\mathbf{X}\right)\right\}=E\left(Y_1\right)-E\left(Y_0\right)\equiv \mu$ corresponds to the ATE.

Let $Q(\mathbf{X};\mathit{\beta })=\int Y\left\{f_Y(Y\mid Z=1,\mathbf{X};\mathit{\beta })-f_Y(Y\mid Z=0,\mathbf{X};\mathit{\beta })\right\}dY$. Given the estimate $\widehat{\mathit{\beta }}$ for the full outcome regression parameter obtained in Section 2.3, we can estimate $E(Y\mid Z,\mathbf{X})$ by $\int Y{f}_Y(Y\mid Z,\mathbf{X};\widehat{\mathit{\beta }})dY$, the CATE $Q\left(\mathbf{X}\right)$ at $\mathbf{X}$ by $Q(\mathbf{X};\widehat{\mathit{\beta }})$, and the ATE by the reference sample data ${\left\{{\mathbf{X}}_i\right\}}_{i=1}^n$ via

$$\widehat{\mu }={\displaystyle \frac{1}{n}}\sum _{i=1}^{n}Q\left({\mathbf{X}}_i;\widehat{\mathit{\beta }}\right).$$

The consistency of the proposed ATE estimator $\widehat{\mu }$ follows directly from the consistency of $\widehat{\mathit{\beta }}$. Also, using the delta method, we can obtain the asymptotic normal distribution of $\widehat{\mu }$; details are provided in Appendix B.

3. The Simulation Studies

In this section, we conduct simulations to assess the performance of the proposed estimators for the regression coefficients of the underlying outcome and treatment distributions, as well as for the ATE. We specifically report results for bias, standard deviation (SD), mean of estimated standard errors (ESE), and coverage probability 1 and 2 (CP1 and CP2) of the $95\%$ and $90\%$ Wald-type confidence intervals, respectively, across 2000 simulation replications.

3.1. The Simulation Setting

We consider a simple simulation setting where there exist independent datasets $D_1$ on $\left\{Y,Z,X_1\right\},D_2$ on $\left\{Y,Z,X_2\right\}$, and $D_3$ on $\left\{X_1,X_2\right\}$, where Y is the common outcome of interest, Z is the common treatment indicator variable, $X_1$ is the covariate observed in dataset $D_1,X_2$ is the covariate observed in dataset $D_2$, and $\left(X_1,X_2\right)$ is the set of covariates in the underlying full regression models for the outcome and the treatment (i.e., the PS). We set the sample sizes of the datasets $D_1$, $D_2$, and $D_3$ to $n_1=500,n_2=500,n=50$, 100 or 200, respectively.

The covariates $\left(X_1,X_2\right)$ in $D_1,D_2$, and $D_3$ are generated by independent standard normal random variables, and given $\left(X_1,X_2\right)$, the treatment assignment Z in both $D_1$ and $D_2$ is generated using the linear logistic regression model: $P\left(Z=1\mid X_1,X_2\right)=expit\left({\alpha }_0+{\alpha }_1{X}_1+{\alpha }_2{X}_2\right)$ with $expit\left(x\right)=exp\left(x\right)/(1+$$exp\left(x\right))$. The outcome Y in both $D_1$ and $D_2$ is binary and generated by the linear logistic regression model for given treatment assignment Z and covariates $\left(X_1,X_2\right)$: $P\left(Y=1\mid Z,X_1,X_2\right)=expit\left({\beta }_0+{\beta }_{z}Z+{\beta }_1{X}_1+{\beta }_2{X}_2\right)$. The coefficients for the covariates, $\left({\alpha }_1,{\alpha }_2\right)$ and $\left({\beta }_1,{\beta }_2\right)$, in the PS and outcome models are both fixed at $(log(1.2),log(1.5\left)\right)$, and the coefficient ${\beta }_Z$ for the treatment indicator is set to 0, $log\left(1.2\right),log\left(1.5\right)$, or $log\left(2\right)$ when generating both datasets $D_1$ and $D_2$. On the other hand, the intercept parameters $\left({\alpha }_0,{\beta }_0\right)$ are set to $(-1,-1)$ and $(0,-1.5)$ when generating datasets $D_1$ and $D_2$, respectively, to reflect different baselines in the underlying distributions of the two datasets. The reduced PS and outcome regression coefficient estimates are obtained from both $D_1$ and $D_2$; in $D_1$, the working PS and outcome regression models are, respectively, the linear logistic models for $\left(Z\mid X_1\right)$ and $\left(Y\mid Z,X_1\right)$, while in $D_2$ they are, respectively, the linear logistic models for $\left(Z\mid X_2\right)$ and $\left(Y\mid Z,X_2\right)$. Such reduced-model regression coefficient estimates are used in the later estimation procedure while individual data in $D_1$ and $D_2$ are not. Individual data in $D_3$ on the complete covariate set $\left(X_1,X_2\right)$ are used as the reference data for the proposed estimation.

The method in Section 2.3 is applied to the regression coefficient estimates from $D_1$ and $D_2$ and the individual data from $D_3$ to obtain the estimates of the parameters in the full PS and outcome regression models, which are correctly specified in the estimation as the data generating models mentioned above. Further, the method in Section 2.4 is applied to estimate the ATE of Z on Y.

We also perform an extended simulation with $K=4$ existing datasets over 1000 simulation replications, and the distribution of the outcome Y is binomial, normal, or Poisson. The complete covariate set is $\mathbf{X}=\left(X_1,X_2,X_3,X_4\right)$, where $\left(X_1,X_2\right)$ are generated from standard normal with a correlation coefficient of $0.5$, and $\left(X_3,X_4\right)$ are Bernoulli random variables with success probabilities expit$(X_1+E_1)$ and expit$(X_2+E_2)$, where $X_1$, $X_2$ are the covariates mentioned above, $E_1$, $E_2$ are independent standard normal, and expit$\left(x\right)=exp\left(x\right)/\{1+exp(x\left)\right\}$. Accordingly, $X_3$ is correlated with $X_1$, and $X_4$ is correlated with $X_2$ and hence also $X_1$ since $X_1$ and $X_2$ are correlated. The four existing datasets contain data on the common outcome variable Y and treatment variable Z, and data on the different covariate sets $\left(X_1,X_2\right)$, $\left(X_2,X_3\right)$, $\left(X_3,X_4\right)$, and $\left(X_1,X_4\right)$, respectively. The treatment variable Z in the kth dataset is generated using the linear logistic regression model: $P_k\left(Z=1\mid X_1,X_2,X_3,X_4\right)=expit\left({\alpha }_{0k}+{\alpha }_1{X}_1+{\alpha }_2{X}_2+{\alpha }_3{X}_3+{\alpha }_4{X}_4\right)$, $\left({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4\right)=\left(log\left(1.2\right),log\left(1.5\right),-log\left(1.2\right),-log\left(1.5\right)\right)$, and ${\alpha }_{0k}=-1,-1.5,-2,-2.5$ for $k=1,2,3,4$. The outcome Y in the kth dataset is generated using the linear logistic regression model: $P_k\left(Y=1\mid Z,X_1,X_2,X_3,X_4\right)=g\left({\beta }_{0k}+{\beta }_{z}Z+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_3{X}_3+{\beta }_4{X}_4\right)$ with g the expit, identity, or exponential function when the distribution of Y is binomial, normal, or Poisson (i.e., the inverse canonical link function), respectively. The true parameter values $\left({\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4\right)=\left(log\left(1.2\right),log\left(1.5\right),-log\left(1.5\right),-log\left(1.2\right)\right)$, ${\beta }_{0k}=-1,-1.5,-2,-2.5$ for $k=1,2,3,4$, and ${\beta }_z=0$ or $log\left(1.5\right)$. The reduced propensity (PS) and outcome regression coefficient estimates are obtained from the four existing datasets, using the models having the same link functions as in the data-generating models but the reduced covariate sets observed in these datasets, as mentioned above. Such reduced regression coefficient estimates are used in the estimation procedure while individual data in the existing datasets are not. Individual data in the reference sample, which contain data on the complete covariate set $\mathbf{X}=\left(X_1,X_2,X_3,X_4\right)$, are also used in the estimation procedure as proposed in Section 2.3. The full PS and outcome regression models in the estimation are correctly specified as the data-generating models mentioned above. Further, the method in Section 2.4 is applied to estimate the average treatment effect (ATE) of Z on Y in the population of each dataset. All four existing studies and the reference sample have the same sample size, set to 200, 500, or 1000.

3.2. The Simulation Results

Table 1. Simulation results (multiplied by 100) for the estimates of the coefficient ${\beta }_Z$ of the treatment variable in the full outcome model and the ATEs with true parameter values ${\beta }_0=-1$ in Study 1, and $-1.5$ in study $2,\left({\beta }_{X_1},{\beta }_{X_2}\right)=(log\left(1.2\right),log\left(1.5\right))$, and different ${\beta }_Z$ and ATE values (in parenthesis), $\mathit{\alpha }=\left({\alpha }_0,{\alpha }_1,{\alpha }_2\right)=(-1,log\left(1.2\right),log\left(1.5\right))$ in Study 1, and $(0,log(1.2),log(1.5\left)\right)$ in Study 2.

	$\mathit{n}=50$					$\mathit{n}=100$					$\mathit{n}=200$
	Bias	SD	ESE	CP1	CP2	Bias	SD	ESE	CP1	CP2	Bias	SD	ESE	CP1	CP2
${\beta }_Z$ (0)	$-0.19$	$17.0$	$17.7$	$96.5$	$91.6$	$-0.32$	$17.7$	$17.5$	$95.1$	$89.8$	$-0.94$	$17.3$	$17.4$	$95.1$	$89.7$
ATE1 $\left(0\right)$	$0.03$	$3.23$	$3.37$	$96.4$	$91.5$	$0.01$	$3.37$	$3.34$	$95.0$	$89.8$	$-0.11$	$3.29$	$3.32$	$94.8$	$89.5$
ATE2 $\left(0\right)$	$-0.01$	$2.53$	$2.64$	$96.5$	$91.5$	$-0.04$	$2.62$	$2.60$	$95.1$	$89.8$	$-0.14$	$2.58$	$2.59$	$95.0$	$89.6$
${\beta }_Z$ $(log(1.2\left)\right)$	$-0.08$	$16.5$	$17.2$	$96.2$	$91.7$	$-0.21$	$17.2$	$17.0$	$95.1$	$89.2$	$-0.53$	$17.0$	$16.9$	$94.7$	$89.8$
ATE1 $\left(0.04\right)$	$-0.00$	$3.33$	$3.47$	$96.3$	$91.4$	$-0.00$	$3.48$	$3.44$	$94.9$	$89.4$	$-0.07$	$3.41$	$3.42$	$94.6$	$89.6$
ATE2 $\left(0.03\right)$	$-0.02$	$2.61$	$2.72$	$95.9$	$91.8$	$-0.05$	$2.70$	$2.68$	$95.3$	$89.3$	$-0.10$	$2.67$	$2.67$	$94.8$	$89.9$
${\beta }_Z$ $(log(1.5\left)\right)$	$-0.13$	$16.2$	$16.8$	$96.1$	$91.5$	$-0.14$	$16.4$	$16.6$	$95.4$	$90.2$	$-0.35$	$16.4$	$16.5$	$94.8$	$90.0$
ATE1 $\left(0.08\right)$	$-0.07$	$3.46$	$3.59$	$96.4$	$91.4$	$-0.04$	$3.52$	$3.55$	$95.5$	$90.2$	$-0.09$	$3.49$	$3.52$	$94.5$	$89.8$
ATE2 $\left(0.07\right)$	$-0.06$	$2.74$	$2.82$	$96.1$	$90.9$	$-0.08$	$2.75$	$2.78$	$95.6$	$91.0$	$-0.10$	$2.75$	$2.77$	$95.3$	$89.6$
${\beta }_Z$ $(log(2\left)\right)$	$0.04$	$15.8$	$16.3$	$96.2$	$91.5$	$0.15$	$15.9$	$16.2$	$95.5$	$90.1$	$-0.17$	$15.9$	$16.1$	$95.6$	$90.3$
ATE1 $\left(0.15\right)$	$-0.12$	$3.56$	$3.69$	$96.2$	$91.2$	$-0.03$	$3.59$	$3.65$	$95.0$	$90.0$	$-0.12$	$3.59$	$3.62$	$95.1$	$90.0$
ATE2 $\left(0.12\right)$	$-0.09$	$2.87$	$2.96$	$95.9$	$90.6$	$-0.08$	$2.85$	$2.91$	$95.4$	$90.2$	$-0.12$	$2.85$	$2.89$	$95.6$	$90.1$

n, size of reference data; SD, standard deviation; ESE, estimated standard error; CP1 and CP2, coverage probability of $95\%$ and $90\%$ confidence intervals, respectively.

presents the simulation results for the proposed estimation of the regression coefficient for the treatment variable in the full outcome model and the ATEs for the populations of the two existing datasets. The findings in Table 1 suggest that the proposed estimation for the full outcome model parameter and the ATE is essentially unbiased, with the absolute bias of the proposed estimates for both parameters being less than 1%. Also, the estimated standard error (ESE) based on the asymptotic theory is close to the simulation standard deviation of the estimator with the absolute difference less than 1%, and the coverage probability of the $95\%$ and $90\%$ Wald-type confidence intervals based on the asymptotic normality of the estimator is close to the nominal levels $0.95$ and $0.90$, respectively. These results reveal that the proposed estimators perform well in finite samples.

Table 2. Simulation results (multiplied by 100) for the estimates of the coefficients $\left({\alpha }_{X_1},{\alpha }_{X_2}\right)$ of the covariate variables in the full propensity score model with true parameter values ${\beta }_0=-1$ in Study 1, and $-1.5$ in study $2,\left({\beta }_{X_1},{\beta }_{X_2}\right)=(log\left(1.2\right)$, $log\left(1.5\right))$, and different ${\beta }_Z$ values, $\mathit{\alpha }=\left({\alpha }_0,{\alpha }_{X_1},{\alpha }_{X_2}\right)=(-1,log\left(1.2\right),log\left(1.5\right))$ in Study 1, and $(0,log(1.2),log(1.5\left)\right)$ in Study 2.

	$\mathit{n}=50$					$\mathit{n}=100$					$\mathit{n}=200$
	Bias	SD	ESE	CP1	CP2	Bias	SD	ESE	CP1	CP2	Bias	SD	ESE	CP1	CP2
${\beta }_Z=0$
${\alpha }_{{\mathrm{X}}_1}$	$0.99$	$13.0$	$12.8$	$96.0$	$90.3$	$1.19$	$12.0$	$11.7$	$95.7$	$89.5$	$0.32$	$11.0$	$11.1$	$95.2$	$90.2$
${\alpha }_{X_2}$	$1.29$	$10.8$	$10.9$	$96.4$	$91.6$	$0.91$	$10.1$	$10.2$	$95.6$	$90.9$	$0.44$	$9.82$	$9.84$	$95.5$	$90.3$
${\beta }_Z=log\left(1.2\right)$
${\alpha }_{X_1}$	$0.99$	$13.0$	$12.8$	$96.2$	$90.4$	$1.20$	$12.0$	$11.7$	$95.6$	$89.6$	$0.33$	$11.0$	$11.1$	$95.2$	$90.3$
${\alpha }_{X_2}$	$1.29$	$10.8$	$10.9$	$96.3$	$91.7$	$0.91$	$10.1$	$10.2$	$95.6$	$90.9$	$0.44$	$9.82$	$9.84$	$95.6$	$90.5$
${\beta }_Z=log\left(1.5\right)$
${\alpha }_{X_1}$	$1.00$	$13.0$	$12.8$	$96.2$	$90.2$	$1.20$	$12.0$	$11.7$	$95.5$	$89.5$	$0.33$	$11.0$	$11.1$	$95.3$	$90.3$
${\alpha }_{X_2}$	$1.30$	$10.8$	$10.9$	$96.3$	$91.9$	$0.93$	$10.1$	$10.2$	$95.9$	$91.1$	$0.44$	$9.82$	$9.84$	$95.4$	$90.3$
${\beta }_Z=log\left(2\right)$
${\alpha }_{X_1}$	$1.00$	$13.0$	$12.8$	$96.0$	$90.4$	$1.20$	$12.0$	$11.7$	$95.7$	$89.5$	$0.33$	$11.0$	$11.1$	$95.3$	$90.3$
${\alpha }_{X_2}$	$1.32$	$10.8$	$10.9$	$96.1$	$91.7$	$0.94$	$10.1$	$10.2$	$95.9$	$91.1$	$0.46$	$9.82$	$9.83$	$95.3$	$90.3$

n, size of reference data; SD, standard deviation; ESE, estimated standard error; CP1 and CP2, coverage probability of $95\%$ and $90\%$ confidence intervals, respectively.

presents the simulation results for the proposed estimation of the regression coefficients in the full PS model under the same simulation settings as in Table 1. The findings demonstrate that our proposed method performs well in the estimation of the full PS model.

The extended simulations, performed under the settings with $K=4$ studies, correlated covariates, and binomial, normal, and Poisson distributed outcome variables (see Section 3.1 for detail), still reveal satisfactory performances of the proposed estimation method; the results are shown in the Appendix C.

The performance of the proposed method in terms of computation time is summarized as follows. When $K=2$, the computation time for a simulation case with $n=200$ is $0.5$ s in a desktop computer with i7-9700 CPU, and the time increases to $1.2$ s when $K=4$ and $n=200$. Essentially, the computation time increases linearly with K. Also, when $K=4$, the computation time for running a case in the simulation study is $5.0$ s when the size of the reference sample increases to $n=1000$.

4. A Real Data Application

In this section, we apply the proposed method to analyze the impact of waist circumference (WC) on hypertension (HT) risk, while controlling for covariates such as age, sex, body mass index (BMI), smoking status (SMK), drinking status (DRK), body fat percentage (BFP), heart rate (HR), and hip circumference (HC) among working-age individuals. The analysis leverages regression analyses from two existing studies on the relationship between WC and HT and a reference dataset encompassing the complete set of covariates.

4.1. Two Existing Studies on the Effect of Waist Circumference on Hypertension

The WC reflects the size of the visceral fat depot and is an effective clinical tool for assessing the risk of diabetes and cardiovascular diseases [19]. Guagnano et al. [20] indicated that WC seems to have a strong association with the risk of hypertension. In recent years, Ren et al. [14] investigated the cut-off values for the obesity indices that represent the elevated incidence of hypertension in Chinese adults aged between 18 and 65. Hu et al. [15] indicated that a combination of WC and BMI was superior to individual indices for identifying hypertension. Data from Ren et al. [14] and Hu et al. [15] were publicly provided (https://dx.doi.org/10.6084/m9.figshare.2151271.v1 (accessed on 15 February 2016), https://doi.org/10.1371/journal.pone.0170238.s001 (accessed on 5 January 2017)). The study of Hu et al. [15], termed Study 1, contains data on the covariate set ${\mathbf{X}}_{\left(1\right)}$ including age, sex, BMI, SMK, DRK, BFP, and HR, while the study of Ren et al. [14], termed Study 2, contains data on the covariate set ${\mathbf{X}}_{\left(2\right)}$ including age, sex, BMI, SMK, DRK, and HC. We focus on working-age (20–65 years old) people and the subsamples from the two studies meeting this criterion are of sizes $n_1=9926$ and $n_2=2970$, respectively, (after removing missing observations).

4.2. The Reference Dataset with Complete Covariates

The Taiwan BioBank (TWB) database, created by Academia Sinica, comprises a community-based cohort of over 200,000 study participants. It includes comprehensive data on demographics, health behaviors, environmental factors, and biomarkers collected through meticulously conducted questionnaires and thorough examinations. Details about the TWB data can be found at https://www.twbiobank.org.tw/ (accessed on 1 June 2024). The reference dataset we employ in the current analysis is based on the released subsample of the TWB cohort consisting of 4575 randomly sampled study subjects aged 20–65 years. The reference dataset contains data on the complete covariate set $\mathbf{X}$ including age, sex, BMI, SMK, DRK, BFP, HR, and HC, but contains no data on either the treatment (WC) or the outcome (HT).

4.3. The Proposed Analysis

In the following analysis, both WC and HT are defined as binary variables, classified according to $\mathrm{WC}=\mathrm{I}(\mathrm{waistcircumference}>80\mathrm{cm})$ and $\mathrm{HT}=\mathrm{I}(\mathrm{SBP}\ge 140$ mmHg or $\mathrm{DBP}\ge 90\mathrm{mmHg})$, where I(.) is the indicator function and SBP and DBP denote systolic and diastolic blood pressures, respectively; the classification rules follow those in Lean et al. (1998) [19]. The covariates age (years), BMI (kg/m${ }^2$), BFP (%), HR (beats/min), and HC (cm) are continuous variables, while the covariates sex (female vs. male), SMK (yes vs. no), and DRK (yes vs. no) are binary.

We apply the proposed methods in Section 2.3 and Section 2.4 to assess the treatment effect of WC on the risk of HT controlling for the covariates age, sex, BMI, SMK, DRK, BFP, HR, and HC, which are regarded as the complete covariate set among working-age people. Specifically, the proposed analysis uses the results of the regression analyses from Study 1 (Hu et al. [15]) and Study 2 (Ren et al. [14]), as well as the reference dataset from the TWB database. Both the analyses in Studies 1 and 2 employ linear logistic regressions to examine the association between WC and HT by adjusting the covariate sets ${\mathbf{X}}_{\left(1\right)}$ and ${\mathbf{X}}_{\left(2\right)}$, respectively; see Section 4.1 for detail. Also, both PS analyses in the two studies are based on the linear logistic regressions for $\mathrm{WC}$ with the covariate sets ${\mathbf{X}}_{\left(1\right)}$ and ${\mathbf{X}}_{\left(2\right)}$, respectively.

In the proposed analysis, only the outcome and regression parameter estimates from Studies 1 and 2 are employed, while individual data are not. The full outcome (HT) model is specified by the linear logistic regression model for HT given the treatment (WC) and the complete covariate set $\mathbf{X}$, and the full PS model is specified as the linear logistic regression model for WC given $\mathbf{X}$; in these logistic regression models only the main effects of the treatment and the covariate variables are considered. To account for possible differences between the baselines of Studies 1 and 2, the intercept parameters of the full models, including the outcome and the PS models, across the studies are treated as different.

4.4. The Analysis Results

The results for the proposed estimation of the logistic regression models for the PS (treatment, WC) and the outcome (HT) adjusting for the complete covariate set are provided in

Table 3. Results of the real-data analysis. The propensity score (PS) for waist circumference (WC) with the covariates age, sex, body mass index (BMI), smoking status (SMK), drinking status (DRK), body fat percentage (BFP), heart rate (HR), and hip circumference (HC) based on results from two studies and a reference sample.

	PS Model
	Study 1	Study 2	Proposed
	${\mathit{n}}_{\mathbf{1}}\mathbf{=}\mathbf{9926}$	${\mathit{n}}_{\mathbf{1}}\mathbf{=}\mathbf{2970}$	$\mathit{n}\mathbf{=}\mathbf{4575}$
	Est (SE)	Est (SE)	Est (SE)
Age	0.025 (0.002) *	0.041 (0.006) *	0.046 (0.003) *
Sex	−1.234 (0.074) *	−0.864 (0.150) *	−1.237 (0.157) *
BMI	0.597 (0.013) *	0.301 (0.021) *	0.310 (0.028) *
SMK	0.009 (0.073)	0.055 (0.150)	0.052 (0.073)
DRK	0.065 (0.067)	0.229 (0.129)	0.118 (0.097)
BFP	0.035 (0.004) *	-	0.025 (0.013)
HR	0.010 (0.002) *	-	0.013 (0.004) *
HC	-	0.217 (0.011) *	0.217 (0.017) *

Est: parameter estimate; SE: standard error; *: p-value $<0.05$.

and

Table 4. Results of the real-data analysis. The risk of hypertension (HT) with treatment of waist circumference (WC) adjusting for age, sex, body mass index (BMI), smoking status (SMK), drinking status (DRK), body fat percentage (BFP), heart rate (HR), and hip circumference (HC) based on results from two studies and a reference sample.

	Outcome Model
	Study 1	Study 2	Proposed
	${\mathit{n}}_{\mathbf{1}}\mathbf{=}\mathbf{9926}$	${\mathit{n}}_{\mathbf{1}}\mathbf{=}\mathbf{2970}$	$\mathit{n}\mathbf{=}\mathbf{4575}$
	Est (SE)	Est (SE)	Est (SE)
WC	0.418 (0.067) *	0.047 (0.131)	0.323 (0.067) *
Age	0.085 (0.003) *	0.061 (0.006) *	0.080 (0.003) *
Sex	−0.196 (0.074) *	−0.382 (0.151) *	−0.251 (0.126) *
BMI	0.101 (0.010) *	0.163 (0.020) *	0.109 (0.024) *
SMK	0.027 (0.071)	0.115 (0.146)	0.044 (0.064)
DRK	−0.060 (0.066)	−0.035 (0.131)	−0.053 (0.084)
BFP	0.005 (0.004)	-	0.005 (0.010)
HR	0.021 (0.002) *	-	0.021 (0.003) *
HC	-	0.002 (0.010)	0.003 (0.013)

Est: parameter estimate; SE: standard error; *: p-value $<0.05$.

, respectively. We can see from Table 3 that, older age, male, larger BMI, higher heart rate, and larger hip circumference tend to have a waist circumference greater than 80 cm (treatment group), and the estimation result seems to nicely summarize, synthesize, and complement the results from the two existing studies. Also, we see from Table 4 that, after adjusting for the covariates age, sex, body mass index (BMI), smoking (SMK), drinking (DRK), body fat percentage (BFP), heart rate (HR), and hip circumference (HC), the effect of waist circumference (WC) on the risk of hypertension (HT) is strongly significant; the odds for hypertension in working-age people with waist circumference greater than 80 cm is 1.4 (≈$exp\left(0.323\right)$) times as high as those with waist circumference no greater than 80 cm (p-value $<0.001$). In contrast, the effect of WC on the risk of HT obtained by adjusting an incomplete covariate set can be somewhat higher (when adjusting only for age, sex, BMI, SMK, DRK, BFP, and HR in Study 1) or lower (when adjusting only for age, sex, BMI, SMK, DRK, and HC in Study 2). Since HR has its own significant effects on both WC and HT, the lower effect of WC on HT without adjusting for HR obtained in Study 2 is likely to be biased.

The average treatment effect of WC on HT, averaged over the covariate distribution, is obtained as 0.044 with SE (standard error) $=0.009$, p-value $<0.001$, in Study 1 population, and as 0.042 with SE $=0.009$, p-value $<0.001$, in Study 2 population. That is overall working-age people with waist circumference larger than 80 cm can have 44 (42) additional cases of HT per 1000 people in the Study 1 (2) population, 95% confidence interval 26–62 (24–60), compared to those who with waist circumference no larger than 80 cm.

From the results mentioned above and shown in Table 3, we conclude that the proposed integrated analysis, using information from both studies and the reference data with complete covariates, can lead to less biased and possibly more efficient analysis results than those from the original individual studies.

5. Discussion and Conclusions

In this study, we propose a new inference framework that integrates the results of the outcome and the treatment (i.e., the PS) regression analyses across multiple existing databases. These databases may vary in their coverage of covariate variables and may contain incomplete data, potentially introducing bias in individual database analyses. Moreover, access to individual-level data from these databases may be restricted. Our proposal integrates the existing PS and outcome regression analyses through a reference sample, which contains only data on the complete covariate set. We obtain consistent estimates for the parameters of the full PS and outcome regression models and for the ATE. The new proposal extends the original generalized meta-analysis method of Kundu et al. [16] by further considering the treatment (propensity score) regression in addition to the outcome regression. Also, the new proposal can apply with a general outcome variable, such as one following a generalized linear model, and hence is more flexible than the work of Li et al. [21], which considers the setting similar to ours but is restricted to normality outcome and linear regression model.

Our approach necessitates a dataset with comprehensive covariate information, which acts as a benchmark for the underlying covariate distribution [16]. Such a reference dataset could be sourced from a large-scale database like the Taiwan Biobank, as outlined in Section 4 of our application. Alternatively, a reference sample might be gathered through a smaller validation study, a method commonly discussed in the epidemiological literature (e.g., Stümer et al. [22]).

As in the existing methods, such as Kundu et al. [16], for integrating common information from different studies, we require the underlying treatment and outcome distributions to be the same across various studies. When this assumption is not satisfied, we should interpret the parameter estimates from the proposed method with caution, since they no longer represent consistent estimates for some common parameters, but instead represent the estimates for some “average effects” over different studies.

In summary, our proposal is the best applicable in the following two scenarios: (1) A multi-center study where individual data from each of the centers are not accessible except for the derived summary statistics (e.g., regression coefficient estimates), and an independent reference sample of complete covariate data is available. (2) A meta-analysis where results for both the outcome and the treatment (propensity score) regression analyses are available for various studies, together with a reference sample of complete covariate data. In both scenarios, our approach optimally integrates analysis results from diverse data sources to yield valid inferences on treatment effects using summarized and synthesized information.