Everything You Wanted to Know but Could Never Find from the Cochran-Mantel-Haenszel Test

José Moral de la Rubia¹, Adrián Valle de la O²

¹School of Psychology, Universidad Autónoma de Nuevo León, Nuevo León, Monterrey, Mexico.
²School of Medicine and Health Science, Tecnológico de Monterrey, Monterrey, Mexico.

DOI: 10.4236/jdaip.2023.113016                

Abstract

The Cochran-Mantel-Haenszel (CMH) test, developed in the 1950s, is a classic in health research, especially in epidemiology and other fields in which dichotomous and polytomous variables are frequent. This nonparametric test makes it possible to measure and check the effect of an antecedent variable X on a health outcome Y, statistically controlling the effect of a third variable Z that acts as a confounding variable in the relationship between X and Y. Both X and Y are measured on a dichotomous qualitative scale and Z on a polytomous-qualitative or ordinal scale. It is assumed that the effect of X on Y is homogeneous between the k strata of Z, which is usually tested by the Breslow-Day test with the Tarone’s correction or the Woolf’s test. The main statistical programs have the CMH test together with a test to verify the assumption of a homogeneous effect across the strata, so that it is easy to apply. However, its fundamentals and details of calculations are a mystery to most researchers, and even difficult to find or understand. The aim of this article is to present these details in a clear and concise way, including the assumptions and alternatives to non-compliance. This technical knowledge is applied to a simulated realistic example of the area of epidemiology in health and, finally, an interpretive synthesis of the analyses is given. In addition, some suggestions for the test report are made.

Keywords

Odds Ratio, Effect Size, Statistical Control, Qualitative Variables, Nonparametric Statistics

1. Introduction

The Cochran-Mantel-Haenszel (CMH) test [1] [2] is a statistical test used to examine the association between two categorical variables while controlling for the effects of a third variable. It is particularly useful when analyzing data with a 2 × 2 contingency table and a potential confounding variable. By calculating a test statistic that follows a chi-squared distribution, the CMH test determines whether the association between the variables of interest is statistically significant after accounting for the stratifying variable. This test, named after its developers, Cochran, Mantel, and Haenszel, is commonly employed in epidemiology and biostatistics to understand relationships in observational studies and clinical trials while adjusting for potential confounders.

This manuscript, focused on statistical aspects of data analysis, aims to present the technical details of the Cochran [1] and Mantel and Haenszel [2] test in order to achieve an adequate understanding of it. First, the hypotheses and assumptions of the test are stated. Next, the point estimation and confidence interval of the common odds ratio and the check of its significance with a Z test are shown; afterward, the computation of the Cochran’s Q statistic and its sampling distribution is examined, as well as the statistical decision making. The issue of effect size for this nonparametric test is also addressed. The verification of the assumption of a homogeneous effect of X on Y in each stratum is deepened through two tests, namely: the Woolf’s test [3] and the Breslow-Day test [4] with Tarone’s correction [5] . The theoretical part ends by indicating how to proceed in case of unfulfillment of the aforementioned homogeneity assumption and what other, more generalized, alternative exists to the CMH test. Finally, an example of application of this test is presented. The effect of tobacco smoke, smoked at home, on bronchial asthma in pediatric minors is analyzed, statistically controlling the effect of air quality in the place of residence. For this purpose, a simulated sample composed of 2000 participants, which is composed of k sub-samples randomly collected from the general population of k countries, was generated. This example is a realistic simulation aimed at showing technical and interpretive aspects of the CMH test.

Due to the current emphasis on practical applications of statistical software in textbooks and teaching in empirical sciences, the mathematical foundations of statistical tests are often overlooked. However, curiosity often arises, especially among educators and researchers, regarding how these seemingly magical results generated by software are obtained, and there is an underlying question of whether these quick and easy results are indeed valid. The motivation of this article is to address this need by providing a comprehensible and practical approach to the CMH test.

With the CMH test, as well as the verification of its assumptions and alternatives, all calculations can be performed using statistical software. The article does not cover the need to know how to conduct the calculations semi-automatically, but focuses on unraveling the mystery of the mathematical foundation of the calculations, a mystery because the information is lost when it is systematically overlooked.

2. Statistical Hypothesis of the CMH Test

This non-parametric test, initially devised by Cochran [1] and deepened by Mantel and Haenszel [2] , was developed to test the relationship between two dichotomous qualitative variables, X and Y, while statistically controlling the effect of a third variable Z with k qualitative or ordinal categories; these categories receive the name of strata. Variable X is a factor that affects health; it is evaluated using two categories: 0 = absent or not exposed to X, and 1 = present or exposed to X. On the other hand, variable Y is the clinical outcome or diagnostic status (case), and it is also evaluated using two categories: 0 = non-case and 1 = case [1] [2] . The information of these three variables is arranged in k 2 × 2 tables of a k × 4 table. The null hypothesis holds that variable X has no effect on variable Y, which implies that the odds ratio (OR) of exposure between cases and non-cases is unitary. The alternative hypothesis states that X has an effect on Y, so that the odds ratio is not unitary. If OR is lower than 1, then X acts as a protective factor for being a case; on the other hand, if OR is higher than 1, then X acts as a risk factor for being a case [6] .

The odds of exposure for an event A are the ratio between the probability of being exposed and the probability of not being exposed: Odds(A) = P(A)/[1 − P(A)]. These oddsare calculated both between cases and non-cases. The odds ratio (OR) is the quotient between the odds of exposure in cases (numerator) and the odds of exposure in non-cases (denominator): OR = Odds(A|case)/Odds(A|non-case). An OR higher than 1 informs how many times the presence of the risk factor is more likely than its absence among cases compared to non-cases. On the other hand, an OR lower than 1 informs how many times is more likely the presence of the protective factor than its absence among non-cases compared to cases. A unitary OR shows that the exposure factor has no effect on being a case or not.

Null hypothesis (H₀): OR = 1

Alternative hypothesis (H₁): OR ≠ 1

3. CMH Test Assumptions

The test requires two dichotomous qualitative variables X and Y. The variable X is supposed to be a risk factor or a protection factor for health, and is measured or evaluated through two categories: exposed and unexposed. Variable Y is a health outcome (be it a disorder, syndrome, or disease), also evaluated through two categories (case vs non-case) assumed to be affected by or dependent on X. There is also a polytomous variable Z, with k qualitative or ordinal categories, that acts as a confounding variable in the relationship between X and Y. Thus, k independent random samples of paired data of X and Y, drawn from the strata or population groups defined by the variable Z, are required. Furthermore, assuming homogeneity across the ORs of the k strata, the effect of X on Y is estimated by the common odds ratio across the k strata.

In order to apply the CMH test, the raw OR of each table, the common OR (the adjusted OR for the stratification variable Z), and the Cochran’s Q test statistic, which is a generalization of the McNemar’s test statistic, are calculated for k 2 × 2 tables [2] . Symbology is expressed in Table 1. Finally, it is necessary to verify the homogeneity of effect between the k strata or equivalence of the OR of the k 2 × 2 tables. The Woolf’s test and the Breslow-Day test with the Tarone’s correction are the most commonly used statistics to verify this assumption [6] [7] .

4. Odds Ratio across Strata

From Table 1, the odds ratio between X and Y in stratum i of the confounding variable Z, OR_i, is calculated using the following probability quotient:

${\widehat{OR}}_i=\frac{Odds\left(X=1|Y=1\right)}{Odds\left(X=1|Y=0\right)}=\frac{\frac{P\left(X=1|Y=1\right)}{1-P\left(X=1|Y=1\right)}}{\frac{P\left(X=1|Y=0\right)}{1-P\left(X=1|Y=0\right)}}=\frac{\frac{a_i}{m_{1i}}/\frac{c_i}{m_{1i}}}{\frac{b_i}{m_{0i}}/\frac{d_i}{m_{0i}}}=\frac{a_i/c_i}{b_i/d_i}=\frac{a_i{d}_i}{b_i{c}_i}$ (1)

If the random sample drawn from a population stratum is large, the significance of each OR_i can be tested and a confidence interval can be estimated from its distributional convergence to a standard normal distribution. Nevertheless, in order to facilitate this convergence, a natural logarithmic transformation is applied [8] .

Statistical hypotheses:

${\text{H}}_0:O{R}_i=1$

${\text{H}}_1:O{R}_i\ne 1$

Test statistic and its asymptotic sampling distribution (Table 1):

$\ln \left({\widehat{OR}}_i\right)=\ln \left(\frac{a_i{d}_i}{b_i{c}_i}\right)$ (2)

${\hat{\sigma }}_{\ln \left(O{R}_i\right)}=\sqrt{\frac{1}{a_i}+\frac{1}{b_i}+\frac{1}{c_i}+\frac{1}{d_i}}$ (3)

$Z=\frac{\ln \left({\widehat{OR}}_i\right)}{{\hat{\sigma }}_{\ln \left(O{R}_i\right)}}=z~N\left(0,1\right)$ (4)

Decision based on the critical value (z_1−α/2) or the critical level (p-value) in a two-tailed test:

$\text{If}z\le z_{1-\frac{\alpha }{2}}\text{or}p\text{-value}=2\times \left[1-P\left(Z\le \left|z\right|\right)\right]\ge \alpha ,{\text{H}}_0\text{isaccepted}$

$\text{If}z>z_{1-\frac{\alpha }{2}}\text{or}p\text{-value}=2\times \left[1-P\left(Z\le \left|z\right|\right)\right]<\alpha ,{\text{H}}_0\text{is}\text{rejected}$

z_1−α/2 = critical value or quantile of order 1 − α/2 of a standard normal distribution N(0, 1). If α = 0.05 (conventional significance level), z_0.975 = 1.96.

P(Z ≤ z) = F_Z(z) = value of the cumulative distribution function or cumulative probability of a standard normal distribution.

Table 1. Joint and marginal frequencies of X and Y in group or stratum i of Z.

Note. X ={0 = not exposed, 1 = exposed} = factor that supposedly generates an effect on health, Y ={0 = no case, 1 = case} = health status or outcome, Z = {1, 2, …, k} = confounding variable in the relationship between X and Y, which defines the k strata. Joint absolute frequencies in stratum i: n_11i = a_i = exposed cases, n_10i = b_i = exposed non-cases, n_01i = c_i = non-exposed cases, and n_00i = d_i = non-exposed non-cases. Absolute marginal frequencies in stratum i: n_1•i = n_1i = exposed, n_0•i = n_0i = non-exposed, n_•1i = m_1i = cases, and n_•0i = m_0i = non-cases. Total stratum of i: n_••i = T_i.

Confidence interval at (1 − α) × 100 of OR_i:

$P\left(\ln \left({\widehat{OR}}_i\right)-z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_i\right)}\le \ln \left(O{R}_i\right)\le \ln \left({\widehat{OR}}_i\right)+z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_i\right)}\right)=1-\alpha$ (5)

$P\left({\text{e}}^{\ln \left({\widehat{OR}}_i\right)-z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_i\right)}}\le O{R}_i\le {\text{e}}^{\ln \left({\widehat{OR}}_i\right)+z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_i\right)}}\right)=1-\alpha$ (6)

$1\in \left[{\text{e}}^{\ln \left({\widehat{OR}}_i\right)-z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_i\right)}},{\text{e}}^{\ln \left({\widehat{OR}}_i\right)+z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_i\right)}}\right]\Rightarrow H_0:O{R}_i=1$

$1\notin \left[{\text{e}}^{\ln \left({\widehat{OR}}_i\right)-z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_i\right)}},{\text{e}}^{\ln \left({\widehat{OR}}_i\right)+z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_i\right)}}\right]\Rightarrow H_1:O{R}_i\ne 1$

5. Odds Ratio across Strata

