Stratified analysis by Mantel-Haenszel method

## http://yihui.name/knitr/options#chunk_options
opts_chunk$set(comment = "", warning = FALSE, message = FALSE, echo = TRUE, 
    tidy = FALSE, fig.width = 8, fig.height = 7)
options(width = 116)

Create data.

If creating tables from raw data,

do xtabs(data = raw.data, formula = ~ diabetes2 +death2 +agecat)

## Data in array-format
tab1.3.array <- array(c( 2, 90, 2, 660,
                        16,360,13,1250,
                        41,557, 7, 774,
                        35,449, 5, 173),
                      dim = c(2,2,4),
                      list(diabetes2 = c("1 Diabetic", "2 Not diabetic"),
                           death2    = c("1 Yes", "2 No"),
                           agecat    = c(1:4)))

Crude (top) and stratified tables (second to fifth)

addmargins(tab1.3.array)[,,c(5,1:4)]
, , agecat = Sum

                death2
diabetes2        1 Yes 2 No  Sum
  1 Diabetic        94   27  121
  2 Not diabetic  1456 2857 4313
  Sum             1550 2884 4434

, , agecat = 1

                death2
diabetes2        1 Yes 2 No Sum
  1 Diabetic         2    2   4
  2 Not diabetic    90  660 750
  Sum               92  662 754

, , agecat = 2

                death2
diabetes2        1 Yes 2 No  Sum
  1 Diabetic        16   13   29
  2 Not diabetic   360 1250 1610
  Sum              376 1263 1639

, , agecat = 3

                death2
diabetes2        1 Yes 2 No  Sum
  1 Diabetic        41    7   48
  2 Not diabetic   557  774 1331
  Sum              598  781 1379

, , agecat = 4

                death2
diabetes2        1 Yes 2 No Sum
  1 Diabetic        35    5  40
  2 Not diabetic   449  173 622
  Sum              484  178 662

Create a vector for each matrix cell (ai,bi,ci,di) for manual calculation.

ai <- tab1.3.array[1,1,]
bi <- tab1.3.array[1,2,]
ci <- tab1.3.array[2,1,]
di <- tab1.3.array[2,2,]

data.frame(ai, bi, ci, di)
  ai bi  ci   di
1  2  2  90  660
2 16 13 360 1250
3 41  7 557  774
4 35  5 449  173

The overall odds ratio can be calculated as follows:

(sum(ai) * sum(di)) / (sum(bi) * sum(ci))
[1] 6.831

The odds ratios for the age category strata are:

(ai * di) / (bi * ci)
    1     2     3     4 
7.333 4.274 8.139 2.697 

Mantel-Haenszel method odds ratio:

\[ OR_{MH} = \frac {\sum_i \frac {a_i d_i} {T_i}} {\sum_i \frac {b_i c_i} {T_i}} \]

Ti <- (ai + bi + ci + di)               # 1. Total for each level separately
numerator   <- sum((ai * di) / Ti)      # 2. (ai * di) / Ti for each level, then sum
denominator <- sum((bi * ci) / Ti)      # 3. (bi * ci) / Ti for each level, then sum
OR.mh <- numerator / denominator        # 4. Divide 2 by 3
OR.mh
[1] 4.951

library(epiR) has a useful function called epi.2by2()

library(epiR)
epi.2by2(tab1.3.array, units = 1)
             Disease +    Disease -      Total        Inc risk *        Odds
Exposed +           94           27        121             0.777       3.481
Exposed -         1456         2857       4313             0.338       0.510
Total             1550         2884       4434             0.350       0.537


Point estimates and 95 % CIs:
---------------------------------------------------------
Inc risk ratio (crude)                    2.3 (2.07, 2.55)
Inc risk ratio (M-H)                      1.69 (1.54, 1.85)
Inc risk ratio (crude:M-H)                1.36
Odds ratio (crude)                        6.83 (4.43, 10.53)
Odds ratio (M-H)                          4.95 (2.5, 9.79)
Odds ratio (crude:M-H)                    1.38
Attrib risk (crude) *                     0.44 (0.36, 0.51)
Attrib risk (M-H) *                       0.32 (-0.42, 1.05)
Attrib risk (crude:M-H)                   1.39
---------------------------------------------------------
 * Cases per population unit 

Odds ratio, relative risk, risk difference (attributable risk) for each strata.

res.strata <- epi.2by2(tab1.3.array, units = 1, verbose = TRUE)
res.strata[c("OR.strata","RR.strata","AR.strata")]
$OR.strata
    est    se weight lower  upper
1 7.333 2.735 0.3633 1.020 52.706
2 4.274 1.460 0.8668 2.037  8.967
3 8.139 1.511 0.8434 3.625 18.276
4 2.697 1.626 0.7893 1.040  6.997

$RR.strata
    est    se weight lower  upper
1 4.167 1.665 0.7712 1.534 11.314
2 2.467 1.190 0.9703 1.755  3.468
3 2.041 1.070 0.9954 1.787  2.331
4 1.212 1.067 0.9958 1.068  1.376

$AR.strata
     est      se weight    lower  upper
1 0.3800 0.25028  15.96 -0.11054 0.8705
2 0.3281 0.09293 115.79  0.14598 0.5103
3 0.4357 0.05271 359.98  0.33238 0.5390
4 0.1531 0.05529 327.10  0.04477 0.2615

library(meta) is a meta-analysis package that can be used for visualization

library(meta)

meta.res <-
        meta::metabin(event.e = ai, n.e = ai + bi,
                      event.c = ci, n.c = ci + di,
                      sm = "OR", studlab = paste("agecat", c(1:4), sep = ""),
                      method="MH",
                      label.e="Exposed", label.c="Non-Exposed", outclab = "Death",
                      comb.fixed = TRUE, comb.random = FALSE)
meta.res
Outcome:    Death

           OR           95%-CI %W(fixed)
agecat1 7.333  [1.020; 52.706]      2.56
agecat2 4.274  [2.037;  8.967]     30.66
agecat3 8.139  [3.625; 18.276]     30.36
agecat4 2.697  [1.040;  6.997]     36.41

Number of studies combined: k=4

                      OR          95%-CI     z  p.value
Fixed effect model 4.952  [3.086; 7.944] 6.632 < 0.0001

Quantifying heterogeneity:
tau^2 = 0.0253; H = 1.05 [1; 2.69]; I^2 = 9.5% [0%; 86.1%]

Test of heterogeneity:
    Q d.f.  p.value
 3.31    3   0.3457

Details on meta-analytical method:
- Mantel-Haenszel method
meta::forest.meta(meta.res)

plot of chunk unnamed-chunk-11