epi.RdThe mhgr function computes the Cochran-Mantel-Haenszel stratified
risk ratio and its confidence limits using the Greenland-Robins variance
estimator.
The lrcum function takes the results of a series of 2x2 tables
representing the relationship between test positivity and diagnosis and
computes positive and negative likelihood ratios (with all their
deficiencies) and the variance of
their logarithms. Cumulative likelihood ratios and their confidence
intervals (assuming independence of tests) are computed, assuming a
string of all positive tests or a string of all negative tests. The
method of Simel et al as described in Altman et al is used.
a binary response variable
a variable with two unique values specifying comparison groups
the stratification variable
confidence level
an object created by mhgr or lrcum
frequency of true positive tests
frequency of false positive tests
frequency of false negative tests
frequency of true negative tests
number of places to the right of the decimal to print for
lrcum
addtitional arguments to be passed to other print functions
Uses equations 4 and 13 from Greenland and Robins.
a list of class "mhgr" or of class "lrcum".
Greenland S, Robins JM (1985): Estimation of a common effect parameter from sparse follow-up data. Biometrics 41:55-68.
Altman DG, Machin D, Bryant TN, Gardner MJ, Eds. (2000): Statistics with Confidence, 2nd Ed. Bristol: BMJ Books, 105-110.
Simel DL, Samsa GP, Matchar DB (1991): Likelihood ratios with confidence: sample size estimation for diagnostic test studies. J Clin Epi 44:763-770.
# Greate Migraine dataset used in Example 28.6 in the SAS PROC FREQ guide
d <- expand.grid(response=c('Better','Same'),
treatment=c('Active','Placebo'),
sex=c('female','male'))
d$count <- c(16, 11, 5, 20, 12, 16, 7, 19)
d
#> response treatment sex count
#> 1 Better Active female 16
#> 2 Same Active female 11
#> 3 Better Placebo female 5
#> 4 Same Placebo female 20
#> 5 Better Active male 12
#> 6 Same Active male 16
#> 7 Better Placebo male 7
#> 8 Same Placebo male 19
# Expand data frame to represent raw data
r <- rep(1:8, d$count)
d <- d[r,]
with(d, mhgr(response=='Better', treatment, sex))
#> Mantel-Haenszel Risk Ratio and 0.95 Greenland-Robins Confidence Interval
#>
#> Common Relative Risk: 2.163597 CI: 1.233568 3.794806
#>
#> N in Each Group
#>
#> group
#> Active Placebo
#> 55 51
# Discrete survival time example, to get Cox-Mantel relative risk and CL
# From Stokes ME, Davis CS, Koch GG, Categorical Data Analysis Using the
# SAS System, 2nd Edition, Sectino 17.3, p. 596-599
#
# Input data in Table 17.5
d <- expand.grid(treatment=c('A','P'), center=1:3)
d$healed2w <- c(15,15,17,12, 7, 3)
d$healed4w <- c(17,17,17,13,17,17)
d$notHealed4w <- c( 2, 7,10,15,16,18)
d
#> treatment center healed2w healed4w notHealed4w
#> 1 A 1 15 17 2
#> 2 P 1 15 17 7
#> 3 A 2 17 17 10
#> 4 P 2 12 13 15
#> 5 A 3 7 17 16
#> 6 P 3 3 17 18
# Reformat to the way most people would collect raw data
d1 <- d[rep(1:6, d$healed2w),]
d1$time <- '2'
d1$y <- 1
d2 <- d[rep(1:6, d$healed4w),]
d2$time <- '4'
d2$y <- 1
d3 <- d[rep(1:6, d$notHealed4w),]
d3$time <- '4'
d3$y <- 0
d <- rbind(d1, d2, d3)
d$healed2w <- d$healed4w <- d$notHealed4w <- NULL
d
#> treatment center time y
#> 1 A 1 2 1
#> 1.1 A 1 2 1
#> 1.2 A 1 2 1
#> 1.3 A 1 2 1
#> 1.4 A 1 2 1
#> 1.5 A 1 2 1
#> 1.6 A 1 2 1
#> 1.7 A 1 2 1
#> 1.8 A 1 2 1
#> 1.9 A 1 2 1
#> 1.10 A 1 2 1
#> 1.11 A 1 2 1
#> 1.12 A 1 2 1
#> 1.13 A 1 2 1
#> 1.14 A 1 2 1
#> 2 P 1 2 1
#> 2.1 P 1 2 1
#> 2.2 P 1 2 1
#> 2.3 P 1 2 1
#> 2.4 P 1 2 1
#> 2.5 P 1 2 1
#> 2.6 P 1 2 1
#> 2.7 P 1 2 1
#> 2.8 P 1 2 1
#> 2.9 P 1 2 1
#> 2.10 P 1 2 1
#> 2.11 P 1 2 1
#> 2.12 P 1 2 1
#> 2.13 P 1 2 1
#> 2.14 P 1 2 1
#> 3 A 2 2 1
#> 3.1 A 2 2 1
#> 3.2 A 2 2 1
#> 3.3 A 2 2 1
#> 3.4 A 2 2 1
#> 3.5 A 2 2 1
#> 3.6 A 2 2 1
#> 3.7 A 2 2 1
#> 3.8 A 2 2 1
#> 3.9 A 2 2 1
#> 3.10 A 2 2 1
#> 3.11 A 2 2 1
#> 3.12 A 2 2 1
#> 3.13 A 2 2 1
#> 3.14 A 2 2 1
#> 3.15 A 2 2 1
#> 3.16 A 2 2 1
#> 4 P 2 2 1
#> 4.1 P 2 2 1
#> 4.2 P 2 2 1
#> 4.3 P 2 2 1
#> 4.4 P 2 2 1
#> 4.5 P 2 2 1
#> 4.6 P 2 2 1
#> 4.7 P 2 2 1
#> 4.8 P 2 2 1
#> 4.9 P 2 2 1
#> 4.10 P 2 2 1
#> 4.11 P 2 2 1
#> 5 A 3 2 1
#> 5.1 A 3 2 1
#> 5.2 A 3 2 1
#> 5.3 A 3 2 1
#> 5.4 A 3 2 1
#> 5.5 A 3 2 1
#> 5.6 A 3 2 1
#> 6 P 3 2 1
#> 6.1 P 3 2 1
#> 6.2 P 3 2 1
#> 11 A 1 4 1
#> 1.17 A 1 4 1
#> 1.21 A 1 4 1
#> 1.31 A 1 4 1
#> 1.41 A 1 4 1
#> 1.51 A 1 4 1
#> 1.61 A 1 4 1
#> 1.71 A 1 4 1
#> 1.81 A 1 4 1
#> 1.91 A 1 4 1
#> 1.101 A 1 4 1
#> 1.111 A 1 4 1
#> 1.121 A 1 4 1
#> 1.131 A 1 4 1
#> 1.141 A 1 4 1
#> 1.15 A 1 4 1
#> 1.16 A 1 4 1
#> 21 P 1 4 1
#> 2.17 P 1 4 1
#> 2.21 P 1 4 1
#> 2.31 P 1 4 1
#> 2.41 P 1 4 1
#> 2.51 P 1 4 1
#> 2.61 P 1 4 1
#> 2.71 P 1 4 1
#> 2.81 P 1 4 1
#> 2.91 P 1 4 1
#> 2.101 P 1 4 1
#> 2.111 P 1 4 1
#> 2.121 P 1 4 1
#> 2.131 P 1 4 1
#> 2.141 P 1 4 1
#> 2.15 P 1 4 1
#> 2.16 P 1 4 1
#> 31 A 2 4 1
#> 3.17 A 2 4 1
#> 3.21 A 2 4 1
#> 3.31 A 2 4 1
#> 3.41 A 2 4 1
#> 3.51 A 2 4 1
#> 3.61 A 2 4 1
#> 3.71 A 2 4 1
#> 3.81 A 2 4 1
#> 3.91 A 2 4 1
#> 3.101 A 2 4 1
#> 3.111 A 2 4 1
#> 3.121 A 2 4 1
#> 3.131 A 2 4 1
#> 3.141 A 2 4 1
#> 3.151 A 2 4 1
#> 3.161 A 2 4 1
#> 41 P 2 4 1
#> 4.15 P 2 4 1
#> 4.21 P 2 4 1
#> 4.31 P 2 4 1
#> 4.41 P 2 4 1
#> 4.51 P 2 4 1
#> 4.61 P 2 4 1
#> 4.71 P 2 4 1
#> 4.81 P 2 4 1
#> 4.91 P 2 4 1
#> 4.101 P 2 4 1
#> 4.111 P 2 4 1
#> 4.12 P 2 4 1
#> 51 A 3 4 1
#> 5.17 A 3 4 1
#> 5.21 A 3 4 1
#> 5.31 A 3 4 1
#> 5.41 A 3 4 1
#> 5.51 A 3 4 1
#> 5.61 A 3 4 1
#> 5.7 A 3 4 1
#> 5.8 A 3 4 1
#> 5.9 A 3 4 1
#> 5.10 A 3 4 1
#> 5.11 A 3 4 1
#> 5.12 A 3 4 1
#> 5.13 A 3 4 1
#> 5.14 A 3 4 1
#> 5.15 A 3 4 1
#> 5.16 A 3 4 1
#> 61 P 3 4 1
#> 6.18 P 3 4 1
#> 6.21 P 3 4 1
#> 6.3 P 3 4 1
#> 6.4 P 3 4 1
#> 6.5 P 3 4 1
#> 6.6 P 3 4 1
#> 6.7 P 3 4 1
#> 6.8 P 3 4 1
#> 6.9 P 3 4 1
#> 6.10 P 3 4 1
#> 6.11 P 3 4 1
#> 6.12 P 3 4 1
#> 6.13 P 3 4 1
#> 6.14 P 3 4 1
#> 6.15 P 3 4 1
#> 6.16 P 3 4 1
#> 12 A 1 4 0
#> 1.18 A 1 4 0
#> 22 P 1 4 0
#> 2.18 P 1 4 0
#> 2.22 P 1 4 0
#> 2.32 P 1 4 0
#> 2.42 P 1 4 0
#> 2.52 P 1 4 0
#> 2.62 P 1 4 0
#> 32 A 2 4 0
#> 3.18 A 2 4 0
#> 3.22 A 2 4 0
#> 3.32 A 2 4 0
#> 3.42 A 2 4 0
#> 3.52 A 2 4 0
#> 3.62 A 2 4 0
#> 3.72 A 2 4 0
#> 3.82 A 2 4 0
#> 3.92 A 2 4 0
#> 42 P 2 4 0
#> 4.16 P 2 4 0
#> 4.22 P 2 4 0
#> 4.32 P 2 4 0
#> 4.42 P 2 4 0
#> 4.52 P 2 4 0
#> 4.62 P 2 4 0
#> 4.72 P 2 4 0
#> 4.82 P 2 4 0
#> 4.92 P 2 4 0
#> 4.102 P 2 4 0
#> 4.112 P 2 4 0
#> 4.121 P 2 4 0
#> 4.13 P 2 4 0
#> 4.14 P 2 4 0
#> 52 A 3 4 0
#> 5.18 A 3 4 0
#> 5.22 A 3 4 0
#> 5.32 A 3 4 0
#> 5.42 A 3 4 0
#> 5.52 A 3 4 0
#> 5.62 A 3 4 0
#> 5.71 A 3 4 0
#> 5.81 A 3 4 0
#> 5.91 A 3 4 0
#> 5.101 A 3 4 0
#> 5.111 A 3 4 0
#> 5.121 A 3 4 0
#> 5.131 A 3 4 0
#> 5.141 A 3 4 0
#> 5.151 A 3 4 0
#> 62 P 3 4 0
#> 6.19 P 3 4 0
#> 6.22 P 3 4 0
#> 6.31 P 3 4 0
#> 6.41 P 3 4 0
#> 6.51 P 3 4 0
#> 6.61 P 3 4 0
#> 6.71 P 3 4 0
#> 6.81 P 3 4 0
#> 6.91 P 3 4 0
#> 6.101 P 3 4 0
#> 6.111 P 3 4 0
#> 6.121 P 3 4 0
#> 6.131 P 3 4 0
#> 6.141 P 3 4 0
#> 6.151 P 3 4 0
#> 6.161 P 3 4 0
#> 6.17 P 3 4 0
# Finally, duplicate appropriate observations to create 2 and 4-week
# risk sets. Healed and not healed at 4w need to be in the 2-week
# risk set as not healed
d2w <- subset(d, time=='4')
d2w$time <- '2'
d2w$y <- 0
d24 <- rbind(d, d2w)
with(d24, table(y, treatment, time, center))
#> , , time = 2, center = 1
#>
#> treatment
#> y A P
#> 0 19 24
#> 1 15 15
#>
#> , , time = 4, center = 1
#>
#> treatment
#> y A P
#> 0 2 7
#> 1 17 17
#>
#> , , time = 2, center = 2
#>
#> treatment
#> y A P
#> 0 27 28
#> 1 17 12
#>
#> , , time = 4, center = 2
#>
#> treatment
#> y A P
#> 0 10 15
#> 1 17 13
#>
#> , , time = 2, center = 3
#>
#> treatment
#> y A P
#> 0 33 35
#> 1 7 3
#>
#> , , time = 4, center = 3
#>
#> treatment
#> y A P
#> 0 16 18
#> 1 17 17
#>
# Matches Table 17.6
with(d24, mhgr(y, treatment, interaction(center, time, sep=';')))
#> Mantel-Haenszel Risk Ratio and 0.95 Greenland-Robins Confidence Interval
#>
#> Common Relative Risk: 1.255945 CI: 1.012105 1.558532
#>
#> N in Each Group
#>
#> group
#> A P
#> 197 204
# Get cumulative likelihood ratios and their 0.95 confidence intervals
# based on the following two tables
#
# Disease Disease
# + - + -
# Test + 39 3 20 5
# Test - 21 17 22 15
lrcum(c(39,20), c(3,5), c(21,22), c(17,15))
#> LR+ Lower 0.95 Upper 0.95 Cum. LR+ Lower 0.95 Upper 0.95
#> 4.333 1.502 12.503 4.333 1.502 12.503
#> 1.905 0.837 4.336 8.254 2.158 31.571
#>
#> LR- Lower 0.95 Upper 0.95 Cum. LR- Lower 0.95 Upper 0.95
#> 0.412 0.279 0.609 0.412 0.279 0.609
#> 0.698 0.476 1.025 0.288 0.166 0.497