The Mantel-Haenszel common odds ratio [2] can be calculated through the following ratio of sums of proportions (Table 1):

${\widehat{OR}}_{MH}=\frac{{\displaystyle {\sum }_{i=1}^k\frac{a_i{d}_i}{T_i}}}{{\displaystyle {\sum }_{i=1}^k\frac{b_i{c}_i}{T_i}}}$ (7)

The significance of the Mantel-Haenszel common odds ratio can be tested using a Z test, as the SPSS program does [9] .

Statistical hypotheses. The null hypothesis holds that population OR_MH is equal to 1, that is, variable X has no effect on variable Y across the homogenous k strata. The alternative hypothesis states that population OR_MH is not equal to 1 and, consequently, variable X has effect on variable Y.

${\text{H}}_0:O{R}_{MH}=1$

${\text{H}}_1:O{R}_{MH}\ne 1$

Assumptions: 1) X is a dichotomous qualitative variable that plays the role of predictor or independent variable; 2) Y, also a dichotomous qualitative variable, serves the role of dependent variable; 3) there are k independent strata defined by the qualitative or ordinal categories of a confounding variable Z; 4) there exists a random sample of paired data from X and Y in each stratum; and 5) the odds ratios are homogeneous across the k strata.

The Z-test statistic and its asymptotic sampling distribution: A natural logarithmic transformation is applied to the Mantel-Haenszel common odds ratio.

$\ln \left({\widehat{OR}}_{MH}\right)$ (8)

The asymptotic standard error of the log-transformed Mantel-Haenszel common odds ratio is estimated through the formula given by Robins, Breslow, and Greenland [10] .

${\hat{\sigma }}_{\ln \left(O{R}_{MH}\right)}=\sqrt{\frac{1}{2}\left(\frac{{\displaystyle {\sum }_{i=1}^k{s}_i{u}_i}}{u^2}+\frac{{\displaystyle {\sum }_{i=1}^k\left(s_i{v}_i+t_i{u}_i\right)}}{uv}+\frac{{\displaystyle {\sum }_{i=1}^k{t}_i{v}_i}}{v^2}\right)}$ (9)

$s_i=\frac{a_i+d_i}{T_i}$ $t_i=\frac{b_i+c_i}{T_i}$ $u_i=\frac{a_i\times d_i}{T_i}$ $v_i=\frac{b_i\times c_i}{T_i}$

${\displaystyle \underset{i=1}{\overset{k}{\sum }}{s}_i{u}_i}={\displaystyle \underset{i=1}{\overset{k}{\sum }}\left(\frac{a_i+d_i}{T_i}\times \frac{a_i\times d_i}{T_i}\right)}$

${\displaystyle \underset{i=1}{\overset{k}{\sum }}{t}_i{v}_i}={\displaystyle \underset{i=1}{\overset{k}{\sum }}\left(\frac{b_i+c_i}{T_i}\times \frac{b_i\times c_i}{T_i}\right)}$

${\displaystyle \underset{i=1}{\overset{k}{\sum }}\left(s_i{v}_i+t_i{u}_i\right)}={\displaystyle \underset{i=1}{\overset{k}{\sum }}\left(\frac{a_i+d_i}{T_i}\times \frac{b_i\times c_i}{T_i}+\frac{b_i+c_i}{T_i}\times \frac{a_i\times d_i}{T_i}\right)}$

$u={\displaystyle \underset{i=1}{\overset{k}{\sum }}{u}_i}={\displaystyle \underset{i=1}{\overset{k}{\sum }}\frac{a_i\times d_i}{T_i}}$

$v={\displaystyle \underset{i=1}{\overset{k}{\sum }}{v}_i}={\displaystyle \underset{i=1}{\overset{k}{\sum }}\frac{b_i\times c_i}{T_i}}$

${\widehat{OR}}_{MH}=u/v$

The quotient between the log-transformed Mantel-Haenszel common odds ratio (numerator) and its standard deviation or error (denominator) constitutes the Z test statistic that converges, in distribution, to a standard normal distribution.

$Z=\frac{\ln \left({\widehat{OR}}_{MH}\right)}{{\hat{\sigma }}_{\ln \left(O{R}_{MH}\right)}}=z~N\left(0,1\right)$ (10)

Decision and confidence interval. The decision on H₀ can be taken based on the critical value (z_1−α/2) or the critical level (p-value) in a two-tailed test.

$\text{If}z\le z_{1-\frac{\alpha }{2}}\text{or}p\text{-value}=2\times \left[1-P\left(Z\le \left|z\right|\right)\right]\ge \alpha ,{\text{H}}_{\text{0}}\text{is}\text{accepted}$

$\text{If}z>z_{1-\frac{\alpha }{2}}\text{or}p\text{-value}=2\times \left[1-P\left(Z\le \left|z\right|\right)\right]<\alpha ,{\text{H}}_0\text{is}\text{rejected}$

For calculating the confidence interval at (1 − α) × 100 for the Mantel-Haenszel common odds ratio, first, the lower and upper limits of the log-transformed OR_MH are obtained, using the Wald’s method [11] ; next, the transformation is undone in order to arrive to the Mantel-Haenszel common odds ratio limits.

$\begin{array}{l}P\left(\ln \left({\widehat{OR}}_{MH}\right)-z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_{MH}\right)}\le \ln \left(O{R}_{MH}\right)\le \ln \left({\widehat{OR}}_{MH}\right)+z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_{MH}\right)}\right)\\ =1-\alpha \end{array}$ (11)

$P\left({\text{e}}^{\ln \left({\widehat{OR}}_{MH}\right)-z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_{MH}\right)}}\le O{R}_{MH}\le {\text{e}}^{\ln \left({\widehat{OR}}_{MH}\right)+z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_{MH}\right)}}\right)=1-\alpha$ (12)

$1\in \left[{\text{e}}^{\ln \left({\widehat{OR}}_{MH}\right)-z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_{MH}\right)}},{\text{e}}^{\ln \left({\widehat{OR}}_{MH}\right)+z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_{MH}\right)}}\right]\Rightarrow O{R}_{HM}=1$

$1\notin \left[{\text{e}}^{\ln \left({\widehat{OR}}_{MH}\right)-z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_{MH}\right)}},{\text{e}}^{\ln \left({\widehat{OR}}_{MH}\right)+z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_{MH}\right)}}\right]\Rightarrow O{R}_{HM}\ne 1$

6. CMH Q Statistic

The Q statistic is a quotient. Its numerator is the square of the difference between the sum of exposed cases in each stratum and the sum of the expected values under the assumption of independence for these frequencies. Its denominator is the sum of the variances of the frequencies in each stratum under the assumption of independence. The higher the numerator is with respect to the denominator, the higher the relationship between the two dichotomous variables [2] [6] .

(13)

The Q test statistic converges to a chi-square distribution with 1 degree of freedom; guarantee of an adequate convergence includes: frequencies expected under the null hypothesis of independence ( $E_{ij}=n_i\times m_j/T_i$ ) higher than or equal to 5, absence of any null observed joint frequency, and totals (T_i) higher than 25 in each stratum. In order to improve the convergence to a chi-square distribution with one degree of freedom, it is recommended to apply the Yates’ correction; the reason for implementing this controversial correction is the need to smooth the transition from a discrete binomial distribution to a continuous chi-square distribution [12] .

${Q^{\prime }}_{CMH}=\frac{{\left(\left|{\displaystyle {\sum }_{i=1}^k{a}_i}-{\displaystyle {\sum }_{i=1}^k\frac{n_{1i}{m}_{1i}}{T_i}}\right|-0.5\right)}^2}{{\displaystyle {\sum }_{i=1}^k\frac{n_{0i}{n}_{1i}{m}_{0i}{m}_{1i}}{T_i^3-T_i^2}}}$ (14)

$Q_{CMH}\text{and}{Q^{\prime }}_{CMH}~{\chi }_1^2$

Decision based on the critical value or critical level in a right-tailed test:

$\text{If}{Q}_{CMH}\le {}_{1-\alpha }\chi {}_1^2\text{or}P\left({\chi }_1^2\ge Q_{CMH}\right)\ge \alpha ,{\text{H}}_0:OR=1\text{is}\text{accepted}$

$\text{If}{Q}_{CMH}>{}_{1-\alpha }\chi {}_1^2\text{or}P\left({\chi }_1^2\ge Q_{CMH}\right)<\alpha ,{\text{H}}_0:OR=1\text{is}\text{rejected}$

${}_{1-\alpha /2}\chi {}_1^2$ = critical value or quantile of order 1 − (α/2) of a chi-square distribution with one degree of freedom.

$P\left({\chi }_1^2\ge Q_{CMH}\right)$ = probability value to the right tail of a chi-square distribution with one degree of freedom.

It should be noted that Mantel and Haenszel [2] , had previously presented a formula to define the confidence interval of the common odds ratio, but it is less efficient and precise, especially with small samples [10] .

$P\left[{\text{e}}^{\left(1-\frac{z_{1-\frac{\alpha }{2}}}{\sqrt{{Q^{\prime }}_{CMH}}}\right)\ln \left({\widehat{OR}}_{MH}\right)}\le O{R}_{MH}\le {\text{e}}^{\left(1+\frac{z_{1-\frac{\alpha }{2}}}{\sqrt{{Q^{\prime }}_{CMH}}}\right)\ln \left({\widehat{OR}}_{MH}\right)}\right]=1-\alpha$ (15)

7. Effect Size

The magnitude of the effect size for the CMH test is the common OR that estimates the effect of X on Y and, for interpretation, it can be transformed into a Cohen’s d with Chinn’s formula [13] . Values of d lower than 0.2 reflect an effect size that is trivial, from 0.2 to 0.49 is small, from 0.5 to 0.79 is large, and higher than or equal to 0.8 is very large (Cohen, 1988).

(16)

8. Tests of Homogeneity of Odds Ratios

In order to avoid overextending this manuscript, only the two most widely used tests for testing the assumption of homogeneity between ORs will be presented, namely: Woolf’s (1957) test and Breslow-Day test (1980) with Tarone’s (1985) correction. Nevertheless, it should be noted that there exist more tests, such as those developed by Bliss; Zelen; Yusuf, Peto, Lewis, Collins, and Sleight; Liang and Self; and Kulinskaya and Dollinger [14] [15] [16] . These last two are especially recommended for small samples [17] .

8.1. Woolf’s Test

Woolf’s test is one of the first tests developed for assessing homogeneity of ORs. Woolf’s method is based on a statistic first introduced by Cochran [1] and it aims to test the assumption of homogeneity of ORs for k independent strata from 2 × 2 tables (effect of a dichotomous variable X on a dichotomous variable Y). The calculation of Woolf’s statistic is simple and its small sample performance is better than that of other tests [7] ; it is available in programs such as Real Statistics Resource Pack for Excel [18] and R software package for statistical analysis [19] .

${\text{H}}_0:O{R}_1=O{R}_2=\cdots =O{R}_k$

${\text{H}}_1:\exists O{R}_i\ne O{R}_j;i\ne j=1,2,\cdots ,k$

Assumptions: 1) a random sample composed of, at least, 25 participants in each stratum; 2) two dichotomous qualitative variables: one is the independent or antecedent variable X and the other is the dependent or consequent variable Y; and 3) a polytomous variable Z that generates the k independent strata in the 2 × 2 tables and that acts as a confounding variable in the relationship between X and Y [20] .

Test statistic:

$Q_{Woolf}=\underset{i=1}{\overset{k}{{\displaystyle \sum }}}\frac{{\left(\ln {\widehat{OR}}_i-\ln {\widehat{OR}}_{comun}\right)}^2}{{\hat{\sigma }}_{\ln \left(O{R}_i\right)}^2}$ (17)

$\ln \left({\widehat{OR}}_i\right)=\ln \left(\frac{a_i{d}_i}{b_i{c}_i}\right)$ (18)

