Can statistical adjustment guided by causal inference improve the accuracy of effect estimation? A simulation and empirical research based on meta-analyses of case–control studies

Background Statistical adjustment is often considered to control confounding bias in observational studies, especially case–control studies. However, different adjustment strategies may affect the estimation of odds ratios (ORs), and in turn affect the results of their pooled analyses. Our study is aimed to investigate how to deal with the statistical adjustment in case–control studies to improve the validity of meta-analyses. Methods Three types of adjustment strategies were evaluated including insufficient adjustment (not all preset confounders were adjusted), full adjustment (all confounders were adjusted under the guidance of causal inference), and improper adjustment (covariates other than confounders were adjusted). We carried out a series of Monte Carlo simulation experiments based on predesigned scenarios, and assessed the accuracy of effect estimations from meta-analyses of case–control studies by combining ORs calculated according to different adjustment strategies. Then we used the data from an empirical review to illustrate the replicability of the simulation results. Results For all scenarios with different strength of causal relations, combining ORs that were comprehensively adjusted for confounders would get the most precise effect estimation. By contrast, combining ORs that were not sufficiently adjusted for confounders or improperly adjusted for mediators or colliders would easily introduce bias in causal interpretation, especially when the true effect of exposure on outcome was weak or none. The findings of the simulation experiments were further verified by the empirical research. Conclusions Statistical adjustment guided by causal inference are recommended for effect estimation. Therefore, when conducting meta-analyses of case–control studies, the causal relationship formulated by exposure, outcome, and covariates should be firstly understood through a directed acyclic graph, and then reasonable original ORs could be extracted and combined by suitable methods.


Supplementary Methods
Method S1. Interpretation of directed acyclic graph in the target population 4 Method S2. Generation of target population 5 Method S3. Generation of case-control studies 8 Method S4. Generation of meta-analyses 10

Supplementary Results
Result S1. Interpretation of effect estimations in scenario Ref 12 Result S2. Interpretation of performances of statistical adjustment strategies 13 Result S3. Construction of directed acyclic graph between passive smoking and breast cancer 15 Supplementary Tables   Table S1. The simulated distributions of exposure, outcome, and covariates in the target population of scenario Ref 17 Table S2. Scenario settings of the simulation 18 Table S3. Performances of statistical adjustment strategies in scenarios with different total effect of exposure on outcome (ORAY) in target population 20 Table S4. Performances of statistical adjustment strategies in scenarios with different independent associations of covariates with exposure (ORUA) in target population 24 Table S5. Performances of statistical adjustment strategies in scenarios with different independent associations of covariates with outcome (ORUY) in target population 28  References 52

Method S1. Interpretation of directed acyclic graph in the target population
Suppose that A is the exposure and Y is the outcome of the target population. The directed acyclic graph (DAG) between A and Y is shown in Figure 1. researchers adjust all covariates that show significance in correlation analyses without regarding any background knowledge, the effect estimations of case-control studies will be challenged. Not to mention that the significant results are also related to the random variability, sample size, collinearity, etc., which will further increase the uncertainty.
Combining these indeterminate results in a meta-analysis is apparently problematic.
Therefore, when doing meta-analyses, researchers should correctly identify the "accurate" causal effects from case-control studies and pool them together by suitable methods; other "inaccurate" original ORs should be analyzed in sensitivity analysis.

Method S2. Generation of target population
Generation of target population was composed of four steps. The first step was to simulate confounder L = [L1, L2, …, L6], which was the ancestor of all variables (as shown in Figure 1). By defining the positive probabilities and correlation matrix of dichotomous variables L1 to L6, a multivariate binomial distribution would be established through the Emrich-Piedmonte algorithm [1], and a population with 10,000 observations would be created. In detail, the multivariate binomial distribution of L was generated using a proxy multivariate normal distribution, the multivariate normal distribution was specified by a mean vector and a variance-covariance matrix, and the components of the variance-covariance matrix were determined by solving equation parameter rUU for all j and k; rUU was specified as 0 in the reference scenario, and as 0.2, 0.5, or 0.8 in scenarios for testing different covariate correlations), and and = 1 − are the positive and negative probabilities of variable , respectively ( was assumed to be equal to 0.2 for all j). The multivariate normal distribution was established from the above-calculated parameters and was then dichotomized to multivariate binominal distribution using the cutoff value of ( ).
The second step was to simulate the children of L, i.e., exposure A and risk factor R by logistic regression models. With the dichotomous variables L1 to L6 generated in step 1, the predicted probabilities of A for 10,000 individuals would be obtained by equation logit( ) = 0 + ∑ 6

