A method of assessing the sensitivity of the Cochran–Mantel–Haenszel test to an unobserved confounder

Observational studies, including the case-control design frequently used in epidemiology, are subject to a number of biases and possible confounding factors. Failure to adjust with them may lead to an erroneous conclusion about the existence of a causal relationship between exposure and disease. The Cochran–Mantel–Haenszel (CMH) test is widely used to measure the strength of the association between an exposure and disease or response, after stratifying on the observed covariates. Thus, observed confounders are accounted for in the analysis. In practice, there may be causal variables that are unknown or difficult to obtain. Hence, they are not incorporated into the analysis. Sensitivity analysis enables investigators to assess the robustness of the findings. A method for assessing the sensitivity of the CMH test to an omitted confounder is presented here. The technique is illustrated by re-examining two datasets: one concerns the effect of maternal hypertension as a risk factor for low birth weight infants and the other focuses on the risk of allopurinol on having a rash. The computer code performing the sensitivity analysis is provided in appendix A.


Introduction
Observational studies are commonly used to evaluate the effect of treatments or chemical exposure, e.g. the association between risk factors and a subsequent event. Unlike randomized experiments, researchers cannot control the assignment of treatments to subjects to ensure that subjects receiving different treatments or exposure are similar with respect to major covariates. The Cochran-Mantel-Haenszel (CMH;Cochran 1954;Mantel & Haenszel 1959) test is widely used to test the strength of the association after controlling for the observed confounders. The data are stratified into multiple 2!2 tables by known confounders, where in each stratum the distributions of the confounders in both groups are similar.
Since observational studies remain subject to unobserved confounding, one should explore how sensitive the results are to a possible omitted variable. Extensive work has been done on how to perform sensitivity analysis for observational data by Cornfield et al. (1959), Rosenbaum & Rubin (1983, 1984, Rosenbaum & Krieger (1990), Gastwirth (1992), Rosenbaum (2002) and others. Most of the literature models the association between an unobserved confounder, denoted by U, and treatment assignment, and then recomputes the test statistics and p values for a range of plausible differences in the prevalences of the confounder in the two groups. If a moderate degree of confounding could change the inference about the treatment effect on the outcome, then one questions the soundness of the conclusions. Gastwirth et al. (1998) proposed a dual and simultaneous sensitivity analysis for matched pair data, which extended a strengthened form of the original Cornfield inequality (Gastwirth 1988, p. 296). Both the relationship between U and treatment assignment and that between U and response are jointly modelled. This simultaneous form allows one to incorporate subject matter knowledge about both relationships into the sensitivity analysis.
This paper extends the 'simultaneous' approach to examine the sensitivity of the CMH test to an unobserved confounder. The rest of this paper is organized as follows. Section 2 reviews the CMH test and describes the stratified model for the simultaneous relationships. Then we apply the sensitivity analysis to a dataset examining maternal hypertension as a risk factor for low birth weight (LBW) infants and another dataset concerned with a side effect (rash) of allopurinol in §3. We discuss the potential use and future development of sensitivity analysis in observational studies in §4.

Sensitivity analysis for the CMH test
Let X be the exposure, e.g. environmental hazard, genetic factor, radiation exposure; and Y be the outcome, e.g. disease, mortality, or in an employment discrimination case, being promoted or passing a pre-employment exam. We are interested in testing the causal relationship between X and Y using the CMH test to control for other known factors related to the response. Usually the data are organized into J 2!2 tables, where in each stratum the covariates have similar values. The data for the jth stratum ( jZ1, ., J ) are shown below.
Let p 1j (p 0j ), jZ1, ., J, be the probability that subjects who are exposed (not exposed) in stratum j has an event. In terms of these parameters, the odds ratio between exposure and outcome in stratum j is q j Z p 1j ð1K p 0j Þ=½p 0j ð1K p 1j Þ; j Z 1; .; J : outcome Y exposure group X event (1) non-event (0) total exposed (1) a j b j m 1j not exposed (0) c j d j m 0j total n 1j n 0j N j Furthermore, assume that the odds ratio has a common value q. Let DZ fa j ; b j ; c j ; d j ; j Z 1; .; J g be the observed data for the J 2!2 tables. When the responses for each group in all strata follow independent binomial distributions, the full data have the product binomial likelihood function If all the confounders are observed and controlled by stratification, then the null hypothesis of no-exposure effect is H 0 : p 0j Zp 1j , jZ1, ., J or H 0 : qZ1. Cochran's (1954) statistic for testing H 0 : qZ1 versus H 1 : qO1 is where w j Zm 1j m 0j /N j , p j Z ðm 1j p 1j C m 0jp0j Þ=N j and p 1j Z a j =m 1j , p 0j Z c j =m 0j are the estimates of p 1j and p 0j , respectively. Using the conditional central hypergeometric likelihood rather than the product binomial (2.1), Mantel & Haenszel (1959) proposed a similar test. The two tests are asymptotically equivalent and are often referred to as the CMH test. The CMH test is known to be the optimal method of combining multiple 2!2 tables with common odds ratio (Woolson et al. 1986).
The large sample distribution of (2.2) converges to a standard normal distribution under the null hypothesis that qZ1 (Cochran 1954). More generally, the statistic converges to a standard normal distribution as N j /N (Woolson et al. 1986). Sometimes, information on known risk factors is not available or other risk factors related to the response have not been discovered. Here we assume that there is an unobserved factor U and that the relationship between binary outcome Y and exposure X and confounder U is modelled by where g is the effect of U, which is called the strength parameter (Yu & Gastwirth 2005), and b is the true effect of exposure X. The null hypothesis of no-exposure effect is H 0 : bZ0 and the alternative hypothesis is H 1 : bs0.
The unobserved variable U may also be associated with the exposure. Assuming that U is binary, the prevalence of U by exposure level in each stratum is modelled by Here, the d j , indicating the association between X and U in each stratum, are called the imbalance parameters (Yu & Gastwirth 2005). To simplify the calculation, we assume that the imbalance parameters d j are equal to d. Then, q j Z e l j =ð1C e l j Þ is the prevalence of U in the non-exposed group in stratum j.
When the confounder U is omitted in the analysis, the marginal probability of response at each exposure level x in stratum j becomes The true null hypothesis of no-exposure effect is H 0 : bZ0, which is really the hypothesis of interest. Note that H 0 : p 0j Zp 1j , jZ1, ., J, and H 0 : bZ0 are not equivalent when there is an omitted confounder. If U is not a confounder, i.e. dZ0 or gZ0, then the equalities p 1j Zp 0j , jZ1, ., J, still hold under the null hypothesis H 0 : bZ0, so that the CMH test C 0 is indeed testing the effect of exposure. When U is indeed a confounder, it has been shown that, if both gO0 and dO0, then p 1j Zp 0j , jZ1, ., J, even when the true null hypothesis holds (Yu & Gastwirth 2005) and the CMH test would be biased. This implies that confounder U could induce spurious positive association between exposure and outcome. However, if gO0 and d!0 then the true association between X and Y will be underestimated if the confounder is omitted. If the strength parameter g, the imbalance parameter d and the prevalences q 1 , ., q J of the confounder are known or can be reliably estimated from other sources, the estimates of p 1j and p 0j can be calculated from the likelihood (2.1) under the null hypothesis H 0 : bZ0 (see appendix A). Then the correct CMH statistic is obtained by substituting the estimatesp 1j ;p 0j À Á for p 0j ; p 1j À Á in (2.3) as the effect of U is now incorporated into the response probabilities p xj . In practice these quantities are unknown, so one assumes a possible range for (q 1 , ., q J , g, d); the test statistic and the p value are calculated using (2.3). A table showing the p values with respect to different parameters for confounding can be used to assess the plausibility that the current inference could be altered due to the effect of an unobserved factor.

Application
The proposed sensitivity analysis will be applied to data from two observational studies. In these applications, CZexp(g) measures the effect of U on binary response Y and DZexp(d) measures the increasing odds of U being positive when X increases by one unit.
(a ) Sensitivity analysis of maternal hypertension as a risk factor for LBW infants Data on 500 singleton births in a London Hospital (Hills & De Stavola 2002) were used to examine the sensitivity of the association of maternal hypertension as a risk factor for LBW infants. The dataset provides the birth information of 500 singletons with the following eight variables: identity number for mother and baby, birth weight of baby (indicator for birth weight less than 2500 g), gestation period (indicator for gestation period less than 37 weeks), maternal age, indicator for maternal hypertension and sex of baby (1, male; and 2, female). Ten mothers with missing gestation period information were dropped from the analysis. The exposure of interest (X ) is maternal hypertension and the event of interest is LBW infant (Y ).
'Advanced' maternal age is defined as any expectant mother who will have reached her 35th birthday by delivery date. The mothers are divided into 'young' and 'advanced' age group according to this criterion. The number of infants with normal birth weight (NBW) and LBW by sex, gestation status, maternal age and maternal hypertension are shown in table 1. The CMH test is significant with pZ0.015 and the odds ratio is 2.74 with 95% CI (1.23, 6.14). The Breslow-Day test for homogeneity shows that the odds ratio estimates are homogeneous across different strata ( pZ0.785). Is it safe to draw the conclusion that maternal hypertension is causally related to LBW?
Despite substantial reductions in US infant mortality during the past several decades, black-white disparities in infant mortality rates persist. Important determinants of racial/ethnic differences in infant mortality are LBW, defined as less than 2500 g (US Department of Health and Human Services 2000). The number of LBW infants among 1000 live births for the whites are 5.7-6.5 and for the blacks are 12.7-13.0 from 1980 to 2000 (Centers for Disease Control and Prevention 2002; http://www.cdc.gov/mmwr/preview/mmwrhtml/mm5127a1.htm). Hence, race is an important risk factor with an apparent odds ratio of approximately 2.0-2.4 for LBW (David & Collins 1997). A recent study (Geronimus et al. 2007) also showed that being black increases the odds ratio of hypertension by 2.11-4.04 for women between 15 and 65 years of age. Also, LBW is related to many other risk factors, e.g. the use of assisted reproductive technology, use of alcohol or drugs. Here, we focus on the race as the unobserved confounder (U ) because it is a recognized risk factor. We assume that the mother's race is not related to the sex of their babies, but that the percentage of mothers who are black varies across the maternal age groups as well as gestation term. The possible scenarios for the percentage of black mothers by maternal age and gestation term are shown in table 2. In scenarios 1 and 2, we assume that mothers with preterm babies consist of a higher percentage of black women. In scenarios 3 and 4, we assume that there are more black mothers with normal maternal age. In scenario 5, the black mothers tend to be younger and have more preterm babies than whites.   Table 3 shows the two-sided p values of the adjusted CMH test with respect to different values of C and D for scenarios 1-5. Note that the strength parameter CZexp(g) is the odds ratio of a black mother having a LBW baby compared with a white mother. The imbalance parameter DZexp(g) is the increased odds ratio of a mother being black if she has hypertension. When CZ1 or DZ1, race is not a confounder and the CMH test is unchanged.
Another way to perform the sensitivity analysis is to obtain the threshold values of the sensitivity parameters C and D, which increase the p value of the CMH test to 0.05. Figure 1 shows the plots of sensitivity parameters C and D, which corresponds to pZ0.05. From figure 1, we see that when CZ1.5 and DZ1.1, then the two-sided p value would reach the 0.05 level. Thus, both table 3 and figure 1 show that the p value of the CMH test is sensitive to moderate changes in C and D, indicating that the significant association between maternal hypertension and LBW might have arisen from a hidden bias due to confounding. Hence, other studies, including race and ethnicity information, are needed to reach the conclusion that maternal hypertension causes LBW. If the studies demonstrate that increased incidence of black mothers with hypertension was slight with odds ratio less than 1.5, then the significance of the London study would not be affected by the omission of race.
Although in Rubin's causal model, sex and race are not used as 'treatment' variables in causal inference because it is difficult to ascribe causality to variables that cannot be altered (Holland 1986), here race is used as a proxy for other factors reflecting socio-economic disparities. More recently, race and sex have been used in the propensity score method to form 'similar strata' (Zanutto et al. 2005) but should not be the main potential outcome of a study as they cannot be altered.
(b ) Sensitivity analysis of allopurinol as a cause of rash Rosenbaum (2002) examined the sensitivity of a study indicating that allopurinol can cause rash (Boston Collaborative Drug Project 1972). Allopurinol is commonly used for the treatment and prevention of attacks of gout and certain types of kidney stones. It is also used to treat elevated uric acid levels in the blood and urine, which can occur in patients receiving chemotherapy for the treatment of leukaemia, lymphoma and other types of cancer. Allopurinol is usually well tolerated by most patients; however, some may experience possible side effects of skin rash, hives and itching. Table 4 shows the joint frequencies of drug use and getting a rash case for males and females, respectively. The CMH test is highly significant with p!0.0001. In a randomized trial or in the absence of unobserved confounders, this indicates strong evidence that allopurinol causes rash (Rosenbaum 2002, p. 131). Because the patients could also be taking amoxicillin and ampicillin, which are known to cause rash, is it possible that the rash was caused by an unobserved factor, e.g. use of amoxicillin or other drugs, instead of allopurinol?
Although information about the other drug use and users' allergic history was not collected, one can assess the sensitivity of the CMH test with respect to an unobserved confounder. Let U be the binary indicator for amoxicillin use, e.g. use amoxicillin (1) and no amoxicillin (0). Again, the prevalences of using amoxicillin among the non-allopurinol users are denoted by q j , jZ1, ., J. The parameter d indicates the imbalance of U between allopurinol users and nonusers. The parameter g indicates the increasing probability of having rash among users of amoxicillin. We assume two possible scenarios for the prevalences of amoxicillin use among the non-allopurinol users, i.e. q j Z10% or 20%, where jZ1, 2.
The increasing odds ratio of having rash for amoxicillin users compared with non-users is not known, but it is likely to be below 4. Based on table 5, we can see that the CMH test of allopurinol causing rash is not sensitive to the omitted variable (use of amoxicillin). Even when the strength parameter CZ4 and the imbalance parameter DZ4, the CMH test remains highly significant with pZ0.002. This example also illustrates the additional information provided by the simultaneous sensitivity analysis compared with the earlier approaches by Cornfield et al. (1959) and Rosenbaum (2002), which focus primarily on the imbalance parameter and assume a confounder with infinitely large strength. Rosenbaum (2002, p. 131) used the parameter G, which is similar to D in our analysis, to measure the imbalance of treatment assignment with respect to the unobserved factor. The upper bounds on the p values based on his analysis are 0.036 and 0.30, for GZ2 and 3, respectively, which indicates that the causal relationship between allopurinol use and rash could have arisen from an unobserved confounder that triples the risk of rash. While such a confounding is unlikely to exist, the conclusion based on the simultaneous analysis is stronger. Even when CZN, the test barely reached non-significance for DZ3 and q j Z10%. This sensitivity analysis is more stringent, indicating the conclusion that allopurinol causes rash is very robust to an unobserved confounder. This is consistent with the current medical practice that lists rash as the one common side effect of allopurinol use (http://www.nlm.nih.gov/medlineplus/druginfo/ medmaster/a682673.html).

Discussion
Because the CMH test is widely used to analyse the data obtained in epidemiological studies, legal cases and social sciences, it is important to conduct a sensitivity analysis to answer whether that inference could plausibly be explained by an unobserved confounder before drawing a firm conclusion. The proposed simultaneous sensitivity analysis approach is similar to the primal sensitivity analysis method by Rosenbaum (2002). The major difference is the introduction of the strength parameter. The primal sensitivity analyses closely parallel the theory of randomized experiments and have been developed for many statistical tests and estimators, e.g. the McNemar and Wilcoxon signed-rank tests and the Hodges-Lehmann estimator (Rosenbaum 2002). The simultaneous approach allows one to use more subject matter knowledge. However, it also requires more assumptions about both relationships. The simultaneous approach is especially useful when an observational study is challenged owing to failure to control for a particular unobserved variable, when there are plausible ranges for the values of g and d (Gastwirth et al. 1998). At the design stage, it is also useful to examine the size and power of the corrected CMH test in the presence of omitted confounders in a range of realistic values of the strength and imbalance parameters. The calculation of size and power for the Cochran-Armitage trend test has been derived by Yu & Gastwirth (2005). A similar technique can be used for the CMH test with omitted variable.
In future studies, we plan to extend this methodology of sensitivity analysis to data where the unobserved confounders are continuous and ordinal. In real-world observational studies, there may be more than one omitted confounder; background information on the joint distribution of the multiple confounders can then be used to replace the univariate distribution (2.5) in a sensitivity analysis. B.Y. was supported in part by the Intramural Research Programme of the National Institute on Aging and J.L.G. was supported in part by grant SES-0317956 from the National Science Foundation.
Appendix A. Estimation of (p 1j ,p 0j ) This appendix describes the estimation of (p 1j ,p 0j ) under the product-binomial likelihood (2.1) with fixed values of (q 1 , ., q J , b, g, d) for the allopurinol data in §3b. For stratum j, the prevalence of U in the non-exposed group is q j Zexp(l j )/[1Cexp(l j )]. Then the prevalence of U in the exposed group is r j Z q j e d =ð1K q j C q j e d Þ. Let A j Z expða j Þ; BZ expðbÞ; C Z expðgÞ. Then for fixed values of (b, d, g, q j ), the marginal probabilities of response are The log-likelihood of function (2.1) is log LðqjDÞZ P J jZ1 [ s ðqjDÞ, where [ j ðqjDÞ Z a j log p 1j C b j log ð1K p 1j Þ C c j log p 0j C d j log ð1K p 0j Þ; for stratum s. The maximum-likelihood estimate (MLE) of a j , and ðp 1j ; p 0j Þ could be obtained from SAS PROC NLP. The code for the following 2!2!2 table from the allopurinol study with q 1 Zq 2 Z0.2, bZ0, gZlog (2) and dZlog (2) is shown below. The MLEs of ðp 1j ; p 0j Þ are given in PIEST.