${\hat{\sigma }}_{\ln \left(O{R}_i\right)}^2=\frac{1}{a_i}+\frac{1}{b_i}+\frac{1}{c_i}+\frac{1}{d_i}$ (19)

${\widehat{OR}}_{common}={\text{e}}^{\frac{{\displaystyle {\sum }_{i=1}^k\left(\ln \left(O{R}_i\right)/{\hat{\sigma }}_{\ln \left(O{R}_i\right)}^2\right)}}{{\displaystyle {\sum }_{i=1}^k\left(1/{\hat{\sigma }}_{\ln \left(O{R}_i\right)}^2\right)}}}$ (20)

Sampling distribution and decision based on the critical value or critical level in a right-tailed test:

$Q_{Woolf}~{\chi }_{k-1}^2$

$\text{If}{Q}_{Woolf}\le {}_{1-\alpha }\chi {}_{k-1}^2\text{or}P\left({\chi }_{k-1}^2\ge Q\right)\ge \alpha ,{\text{H}}_{\text{0}}\text{is}\text{accepted}$

8.2. Breslow-Day Test Using Tarone’s Bias Correction

The Breslow-Day test requires large samples in each stratum, and is not adequate to use it with small sample sizes (n < 20). Its test statistic is consistent, but it presents a slight bias, and that is why Tarone [5] developed a correction. Precisely, Breslow [21] recommends the use of this correction, although he points out that, when the sample size is large, the resulting change in the estimate is minimal and only perceptible beyond the third or fourth decimal place of the probability value. The Breslow-Day test using Tarone’s correction is the most powerful when comparing various OR homogeneity tests, all of which are not very powerful with small strata sizes [17] . It is worth noting that the Breslow-Day test using Tarone’s correction is the most frequently OR homogeneity test included in several statistical software packages available in the market and, consequently, it is the most reported OR homogeneity test in research articles [7] .

The omnibus formulation of hypotheses and assumptions of the Breslow-Day test, with or without the Tarone’s correction, are the same as in the Woolf’s test; nevertheless, its asymptotic approximation to the chi-square distribution with k − 1 degrees of freedom is better when Tarone’s correction is used.

Breslow-Day test statistic is a goodness-of-fit test. The expected absolute frequencies are calculated under the assumption of a common OR estimated by the formula of Mantel and Haenszel [2] .

$\begin{array}{c}{Q}_{BD}=\underset{i=1}{\overset{k}{{\displaystyle \sum }}}\frac{{\left(a_i-E\left(a_i|O{R}_{MH}\right)\right)}^2}{E\left(a_i|O{R}_{MH}\right)}+\frac{{\left(b_i-E\left(b_i|O{R}_{MH}\right)\right)}^2}{E\left(b_i|O{R}_{MH}\right)}\\ +\frac{{\left(c_i-E\left(c_i|O{R}_{MH}\right)\right)}^2}{E\left(c_i|O{R}_{MH}\right)}+\frac{{\left(d_i-E\left(d_i|O{R}_{MH}\right)\right)}^2}{E\left(d_i|O{R}_{MH}\right)}\end{array}$ (21)

a_i = n_11i = absolute (joint) frequency of exposed cases in stratum j.

b_i = n_10i = absolute (joint) frequency of exposed non-cases in stratum j.

c_i = n_01i = absolute (joint) frequency of non-exposed cases in stratum j.

d_i = n_00i = absolute (joint) frequency of non-exposed non-cases in stratum j.

OR_MH = common odds ratio of k strata estimated through Mantel-Haenszel method

E(a_i|OR_MH) = expected value for the number of exposed cases in stratum j conditioned on OR_MH. This value is obtained upon solving the following quadratic equation:

$\begin{array}{l}\left(1-{\widehat{OR}}_{MH}\right)\times E^2\left(a_i|O{R}_{MH}\right)+\left[n_{0i}-m_{1i}+{\widehat{OR}}_{MH}\left(n_{1i}+m_{1i}\right)\right]\\ \times E\left(a_i|O{R}_{MH}\right)-m_{1i}{n}_{1i}{\widehat{OR}}_{MH}=0\end{array}$ (22)

$\alpha =1-{\widehat{OR}}_{MH}$

$\gamma =-m_{1i}{n}_{1i}{\widehat{OR}}_{MH}$

$E\left(a_i|O{R}_{MH}\right)=\frac{-\beta \mp \sqrt{{\beta }^2-4\alpha \gamma }}{2\alpha }$ (23)

Every quadratic equation has two solutions. In this case, the solution in the set of real numbers with a value higher than 0 and lower than or equal to the minimum of m_1i and n_1i is chosen. The expected values are not rounded to integers, but are handled with their decimal places.

$0<E\left(a_i|O{R}_{MH}\right)\le \min \left(m_{1i},n_{1i}\right)$

Once the expected frequency of the exposed cases, E(a_i|OR_MH), is estimated, the other expected frequencies can be obtained under the condition of a common odds ratio (OR_MH).

(24)

$E\left(c_i|O{R}_{MH}\right)=m_{1i}-E\left(a_i|O{R}_{MH}\right)$ (25)

$E\left(d_i|O{R}_{MH}\right)=T_i-m_{1i}-n_{1i}+E\left(a_i|O{R}_{MH}\right)$ (26)

$T_i=E\left(a_i|O{R}_{MH}\right)+E\left(b_i|O{R}_{MH}\right)+E\left(c_i|O{R}_{MH}\right)+E\left(d_i|O{R}_{MH}\right)$

Next, taking into account the expected values E(a_i|OR_MH) of each stratum, the Mantel-Haenszel common odds ratio (1959) can be calculated, and this allows checking that there are no errors in the calculation of the expected frequencies, as the OR_MH value is already known.

(27)

The squares of the differences between the observed and expected frequency within each stratum are equal.

$\begin{array}{c}{\left(a_i-E\left(a_i|O{R}_{MH}\right)\right)}^2={\left(b_i-E\left(b_i|O{R}_{MH}\right)\right)}^2={\left(c_i-E\left(c_i|O{R}_{MH}\right)\right)}^2\\ ={\left(d_i-E\left(d_i|O{R}_{MH}\right)\right)}^2\end{array}$

Likewise, the variances of the frequencies conditioned to a common odds ratio OR_MH coincide within each stratum, and they are obtained through the following formula:

$\begin{array}{l}{s}^2\left(a_i|O{R}_{MH}\right)=\frac{1}{\frac{1}{E\left(a_i|O{R}_{MH}\right)}+\frac{1}{E\left(b_i|O{R}_{MH}\right)}+\frac{1}{E\left(c_i|O{R}_{MH}\right)}+\frac{1}{E\left(d_i|O{R}_{MH}\right)}}\\ =\frac{1}{\frac{1}{E\left(a_i|O{R}_{MH}\right)}+\frac{1}{n_{1i}-E\left(a_i|O{R}_{MH}\right)}+\frac{1}{m_{1i}-E\left(a_i|O{R}_{MH}\right)}+\frac{1}{T_i-m_{1i}-n_{1i}+E\left(a_i|O{R}_{MH}\right)}}\end{array}$ (28)

Thus, the Breslow-Day statistic can be simplified to the following formula with k addends:

$Q_{BD}={\displaystyle \underset{i=1}{\overset{k}{\sum }}\frac{{\left(a_i-E\left(a_i|O{R}_{MH}\right)\right)}^2}{s^2\left(a_i|O{R}_{MH}\right)}}$ (29)

Tarone’s correction for the Breslow-Day statistic involves subtracting from the Q_BD a quotient whose numerator is the square of the difference between the sum of the exposed cases and the sum of their expected values conditioned to the Mantel-Haenszel common odds ratio, whereas its denominator is the sum of the variances of the exposed cases conditioned to the Mantel-Haenszel common odds ratio; hence, once the data for the Q_BD calculation are available, its correction is very simple.

$\begin{array}{c}{Q}_{BDT}=Q_{BD}-\frac{{\left(a-E\left(a|O{R}_{MH}\right)\right)}^2}{s^2\left(a|O{R}_{MH}\right)}\\ =\underset{i=1}{\overset{k}{{\displaystyle \sum }}}\frac{{\left(a_i-E\left(a_i|O{R}_{MH}\right)\right)}^2}{s^2\left(a_i|O{R}_{MH}\right)}-\frac{{\left(a-E\left(a|O{R}_{MH}\right)\right)}^2}{s^2\left(a|O{R}_{MH}\right)}\end{array}$ (30)

$a={\displaystyle \underset{i=1}{\overset{k}{\sum }}{a}_i}$

$E\left(a|O{R}_{MH}\right)={\displaystyle \underset{i=1}{\overset{k}{\sum }}E}\left(a_i|O{R}_{MH}\right)$ (31)

$s^2\left(a|O{R}_{MH}\right)={\displaystyle \underset{i=1}{\overset{k}{\sum }}{s}^2}\left(a_i|O{R}_{MH}\right)$ (32)

Sampling distribution and decision based on the critical value or critical level in a right-tailed test:

$Q_{BDT}~{\chi }_{k-1}^2$

$\text{If}{Q}_{BDT}\le {}_{1-\alpha }\chi {}_{k-1}^2\text{or}P\left({\chi }_{k-1}^2\ge Q\right)\ge \alpha ,{\text{H}}_{\text{0}}\text{is}\text{accepted}$

$\text{If}{Q}_{BDT}>{}_{1-\alpha }\chi {}_{k-1}^2\text{or}P\left({\chi }_{k-1}^2\ge Q\right)<\alpha ,{\text{H}}_{\text{0}}\text{is}\text{rejected}$

9. Unfulfillment of the Assumption of Homogeneous Effect in the k Strata

In the event that the null hypothesis of equivalence of the odds ratio between the k strata is rejected, it is necessary to identify which strata are homogeneous and which strata are heterogeneous in order that the CMH test can be applied to the homogenous ones, and the effect size can be reported either by strata groups or individual strata.

Pairwise comparisons can be made with a Z test [7] with the Holm-Bonferroni or Holm-Sidak correction in order to control the increase in type-1 error or with the Benjamin-Hochberg correction in order to control the false discovery rate [22] .

Statistical hypotheses:

$H_0:O{R}_i=O{R}_j;i\ne j=1,2,\cdots ,k$

$H_1:O{R}_i\ne O{R}_j$

Test statistic:

$Z=\frac{\ln \left({\widehat{OR}}_i\right)-\ln \left({\widehat{OR}}_j\right)}{\sqrt{{\hat{\sigma }}_{\ln \left(O{R}_i\right)}^2+{\hat{\sigma }}_{\ln \left(O{R}_j\right)}^2}}=\frac{\ln \left(\frac{a_i{d}_i}{b_i{c}_i}\right)-\ln \left(\frac{a_j{d}_j}{b_j{c}_j}\right)}{\sqrt{\frac{1}{a_i}+\frac{1}{b_i}+\frac{1}{c_i}+\frac{1}{d_i}}\times \sqrt{\frac{1}{a_j}+\frac{1}{b_j}+\frac{1}{c_j}+\frac{1}{d_j}}}~N\left(0,1\right)$ (33)

Decision based on the critical value or critical level in a two-tailed test:

$\text{If}z\le z_{1-\frac{\alpha }{2}}\text{or}2\times \left(1-P\left(Z\le \left|z\right|\right)\right)\ge \alpha ,{\text{H}}_{\text{0}}\text{is}\text{accepted}$

$\text{If}z>z_{1-\frac{\alpha }{2}}\text{or}2\times \left(1-P\left(Z\le \left|z\right|\right)\right)<\alpha ,{\text{H}}_{\text{0}}\text{is}\text{rejected}$

