## 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)
## 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)))
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
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
(sum(ai) * sum(di)) / (sum(bi) * sum(ci))
[1] 6.831
(ai * di) / (bi * ci)
1 2 3 4
7.333 4.274 8.139 2.697
\[ 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)
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
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)
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)