=1
, where denotes the i th individual, intercept 0 = ln(P( = 1| = )) ln(P( = 0| = )) ⁄ reflects the hypothetical prevalence of A in individuals unexposed to L (P( = 1| = ) was assumed to be 0.2), and = ln (OR ) , = 1, … , 6 reflect the independent associations of with A (OR was assumed to be equal to a parameter ORUA for all j; ORUA was specified as 2 in the reference scenario, and as 0. to have no causal effect on R), where 0 = ln(P( = 1| = )) ln(P( = 0| = )) ⁄ reflects the hypothetical prevalence of R in individuals unexposed to L (P ( = 1| = ) was assumed to be 0.2), and = ln (OR ) , = 1, … , 4 reflect the independent associations of L with R (OR was assumed to be equal to ORUA for all j).
The third step was to deal with the mediation effect from exposure A to outcome  The last step was to simulate collider C, which was the descendent of all variables, by equation logit( ) = 0 + 1 + 2 , where 0 = ln(P( = 1| = 0, = 0)) ln(P( = 0| = 0, = 0)) ⁄ reflects the hypothetical prevalence of C in individuals unexposed to A and Y (P( = 1| = 0, = 0) was assumed to be 0.2), 1 = ln(OR ) reflects the effect of A on C (ORAC was assumed to be equal to ORUA), and 2 = ln(OR ) reflects the effect of Y on C (ORYC was assumed to be equal to ORUY).
Through the four-step process, a target population with 10,000 records and 11 variables would be generated. The simulated distributions of exposure, outcome, and covariates were close to what we specified (as shown in Supplementary Table S1).

Method S3. Generation of case-control studies
Observations of case-control studies were randomly sampled from the target population with pre-specified sample size and matching approach. For frequency matched studies, a same number of cases and controls were selected from the strata Y=1 and 0, respectively, using stratified random sampling method (the number of cases was specified as 100 in the reference scenario, and as 20 or 500 in scenarios for testing different sample sizes). For individual matched studies, cases were selected from the stratum Y=1 using simple random sampling method at first, and then 1, 2, or 4 controls were matched with each case by L6.
In every case-control study, 11 original ORs were calculated according to the following adjustment strategies.
Insufficient adjustment strategies:

Method S4. Generation of meta-analyses
Meta-analyses were generated by combining several case-control studies (the number of studies was specified as 20 in the reference scenario, and as 5 or 50 in scenarios for testing different meta scales). In the reference scenario, all original studies came from a single target population and shared a same sampling process, so that the net effect of adjustment strategies on meta-analyses could be evaluated without interference from clinical and methodological heterogeneity. In scenarios for testing different pooling methods, original studies were sampled from different sources of populations and were exposed to different sources of biases other than confounding (a systematic error of original OR that obeyed normal distribution was simulated for each case-control study). Therefore, heterogeneity might play a role in the pooled effect estimation, and should be addressed by suitable statistical models.  The heterogeneity statistic is given by Under the null hypothesis that there are no differences in effect estimation among original studies, Q follows a chi-squared distribution with N−1 degrees of freedom. I 2 is calculated as Random-effects model: The individual effect sizes are assumed to have a distribution ~( , 2 ), where the 2 is estimated by Each original effect sizes are weighted by and thus the summary estimate is The confidence interval for ̂ is estimated by ̂± Although we found some variations in effect estimations by different adjustment strategies, the causal interpretation had not been affected in this specific scenario, for the pooled ORs of meta-analyses ranged from 1.69 to 2.82, which could all be interpreted as medium effect. Similar phenomena were also observed when ORAY was 0.2, 0.5, or 5 (Supplementary Table S3). Hence, when the true effect of exposure A on outcome Y was medium to strong, adjustment strategies hardly influenced the casual effect estimation. However, when the true effect was weak to none, inappropriate adjustment strategies might distort the causal interpretation and need special attention (Supplementary Table S3). Besides, if the covariates U = [L, R, M, C] were continuous rather than categorical, the pooled ORs also showed more sensitivity to the adjustment strategies (Supplementary Figure S2).