Given k independent groups, the number of pairwise comparisons corresponds to combinatorics without repetition: k of 2, C(k, 2), which means m = k(k − 1)/2 comparisons. To make the Holm-Bonferroni correction, the probabilities are ordered ascendingly: an order value (i) equal to 1 is assigned to the smallest probability, an order value (i) equal to 2 is given to the second smallest probability, and so on until assigning a value equal to m to the highest probability; in the event that two or more probability values are repeated, an average order value is assigned (sum of the corresponding orders divided by the number of added orders). Next, the significance level α (usually 0.05) is divided by m + 1 − i, thus obtaining the corrected significance level: ${\alpha }_{HB}=\alpha /\left(m+\text{1}-i\right)$ . In the Holm-Sidak correction, the corrected significance level is: ${\alpha }_{HS}=1-\left(1-\alpha \right)1/\left(m+1-i\right)$ . The Benjamini-Hochberg corrected significance level is: ${\alpha }_{BH}=\left(i/m\right)\times \alpha$ [22] . See an example with four independent groups (Table 2).

$m=\left(\begin{array}{c}k\\ 2\end{array}\right)=\left(\begin{array}{c}4\\ 2\end{array}\right)=\frac{k\left(k-1\right)}{2}=\frac{4\times 3}{2}=\frac{12}{2}=6$

$p_1<p_2<p_3<p_4<p_5<p_6$

$\text{If}{p}_i<{\alpha }_c,{\text{H}}_0\text{is}\text{rejected}$

A more general test than the CMH test is conditional logistic regression. This multivariate technique can include more than one antecedent or predictor variable, whether qualitative (dichotomous, polytomous, ordinal) or quantitative (discrete or continuous). When conditional logistic regression handles a single dichotomous predictor, apart from the confounding variable that defines k strata, the score test statistic corresponds to the CMH Q statistic [20] . Being θ the estimated regression parameter or weight, the score test statistic, denoted by S(θ), is the quotient between the derivative in θ of the natural logarithm of the likelihood function ln(θ|x) and the Fisher information about the parameter θ. The statistic, in this one-parameter situation, follows a chi-square distribution with one degree of freedom, as the CMH Q statistic does.

$S\left(\theta \right)=\frac{\frac{\delta }{\delta \theta }\ln \mathcal{L}\left(\theta |x\right)}{I\left(\theta \right)}=\frac{\frac{\delta }{\delta \theta }\ln \left(\mathcal{L}\left(\theta |x\right)\right)}{-E\left[\frac{{\delta }^2}{\delta {\theta }^2}\ln \left(f_X\left(x|\theta \right)\right)\right]}=Q_{CMH}~{\chi }_1^2$

Conditional logistic regression was developed in 1978 by Breslow, Day, Halvorsen, Prentice, and Sabai [23] . When data for the predicted dichotomous variable and predictors are recorded by stratum, this technique has a clear

Table 2. Holm-Bonferroni, Holm-Sidak, and Benjamin-Hochberg corrections for a significance level α at 0.05 and four compared groups.

advantage over unconditional logistic regression, especially when there are many small strata sizes, such as the extreme case of paired data or n/2 strata of size 2. Logistic regression includes the confounding variable that generates the k strata as one more predictor, and estimate location parameters for each of its k strata. Nevertheless, conditional logistic regression does not require calculating these parameters and achieves more efficiency in its estimates and power, which becomes more evident as the number of strata increases [4] .

Conditional logistic regression is available in the R statistical software using the “clogit” function in the survival package. It is located in this package because the log-likelihood of a conditional logistic model is the same as the log-likelihood of a Cox model with a paired or stratified data structure [24] .

10. Simulated Example of the Application of the CMH Test

In order to see a more realistic example, a feasible study in epidemiology is presented. A sample composed of 2000 Mexican males from 0 to 14 years old (pediatric population) was collected and classified into two groups: 1 = case of bronchial asthma and 0 = no case. The diagnosis of unspecified uncomplicated bronchial asthma (CA23.32) was made by pediatricians from public health centers using the criteria of the eleventh edition of the International Classification of Diseases (ICD-11) of the World Health Organization [25] . The exposure factor was established based on indoor smoking of tobacco at home by parents and other adults who reside with the minor: 1 = daily exposure to tobacco smoke at home and 0 = no exposure to household tobacco smoke. Tobacco is understood as any product processed from the leaves of the plant known as Nicotianatabacum, such as cigarettes, cigars, or pipe tobacco. Of the 2000 participants, 595 lived in an urban location with good air quality most days of the year, 1007 in a location with usually acceptable air quality, and 398 with usually poor air quality, according to the official Mexican standard: NOM-172-SEMARNAT-2019 on Guidelines to Obtain and Communicate the Air Quality and Health Risks Index [26] . The objective of the study is to verify the relation between daily exposure to tobacco smoke at home and the child’s bronchial asthma by controlling the confounding variable (air quality in the urban area of residence) in males aged 0 to 14 years using the CMH test [2] [6] with a significance level at 5%. In case of positive association, the effect size with the Mantel-Haenszelcommon odds ratio [2] is estimated. Also, the assumption of homogeneous odds ratio across the three strata is tested by the Woolf’s test [3] and Breslow-Day test with the Tarone’s correction [5] with a significance level at 5%.

Table 3 shows the three 2 × 2 contingency tables stratified by the confounding variable (air quality in the urban area of residence) of the pediatric participants (0 to 14 years of age). The three strata are good, fair, and poor air quality most days of the year. The two crossed dichotomous variables are: X = daily exposure to tobacco smoke at home = {0 = no, 1 = yes} and Y = bronchial asthma = {0 = no case, 1 = case}.

Table 3. Joint absolute and marginal frequencies of daily exposure to tobacco smoke at home and whether or not the participant is a case of bronchial asthma, stratifying by the confounding variable (air quality in the urban area of residence) of the participants.

Note. X = daily exposure to tobacco smoke at home = {0 = no, 1 = yes}, Y = bronchial asthma = {0 = no case, 1 = case}, Z = air quality in the urban area of residence = {1 = good, 2 = fair, 3 = poor}, ∑_column = sum per column in each stratum, ∑_row = sum per row in each stratum_i.

The first step is to calculate the odds ratio of each stratum both in its point estimate, using Equation (1), and 95% interval confidences, using Equations(2), (3), (5), and (6). In order to facilitate these calculations, Table 4 has been created. The significance of each OR is calculated with the Z-test (Equation 4), and the effect size in each stratum with the OR transformed into a Cohen’s d with Chinn’s formula [13] that appears in the Equation (16).

Let Z be the confounding variable, that is, air quality in the urban area of residence of the participant: Z = {1 = good, 2 = fair, 3 poor}.

The OR for strata 1 (good air quality) and 2 (fair quality) are significant with a significance level at 5%, since they do not include the value 1 when they are estimated with a confidence level at 95%, and the two-tailed probability-values of their Z-test statistics are lower than the α significance level of 0.05. They reflect a small effect size (0.2 ≤ d < 0.5) and a detrimental or health risk effect (OR > 1). The OR of stratum 3 (poor air quality) is not significant at significance level equal to 0.05 and shows a trivial effect size (d < 0.2). Nevertheless, the 95% confidence intervals of the three OR overlap, so that they seem to be equivalent from this first analysis, and the significant, small, and detrimental effect that daily exposure to tobacco smoke at home has on bronchial asthma in pediatric male patients remains unclear.

${\widehat{OR}}_1=\frac{a_1{d}_1}{b_1{c}_1}=\frac{80\times 121}{381\times 13}=\frac{9680}{4953}=1.9544$

$\ln \left({\widehat{OR}}_1\right)=\ln \left(1.9544\right)=0.6701$

${\hat{\sigma }}_{\ln \left(O{R}_1\right)}=\sqrt{\frac{1}{a_1}+\frac{1}{b_1}+\frac{1}{c_1}+\frac{1}{d_1}}=\sqrt{\frac{1}{80}+\frac{1}{381}+\frac{1}{13}+\frac{1}{121}}=0.3167$

$Z=\frac{\ln \left({\widehat{OR}}_1\right)}{{\hat{\sigma }}_{\ln \left(O{R}_1\right)}}=\frac{\ln \left(1.9544\right)}{0.3167}=\frac{0.6701}{0.3167}=2.1156$

$Sig=2\times \left(1-P\left(Z\le 2.1156\right)\right)=0.0344<\alpha =0.05$

Table 4. Joint frequencies between X and Y, and total for each air quality stratum.

Note. a_i = n_11i = exposed cases in each stratum i, b_i = n_10i = exposed non-cases in each stratum i, c_i = n_01i = non-exposed cases in each stratum i, d_i = n_00i = non-exposed non-cases in each stratum i, T_i = n_••i = total participants in each stratum or sum per row (i from 1 to 3), ∑_i = sum by column i from 1 to 3.

$P\left({\text{e}}^{\ln \left({\widehat{OR}}_1\right)-z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_1\right)}}\le O{R}_1\le {\text{e}}^{\ln \left({\widehat{OR}}_1\right)+z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_1\right)}}\right)=1-\alpha$

$P\left({\text{e}}^{\ln \left(1.9544\right)-1.96\times 0.3167}\le O{R}_1\le {\text{e}}^{\ln \left(1.9544\right)+1.96\times 0.3167}\right)=0.95$

$P\left(O{R}_1\in \left[1.0505,3.6358\right]\right)=0.95$

$1\notin \left[1.0505,3.6358\right]\Rightarrow O{R}_1\ne 1,O{R}_1>1$

$d_1=\frac{\ln \left({\widehat{OR}}_1\right)}{\pi /\sqrt{3}}=\frac{\ln \left(1.9544\right)}{\pi /\sqrt{3}}=0.3694$

${\widehat{OR}}_2=\frac{a_2{d}_2}{b_2{c}_2}=\frac{145\times 194}{642\times 26}=\frac{28130}{16692}=1.6852$

$\ln \left({\widehat{OR}}_2\right)=\ln \left(1.6852\right)=0.5219$

${\hat{\sigma }}_{\ln \left(O{R}_2\right)}=\sqrt{\frac{1}{a_2}+\frac{1}{b_2}+\frac{1}{c_2}+\frac{1}{d_2}}=\sqrt{\frac{1}{145}+\frac{1}{642}+\frac{1}{26}+\frac{1}{194}}=0.2282$

$Z=\frac{\ln \left({\widehat{OR}}_2\right)}{{\hat{\sigma }}_{\ln \left(O{R}_2\right)}}=\frac{\ln \left(1.6852\right)}{0.2282}=\frac{0.5219}{0.2282}=2.2872$

$Sig=2\times \left(1-P\left(Z\le 2.2872\right)\right)=0.0222<\alpha =0.05$

$P\left({\text{e}}^{\ln \left({\widehat{OR}}_2\right)-z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_2\right)}}\le O{R}_2\le {\text{e}}^{\ln \left({\widehat{OR}}_2\right)+z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_2\right)}}\right)=1-\alpha$

$P\left({\text{e}}^{\ln \left(1.6852\right)-1.96\times 0.2282}\le O{R}_2\le {\text{e}}^{\ln \left(1.6852\right)+1.96\times 0.2282}\right)=0.95$

$P\left(O{R}_2\in \left[1.0775,2.6357\right]\right)=0.95$

$1\notin \left[1.0775,2.6357\right]\Rightarrow O{R}_2\ne 1,O{R}_2>1$

$d_2=\frac{\ln \left({\widehat{OR}}_2\right)}{\pi /\sqrt{3}}=\frac{\ln \left(1.6852\right)}{\pi /\sqrt{3}}=0.2877$

${\widehat{OR}}_3=\frac{a_3{d}_3}{b_3{c}_3}=\frac{62\times 102}{210\times 24}=\frac{6324}{5040}=1.2548$