Result S2. Interpretation of performances of statistical adjustment strategies
Supplementary  Figure S4). The simulated results confirmed that whether adjusted for these variables would not affect the estimation of causal effect.
In above scenarios, we assumed all confounders L = [L1, L2, …, L6] were independent to each other, and full adjustment strategy was recommended for accurate effect estimation. However, if the confounders were correlated (rUU  0), less than 6 of them would be sufficient to control confounding bias, and the stronger the correlation, the fewer the variables needed to be adjusted (Supplementary To estimate the effect of passive smoking on breast cancer, women's characteristics, including demographics (e.g., age, body mass index), socioeconomics (e.g., education, occupation), and lifestyle behaviors (e.g., physical activity, alcohol intake), were adjusted in almost every case-control study. They affected passive smoking through affecting surrounding people's characteristics as well as smoking status, and affected breast cancer directly or through affecting women's reproductive factors (e.g., menarche, parity, menopause) [2]. Since women's characteristics were the ancestor variables of both passive smoking and breast cancer, they could be considered as confounders and adjusted in original studies; yet not all of them had to be adjusted due to collinearity. Reproductive factors were also adjusted in many studies, but they were risk factors of breast cancer that had no causations with passive smoking in our hypothesis. Hence, whether adjusted for reproductive factors did not influence the causal effect estimations. Family history was clearly a cause of breast cancer [2], and it might change the smoking status of surrounding people and in turn reduce the passive smoking exposure of the woman; i.e., it was a confounder. Benign breast disease was also a cause of breast cancer [2], but it might be induced by passive smoking; i.e., it was a mediator. Except for the above-mentioned variables, a few studies also adjusted for radiation therapy (a cause of breast cancer [3,4]), cardiovascular disease (an effect of passive smoking that showed associations with breast cancer [5]), or aspirin/nonsteroidal anti-inflammatory drug use (an indicator of high C-reactive protein that caused breast cancer and cardiovascular disease [6]). Despite the possibility of reversal causality, radiation therapy and aspirin/nonsteroidal anti-inflammatory drug use would not influence the causal pathway from passive smoking to breast cancer.
Cardiovascular disease, however, might be a collider that need to be treated with more caution. For individual matched studies, cases and controls were also matched by time of recruitment and vital status to balance unmeasured confounders. The corresponding DAG is displayed in Supplementary Figure S5. OR, odds ratio; A, exposure; Y, outcome; U, covariate; L, confounder; R, risk factor; M, mediator; C, collider. * The positive probabilities of L1 to L6 were specified as 20%. † The correlations among confounders L1 to L6 were specified as 0 (rUU = 0). ‡ The associations of L1 to L6, M, and C with A were specified as 2 (ORUA = 2). § The associations of L1 to L6, R, M, and C with Y were specified as 2 (ORUY = 2). ORs were calculated by multiple logistic regression models that adjusted for appropriate variables according to the assumption of causal relationships defined in Figure 1.
The simulated distributions of exposure, outcome, and covariates were close to what the pre-specification.     If the confounders were correlated (rUU  0), less than 6 of them would be sufficient to control confounding bias. Except for unadjustment strategy that greatly overestimated the true effect when rUU was strong, collinearity hardly influence the effect estimations under other adjustment strategies.    If some included studies were obtained from different sources, random-effects model presented significantly better performance than fixed-effects model. However, neither the source of original studies nor the model used in meta-analyses affected the accuracy of full adjustment strategy.  Lee and Hamling from data reported in original case-control studies. ‡ "" means overestimation, and "" means underestimation. Among 29 case-control studies included in the meta-analysis, 8 (27.6%) gave reasonable effect estimations and were included in the primary analysis; 12 (41.4%) were underestimated due to not adjusting for negative confounders of family history (2/12), adjusting for mediators of benign breast disease (9/12), or adjusting for colliders of cardiovascular disease (1/12); 9 (31.0%) were overestimated due to not adjusting for positive confounders of age or body mass index. 47 Figure S1. Flow chart of the simulation study OR, odds ratio. A target population was simulated with pre-determined exposure, outcome, and covariates. Then case-control studies were randomly selected from the population, and a series of original ORs were calculated in each case-control study by different adjustment strategies. Then meta-analyses were conducted to pool these ORs. The process was repeated for 1000 times to obtain the empirical distribution of pooled OR.  In the causal path from passive smoking to breast cancer, women's characteristics and family history of breast cancer were considered to be confounders; reproductive factors, radiation therapy, and aspirin/NSAID use were considered to be risk factors; benign breast disease was considered to be a mediator; cardiovascular disease was considered to be a collider. For individual matched studies, cases and controls were also matched by time of recruitment and vital status to balance unmeasured confounders.