$\ln \left({\widehat{OR}}_3\right)=\ln \left(1.2548\right)=0.2269$

${\hat{\sigma }}_{\ln \left(O{R}_3\right)}=\sqrt{\frac{1}{a_3}+\frac{1}{b_3}+\frac{1}{c_3}+\frac{1}{d_3}}=\sqrt{\frac{1}{62}+\frac{1}{210}+\frac{1}{24}+\frac{1}{102}}=0.2690$

$Z=\frac{\ln \left({\widehat{OR}}_3\right)}{{\hat{\sigma }}_{\ln \left(O{R}_3\right)}}=\frac{\ln \left(1.2548\right)}{0.2690}=\frac{0.2269}{0.2690}=0.8437$

$Sig=2\times \left(1-P\left(Z\le 0.8437\right)\right)=0.3989>\alpha =0.05$

$P\left({\text{e}}^{\ln \left({\widehat{OR}}_3\right)-z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_3\right)}}\le O{R}_3\le {\text{e}}^{\ln \left({\widehat{OR}}_3\right)+z_{1-\frac{\alpha }{2}}{\hat{\sigma }}_{\ln \left(O{R}_3\right)}}\right)=1-\alpha$

$P\left({\text{e}}^{\ln \left(1.2548\right)-1.96\times 0.2690}\le O{R}_3\le {\text{e}}^{\ln \left(1.2548\right)+1.96\times 0.2690}\right)=0.95$

$P\left(O{R}_3\in \left[0.7406,2.1259\right]\right)=0.95$

$1\in \left[0.7406,2.1259\right]\Rightarrow O{R}_3=1$

$d_3=\frac{\ln \left({\widehat{OR}}_3\right)}{\pi /\sqrt{3}}=\frac{\ln \left(1.2548\right)}{\pi /\sqrt{3}}=0.1251$

The second step is to calculate the common odds ratio for the three strata using the Mantel-Haenszel formula [2] . Firstly, Mantel-Haenszel point estimate is calculated using Equation (7).

${\widehat{OR}}_{MH}=\frac{{\displaystyle {\sum }_{i=1}^k\frac{a_i{d}_i}{T_i}}}{{\displaystyle {\sum }_{i=1}^k\frac{b_i{c}_i}{T_i}}}=\frac{\frac{80\times 121}{595}+\frac{145\times 194}{1007}+\frac{62\times 102}{398}}{\frac{381\times 13}{595}+\frac{642\times 26}{1007}+\frac{210\times 24}{398}}=\frac{60.0928}{37.5637}=1.5998$

Secondly, the 95% confidence interval for the Mantel-Haenszel common OR (OR_MH) is calculated by applying a logarithmic transformation upon OR_MH, as shown in the Equation (8); next, the asymptotic standard deviation of the natural logarithm of OR_MH is computed [10] , as in Equation (9); then, the lower and upper limits for the log-transformed OR_MH are obtained using the Wald’s method [11] , as in Equation (11); and lastly, the logarithmic transformation is undone in order to obtain the limits for the OR_MH, as in the Equation (12).The effect size, which turns out to be small, is calculated with Equation (16), which transforms an OR into a Cohen’s d.

$\ln \left(O{R}_{MH}\right)=\ln \left(1.5998\right)=0.4699$

$u={\displaystyle \underset{i=1}{\overset{k}{\sum }}{u}_i}={\displaystyle \underset{i=1}{\overset{k}{\sum }}\frac{a_i\times d_i}{T_i}}=\frac{80\times 121}{595}+\frac{145\times 194}{1007}+\frac{62\times 102}{398}=60.0928$

$\begin{array}{c}{\displaystyle \underset{i=1}{\overset{k}{\sum }}{s}_i{u}_i}={\displaystyle \underset{i=1}{\overset{k}{\sum }}\left(\frac{a_i+d_i}{T_i}\times \frac{a_i\times d_i}{T_i}\right)}\\ =\frac{80+121}{595}\times \frac{80\times 121}{595}+\frac{145+194}{1007}\times \frac{145\times 194}{1007}+\frac{62+102}{398}\times \frac{62\times 102}{398}\\ =21.4472\end{array}$

$v={\displaystyle \underset{i=1}{\overset{k}{\sum }}{v}_i}={\displaystyle \underset{i=1}{\overset{k}{\sum }}\frac{b_i\times c_i}{T_i}}=\frac{381\times 13}{595}+\frac{642\times 26}{1007}+\frac{210\times 24}{398}=37.5637$

$\begin{array}{c}{\displaystyle \underset{i=1}{\overset{k}{\sum }}{t}_i{v}_i}={\displaystyle \underset{i=1}{\overset{k}{\sum }}\left(\frac{b_i+c_i}{T_i}\times \frac{b_i\times c_i}{T_i}\right)}\\ =\frac{381+13}{595}\times \frac{381\times 13}{595}+\frac{642+26}{1007}\times \frac{642\times 26}{1007}+\frac{210+24}{398}\times \frac{210\times 24}{398}\\ =23.9533\end{array}$

$\begin{array}{l}{\displaystyle \underset{i=1}{\overset{k}{\sum }}\left(s_i{v}_i+t_i{u}_i\right)}={\displaystyle \underset{i=1}{\overset{k}{\sum }}{s}_i{v}_i}+{\displaystyle \underset{i=1}{\overset{k}{\sum }}{t}_i{u}_i}\\ =\frac{80+121}{595}\times \frac{381\times 13}{595}+\frac{145+194}{1007}\times \frac{642\times 26}{1007}+\frac{62+102}{398}\times \frac{210\times 42}{398}\\ +\frac{381+13}{595}\times \frac{80\times 121}{595}+\frac{642+26}{1007}\times \frac{145\times 194}{1007}+\frac{210+24}{398}\times \frac{62\times 102}{398}\\ =13.6103+38.6456=52.2559\end{array}$

$\begin{array}{c}{\hat{\sigma }}_{\ln \left(O{R}_{MH}\right)}=\sqrt{\frac{1}{2}\left(\frac{{\displaystyle {\sum }_{i=1}^k{s}_i{u}_i}}{u^2}+\frac{{\displaystyle {\sum }_{i=1}^k\left(s_i{u}_i+t_i{v}_i\right)}}{uv}+\frac{{\displaystyle {\sum }_{i=1}^k{t}_i{v}_i}}{v^2}\right)}\\ =\sqrt{\frac{1}{2}\left(\frac{21.4472}{{60.0928}^2}+\frac{52.2559}{60.0928\times 37.5637}+\frac{23.9533}{{37.5637}^2}\right)}\\ =\sqrt{\frac{0.0461}{2}}=0.1518\end{array}$

$P\left({\text{e}}^{\ln \left({\widehat{OR}}_{MH}\right)-z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_{MH}\right)}}\le O{R}_{MH}\le {\text{e}}^{\ln \left({\widehat{OR}}_{MH}\right)+z_{1-\frac{\alpha }{2}}\times {\hat{\sigma }}_{\ln \left(O{R}_{MH}\right)}}\right)=1-\alpha$

$P\left({\text{e}}^{\ln \left(1.5998\right)-1.96\times 0.1518}\le O{R}_{MH}\le {\text{e}}^{\ln \left(1.5998\right)+1.96\times 0.1518}\right)=0.95$

$P\left(O{R}_{MH}\in \left[1.1882,2.1540\right]\right)=0.95$

$1\notin \left[1.1882,2.1540\right]\Rightarrow O{R}_{MH}\ne 1,O{R}_{MH}>1$

$d=\frac{\ln \left({\widehat{OR}}_{MH}\right)}{\pi /\sqrt{3}}=\frac{\ln \left(1.5998\right)}{\pi /\sqrt{3}}=0.2590$

Using the logarithmic transformation of the OR_MH and its standard deviation, the significance of the effect of X on Y can be tested with a Z test, just as it was done with the odds ratio in each one of the three strata [10] . See Equation (10). In a two-tailed test with a significance level of 5%, the Mantel-Haenszel common OR is different from 1. Thus, exposure acts as a risk factor for being a case.

$Z=\frac{\ln \left({\widehat{OR}}_{MH}\right)}{{\hat{\sigma }}_{MH}}=\frac{\ln \left(1.5998\right)}{0.1518}=\frac{0.4699}{0.1518}=3.0959>z_{0.975}=1.96$

$2\times \left(1-P\left(Z\le 3.0959\right)\right)=0.0020<\alpha =0.05$

Thirdly, the CMH Q statistic is computed using Equation (13). Table 5 shows the frequencies of the exposed cases (a_i), the marginal frequencies

Table 5. Frequency of exposed cases, their expected values and their variances assuming independence between X and Y.

Note. Strata or categories of the confounding variable Z = air quality in the area of residence = {1 = good, 2 = fair, 3 = bad}, X = antecedent variable: daily exposure to tobacco smoke at home = {0 = no, 1 = yes}, Y = outcome variable: case of bronchial asthma = {0 = no, 1 = yes}, a_i = n_11i = exposed cases in stratum i, (cases of bronchial asthma that were exposed to tobacco smoke at home in stratum i), n_1i = n_1•i = total of participants (with and without bronchial asthma) exposed to tobacco smoke at home in stratum i, n_0i = n_0•i = total of participants (with and without bronchial asthma) non-exposed to tobacco in stratum i, m_1i = n_•1i = total cases of bronchial asthma in stratum i, m_0i = m_•0i = total of non-cases (participants without bronchial asthma) in stratum i, T_i = n_••i = total number of participants in stratum I, e_i = (n_1i × m_1i)/T_i = expected number of a_i under the assumption of independence between X and Y, $s_i^2$ = variance of a_i under the assumption of independence between X and Y.

(n_1i, n_0i, m_1i, m_0i) and the total (T_i) in each stratum. With these data, the expected value, e_i, and the variance of the exposed cases in each stratum, $s_i^2$ , are calculated under an assumption of independence between X and Y.

The null hypothesis regarding of unitary OR is rejected with the Q statistic without the Yates’ continuity correction. The latter statistic is calculated using Equation (14).

$\begin{array}{c}{Q}_{CMH}=\frac{{\left({\displaystyle {\sum }_{i=1}^k{a}_i}-{\displaystyle {\sum }_{i=1}^k{e}_i}\right)}^2}{{\displaystyle {\sum }_{i=1}^k{s}_i^2}}=\frac{{\left({\displaystyle {\sum }_{i=1}^k{a}_i}-{\displaystyle {\sum }_{i=1}^k\frac{n_{1i}{m}_{1i}}{T_i}}\right)}^2}{{\displaystyle {\sum }_{i=1}^k\frac{n_{0i}{n}_{1i}{m}_{0i}{m}_{1i}}{T_i^3-T_i^2}}}\\ =\frac{{\left(287-264.4708\right)}^2}{52.6001}=\frac{507.5630}{52.6001}=9.6495>{}_{0.95}\chi {}_1^2=3.8415\end{array}$

$P\left({\chi }_1^2\ge 9.6495\right)=0.0019<\alpha =0.05$

Owing to the use of an asymptotic approximation of a discrete distribution to a chi-square distribution with one degree of freedom, it is advisable to use the Yates’ correction. Nevertheless, considering that the sample size is large, the difference between the statistics without or with continuity correction is minimal, Q_CMH = 9.6495, $P\left({\chi }_1^2\ge 9.6495\right)=0.0019$ and ${Q^{\prime }}_{CMH}=9.2259$ , , respectively.

$\begin{array}{c}{Q^{\prime }}_{CMH}=\frac{{\left(\left|{\displaystyle {\sum }_{i=1}^k{a}_i}-{\displaystyle {\sum }_{i=1}^k\frac{n_{1i}{m}_{1i}}{T_i}}\right|-0.5\right)}^2}{{\displaystyle {\sum }_{i=1}^k\frac{n_{0i}{n}_{1i}{m}_{0i}{m}_{1i}}{T_i^3-T_i^2}}}\\ =\frac{{\left(\left|287-264.4708\right|-0.5\right)}^2}{52.6001}=\frac{485.2838}{52.6001}=9.2259>{}_{0.95}\chi {}_1^2=3.8415\end{array}$

$P\left({\chi }_1^2\ge 9.2259\right)=0.0024<\alpha =0.05$

Likewise, because the sample size is large, when the confidence interval of the common odds ratio is calculated with the Mantel-Haenszel formula [2] , as shown in the Equation (15), 95% CI [1.1814, 2.1663], the result is very similar to that obtained with the formula given by Robins et al. [10] , as shown in the Equation (12), 95% CI [1.1882, 2.1540]. The latter is more efficient and accurate; however, this difference becomes evident especially when the sample is small.

$P\left({\text{e}}^{\left(1-\frac{z_{1-\frac{\alpha }{2}}}{\sqrt{{\widehat{Q^{\prime }}}_{CMH}}}\right)\ln \left({\widehat{OR}}_{MH}\right)}\le O{R}_{MH}\le {\text{e}}^{\left(1+\frac{z_{1-\frac{\alpha }{2}}}{\sqrt{{\widehat{Q^{\prime }}}_{CMH}}}\right)\ln \left({\widehat{OR}}_{MH}\right)}\right)=1-\alpha$

$P\left({\text{e}}^{\left(1-\frac{1.96}{\sqrt{9.2259}}\right)\ln \left(1.5998\right)}\le O{R}_{MH}\le {\text{e}}^{\left(1+\frac{1.96}{\sqrt{9.2259}}\right)\ln \left(1.5998\right)}\right)=0.95$

$P\left(1.1814\le O{R}_{MH}\le 2.1663\right)=0.95$

Next, the assumption of homogeneity of the OR across the three strata is tested using the Woolf’s test [3] , as shown inin Equation (17). The required calculations for this test are presented in Table 6, which includes the odds ratios (OR_i), the logarithmic transformations of these odds ratios (ln(OR_i)) using Equation (18), and the variances of the logarithms of the odds ratio for each stratum (s²[ln(OR_i)]) using Equation (19). The natural base exponential of the quotient between the sum of the three quantities in the fifth column of Table 6 and the sum of the three quantities in the sixth column of the same table yields the common odds ratio, calculated with Equation (20).

${\widehat{OR}}_{common}={\text{e}}^{\frac{{\displaystyle {\sum }_{i=1}^k\left(\ln \left(O{R}_i\right)/{\hat{\sigma }}_{\ln \left(O{R}_i\right)}^2\right)}}{{\displaystyle {\sum }_{i=1}^k\left(1/{\hat{\sigma }}_{\ln \left(O{R}_i\right)}^2\right)}}}={\text{e}}^{\frac{19.8392}{42.9932}}={\text{e}}^{0.4615}=1.5864$

Table 6. Odds ratio by stratum and calculations for Woolf’s test.

Note. Strata or categories of the confounding variable Z = AQ = air quality = {G = good, F = fair, P = poor}, OR_i = odds ratio of stratum i, ln(OR_i) = natural logarithm of the odds ratio of stratum i, s²[ln(OR_i)] = variance of the natural logarithm of the odds ratio of stratum i, ln(OR_i)/s²[ln(OR_i)] = quotient between the natural logarithm of the odds ratio of stratum i and its variance, 1/s²[ln(OR_i)] = inverse of the variance of the natural logarithm of the odds ratio of stratum i, [ln(OR_i) − ln(OR)]²/s²[ln(OR_i)] = quotient between the square of the difference between the natural logarithm of the odds ratio of stratum i and the natural logarithm of the common odds ratio (numerator) and the variance of the natural logarithm of the odds ratio of stratum i (denominator), ∑_i = sum per column (i from G toP).

The sum of the three quantities in the sixth column of Table 6 yields the Woolf’s Q test statistic. Each of these three quantities is the quotient between the square of the difference between the natural logarithm of the odds ratio of each stratum and the natural logarithm of the common odds ratio (numerator) and the variance of the natural logarithm of the odds ratio of each stratum (denominator), as shown in Equation (17).

$Q_{Woolf}=\underset{i=1}{\overset{k}{{\displaystyle \sum }}}\frac{{\left(\ln {\widehat{OR}}_i-\ln {\widehat{OR}}_{common}\right)}^2}{{\hat{\sigma }}_{\ln \left(O{R}_i\right)}^2}=0.4339+0.0702+0.7600=1.2640$

The value of the test statistic Q is small and lower than the quantile of order 0.95 of a chi-square distribution with two degrees of freedom (k − 1 = 3 − 1 = 2); thus, the null hypothesis of equality of effect across the three strata is hold.

$Q_{Woolf}=1.2640<{}_{1-\alpha }\chi {}_{k-1}^2={}_{0.95}\chi {}_2^2=5.991$

$P\left({\chi }_2^2\ge Q_{Woolf}=1.2640\right)=0.5315>\alpha =0.05$

Finally, the assumption of homogeneity of the effect across the three strata is tested by Breslow-Day test with Tarone’s correction [5] , as shown in Equation (21). As in the case of the Woolf’s test, $Q_{Woolf}=1.2640<{}_{0.95}\chi {}_2^2=5.991$ , $P\left({\chi }_2^2\ge 1.2640\right)=0.5315>\alpha =0.05$ , the null hypothesis of equivalent odds ratio across the three strata is fulfilled. Because the samples per stratum are large, the difference of the Breslow-Day test in regard to Woolf’s test is very small, Q_BD = 1.270705, $P\left({\chi }_2^2\ge 1.270705\right)=0.529749$ , and it becomes minimal upon applying the Tarone’s correction, Q_BDT = 1.270382, $P\left({\chi }_2^2\ge 1.270382\right)=0.529834$ .

The computation of the Breslow-Day test statistic starts with the calculation of the expected value, E(a₁|OR_MH), using Equation (23), and the variance of the number of exposed cases in relation to the Mantel-Haenszel common OR in the first stratum, s²(a₁|OR_MH), using Equation (28). In order to obtain the expected value E(a₁|OR_MH), it is necessary to pose a quadratic equation that should be solved, which is shown in Equation (22). The joint and marginal frequencies observed for stratum 1 (good air quality) are reproduced in Table 7 in order to facilitate the follow-up of the calculations.

The expected value of a₁ (exposed cases from the first stratum) conditioned to the common OR estimated by the Mantel and Haenszel formula [2] , upon having solved the quadratic equation and chosen a result within the range of 0 to 93 (minimum value for n₁₁ and m₁₁), is equal to 77.8904. It can be verified that this expected frequency is close to the observed one, which is equal to 80.

$\begin{array}{l}\left(1-{\widehat{OR}}_{MH}\right)\times E^2\left(a_1|O{R}_{MH}\right)+\left[n_{01}-m_{11}+{\widehat{OR}}_{MH}\left(n_{11}+m_{11}\right)\right]\\ \times E\left(a_1|O{R}_{MH}\right)-m_{11}{n}_{11}{\widehat{OR}}_{MH}=0\end{array}$

$\begin{array}{l}\left(1-1.5998\right)\times E^2\left(a_1|O{R}_{MH}\right)+\left[134-93+1.5998\times \left(461+93\right)\right]\\ \times E\left(a_1|O{R}_{MH}\right)-93\times 461\times 1.5998=0\end{array}$

Table 7.Joint and marginal frequencies of stratum 1 (good air quality).

Note. X = exposure to indoor tobacco smoke at home, Y = bronchial asthma in pediatric minors. ∑_j = sum per rows (j = 0 and 1), ∑_i sum per columns (i = 0 and 1).

$-0.5998\times E^2\left(a_1|O{R}_{MH}\right)+927.2668\times E\left(a_1|O{R}_{MH}\right)-68586.4895=0$

$\alpha =-0.5998$

$\beta =927.2668$

$\gamma =-68586.4895$

$\begin{array}{c}E\left(a_1|O{R}_{MH}\right)=\frac{-\beta \mp \sqrt{{\beta }^2-4\alpha \gamma }}{2\alpha }\\ =\frac{-927.2668\mp \sqrt{{927.2668}^2-4\times \left(-0.5998\right)\times \left(-68586.4895\right)}}{2\times \left(-0.5998\right)}\\ =\{\begin{array}{l}1468.1739\\ 77.8904\end{array}\end{array}$

$0<E\left(a_1|O{R}_{MH}\right)=77.8904\le \min \left(m_{11},n_{11}\right)=\min \left(93,461\right)=93$

With E(a₁|OR_MH), the expected values for the other three cells of the 2 × 2 table of the first stratum are computed. The expected frequencies are close to the observed ones, namely: b₁ = 381, c₁ = 13, and d₁ = 121. These calculations correspond to Equations (24), (25), and (26).

$E\left(c_1|O{R}_{MH}\right)=m_{11}-E\left(a_1|O{R}_{MH}\right)=93-77.8904=15.10960$

$\begin{array}{c}E\left(d_1|O{R}_{MH}\right)=T_1-m_{11}-n_{11}+E\left(a_1|O{R}_{MH}\right)\\ =595-461-93-77.8904=118.8904\end{array}$

The sum of the four expected frequencies is equal to the total of 595 participants in the first stratum, and it can be verified that the expected frequency for a₁ allows to correctly obtain the OR_MH. See Equation (27).

$\begin{array}{c}{T}_1=E\left(a_1|O{R}_{MH}\right)+E\left(b_1|O{R}_{MH}\right)+E\left(c_1|O{R}_{MH}\right)+E\left(d_1|O{R}_{MH}\right)\\ =77.8904+383.1096+15.10960+118.8904=595\end{array}$

$\begin{array}{c}{\widehat{OR}}_{MH}=\frac{E\left(a_1|O{R}_{MH}\right)\left(n_{01}-m_{11}+E\left(a_1|O{R}_{MH}\right)\right)}{\left(m_{11}-E\left(a_1|O{R}_{MH}\right)\right)\left(n_{11}-E\left(a_1|O{R}_{MH}\right)\right)}\\ =\frac{77.8904\times \left(134-94+77.8904\right)}{\left(94-77.8904\right)\left(461-77.8904\right)}=1.5998\end{array}$

The variance of the exposed cases, conditioned to the common OR estimated by the formula of Mantel and Haenszel formula [2] in the first stratum, is equal to 11.1058. See Equation (28).

$\begin{array}{c}{s}^2\left(a_1|O{R}_{MH}\right)=\frac{1}{\frac{1}{E\left(a_1|O{R}_{MH}\right)}+\frac{1}{E\left(b_1|O{R}_{MH}\right)}+\frac{1}{E\left(c_1|O{R}_{MH}\right)}+\frac{1}{E\left(d_1|O{R}_{MH}\right)}}\\ =\frac{1}{\frac{1}{77.8904}+\frac{1}{383.1096}+\frac{1}{15.1096}+\frac{1}{118.8904}}=11.1058\end{array}$

We proceed in the same way with the calculation of the expected value and the conditioned variance to the common OR of the exposed cases in the second stratum. The observed joint and marginal frequencies of stratum 2 (fair air quality) are shown in Table 8.

The expected value of a₂ (exposed cases from the second stratum) conditioned to the Mantel-Haenszel common OR [2] in Equation (23), upon having solved the quadratic Equation (22) and chosen a result within the range of 0 to 171 (minimum of n₁₂ and m₁₂), is equal to 143.9868. As it happened in the first stratum, this expected frequency is close to the observed one, which is equal to 145.

$\begin{array}{l}\left(1-{\widehat{OR}}_{MH}\right)\times E^2\left(a_2|O{R}_{MH}\right)+\left[n_{02}-m_{12}+{\widehat{OR}}_{MH}\left(n_{12}+m_{12}\right)\right]\\ \times E\left(a_2|O{R}_{MH}\right)-m_{12}{n}_{12}{\widehat{OR}}_{MH}=0\end{array}$

$\begin{array}{l}\left(1-1.5998\right)\times E^2\left(a_2|O{R}_{MH}\right)+\left[220-171+1.5998\times \left(787+171\right)\right]\\ \times E\left(a_2|O{R}_{MH}\right)-171\times 787\times 1.5998=0\end{array}$

$-0.5998\times E^2\left(a_2|O{R}_{MH}\right)+1581.5696\times E\left(a_2|O{R}_{MH}\right)-215290.8355=0$

$\alpha =-0.5998$

$\beta =1581.5696$

$\gamma =-215290.8355$

$\begin{array}{c}E\left(a_2|O{R}_{MH}\right)=\frac{-\beta \mp \sqrt{{\beta }^2-4\alpha \gamma }}{2\alpha }\\ =\frac{-1581.5696\mp \sqrt{{1581.5696}^2-4\times \left(-0.5998\right)\times \left(-215290.8355\right)}}{2\times \left(-0.5998\right)}\\ =\{\begin{array}{l}2493.0195\\ 143.9868\end{array}\end{array}$

$0<E\left(a_2|O{R}_{MH}\right)=143.9868\le \min \left(m_{12},n_{12}\right)=\min \left(171,787\right)=171$

The expected frequencies for b₂, c₂ and d₂ are also close to the observed ones, which are equal to 642, 26, and 194, respectively. These calculations correspond to Equations (24), (25), and (26).

$E\left(b_2|O{R}_{MH}\right)=n_{12}-E\left(a_2|O{R}_{MH}\right)=787-143.9868=643.0132$

$E\left(c_2|O{R}_{MH}\right)=m_{12}-E\left(a_2|O{R}_{MH}\right)=171-143.9868=27.0132$

$\begin{array}{c}E\left(d_2|O{R}_{MH}\right)=T_2-m_{12}-n_{12}+E\left(a_2|O{R}_{MH}\right)\\ =1007-171-787-143.9868=192.9868\end{array}$

Table 8. Joint and marginal frequencies of stratum 2 (fair air quality).

Note. X = exposure to indoor tobacco smoke at home, Y = bronchial asthma in pediatric minors. ∑_j = sum per rows (j = 0 and 1), ∑_i sum per columns (i = 0 and 1).

The sum of the four expected frequencies is equal to the total of 1007 participants in the second stratum, and it can be verified that the expected value of a₂ allows to correctly obtain the OR_MH. See Equation (27).

$\begin{array}{c}{T}_2=E\left(a_2|O{R}_{MH}\right)+E\left(b_2|O{R}_{MH}\right)+E\left(c_2|O{R}_{MH}\right)+E\left(d_2|O{R}_{MH}\right)\\ =143.9868+643.0132+27.0132+192.9868=1007\end{array}$

$\begin{array}{c}{\widehat{OR}}_{MH}=\frac{E\left(a_2|O{R}_{MH}\right)\left(n_{02}-m_{12}+E\left(a_2|O{R}_{MH}\right)\right)}{\left(m_{12}-E\left(a_2|O{R}_{MH}\right)\right)\left(n_{12}-E\left(a_2|O{R}_{MH}\right)\right)}\\ =\frac{143.9868\times \left(787-171+143.9868\right)}{\left(171-143.9868\right)\left(787-143.9868\right)}=1.5998\end{array}$

The variance of the exposed cases, a₂, conditioned to the common OR estimated by the Mantel and Haenszel formula [2] in the second stratum, is equal to 19.7235. See Equation (28).

$\begin{array}{c}{s}^2\left(a_2|O{R}_{MH}\right)=\frac{1}{\frac{1}{E\left(a_2|O{R}_{MH}\right)}+\frac{1}{E\left(b_2|O{R}_{MH}\right)}+\frac{1}{E\left(c_2|O{R}_{MH}\right)}+\frac{1}{E\left(d_2|O{R}_{MH}\right)}}\\ =\frac{1}{\frac{1}{143.9868}+\frac{1}{643.0132}+\frac{1}{27.0132}+\frac{1}{192.9868}}=19.7235\end{array}$

Finally, the expected values and the variance conditioned to the common OR of the cases exposed in the third stratum are calculated. The observed joint and marginal frequencies of stratum 3 (poor air quality) are shown in Table 9.

The expected value of a₃ (exposed cases from the third stratum) conditioned to the Mantel-Haenszel common OR [2] in Equation (23), upon having solved the quadratic Equation (22) and chosen a result within the range of 0 to 86 (minimum of n₁₃ and m₁₃), is equal to 65.2416. As in the first and second strata, this expected frequency is close to the observed one, which is equal to 62.

$\begin{array}{l}\left(1-{\widehat{OR}}_{MH}\right)\times E^2\left(a_3|O{R}_{MH}\right)+\left[n_{03}-m_{13}+{\widehat{OR}}_{MH}\left(n_{13}+m_{13}\right)\right]\\ \times E\left(a_3|O{R}_{MH}\right)-m_{13}{n}_{13}{\widehat{OR}}_{MH}=0\end{array}$

$\begin{array}{l}\left(1-1.5998\right)\times E^2\left(a_3|O{R}_{MH}\right)+\left[126-86+1.5998\times \left(272+86\right)\right]\\ \times E\left(a_3|O{R}_{MH}\right)-86\times 272\times 1.5998=0\end{array}$

Table 9. Joint and marginal frequencies of stratum 3 (poor quality).

Note. X = exposure to indoor tobacco smoke at home, Y = bronchial asthma in pediatric minors. ∑_j = sum per rows (j = 0 and 1), ∑_i sum per columns (i = 0 and 1).

$-0.5998\times E^2\left(a_3|O{R}_{MH}\right)+612.7139\times E\left(a_3|O{R}_{MH}\right)-37421.5744=0$

$\alpha =-0.5998$

$\beta =612.7139$

$\gamma =-37421.5744$

$\begin{array}{c}E\left(a_2|O{R}_{MH}\right)=\frac{-\beta \mp \sqrt{{\beta }^2-4\alpha \gamma }}{2\alpha }\\ =\frac{-612.7139\mp \sqrt{{612.7139}^2-4\times \left(-0.5998\right)\times \left(-37421.5744\right)}}{2\times \left(-0.5998\right)}\\ =\{\begin{array}{l}956.3577\\ 65.2416\end{array}\end{array}$

$0<E\left(a_3|O{R}_{MH}\right)=65.2416\le \min \left(m_{13},n_{13}\right)=\min \left(86,272\right)=86$

The expected frequencies for b₃, c₃, andd₃ are also close to the observed ones, which are equal to 210, 24, and 102, respectively. These calculations correspond to Equations (24), (25), and (26).

$E\left(c_3|O{R}_{MH}\right)=m_{13}-E\left(a_3|O{R}_{MH}\right)=86-65.2416=20.7584$

$\begin{array}{c}E\left(d_3|O{R}_{MH}\right)=T_3-m_{13}-n_{13}+E\left(a_3|O{R}_{MH}\right)\\ =398-86-272-65.2416=105.2416\end{array}$

The sum of the four expected frequencies is equal to the total of 398 participants in the third stratum, and it can be verified that the expected value of a₃ allows to correctly obtain the OR_MH. See Equation (27).

$\begin{array}{c}{s}^2\left(a_3|O{R}_{MH}\right)=\frac{1}{\frac{1}{E\left(a_3|O{R}_{MH}\right)}+\frac{1}{E\left(b_3|O{R}_{MH}\right)}+\frac{1}{E\left(c_3|O{R}_{MH}\right)}+\frac{1}{E\left(d_3|O{R}_{MH}\right)}}\\ =\frac{1}{\frac{1}{65.2416}+\frac{1}{206.7584}+\frac{1}{20.7584}+\frac{1}{105.2416}}=12.8470\end{array}$

Once the expected values of the four frequencies of each stratum and the variances of the exposed cases with regard to the common OR of each stratum are available, the Breslow-Day Q test statistic [4] can be calculated, as shown in Equation (21). This statistic can be calculated using two formulas. The first formula is a chi-square goodness-of-fit statistic and it needs the expected values for the four frequencies of each stratum. See Table 10. See Equation (28).

Table 10. Joint frequencies observed in the three strata, expected frequencies, and variances conditioned upon the assumption of a common OR estimated by the Mantel-Haenszel formula.

Note. Z_k = confounding variable with three ordinal categories = air quality in the place of residence (k: 1 = good, 2 = fair, and 3 = poor), X_i = exposure with two categories (i: 0 = no and 1 = yes), Y_j = bronchial asthma in pediatric minors with two categories (j: 0 = non case and 1 = case), O_ijk = joint absolute frequency of category i of X and category j of Y in stratum k of Z, E_ijk = expected absolute frequency conditioned to a common odds ratio estimated by the Mantel-Haenszel formula, (O_ijk − E_ijk)² = square of the difference between the observed and expected probability, (O_ijk − E_ijk)²/E_ijk probability divided by the expected frequency, s²(O_ijk|) = Var(O_ijk|OR_MH) = variance of the observed frequencies conditioned to a common odds ratio estimated by the Mantel-Haenszel formula, (O − E)²/s² = (O_11k − E_11k)²/Var(O_11k|OR_MH) = quotient between the square of the difference between the observed frequency and the expected frequency of the exposed cases (numerator) and the variance of the exposed cases conditioned to a common odds ratio estimated by the Mantel-Haenszel formula (denominator), ∑ = sum per column.

$\begin{array}{c}{Q}_{BD}=\underset{i=1}{\overset{k}{{\displaystyle \sum }}}\frac{{\left(a_i-E\left(a_i|O{R}_{MH}\right)\right)}^2}{E\left(a_i|O{R}_{MH}\right)}+\frac{{\left(b_i-E\left(b_i|O{R}_{MH}\right)\right)}^2}{E\left(b_i|O{R}_{MH}\right)}\\ +\frac{{\left(c_i-E\left(c_i|O{R}_{MH}\right)\right)}^2}{E\left(c_i|O{R}_{MH}\right)}+\frac{{\left(d_i-E\left(d_i|O{R}_{MH}\right)\right)}^2}{E\left(d_i|O{R}_{MH}\right)}\end{array}$

$\begin{array}{c}{Q}_{BD}=0.0571+0.0116+0.2945+0.0374+0.0071+0.0016\\ +0.0380+0.0053+0.1611+0.0508+0.5062+0.0998\\ =1.2707\end{array}$

The second formula, given in Equation (30), is a simplified form of the first one, shown in Equation (21), considering that the squares of the differences between the observed and expected frequency for each stratum are equal, as well as the variances conditioned to a common odds ratio. This second formula needs the expected frequencies of the exposed cases of the three strata and their variances conditioned to the Mantel-Haenszel common OR. This mathematical expectation is calculated with Equation (31) and the variance is calculated with Equation (32). See Table 10.

$\begin{array}{c}{Q}_{BD}=\underset{i=1}{\overset{k}{{\displaystyle \sum }}}\frac{{\left(a_i-E\left(a_i|O{R}_{MH}\right)\right)}^2}{s^2\left(a_i|O{R}_{MH}\right)}\\ =\frac{{\left(80-77.8904\right)}^2}{11.1058}+\frac{{\left(145-143.9868\right)}^2}{19.7235}+\frac{{\left(62-65.2416\right)}^2}{12.8470}\\ =0.400729+0.052047+0.817930=1.270705<{}_{1-\alpha }\chi {}_{k-1}^2={}_{0.95}\chi {}_2^2=5.991\end{array}$

$P\left({\chi }_2^2\ge 1.270705\right)=0.529749>\alpha =0.05$

The Breslow-Day statistic with Tarone’s correction [5] undergoes a minimal change than without correction. See Equation (30), (31), and (32).

$a={\displaystyle \underset{i=1}{\overset{k}{\sum }}{a}_i}=80+145+62=287$

$E\left(a|O{R}_{MH}\right)={\displaystyle \underset{i=1}{\overset{k}{\sum }}E\left(a_i|m_{1i},O{R}_{MH}\right)}=77.8904+143.9868+65.2416=287.1188$

$s^2\left(a|O{R}_{MH}\right)={\displaystyle \underset{i=1}{\overset{k}{\sum }}Var\left(a_i|m_{1i},O{R}_{MH}\right)}=11.1058+19.7235+12.8470=43.6763$

$\begin{array}{c}{Q}_{BDT}=Q_{BD}-\frac{{\left(a-E\left(a|O{R}_{MH}\right)\right)}^2}{s^2\left(a|O{R}_{MH}\right)}\\ ={\displaystyle \underset{i=1}{\overset{k}{\sum }}\frac{a_i-E\left(a_i|m_{1i},O{R}_{MH}\right)}{s^2\left(a_i|m_{1i},O{R}_{MH}\right)}}-\frac{{\left(a-E\left(a|O{R}_{MH}\right)\right)}^2}{s^2\left(a|O{R}_{MH}\right)}\\ =1.270705-\frac{{\left(287-287.1188\right)}^2}{43.6763}\\ =1.270705-0.000323=1.270382\end{array}$

$P\left({\chi }_2^2\ge 1.270382\right)=0.529834>\alpha =0.05$

11. Statistical Power and Type-2 Error

When rejecting the null hypothesis, it is important to also report the statistical power or probability of rejecting the null hypothesis conditional on it being false (ϕ). One way to calculate power is by the complement of the cumulative distribution function of the noncentral chi-square distribution or right-tailed probability. The degrees of freedom are 1 (gl = 1). The non-centrality parameter is the value of the test statistic: PNC = Q_CMH. This function is evaluated at the critical point: ${}_{1-\alpha }\chi {}_1^2$ . If the null hypothesis is accepted, the type-2 error or probability of maintaining the null hypothesis conditional on it being false (β) is calculated, which is the probability at the left tail of the previously defined no-central distribution.

Returning to the example in which the null hypothesis of no effect was rejected and using the Q_CMH statistic with the Yates’ correction, the statistical power is high:

$\varphi =1-\beta =1-NC{\chi }_{gl=1,PNC=9.2259}^2\left({}_{0.95}\chi {}_1^2=3.8415\right)=1-0.1406=0.8594$

12. Conclusions

Based on the analyses conducted on the simulated and realistic data presented in this study, it can be concluded that exposure to indoor tobacco smoke at home has a detrimental effect on bronchial asthma in pediatric males, when statistically controlling for the confounding variable of air quality in the urban area where the participants reside. The effect size of tobacco smoke is small, OR_MH = 1.60, d = 0.26, and its effect is homogeneous across the three strata: good, fair, and poor air quality. Exposure is 1.6 times more likely than non-exposure in males 0 - 14 years of age who suffer bronchial asthma compared with age-matched males without bronchial asthma. The implication that would follow from these results is that community preventive campaigns and pediatric interventions should focus on or take into account this health risk, namely smoking by parents or adults in the home. On the other hand, the confounding variable also has a significant effect, making air quality another health issue to be addressed by health authorities. These conclusions are similar to those obtained by Paciência, Cavaleiro-Rufo, and Moreira in a review study [27] .

As has been done in this study, it is important to check the assumptions of the test, calculate the statistical power in case the null hypothesis of no effect is rejected or the type-2 error in case it holds, report the effect size in case the Mantel-Haenszel common OR is non-zero, and adjust for the limitations involved in using qualitative variables, one dichotomous as antecedent, another dichotomous as an outcome, and a third statistically controlled polychotomous variable. A broader alternative to the CMH test is conditional logistic regression. This latter statistical technique is used to analyze the relationship between a dependent variable and independent variables in matched or clustered data. It accounts for within-set dependencies and allows for controlling confounding factors. It estimates the odds ratio and uses maximum likelihood estimation to estimate regression coefficients, providing a valuable tool for analyzing matched or clustered data while controlling for confounders.

Acknowledgements

The authors thank the reviewers and editor for their helpful comments.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Cochran, W.G. (1954) On the Methods for Strengthening the Common χ2 Tests. Biometrics, 10, 417-451. https://doi.org/10.2307/3001616
[2]	Mantel, N. and Haenszel, W. (1959) Statistical Aspects of the Analysis of Data from Retrospective Studies of Disease. Journal of the National Cancer Institute, 22, 719-748.
[3]	Woolf, B. (1955) On Estimating the Relation between Blood Group and Disease. Annals of Human Genetics, 19, 251-253. https://doi.org/10.1111/j.1469-1809.1955.tb01348.x
[4]	Breslow, N.E. and Day, N.E. (1980) Statistical Methods in Cancer Research: The Analysis of Case-Control Studies. International Agency for Research on Cancer Scientific, Lyon.
[5]	Tarone, R.E. (1985) On Heterogeneity Tests Based on Efficient Scores. Biometrika, 72, 91-95. https://doi.org/10.1093/biomet/72.1.91
[6]	Agresti, A. (2019) An Introduction to the Categorical Data Analysis. 3th Edition, John Wiley and Sons, Hoboken.
[7]	Yilmaz, A.E. and Altunay, S.A. (2022) Post-Hoc Comparison Tests for Odds Ratios. Electronic Journal of Applied Statistical Analysis, 15, 75-94.
[8]	Schwarzer, G. and Rücker, G. (2021) Meta-Analysis of Proportions. In: Evangelou, E. and Veroniki, A.A., Ed., Meta-Research: Methods and Protocols, Springer, New York, 159-172. https://doi.org/10.1007/978-1-0716-1566-9_10
[9]	International Business Machines Corporation (2022) IBM SPSS Statistics Base 28. https://www.ibm.com/docs/en/SSLVMB_28.0.0/pdf/IBM_SPSS_Statistics_Base.pdf
[10]	Robins, J.M., Breslow, N. and Greenland, S. (1986) Estimators of the Mantel-Haenszel Variance Consistent in Both Sparse Data and Large-Strata Limiting Models. Biometrics, 42, 311-323. https://doi.org/10.2307/2531052
[11]	Böhning, D., Sangnawakij, P. and Holling, H. (2022) Confidence Interval Estimation for The Mantel—Haenszel Estimator of the Risk Ratio and Risk Difference in Rare Event Meta-Analysis with Emphasis on the Bootstrap. Journal of Statistical Computation and Simulation, 92, 1267-1291. https://doi.org/10.1080/00949655.2021.1991347
[12]	Cucchetti, A. (2022) Analyzing Categorical Variable: Descriptive Statistics and Comparisons. In: Ceresoli, M., Abu-Zidan, F.M., Staudenmayer, K.L., Catena, F. and Coccolini, F., Eds., Statistics and Research Methods for Acute Care and General Surgeons. Hot Topics in Acute Care Surgery and Trauma, Springer, Cham, 67-73. https://doi.org/10.1007/978-3-031-13818-8_6
[13]	Deeks, J.J., Higgins, J.P. and Altman, D.G. (2019) Analysing Data and Undertaking Meta-Analyses. In: Higgins, J.P.T., Thomas, J., Chandler, J., Cumpston, M., Li, T., Page, M.J. and Welch, V.A., Eds., Cochrane Handbook for Systematic Reviews of Interventions, Wiley Online Library, New York, 1-48. https://training.cochrane.org/handbook/current/chapter-10
[14]	Kulinskaya, E. and Dollinger, M.B. (2015) An Accurate Test for Homogeneity of Odds Ratios Based on Cochran’s Q-Statistic. BMC Medical Research Methodology, 15, Article No. 49. https://doi.org/10.1186/s12874-015-0034-x
[15]	Zhang, C., Wang, X., Chen, M. and Wang, T.A. (2021) Comparison of Hypothesis Tests for Homogeneity in Meta-Analysis with Focus on Rare Binary Events. Research Synthesis Methods, 12, 408-428. https://doi.org/10.1002/jrsm.1484
[16]	Shim, S.R., Kim, S.J. and Lee, J. (2019) Diagnostic Test Accuracy: Application and Practice Using R Software. Epidemiology and Health, 41, e2019007. https://doi.org/10.4178/epih.e2019007
[17]	Almalik, O. and van den Heuvel, E.R. (2018) Testing Homogeneity of Effect Sizes in Pooling 2x2 Contingency Tables from Multiple Studies: A Comparison of Methods. Cogent Mathematics & Statistics, 5, Article ID: 1478698. https://doi.org/10.1080/25742558.2018.1478698
[18]	Zaiontz, C. (2022) Cochran-Mantel-Haenszel (CMH) Test. In Real Statistics Using Excel. https://www.real-statistics.com/chi-square-and-f-distributions/cochran-mantel-haenszel/
[19]	Meyer, D., Zeileis, A., Hornik, K. and Friendly, M. (2022) Woolf Test for Homogeneity in 2 × 2 × k Tables. https://rdrr.io/cran/DescTools/man/WoolfTest.html
[20]	Vierra, A., Razzaq, A. and Andreadis, A. (2023) Categorical Variable Analyses: Chi-Square, Fisher Exact, and Mantel-Haenszel. In: Eltorai, A.E.M., Bakal, J.A., Newell, P.C. and Osband, A.J., Eds., Translational Surgery. Handbook for Designing and Conducting Clinical and Translational Research, Elsevier Science, Amsterdam, 171-175. https://doi.org/10.1016/B978-0-323-90300-4.00095-1
[21]	Breslow, N.E. (1996) Statistics in Epidemiology: The Case-Control Study. Journal of the American Statistical Association, 91, 14-28. https://doi.org/10.1080/01621459.1996.10476660
[22]	Korthauer, K., Kimes, P.K., Duvallet, C., Reyes, A., Subramanian, A., Teng, M., Shukla, C., Alm, E.J. and Hicks, S.C. (2019) A Practical Guide to Methods Controlling False Discoveries in Computational Biology. Genome Biology, 20, Article No. 118. https://doi.org/10.1186/s13059-019-1716-1
[23]	Breslow, N.E., Day, N.E., Halvorsen, K.T., Prentice, R.L. and Sabai, C. (1978) Estimation of Multiple Relative Risk Functions in Matched Case-Control Studies. American Journal of Epidemiology, 108, 299-307. https://doi.org/10.1093/oxfordjournals.aje.a112623
[24]	Avalos, M., Pouyes, H., Kwemou, M. and Xu, B. (2022) Package ‘clogitLasso’. https://cran.r-project.org/web/packages/clogitLasso/clogitLasso.pdf
[25]	World Health Association (2020) Eleventh Revision of the International Classification of Diseases (ICD-11). Geneva. https://www.who.int/classifications/icd/en/
[26]	Secretaria de Gobernación (2019) Norma Oficial Mexicana NOM-172-SEMARNAT-2019. Lineamientos para la Obtención y Comunicación del índice de Calidad del Aire y Riesgos a la Salud [Mexican Official Standard NOM-172-SEMARNAT-2019. Guidelines for Obtaining and Communicating the Air Quality Index and Health Risks]. Diario Oficial de la Federación. https://www.dof.gob.mx/nota_detalle.php?codigo=5579387&fecha=20/11/2019#gsc.tab=0
[27]	Paciência, I., Cavaleiro-Rufo, J. and Moreira, A. (2022) Environmental Inequality: Air Pollution and Asthma in Children. Pediatric Allergy and Immunology, 33, e13818. https://doi.org/10.1111/pai.13818