Currently the functions in the epiR package either print the results (verbose = FALSE) or return the unformatted values (verbose = TRUE). This can be improved using the S3 methods. Here, I created a rudimentary version of S3 methods for the epi.2by2() function.
What this version does)
The epi.2by2.new() function does not have the verbose option. When executed it returns an object of class “epi.2by2” instead of printing out the results. The returned object is then automatically formatted and printed by the print.epi.2by2() function, which is called by the generic function print().
The verbose = TRUE equivalent is implemented as the summary method. summary(epi.2by2_class_object) will return what we can get via verbose = TRUE.
Why bother?) This way the result can be saved in an object, and the result can be printed in the correct format at a later time, sparing the need for running epi.2by2() for the same data.
Internal workings) Internally, all measures created at the beginning of the code are saved in a list named “elements”. The “RESULTS” variables are saved in a list named “res”. The verbose equivalent values are saved in a list named “rval”. The table is saved as a formatted data frame named “tab”.
All the above, the method argument, and the conf.level arguments are passed to the returned object.
The printing part is implemented in an S3 method print.epi.2by2(). Depending on the original method argument, the results are formatted and printed.
## Load epiR for comparison
library(epiR)
Loading required package: survival
Loading required package: splines
Package epiR 0.9-51 is loaded
Type help(epi.about) for summary information
### Define new function, class, and methods
## New function using S3 methods
source("./epi.2by2.r")
## Print method
source("./print.epi.2by2.r")
## Summary method to print the "verbose" part
source("./summary.epi.2by2.r")
## Show available methods
methods(class = "epi.2by2")
[1] print.epi.2by2 summary.epi.2by2
### cohort.count
## Matrix example
mat <- matrix(c(80,3,20,47), ncol = 2)
epi.2by2(dat = mat, method = c("cohort.count","cohort.time","case.control","cross.sectional")[1],
units = 1, homogeneity = c("breslow.day","woolf")[1], verbose = FALSE
)
Disease + Disease - Total Inc risk * Odds
Exposed + 80 20 100 0.800 4.0000
Exposed - 3 47 50 0.060 0.0638
Total 83 67 150 0.553 1.2388
Point estimates and 95 % CIs:
---------------------------------------------------------
Inc risk ratio 13.33 (4.43, 40.11)
Odds ratio 60.37 (16.88, 336.62)
Attrib risk * 0.74 (0.64, 0.84)
Attrib risk in population * 0.49 (0.39, 0.6)
Attrib fraction in exposed (%) 92.5 (77.44, 97.51)
Attrib fraction in population (%) 89.16 (68.5, 96.27)
---------------------------------------------------------
* Cases per population unit
verbose <- epi.2by2(dat = mat, method = c("cohort.count","cohort.time","case.control","cross.sectional")[1],
units = 1, homogeneity = c("breslow.day","woolf")[1], verbose = TRUE
)
out <- epi.2by2.new(dat = mat, method = c("cohort.count","cohort.time","case.control","cross.sectional")[1],
units = 1, homogeneity = c("breslow.day","woolf")[1]
)
out
Disease + Disease - Total Inc risk * Odds
Exposed + 80 20 100 0.800 4.0000
Exposed - 3 47 50 0.060 0.0638
Total 83 67 150 0.553 1.2388
Point estimates and 95 % CIs:
---------------------------------------------------------
Inc risk ratio 13.33 (4.43, 40.11)
Odds ratio 60.37 (16.88, 336.62)
Attrib risk * 0.74 (0.64, 0.84)
Attrib risk in population * 0.49 (0.39, 0.60)
Attrib fraction in exposed (%) 92.50 (77.44, 97.51)
Attrib fraction in population (%) 89.16 (68.50, 96.27)
---------------------------------------------------------
* Cases per population unit
identical(verbose, summary(out))
[1] TRUE
## Array example
arr <- array(c(80,3,20,47, 50,30,20,10), dim = c(2,2,2))
epi.2by2(dat = arr, method = c("cohort.count","cohort.time","case.control","cross.sectional")[1],
units = 1, homogeneity = c("breslow.day","woolf")[1], verbose = FALSE
)
Disease + Disease - Total Inc risk * Odds
Exposed + 130 40 170 0.765 3.250
Exposed - 33 57 90 0.367 0.579
Total 163 97 260 0.627 1.680
Point estimates and 95 % CIs:
---------------------------------------------------------
Inc risk ratio (crude) 2.09 (1.57, 2.77)
Inc risk ratio (M-H) 2.13 (1.58, 2.86)
Inc risk ratio (crude:M-H) 0.98
Odds ratio (crude) 5.57 (3.1, 10.18)
Odds ratio (M-H) 5.06 (2.93, 8.72)
Odds ratio (crude:M-H) 1.1
Attrib risk (crude) * 0.4 (0.28, 0.52)
Attrib risk (M-H) * 0.4 (0.29, 0.52)
Attrib risk (crude:M-H) 0.98
---------------------------------------------------------
* Cases per population unit
verbose <- epi.2by2(dat = arr, method = c("cohort.count","cohort.time","case.control","cross.sectional")[1],
units = 1, homogeneity = c("breslow.day","woolf")[1], verbose = TRUE
)
out <- epi.2by2.new(dat = arr, method = c("cohort.count","cohort.time","case.control","cross.sectional")[1],
units = 1, homogeneity = c("breslow.day","woolf")[1]
)
out
Disease + Disease - Total Inc risk * Odds
Exposed + 130 40 170 0.765 3.250
Exposed - 33 57 90 0.367 0.579
Total 163 97 260 0.627 1.680
Point estimates and 95 % CIs:
---------------------------------------------------------
Inc risk ratio (crude) 2.09 (1.57, 2.77)
Inc risk ratio (M-H) 2.13 (1.58, 2.86)
Inc risk ratio (crude:M-H) 0.98
Odds ratio (crude) 5.57 (3.10, 10.18)
Odds ratio (M-H) 5.06 (2.93, 8.72)
Odds ratio (crude:M-H) 1.10
Attrib risk (crude) * 0.40 (0.28, 0.52)
Attrib risk (M-H) * 0.40 (0.29, 0.52)
Attrib risk (crude:M-H) 0.98
---------------------------------------------------------
* Cases per population unit
identical(verbose, summary(out))
[1] TRUE
### cohort.time
## Matrix example
mat <- matrix(c(80,30,200,500), ncol = 2)
epi.2by2(dat = mat, method = c("cohort.time","cohort.time","case.control","cross.sectional")[1],
units = 1, homogeneity = c("breslow.day","woolf")[1], verbose = FALSE
)
Disease + Time at risk Inc rate *
Exposed + 80 200 0.400
Exposed - 30 500 0.060
Total 110 700 0.157
Point estimates and 95 % CIs:
---------------------------------------------------------
Inc rate ratio 6.67 (4.33, 10.51)
Attrib rate * 0.34 (0.25, 0.43)
Attrib rate in population * 0.1 (0.06, 0.13)
Attrib fraction in exposed (%) 85 (76.92, 90.48)
Attrib fraction in population (%) 61.82 (54.78, 67.71)
---------------------------------------------------------
* Cases per unit of population time at risk
verbose <- epi.2by2(dat = mat, method = c("cohort.time","cohort.time","case.control","cross.sectional")[1],
units = 1, homogeneity = c("breslow.day","woolf")[1], verbose = TRUE
)
out <- epi.2by2.new(dat = mat, method = c("cohort.time","cohort.time","case.control","cross.sectional")[1],
units = 1, homogeneity = c("breslow.day","woolf")[1]
)
out
Disease + Time at risk Inc rate *
Exposed + 80 200 0.400
Exposed - 30 500 0.060
Total 110 700 0.157
Point estimates and 95 % CIs:
---------------------------------------------------------
Inc rate ratio 6.67 (4.33, 10.51)
Attrib rate * 0.34 (0.25, 0.43)
Attrib rate in population * 0.10 (0.06, 0.13)
Attrib fraction in exposed (%) 85.00 (76.92, 90.48)
Attrib fraction in population (%) 61.82 (54.78, 67.71)
---------------------------------------------------------
* Cases per unit of population time at risk
identical(verbose, summary(out))
[1] TRUE
## Array example
arr <- array(c(80,30,200,500, 50,30,200,100), dim = c(2,2,2))
epi.2by2(dat = arr, method = c("cohort.time","cohort.time","case.control","cross.sectional")[1],
units = 1, homogeneity = c("breslow.day","woolf")[1], verbose = FALSE
)
Disease + Time at risk Inc rate *
Exposed + 130 400 0.325
Exposed - 60 600 0.100
Total 190 1000 0.190
Point estimates and 95 % CIs:
---------------------------------------------------------
Inc rate ratio (crude) 3.25 (2.38, 4.49)
Inc rate ratio (M-H) 2.58 (1.97, 3.39)
Inc rate ratio (crude:M-H) 1.26
Attrib rate (crude) * 0.22 (0.16, 0.29)
Attrib rate (M-H) * 0.22 (0.14, 0.29)
Attrib rate (crude:M-H) 1.04
---------------------------------------------------------
* Cases per unit of population time at risk
verbose <- epi.2by2(dat = arr, method = c("cohort.time","cohort.time","case.control","cross.sectional")[1],
units = 1, homogeneity = c("breslow.day","woolf")[1], verbose = TRUE
)
out <- epi.2by2.new(dat = arr, method = c("cohort.time","cohort.time","case.control","cross.sectional")[1],
units = 1, homogeneity = c("breslow.day","woolf")[1]
)
out
Disease + Time at risk Inc rate *
Exposed + 130 400 0.325
Exposed - 60 600 0.100
Total 190 1000 0.190
Point estimates and 95 % CIs:
---------------------------------------------------------
Inc rate ratio (crude) 3.25 (2.38, 4.49)
Inc rate ratio (M-H) 2.58 (1.97, 3.39)
Inc rate ratio (crude:M-H) 1.26
Attrib rate (crude) * 0.23 (0.16, 0.29)
Attrib rate (M-H) * 0.22 (0.14, 0.29)
Attrib rate (crude:M-H) 1.04
---------------------------------------------------------
* Cases per unit of population time at risk
identical(verbose, summary(out))
[1] TRUE
### case.control
## Matrix example
mat <- matrix(c(80,30,200,500), ncol = 2)
epi.2by2(dat = mat, method = c("case.control"),
units = 1, homogeneity = c("breslow.day","woolf")[1], verbose = FALSE
)
Disease + Disease - Total Prevalence * Odds
Exposed + 80 200 280 0.2857 0.400
Exposed - 30 500 530 0.0566 0.060
Total 110 700 810 0.1358 0.157
Point estimates and 95 % CIs:
---------------------------------------------------------
Odds ratio 6.65 (4.17, 10.83)
Attrib prevalence * 0.23 (0.17, 0.29)
Attrib prevalence in population * 0.08 (0.05, 0.11)
Attrib fraction (est) in exposed (%) 84.96 (76.03, 90.76)
Attrib fraction (est) in population (%) 61.82 (48.01, 71.96)
---------------------------------------------------------
* Cases per population unit
verbose <- epi.2by2(dat = mat, method = c("case.control"),
units = 1, homogeneity = c("breslow.day","woolf")[1], verbose = TRUE
)
out <- epi.2by2.new(dat = mat, method = c("case.control"),
units = 1, homogeneity = c("breslow.day","woolf")[1]
)
out
Disease + Disease - Total Prevalence * Odds
Exposed + 80 200 280 0.2857 0.400
Exposed - 30 500 530 0.0566 0.060
Total 110 700 810 0.1358 0.157
Point estimates and 95 % CIs:
---------------------------------------------------------
Odds ratio 6.65 (4.17, 10.83)
Attrib prevalence * 0.23 (0.17, 0.29)
Attrib prevalence in population * 0.08 (0.05, 0.11)
Attrib fraction (est) in exposed (%) 84.96 (76.03, 90.76)
Attrib fraction (est) in population (%) 61.82 (48.01, 71.96)
---------------------------------------------------------
* Cases per population unit
identical(verbose, summary(out))
[1] TRUE
## Array example
arr <- array(c(80,30,200,500, 50,30,200,100), dim = c(2,2,2))
epi.2by2(dat = arr, method = c("case.control"),
units = 1, homogeneity = c("breslow.day","woolf")[1], verbose = FALSE
)
Disease + Disease - Total Prevalence * Odds
Exposed + 130 400 530 0.2453 0.325
Exposed - 60 600 660 0.0909 0.100
Total 190 1000 1190 0.1597 0.190
Point estimates and 95 % CIs:
---------------------------------------------------------
Odds ratio (crude) 3.25 (2.31, 4.6)
Odds ratio (M-H) 2.7 (1.95, 3.73)
Odds ratio (crude:M-H) 1.2
Attrib prevalence (crude) * 0.15 (0.11, 0.2)
Attrib prevalence (M-H) * 0.15 (0.08, 0.21)
Attrib prevalence (crude:M-H) 1.04
---------------------------------------------------------
* Cases per population unit
verbose <- epi.2by2(dat = arr, method = c("case.control"),
units = 1, homogeneity = c("breslow.day","woolf")[1], verbose = TRUE
)
out <- epi.2by2.new(dat = arr, method = c("case.control"),
units = 1, homogeneity = c("breslow.day","woolf")[1]
)
out
Disease + Disease - Total Prevalence * Odds
Exposed + 130 400 530 0.2453 0.325
Exposed - 60 600 660 0.0909 0.100
Total 190 1000 1190 0.1597 0.190
Point estimates and 95 % CIs:
---------------------------------------------------------
Odds ratio (crude) 3.25 (2.31, 4.60)
Odds ratio (M-H) 2.70 (1.95, 3.73)
Odds ratio (crude:M-H) 1.20
Attrib prevalence (crude) * 0.15 (0.11, 0.20)
Attrib prevalence (M-H) * 0.15 (0.08, 0.21)
Attrib prevalence (crude:M-H) 1.05
---------------------------------------------------------
* Cases per population unit
identical(verbose, summary(out))
[1] TRUE
### cross.sectional
## Matrix example
mat <- matrix(c(80,30,200,500), ncol = 2)
epi.2by2(dat = mat, method = c("cross.sectional"),
units = 1, homogeneity = c("breslow.day","woolf")[1], verbose = FALSE
)
Disease + Disease - Total Prevalence * Odds
Exposed + 80 200 280 0.2857 0.400
Exposed - 30 500 530 0.0566 0.060
Total 110 700 810 0.1358 0.157
Point estimates and 95 % CIs:
---------------------------------------------------------
Prevalence ratio 5.05 (3.4, 7.48)
Odds ratio 6.65 (4.17, 10.83)
Attrib prevalence * 0.23 (0.17, 0.29)
Attrib prevalence in population * 0.08 (0.05, 0.11)
Attrib fraction in exposed (%) 80.19 (70.63, 86.64)
Attrib fraction in population (%) 58.32 (44.2, 68.86)
---------------------------------------------------------
* Cases per population unit
verbose <- epi.2by2(dat = mat, method = c("cross.sectional"),
units = 1, homogeneity = c("breslow.day","woolf")[1], verbose = TRUE
)
out <- epi.2by2.new(dat = mat, method = c("cross.sectional"),
units = 1, homogeneity = c("breslow.day","woolf")[1]
)
out
Disease + Disease - Total Prevalence * Odds
Exposed + 80 200 280 0.2857 0.400
Exposed - 30 500 530 0.0566 0.060
Total 110 700 810 0.1358 0.157
Point estimates and 95 % CIs:
---------------------------------------------------------
Prevalence ratio 5.05 (3.40, 7.48)
Odds ratio 6.65 (4.17, 10.83)
Attrib prevalence * 0.23 (0.17, 0.29)
Attrib prevalence in population * 0.08 (0.05, 0.11)
Attrib fraction in exposed (%) 80.19 (70.63, 86.64)
Attrib fraction in population (%) 58.32 (44.20, 68.86)
---------------------------------------------------------
* Cases per population unit
identical(verbose, summary(out))
[1] TRUE
## Array example
arr <- array(c(80,30,200,500, 50,30,200,100), dim = c(2,2,2))
epi.2by2(dat = arr, method = c("cross.sectional"),
units = 1, homogeneity = c("breslow.day","woolf")[1], verbose = FALSE
)
Disease + Disease - Total Prevalence * Odds
Exposed + 130 400 530 0.2453 0.325
Exposed - 60 600 660 0.0909 0.100
Total 190 1000 1190 0.1597 0.190
Point estimates and 95 % CIs:
---------------------------------------------------------
Prevalence ratio (crude) 2.7 (2.03, 3.58)
Prevalence ratio (M-H) 2.31 (1.78, 2.98)
Prevalence ratio (crude:M-H) 1.17
Odds ratio (crude) 3.25 (2.31, 4.6)
Odds ratio (M-H) 2.7 (1.95, 3.73)
Odds ratio (crude:M-H) 1.2
Atributable prevalence (crude) * 0.15 (0.11, 0.2)
Atributable prevalence (M-H) * 0.15 (0.08, 0.21)
Atributable prevalence (crude:M-H) 1.05
---------------------------------------------------------
* Cases per population unit
verbose <- epi.2by2(dat = arr, method = c("cross.sectional"),
units = 1, homogeneity = c("breslow.day","woolf")[1], verbose = TRUE
)
out <- epi.2by2.new(dat = arr, method = c("cross.sectional"),
units = 1, homogeneity = c("breslow.day","woolf")[1]
)
out
Disease + Disease - Total Prevalence * Odds
Exposed + 130 400 530 0.2453 0.325
Exposed - 60 600 660 0.0909 0.100
Total 190 1000 1190 0.1597 0.190
Point estimates and 95 % CIs:
---------------------------------------------------------
Prevalence ratio (crude) 2.70 (2.03, 3.58)
Prevalence ratio (M-H) 2.31 (1.78, 2.98)
Prevalence ratio (crude:M-H) 1.17
Odds ratio (crude) 3.25 (2.31, 4.60)
Odds ratio (M-H) 2.70 (1.95, 3.73)
Odds ratio (crude:M-H) 1.20
Atributable prevalence (crude) * 0.15 (0.11, 0.20)
Atributable prevalence (M-H) * 0.15 (0.08, 0.21)
Atributable prevalence (crude:M-H) 1.05
---------------------------------------------------------
* Cases per population unit
identical(verbose, summary(out))
[1] TRUE
## The main function
epi.2by2.new
function (dat, method = "cohort.count", conf.level = 0.95, units = 100,
homogeneity = "breslow.day")
{
elements <- list()
elements <- within(elements, {
if (length(dim(dat)) == 2) {
a <- dat[1]
A <- a
b <- dat[3]
B <- b
c <- dat[2]
C <- c
d <- dat[4]
D <- d
}
if (length(dim(dat)) > 2) {
a <- dat[1, 1, ]
A <- a
b <- dat[1, 2, ]
B <- b
c <- dat[2, 1, ]
C <- c
d <- dat[2, 2, ]
D <- d
}
for (i in 1:length(a)) {
if (a[i] < 1 | b[i] < 1 | c[i] < 1 | d[i] < 1) {
a[i] <- a[i] + 0.5
b[i] <- b[i] + 0.5
c[i] <- c[i] + 0.5
d[i] <- d[i] + 0.5
}
}
.funincrisk <- function(dat, conf.level) {
N. <- 1 - ((1 - conf.level)/2)
a <- dat[, 1]
n <- dat[, 2]
b <- n - a
p <- a/n
a. <- ifelse(a == 0, a + 1, a)
b. <- ifelse(b == 0, b + 1, b)
low <- a./(a. + (b. + 1) * (1/qf(1 - N., 2 * a.,
2 * b. + 2)))
up <- (a. + 1)/(a. + 1 + b./(1/qf(1 - N., 2 * b.,
2 * a. + 2)))
low <- ifelse(a == 0, 0, low)
up <- ifelse(a == n, 1, up)
rval <- data.frame(p, low, up)
names(rval) <- c("est", "lower", "upper")
rval
}
.funincrate <- function(dat, conf.level) {
N. <- 1 - ((1 - conf.level)/2)
a <- dat[, 1]
n <- dat[, 2]
p <- a/n
low <- 0.5 * qchisq(p = N., df = 2 * a + 2, lower.tail = FALSE)/n
up <- 0.5 * qchisq(p = 1 - N., df = 2 * a + 2, lower.tail = FALSE)/n
rval <- data.frame(p, low, up)
names(rval) <- c("est", "lower", "upper")
rval
}
N. <- 1 - ((1 - conf.level)/2)
z <- qnorm(N., mean = 0, sd = 1)
a <- as.numeric(a)
A <- as.numeric(A)
b <- as.numeric(b)
B <- as.numeric(B)
c <- as.numeric(c)
C <- as.numeric(C)
d <- as.numeric(d)
D <- as.numeric(D)
M1 <- a + c
M0 <- b + d
N1 <- a + b
N0 <- c + d
total <- a + b + c + d
n.strata <- length(a)
if (sum(A) > 0 & sum(B) > 0 & sum(C) > 0 & sum(D) > 0) {
sa <- sum(A)
sb <- sum(B)
sc <- sum(C)
sd <- sum(D)
}
if (sum(A) == 0 | sum(B) == 0 | sum(C) == 0 | sum(D) ==
0) {
sa <- sum(a)
sb <- sum(b)
sc <- sum(c)
sd <- sum(d)
}
sM1 <- sa + sc
sM0 <- sb + sd
sN1 <- sa + sb
sN0 <- sc + sd
stotal <- sa + sb + sc + sd
tmp <- .funincrisk(as.matrix(cbind(a, N1)), conf.level = conf.level)
IRiske.p <- as.numeric(tmp[, 1]) * units
IRiske.l <- as.numeric(tmp[, 2]) * units
IRiske.u <- as.numeric(tmp[, 3]) * units
tmp <- .funincrisk(as.matrix(cbind(c, N0)), conf.level = conf.level)
IRisko.p <- as.numeric(tmp[, 1]) * units
IRisko.l <- as.numeric(tmp[, 2]) * units
IRisko.u <- as.numeric(tmp[, 3]) * units
tmp <- .funincrisk(as.matrix(cbind(M1, total)), conf.level = conf.level)
IRiskpop.p <- as.numeric(tmp[, 1]) * units
IRiskpop.l <- as.numeric(tmp[, 2]) * units
IRiskpop.u <- as.numeric(tmp[, 3]) * units
tmp <- .funincrate(as.matrix(cbind(a, b)), conf.level = conf.level)
IRatee.p <- as.numeric(tmp[, 1]) * units
IRatee.l <- as.numeric(tmp[, 2]) * units
IRatee.u <- as.numeric(tmp[, 3]) * units
tmp <- .funincrate(as.matrix(cbind(c, d)), conf.level = conf.level)
IRateo.p <- as.numeric(tmp[, 1]) * units
IRateo.l <- as.numeric(tmp[, 2]) * units
IRateo.u <- as.numeric(tmp[, 3]) * units
tmp <- .funincrate(as.matrix(cbind(M1, M0)), conf.level = conf.level)
IRatepop.p <- as.numeric(tmp[, 1]) * units
IRatepop.l <- as.numeric(tmp[, 2]) * units
IRatepop.u <- as.numeric(tmp[, 3]) * units
Al <- (qbinom(1 - N., size = a + b, prob = (a/(a + b))))/(a +
b)
Au <- (qbinom(N., size = a + b, prob = (a/(a + b))))/(a +
b)
Oe.p <- (a/b)
Oe.l <- (Al/(1 - Al))
Oe.u <- (Au/(1 - Au))
Al <- (qbinom(1 - N., size = c + d, prob = (c/(c + d))))/(c +
d)
Au <- (qbinom(N., size = c + d, prob = (c/(c + d))))/(c +
d)
Oo.p <- (c/d)
Oo.l <- (Al/(1 - Al))
Oo.u <- (Au/(1 - Au))
Al <- (qbinom(1 - N., size = M1 + M0, prob = (M1/(M1 +
M0))))/(M1 + M0)
Au <- (qbinom(N., size = M1 + M0, prob = (M1/(M1 + M0))))/(M1 +
M0)
Opop.p <- (M1/M0)
Opop.l <- (Al/(1 - Al))
Opop.u <- (Au/(1 - Au))
tmp <- .funincrisk(as.matrix(cbind(sa, sN1)), conf.level = conf.level)
cIRiske.p <- as.numeric(tmp[, 1]) * units
cIRiske.l <- as.numeric(tmp[, 2]) * units
cIRiske.u <- as.numeric(tmp[, 3]) * units
tmp <- .funincrisk(as.matrix(cbind(sc, sN0)), conf.level = conf.level)
cIRisko.p <- as.numeric(tmp[, 1]) * units
cIRisko.l <- as.numeric(tmp[, 2]) * units
cIRisko.u <- as.numeric(tmp[, 3]) * units
tmp <- .funincrisk(as.matrix(cbind(sM1, stotal)), conf.level = conf.level)
cIRiskpop.p <- as.numeric(tmp[, 1]) * units
cIRiskpop.l <- as.numeric(tmp[, 2]) * units
cIRiskpop.u <- as.numeric(tmp[, 3]) * units
tmp <- .funincrate(as.matrix(cbind(sa, sb)), conf.level = conf.level)
cIRatee.p <- as.numeric(tmp[, 1]) * units
cIRatee.l <- as.numeric(tmp[, 2]) * units
cIRatee.u <- as.numeric(tmp[, 3]) * units
tmp <- .funincrate(as.matrix(cbind(sc, sd)), conf.level = conf.level)
cIRateo.p <- as.numeric(tmp[, 1]) * units
cIRateo.l <- as.numeric(tmp[, 2]) * units
cIRateo.u <- as.numeric(tmp[, 3]) * units
tmp <- .funincrate(as.matrix(cbind(sM1, sM0)), conf.level = conf.level)
cIRatepop.p <- as.numeric(tmp[, 1]) * units
cIRatepop.l <- as.numeric(tmp[, 2]) * units
cIRatepop.u <- as.numeric(tmp[, 3]) * units
Al <- (qbinom(1 - N., size = sa + sb, prob = (sa/(sa +
sb))))/(sa + sb)
u <- (qbinom(N., size = sa + sb, prob = (sa/(sa + sb))))/(sa +
sb)
cOe.p <- sa/sb
cOe.l <- Al/(1 - Al)
cOe.u <- Au/(1 - Au)
Al <- (qbinom(1 - N., size = sc + sd, prob = (sc/(sc +
sd))))/(sc + sd)
u <- (qbinom(N., size = sc + sd, prob = (sc/(sc + sd))))/(sc +
sd)
cOo.p <- sc/sd
cOo.l <- Al/(1 - Al)
cOo.u <- Au/(1 - Au)
Al <- (qbinom(1 - N., size = sM1 + sM0, prob = (sM1/(sM1 +
sM0))))/(sM1 + sM0)
u <- (qbinom(N., size = sM1 + sM0, prob = (sM1/(sM1 +
sM0))))/(sM1 + sM0)
cOpop.p <- sM1/sM0
cOpop.l <- Al/(1 - Al)
cOpop.u <- Au/(1 - Au)
RR.p <- (a/N1)/(c/N0)
lnRR <- log(RR.p)
lnRR.var <- (1/a) - (1/N1) + (1/c) - (1/N0)
lnRR.se <- sqrt((1/a) - (1/N1) + (1/c) - (1/N0))
RR.se <- exp(lnRR.se)
RR.l <- exp(lnRR - (z * lnRR.se))
RR.u <- exp(lnRR + (z * lnRR.se))
RR.w <- 1/(exp(lnRR.var))
IRR.p <- (a/b)/(c/d)
lnIRR <- log(IRR.p)
lnIRR.var <- (1/a) + (1/c)
lnIRR.se <- sqrt((1/a) + (1/c))
IRR.se <- exp(lnIRR.se)
pl <- a/(a + (c + 1) * (1/qf(1 - N., 2 * a, 2 * c + 2)))
ph <- (a + 1)/(a + 1 + c/(1/qf(1 - N., 2 * c, 2 * a +
2)))
IRR.l <- pl * d/((1 - pl) * b)
IRR.u <- ph * d/((1 - ph) * b)
IRR.w <- 1/(exp(lnIRR.var))
OR.p <- c()
OR.l <- c()
OR.u <- c()
if (length(dim(dat)) == 3) {
for (i in 1:dim(dat)[3]) {
OR.tmp <- fisher.test(dat[, , i], conf.int = TRUE,
conf.level = conf.level)
tOR.p <- as.numeric(OR.tmp$estimate)
OR.p <- c(OR.p, tOR.p)
tOR.l <- as.numeric(OR.tmp$conf.int)[1]
OR.l <- c(OR.l, tOR.l)
tOR.u <- as.numeric(OR.tmp$conf.int)[2]
OR.u <- c(OR.u, tOR.u)
}
}
if (length(dim(dat)) == 2) {
OR.tmp <- fisher.test(dat, conf.int = TRUE, conf.level = conf.level)
tOR.p <- as.numeric(OR.tmp$estimate)
OR.p <- c(OR.p, tOR.p)
tOR.l <- as.numeric(OR.tmp$conf.int)[1]
OR.l <- c(OR.l, tOR.l)
tOR.u <- as.numeric(OR.tmp$conf.int)[2]
OR.u <- c(OR.u, tOR.u)
}
lnOR <- log(OR.p)
lnOR.var <- 1/a + 1/b + 1/c + 1/d
lnOR.se <- sqrt(1/a + 1/b + 1/c + 1/d)
lnOR.l <- lnOR - (z * lnOR.se)
lnOR.u <- lnOR + (z * lnOR.se)
OR.se <- exp(lnOR.se)
OR.w <- 1/(exp(lnOR.var))
cRR.p <- OR.p/((1 - N0) + (N0 * OR.p))
cRR.l <- OR.l/((1 - N0) + (N0 * OR.l))
cRR.u <- OR.u/((1 - N0) + (N0 * OR.u))
ARisk.p <- ((a/N1) - (c/N0)) * units
ARisk.se <- (sqrt(((a * (N1 - a))/N1^3) + ((c * (N0 -
c))/N0^3))) * units
ARisk.l <- (ARisk.p - (z * ARisk.se))
ARisk.u <- (ARisk.p + (z * ARisk.se))
ARisk.w <- 1/(ARisk.se/units)^2
ARate.p <- ((a/b) - (c/d)) * units
ARate.var <- (a/b^2) + (c/d^2)
ARate.se <- (sqrt((a/b^2) + (c/d^2))) * units
ARate.l <- ARate.p - (z * ARate.se)
ARate.u <- ARate.p + (z * ARate.se)
ARate.w <- 1/(ARate.var)
AFRisk.p <- ((RR.p - 1)/RR.p)
AFRisk.l <- (RR.l - 1)/RR.l
AFRisk.u <- (RR.u - 1)/RR.u
AFRate.p <- (IRR.p - 1)/IRR.p
AFRate.l <- (IRR.l - 1)/IRR.l
AFRate.u <- (IRR.u - 1)/IRR.u
AFest.p <- (OR.p - 1)/OR.p
AFest.l <- (OR.l - 1)/OR.l
AFest.u <- (OR.u - 1)/OR.u
PARisk.p <- ((M1/total) - (c/N0)) * units
PARisk.se <- (sqrt(((M1 * (total - M1))/total^3) + ((c *
(N0 - c))/N0^3))) * units
PARisk.l <- PARisk.p - (z * PARisk.se)
PARisk.u <- PARisk.p + (z * PARisk.se)
PARate.p <- ((M1/M0) - (c/d)) * units
PARate.se <- (sqrt((M1/M0^2) + (c/d^2))) * units
PARate.l <- PARate.p - (z * PARate.se)
PARate.u <- PARate.p + (z * PARate.se)
PAFRisk.p <- ((a * d) - (b * c))/((a + c) * (c + d))
PAFRisk.var <- (b + (PAFRisk.p * (a + d)))/(total * c)
PAFRisk.l <- 1 - exp(log(1 - PAFRisk.p) + (z * sqrt(PAFRisk.var)))
PAFRisk.u <- 1 - exp(log(1 - PAFRisk.p) - (z * sqrt(PAFRisk.var)))
PAFRate.p <- (IRatepop.p - IRateo.p)/IRatepop.p
PAFRate.l <- min((IRatepop.l - IRateo.l)/IRatepop.l,
(IRatepop.u - IRateo.u)/IRatepop.u)
PAFRate.u <- max((IRatepop.l - IRateo.l)/IRatepop.l,
(IRatepop.u - IRateo.u)/IRatepop.u)
PAFest.p <- ((a * d) - (b * c))/(d * (a + c))
PAFest.var <- (a/(c * (a + c))) + (b/(d * (b + d)))
PAFest.l <- 1 - exp(log(1 - PAFest.p) + (z * sqrt(PAFest.var)))
PAFest.u <- 1 - exp(log(1 - PAFest.p) - (z * sqrt(PAFest.var)))
cRR.p <- (sa/sN1)/(sc/sN0)
clnRR <- log(cRR.p)
clnRR.var <- (1/sa) - (1/sN1) + (1/sc) - (1/sN0)
clnRR.se <- sqrt((1/sa) - (1/sN1) + (1/sc) - (1/sN0))
clnRR.l <- clnRR - (z * clnRR.se)
clnRR.u <- clnRR + (z * clnRR.se)
cRR.se <- exp(clnRR.se)
cRR.l <- exp(clnRR.l)
cRR.u <- exp(clnRR.u)
cIRR.p <- (sa/sb)/(sc/sd)
clnIRR <- log(cIRR.p)
clnIRR.se <- sqrt((1/sa) + (1/sc))
cIRR.se <- exp(clnIRR.se)
pl <- sa/(sa + (sc + 1) * (1/qf(1 - N., 2 * sa, 2 * sc +
2)))
ph <- (sa + 1)/(sa + 1 + sc/(1/qf(1 - N., 2 * sc, 2 *
sa + 2)))
cIRR.l <- pl * sd/((1 - pl) * sb)
cIRR.u <- ph * sd/((1 - ph) * sb)
cOR.tmp <- fisher.test(apply(dat, MARGIN = c(1, 2), FUN = sum),
conf.int = TRUE, conf.level = conf.level)
cOR.p <- as.numeric(cOR.tmp$estimate)
cOR.l <- as.numeric(cOR.tmp$conf.int)[1]
cOR.u <- as.numeric(cOR.tmp$conf.int)[2]
clnOR <- log(cOR.p)
clnOR.se <- sqrt(1/sa + 1/sb + 1/sc + 1/sd)
clnOR.l <- clnOR - (z * clnOR.se)
clnOR.u <- clnOR + (z * clnOR.se)
cOR.se <- exp(clnOR.se)
cARisk.p <- ((sa/sN1) - (sc/sN0)) * units
cARisk.se <- (sqrt(((sa * (sN1 - sa))/sN1^3) + ((sc *
(sN0 - sc))/sN0^3))) * units
cARisk.l <- cARisk.p - (z * cARisk.se)
cARisk.u <- cARisk.p + (z * cARisk.se)
cARate.p <- ((sa/sb) - (sc/sd)) * units
cARate.se <- (sqrt((sa/sb^2) + (sc/sd^2))) * units
cARate.l <- cARate.p - (z * cARate.se)
cARate.u <- cARate.p + (z * cARate.se)
cAFRisk.p <- (cRR.p - 1)/cRR.p
cAFRisk.l <- min((cRR.l - 1)/cRR.l, (cRR.u - 1)/cRR.u)
cAFRisk.u <- max((cRR.l - 1)/cRR.l, (cRR.u - 1)/cRR.u)
cAFRate.p <- (cIRR.p - 1)/cIRR.p
cAFRate.l <- min((cIRR.l - 1)/cIRR.l, (cIRR.u - 1)/cIRR.u)
cAFRate.u <- max((cIRR.l - 1)/cIRR.l, (cIRR.u - 1)/cIRR.u)
cAFest.p <- (cOR.p - 1)/cOR.p
cAFest.l <- min((cOR.l - 1)/cOR.l, (cOR.u - 1)/cOR.u)
cAFest.u <- max((cOR.l - 1)/cOR.l, (cOR.u - 1)/cOR.u)
cPARisk.p <- ((sM1/stotal) - (sc/sN0)) * units
cPARisk.se <- (sqrt(((sM1 * (stotal - sM1))/stotal^3) +
((sc * (sN0 - sc))/sN0^3))) * units
cPARisk.l <- cPARisk.p - (z * cPARisk.se)
cPARisk.u <- cPARisk.p + (z * cPARisk.se)
cPARate.p <- ((sM1/sM0) - (sc/sd)) * units
cPARate.se <- (sqrt((sM1/sM0^2) + (sc/sd^2))) * units
cPARate.l <- cPARate.p - (z * cPARate.se)
cPARate.u <- cPARate.p + (z * cPARate.se)
cPAFRisk.p <- (cIRiskpop.p - cIRisko.p)/cIRiskpop.p
cPAFRisk.l <- min((cIRiskpop.l - cIRisko.l)/cIRiskpop.l,
(cIRiskpop.u - cIRisko.u)/cIRiskpop.u)
cPAFRisk.u <- max((cIRiskpop.l - cIRisko.l)/cIRiskpop.l,
(cIRiskpop.u - cIRisko.u)/cIRiskpop.u)
cPAFRate.p <- ((cIRR.p - 1)/cIRR.p) * (sa/sM1)
cPAFRate.l <- ((cIRR.p - 1)/cIRR.p) * (sa/sM1)
cPAFRate.u <- ((cIRR.p - 1)/cIRR.p) * (sa/sM1)
cPAFRate.p <- (cIRatepop.p - cIRateo.p)/cIRatepop.p
cPAFRate.l <- min((cIRatepop.l - cIRateo.l)/cIRatepop.l,
(cIRatepop.u - cIRateo.u)/cIRatepop.u)
cPAFRate.u <- max((cIRatepop.l - cIRateo.l)/cIRatepop.l,
(cIRatepop.u - cIRateo.u)/cIRatepop.u)
cPAFest.p <- (cOpop.p - cOo.p)/cOpop.p
cPAFest.l <- min((cOpop.l - cOo.l)/cOpop.l, (cOpop.u -
cOo.u)/cOpop.u)
cPAFest.u <- max((cOpop.l - cOo.l)/cOpop.l, (cOpop.u -
cOo.u)/cOpop.u)
exp.a <- (N1 * M1)/total
exp.b <- (N1 * M0)/total
exp.c <- (N0 * M1)/total
exp.d <- (N0 * M0)/total
chi2 <- (((a - exp.a)^2)/exp.a) + (((b - exp.b)^2)/exp.b) +
(((c - exp.c)^2)/exp.c) + (((d - exp.d)^2)/exp.d)
p.chi2 <- 1 - pchisq(chi2, df = 1)
exp.sa <- (sN1 * sM1)/stotal
exp.sb <- (sN1 * sM0)/stotal
exp.sc <- (sN0 * sM1)/stotal
exp.sd <- (sN0 * sM0)/stotal
chi2s <- (((sa - exp.sa)^2)/exp.sa) + (((sb - exp.sb)^2)/exp.sb) +
(((sc - exp.sc)^2)/exp.sc) + (((sd - exp.sd)^2)/exp.sd)
p.chi2s <- 1 - pchisq(chi2s, df = 1)
if (length(a) > 1) {
chi2m <- as.numeric(mantelhaen.test(dat)$statistic)
p.chi2m <- as.numeric(mantelhaen.test(dat)$p.value)
}
sRR.p <- sum((a * N0/total))/sum((c * N1/total))
varLNRR.s <- sum(((M1 * N1 * N0)/total^2) - ((a * c)/total))/(sum((a *
N0)/total) * sum((c * N1)/total))
lnRR.s <- log(sRR.p)
sRR.se <- (sqrt(varLNRR.s))
sRR.l <- exp(lnRR.s - (z * sqrt(varLNRR.s)))
sRR.u <- exp(lnRR.s + (z * sqrt(varLNRR.s)))
sIRR.p <- sum((a * d)/M0)/sum((c * b)/M0)
lnIRR.s <- log(sIRR.p)
varLNIRR.s <- (sum((M1 * b * d)/M0^2))/(sum((a * d)/M0) *
sum((c * b)/M0))
sIRR.se <- sqrt(varLNIRR.s)
sIRR.l <- exp(lnIRR.s - (z * sqrt(varLNIRR.s)))
sIRR.u <- exp(lnIRR.s + (z * sqrt(varLNIRR.s)))
sOR.p <- sum((a * d/total))/sum((b * c/total))
G <- a * d/total
H <- b * c/total
P <- (a + d)/total
Q <- (b + c)/total
GQ.HP <- G * Q + H * P
sumG <- sum(G)
sumH <- sum(H)
sumGP <- sum(G * P)
sumGH <- sum(G * H)
sumHQ <- sum(H * Q)
sumGQ <- sum(G * Q)
sumGQ.HP <- sum(GQ.HP)
varLNOR.s <- sumGP/(2 * sumG^2) + sumGQ.HP/(2 * sumG *
sumH) + sumHQ/(2 * sumH^2)
lnOR.s <- log(sOR.p)
sOR.se <- sqrt(varLNOR.s)
sOR.l <- exp(lnOR.s - z * sqrt(varLNOR.s))
sOR.u <- exp(lnOR.s + z * sqrt(varLNOR.s))
sARisk.p <- (sum(((a * N0) - (c * N1))/total)/sum((N1 *
N0)/total)) * units
w <- (N1 * N0)/total
var.p1 <- (((a * d)/(N1^2 * (N1 - 1))) + ((c * b)/(N0^2 *
(N0 - 1))))
var.p1[N0 == 1] <- 0
var.p1[N1 == 1] <- 0
varARisk.s <- sum(w^2 * var.p1)/sum(w)^2
sARisk.se <- (sqrt(varARisk.s)) * units
sARisk.l <- sARisk.p - (z * sARisk.se)
sARisk.u <- sARisk.p + (z * sARisk.se)
sARate.p <- sum(((a * d) - (c * b))/M0)/sum((b * d)/M0) *
units
varARate.s <- sum(((b * d)/M0)^2 * ((a/b^2) + (c/d^2)))/sum((b *
d)/M0)^2
sARate.se <- sqrt(varARate.s) * units
sARate.l <- sARate.p - (z * sARate.se)
sARate.u <- sARate.p + (z * sARate.se)
RR.conf.p <- (cRR.p/sRR.p)
RR.conf.l <- (cRR.l/sRR.l)
RR.conf.u <- (cRR.u/sRR.u)
IRR.conf.p <- (cIRR.p/sIRR.p)
IRR.conf.l <- (cIRR.l/sIRR.l)
IRR.conf.u <- (cIRR.u/sIRR.u)
OR.conf.p <- (cOR.p/sOR.p)
OR.conf.l <- (cOR.l/sOR.l)
OR.conf.u <- (cOR.u/sOR.u)
ARisk.conf.p <- (cARisk.p/sARisk.p)
ARisk.conf.l <- (cARisk.l/sARisk.l)
ARisk.conf.u <- (cARisk.u/sARisk.u)
ARate.conf.p <- (cARate.p/sARate.p)
ARate.conf.l <- (cARate.l/sARate.l)
ARate.conf.u <- (cARate.u/sARate.u)
if (length(a) > 1) {
if (homogeneity == "woolf") {
lnRR. <- log((a/(a + b))/(c/(c + d)))
lnRR.var. <- (b/(a * (a + b))) + (d/(c * (c +
d)))
wRR. <- 1/lnRR.var.
lnRR.s. <- sum(wRR. * lnRR.)/sum(wRR.)
RR.homogeneity <- sum(wRR. * (lnRR. - lnRR.s.)^2)
RR.homogeneity.p <- 1 - pchisq(RR.homogeneity,
df = n.strata - 1)
RR.homog <- data.frame(test.statistic = RR.homogeneity,
df = n.strata - 1, p.value = RR.homogeneity.p)
lnOR. <- log(((a + 0.5) * (d + 0.5))/((b + 0.5) *
(c + 0.5)))
lnOR.var. <- (1/(a + 0.5)) + (1/(b + 0.5)) +
(1/(c + 0.5)) + (1/(d + 0.5))
wOR. <- 1/lnOR.var.
lnOR.s. <- sum((wOR. * lnOR.))/sum(wOR.)
OR.homogeneity <- sum(wOR. * (lnOR. - lnOR.s.)^2)
OR.homogeneity.p <- 1 - pchisq(OR.homogeneity,
df = n.strata - 1)
OR.homog <- data.frame(test.statistic = OR.homogeneity,
df = n.strata - 1, p.value = OR.homogeneity.p)
}
if (homogeneity == "breslow.day") {
n11k = dat[1, 1, ]
n21k = dat[2, 1, ]
n12k = dat[1, 2, ]
n22k = dat[2, 2, ]
row1sums = n11k + n12k
row2sums = n21k + n22k
col1sums = n11k + n21k
Amax = apply(cbind(row1sums, col1sums), 1, min)
bb = row2sums + row1sums * sRR.p - col1sums *
(1 - sRR.p)
determ = sqrt(bb^2 + 4 * (1 - sRR.p) * sRR.p *
row1sums * col1sums)
Astar = (-bb + cbind(-determ, determ))/(2 - 2 *
sRR.p)
Astar = ifelse(Astar[, 1] <= Amax & Astar[, 1] >=
0, Astar[, 1], Astar[, 2])
Bstar = row1sums - Astar
Cstar = col1sums - Astar
Dstar = row2sums - col1sums + Astar
Var = apply(1/cbind(Astar, Bstar, Cstar, Dstar),
1, sum)^(-1)
RR.homogeneity = sum((dat[1, 1, ] - Astar)^2/Var)
RR.homogeneity.p = 1 - pchisq(RR.homogeneity,
df = n.strata - 1)
RR.homog <- data.frame(test.statistic = RR.homogeneity,
df = n.strata - 1, p.value = RR.homogeneity.p)
bb = row2sums + row1sums * sOR.p - col1sums *
(1 - sOR.p)
determ = sqrt(bb^2 + 4 * (1 - sOR.p) * sOR.p *
row1sums * col1sums)
Astar = (-bb + cbind(-determ, determ))/(2 - 2 *
sOR.p)
Astar = ifelse(Astar[, 1] <= Amax & Astar[, 1] >=
0, Astar[, 1], Astar[, 2])
Bstar = row1sums - Astar
Cstar = col1sums - Astar
Dstar = row2sums - col1sums + Astar
Var = apply(1/cbind(Astar, Bstar, Cstar, Dstar),
1, sum)^(-1)
OR.homogeneity = sum((dat[1, 1, ] - Astar)^2/Var)
OR.homogeneity.p = 1 - pchisq(OR.homogeneity,
df = n.strata - 1)
OR.homog <- data.frame(test.statistic = OR.homogeneity,
df = n.strata - 1, p.value = OR.homogeneity.p)
}
}
})
res <- list(RR.strata = with(elements, data.frame(est = RR.p,
se = RR.se, weight = RR.w, lower = RR.l, upper = RR.u)),
IRR.strata = with(elements, data.frame(est = IRR.p, se = IRR.se,
weight = IRR.w, lower = IRR.l, upper = IRR.u)), OR.strata = with(elements,
data.frame(est = OR.p, se = OR.se, weight = OR.w,
lower = OR.l, upper = OR.u)), cRR.strata = with(elements,
data.frame(est = cRR.p, lower = cRR.l, upper = cRR.u)),
ARisk.strata = with(elements, data.frame(est = ARisk.p,
se = ARisk.se, weight = ARisk.w, lower = ARisk.l,
upper = ARisk.u)), ARate.strata = with(elements,
data.frame(est = ARate.p, se = ARate.se, lower = ARate.l,
upper = ARate.u)), AFRisk.strata = with(elements,
data.frame(est = AFRisk.p, lower = AFRisk.l, upper = AFRisk.u)),
AFRate.strata = with(elements, data.frame(est = AFRate.p,
lower = AFRate.l, upper = AFRate.u)), AFest.strata = with(elements,
data.frame(est = AFest.p, lower = AFest.l, upper = AFest.u)),
PARisk.strata = with(elements, data.frame(est = PARisk.p,
se = PARisk.se, lower = PARisk.l, upper = PARisk.u)),
PARate.strata = with(elements, data.frame(est = PARate.p,
se = PARate.se, lower = PARate.l, upper = PARate.u)),
PAFRisk.strata = with(elements, data.frame(est = PAFRisk.p,
lower = PAFRisk.l, upper = PAFRisk.u)), PAFRate.strata = with(elements,
data.frame(est = PAFRate.p, lower = PAFRate.l, upper = PAFRate.u)),
PAFest.strata = with(elements, data.frame(est = PAFest.p,
lower = PAFest.l, upper = PAFest.u)), RR.crude = with(elements,
data.frame(est = cRR.p, se = cRR.se, lower = cRR.l,
upper = cRR.u)), IRR.crude = with(elements, data.frame(est = cIRR.p,
se = cIRR.se, lower = cIRR.l, upper = cIRR.u)), OR.crude = with(elements,
data.frame(est = cOR.p, se = cOR.se, lower = cOR.l,
upper = cOR.u)), ARisk.crude = with(elements,
data.frame(est = cARisk.p, se = cARisk.se, lower = cARisk.l,
upper = cARisk.u)), ARate.crude = with(elements,
data.frame(est = cARate.p, se = cARate.se, lower = cARate.l,
upper = cARate.u)), AFRisk.crude = with(elements,
data.frame(est = cAFRisk.p, lower = cAFRisk.l, upper = cAFRisk.u)),
AFRate.crude = with(elements, data.frame(est = cAFRate.p,
lower = cAFRate.l, upper = cAFRate.u)), AFest.crude = with(elements,
data.frame(est = cAFest.p, lower = cAFest.l, upper = cAFest.u)),
PARisk.crude = with(elements, data.frame(est = cPARisk.p,
se = cPARisk.se, lower = cPARisk.l, upper = cPARisk.u)),
PARate.crude = with(elements, data.frame(est = cPARate.p,
se = cPARate.se, lower = cPARate.l, upper = cPARate.u)),
PAFRisk.crude = with(elements, data.frame(est = cPAFRisk.p,
lower = cPAFRisk.l, upper = cPAFRisk.u)), PAFRate.crude = with(elements,
data.frame(est = cPAFRate.p, lower = cPAFRate.l,
upper = cPAFRate.u)), PAFest.crude = with(elements,
data.frame(est = cPAFest.p, lower = cPAFest.l, upper = cPAFest.u)),
RR.mh = with(elements, data.frame(est = sRR.p, se = sRR.se,
lower = sRR.l, upper = sRR.u)), IRR.mh = with(elements,
data.frame(est = sIRR.p, se = sIRR.se, lower = sIRR.l,
upper = sIRR.u)), OR.mh = with(elements, data.frame(est = sOR.p,
se = sOR.se, lower = sOR.l, upper = sOR.u)), ARisk.mh = with(elements,
data.frame(est = sARisk.p, se = sARisk.se, lower = sARisk.l,
upper = sARisk.u)), ARate.mh = with(elements,
data.frame(est = sARate.p, se = sARate.se, lower = sARate.l,
upper = sARate.u)), RR.conf = with(elements,
data.frame(est = RR.conf.p, lower = RR.conf.l, upper = RR.conf.u)),
IRR.conf = with(elements, data.frame(est = IRR.conf.p,
lower = IRR.conf.l, upper = IRR.conf.u)), OR.conf = with(elements,
data.frame(est = OR.conf.p, lower = OR.conf.l, upper = OR.conf.u)),
ARisk.conf = with(elements, data.frame(est = ARisk.conf.p,
lower = ARisk.conf.l, upper = ARisk.conf.u)), ARate.conf = with(elements,
data.frame(est = ARate.conf.p, lower = ARate.conf.l,
upper = ARate.conf.u)), chisq.strata = with(elements,
data.frame(test.statistic = chi2, df = 1, p.value = p.chi2)),
chisq.crude = with(elements, data.frame(test.statistic = chi2s,
df = 1, p.value = p.chi2s)), count.units = ifelse(units ==
1, "Cases per population unit", paste("Cases per ",
units, " population units", sep = "")), time.units = ifelse(units ==
1, "Cases per unit of population time at risk", paste("Cases per ",
units, " units of population time at risk", sep = "")))
if (length(dim(dat)) > 2) {
res$chisq.mh <- with(elements, data.frame(test.statistic = chi2m,
df = 1, p.value = p.chi2m))
}
if (method == "cohort.count" & length(elements$a) == 1) {
rval <- list(RR = res$RR.strata, OR = res$OR.strata,
AR = res$ARisk.strata, ARp = res$PARisk.strata, AFe = res$AFRisk.strata,
AFp = res$PAFRisk.strata, chisq = res$chisq.strata)
r1 <- with(elements, c(a, b, N1, cIRiske.p, cOe.p))
r2 <- with(elements, c(c, d, N0, cIRisko.p, cOo.p))
r3 <- with(elements, c(M1, M0, M0 + M1, cIRiskpop.p,
cOpop.p))
tab <- as.data.frame(rbind(r1, r2, r3))
colnames(tab) <- c(" Disease +", " Disease -", " Total",
" Inc risk *", " Odds")
rownames(tab) <- c("Exposed +", "Exposed -", "Total")
tab <- format.data.frame(tab, digits = 3, justify = "right")
out <- list(method = "cohort.count", conf.level = conf.level,
elements = elements, res = res, rval = rval, tab = tab)
}
if (method == "cohort.count" & length(elements$a) > 1) {
rval <- list(RR.strata = res$RR.strata, RR.crude = res$RR.crude,
RR.mh = res$RR.mh, OR.strata = res$OR.strata, OR.crude = res$OR.crude,
OR.mh = res$OR.mh, AR.strata = res$ARisk.strata,
AR.crude = res$ARisk.crude, AR.mh = res$ARisk.mh,
ARp.strata = res$PARisk.strata, AFe.strata = res$AFRisk.strata,
AFp.strata = res$PAFRisk.strata, chisq.strata = res$chisq.strata,
chisq.crude = res$chisq.crude, chisq.mh = res$chisq.mh,
RR.homog = elements$RR.homog, OR.homog = elements$OR.homog)
r1 <- with(elements, c(sa, sb, sN1, cIRiske.p, cOe.p))
r2 <- with(elements, c(sc, sd, sN0, cIRisko.p, cOo.p))
r3 <- with(elements, c(sM1, sM0, sM0 + sM1, cIRiskpop.p,
cOpop.p))
tab <- as.data.frame(rbind(r1, r2, r3))
colnames(tab) <- c(" Disease +", " Disease -", " Total",
" Inc risk *", " Odds")
rownames(tab) <- c("Exposed +", "Exposed -", "Total")
tab <- format.data.frame(tab, digits = 3, justify = "right")
out <- list(method = "cohort.count", conf.level = conf.level,
elements = elements, res = res, rval = rval, tab = tab)
}
if (method == "cohort.time" & length(elements$a) == 1) {
rval <- list(IRR = res$IRR.strata, AR = res$ARate.strata,
ARp = res$PARate.strata, AFe = res$AFRate.strata,
AFp = res$PAFRate.strata, chisq = res$chisq.strata)
r1 <- with(elements, c(a, b, cIRatee.p))
r2 <- with(elements, c(c, d, cIRateo.p))
r3 <- with(elements, c(M1, M0, cIRatepop.p))
tab <- as.data.frame(rbind(r1, r2, r3))
colnames(tab) <- c(" Disease +", " Time at risk",
" Inc rate *")
rownames(tab) <- c("Exposed +", "Exposed -", "Total")
tab <- format.data.frame(tab, digits = 3, justify = "right")
out <- list(method = "cohort.time", conf.level = conf.level,
elements = elements, res = res, rval = rval, tab = tab)
}
if (method == "cohort.time" & length(elements$a) > 1) {
rval <- list(IRR.strata = res$IRR.strata, IRR.crude = res$IRR.crude,
IRR.mh = res$IRR.mh, AR.strata = res$ARate.strata,
AR.crude = res$ARate.crude, AR.mh = res$ARate.mh,
ARp.strata = res$PARate.strata, AFp.strata = res$PAFRate.strata,
chisq.strata = res$chisq.strata, chisq.crude = res$chisq.crude,
chisq.mh = res$chisq.mh)
r1 <- with(elements, c(sa, sb, cIRatee.p))
r2 <- with(elements, c(sc, sd, cIRateo.p))
r3 <- with(elements, c(sM1, sM0, cIRatepop.p))
tab <- as.data.frame(rbind(r1, r2, r3))
colnames(tab) <- c(" Disease +", " Time at risk",
" Inc rate *")
rownames(tab) <- c("Exposed +", "Exposed -", "Total")
tab <- format.data.frame(tab, digits = 3, justify = "right")
out <- list(method = "cohort.time", conf.level = conf.level,
elements = elements, res = res, rval = rval, tab = tab)
}
if (method == "case.control" & length(elements$a) == 1) {
rval <- list(OR = res$OR.strata, AR = res$ARisk.strata,
ARp = res$PARisk.strata, AFest = res$AFest.strata,
AFp = res$PAFest.strata, chisq = res$chisq.strata)
r1 <- with(elements, c(a, b, N1, cIRiske.p, cOe.p))
r2 <- with(elements, c(c, d, N0, cIRisko.p, cOo.p))
r3 <- with(elements, c(M1, M0, M0 + M1, cIRiskpop.p,
cOpop.p))
tab <- as.data.frame(rbind(r1, r2, r3))
colnames(tab) <- c(" Disease +", " Disease -", " Total",
" Prevalence *", " Odds")
rownames(tab) <- c("Exposed +", "Exposed -", "Total")
tab <- format.data.frame(tab, digits = 3, justify = "right")
out <- list(method = "case.control", conf.level = conf.level,
elements = elements, res = res, rval = rval, tab = tab)
}
if (method == "case.control" & length(elements$a) > 1) {
rval <- list(OR.strata = res$OR.strata, OR.crude = res$OR.crude,
OR.mh = res$OR.mh, AR.strata = res$ARisk.strata,
AR.crude = res$ARisk.crude, AR.mh = res$ARisk.mh,
ARp.strata = res$PARisk.strata, AFest.strata = res$AFest.strata,
AFpest.strata = res$PAFest.strata, chisq.strata = res$chisq.strata,
chisq.crude = res$chisq.crude, chisq.mh = res$chisq.mh,
OR.homog = elements$OR.homog)
r1 <- with(elements, c(sa, sb, sN1, cIRiske.p, cOe.p))
r2 <- with(elements, c(sc, sd, sN0, cIRisko.p, cOo.p))
r3 <- with(elements, c(sM1, sM0, sM0 + sM1, cIRiskpop.p,
cOpop.p))
tab <- as.data.frame(rbind(r1, r2, r3))
colnames(tab) <- c(" Disease +", " Disease -", " Total",
" Prevalence *", " Odds")
rownames(tab) <- c("Exposed +", "Exposed -", "Total")
tab <- format.data.frame(tab, digits = 3, justify = "right")
out <- list(method = "case.control", conf.level = conf.level,
elements = elements, res = res, rval = rval, tab = tab)
}
if (method == "cross.sectional" & length(elements$a) == 1) {
rval <- list(PR = res$RR.strata, OR = res$OR.strata,
AR = res$ARisk.strata, ARp = res$PARisk.strata, AFe = res$AFRisk.strata,
AFp = res$PAFRisk.strata, chisq = res$chisq.strata)
r1 <- with(elements, c(a, b, N1, cIRiske.p, cOe.p))
r2 <- with(elements, c(c, d, N0, cIRisko.p, cOo.p))
r3 <- with(elements, c(M1, M0, M0 + M1, cIRiskpop.p,
cOpop.p))
tab <- as.data.frame(rbind(r1, r2, r3))
colnames(tab) <- c(" Disease +", " Disease -", " Total",
" Prevalence *", " Odds")
rownames(tab) <- c("Exposed +", "Exposed -", "Total")
tab <- format.data.frame(tab, digits = 3, justify = "right")
out <- list(method = "cross.sectional", conf.level = conf.level,
elements = elements, res = res, rval = rval, tab = tab)
}
if (method == "cross.sectional" & length(elements$a) > 1) {
rval <- list(PR.strata = res$RR.strata, PR.crude = res$RR.crude,
PR.mh = res$RR.mh, OR.strata = res$OR.strata, OR.crude = res$OR.crude,
OR.mh = res$OR.mh, AR.strata = res$ARisk.strata,
AR.crude = res$ARisk.crude, AR.mh = res$ARisk.mh,
ARp.strata = res$PARisk.strata, AFe.strata = res$AFRisk.strata,
AFp.strata = res$PAFRisk.strata, chisq.strata = res$chisq.strata,
chisq.crude = res$chisq.crude, chisq.mh = res$chisq.mh,
PR.homog = elements$RR.homog, OR.homog = elements$OR.homog)
r1 <- with(elements, c(sa, sb, sN1, cIRiske.p, cOe.p))
r2 <- with(elements, c(sc, sd, sN0, cIRisko.p, cOo.p))
r3 <- with(elements, c(sM1, sM0, sM1 + sM0, cIRiskpop.p,
cOpop.p))
tab <- as.data.frame(rbind(r1, r2, r3))
colnames(tab) <- c(" Disease +", " Disease -", " Total",
" Prevalence *", " Odds")
rownames(tab) <- c("Exposed +", "Exposed -", "Total")
tab <- format.data.frame(tab, digits = 3, justify = "right")
out <- list(method = "cross.sectional", conf.level = conf.level,
elements = elements, res = res, rval = rval, tab = tab)
}
class(out) <- "epi.2by2"
return(out)
}
## Print method
print.epi.2by2
function (obj)
{
if (obj$method == "cohort.count" & length(obj$elements$a) ==
1) {
print(obj$tab)
cat("\nPoint estimates and", obj$conf.level * 100, "%",
"CIs:")
cat("\n---------------------------------------------------------")
with(obj$res, {
cat(sprintf("\nInc risk ratio %.2f (%.2f, %.2f)",
RR.crude$est, RR.crude$lower, RR.crude$upper))
cat(sprintf("\nOdds ratio %.2f (%.2f, %.2f)",
OR.crude$est, OR.crude$lower, OR.crude$upper))
cat(sprintf("\nAttrib risk * %.2f (%.2f, %.2f)",
ARisk.crude$est, ARisk.crude$lower, ARisk.crude$upper))
cat(sprintf("\nAttrib risk in population * %.2f (%.2f, %.2f)",
PARisk.crude$est, PARisk.crude$lower, PARisk.crude$upper))
cat(sprintf("\nAttrib fraction in exposed (%%) %.2f (%.2f, %.2f)",
AFRisk.crude$est * 100, AFRisk.crude$lower *
100, AFRisk.crude$upper * 100))
cat(sprintf("\nAttrib fraction in population (%%) %.2f (%.2f, %.2f)",
PAFRisk.strata$est * 100, PAFRisk.strata$lower *
100, PAFRisk.strata$upper * 100))
})
cat("\n---------------------------------------------------------")
cat("\n", "*", obj$res$count.units, "\n")
}
if (obj$method == "cohort.count" & length(obj$elements$a) >
1) {
print(obj$tab)
cat("\n")
cat("\nPoint estimates and", obj$conf.level * 100, "%",
"CIs:")
cat("\n---------------------------------------------------------")
with(obj$res, {
cat(sprintf("\nInc risk ratio (crude) %.2f (%.2f, %.2f)",
RR.crude$est, RR.crude$lower, RR.crude$upper))
cat(sprintf("\nInc risk ratio (M-H) %.2f (%.2f, %.2f)",
RR.mh$est, RR.mh$lower, RR.mh$upper))
cat(sprintf("\nInc risk ratio (crude:M-H) %.2f",
RR.conf$est))
cat(sprintf("\nOdds ratio (crude) %.2f (%.2f, %.2f)",
OR.crude$est, OR.crude$lower, OR.crude$upper))
cat(sprintf("\nOdds ratio (M-H) %.2f (%.2f, %.2f)",
OR.mh$est, OR.mh$lower, OR.mh$upper))
cat(sprintf("\nOdds ratio (crude:M-H) %.2f",
OR.conf$est))
cat(sprintf("\nAttrib risk (crude) * %.2f (%.2f, %.2f)",
ARisk.crude$est, ARisk.crude$lower, ARisk.crude$upper))
cat(sprintf("\nAttrib risk (M-H) * %.2f (%.2f, %.2f)",
ARisk.mh$est, ARisk.mh$lower, ARisk.mh$upper))
cat(sprintf("\nAttrib risk (crude:M-H) %.2f",
ARisk.conf$est))
})
cat("\n---------------------------------------------------------")
cat("\n", "*", obj$res$count.units, "\n")
}
if (obj$method == "cohort.time" & length(obj$elements$a) ==
1) {
print(obj$tab)
cat("\nPoint estimates and", obj$conf.level * 100, "%",
"CIs:")
cat("\n---------------------------------------------------------")
with(obj$res, {
cat(sprintf("\nInc rate ratio %.2f (%.2f, %.2f)",
IRR.crude$est, IRR.crude$lower, IRR.crude$upper))
cat(sprintf("\nAttrib rate * %.2f (%.2f, %.2f)",
ARate.crude$est, ARate.crude$lower, ARate.crude$upper))
cat(sprintf("\nAttrib rate in population * %.2f (%.2f, %.2f)",
PARate.crude$est, PARate.crude$lower, PARate.crude$upper))
cat(sprintf("\nAttrib fraction in exposed (%%) %.2f (%.2f, %.2f)",
AFRate.crude$est * 100, AFRate.crude$lower *
100, AFRate.crude$upper * 100))
cat(sprintf("\nAttrib fraction in population (%%) %.2f (%.2f, %.2f)",
PAFRate.crude$est * 100, PAFRate.crude$lower *
100, PAFRate.crude$upper * 100))
})
cat("\n---------------------------------------------------------")
cat("\n", "*", obj$res$time.units, "\n")
}
if (obj$method == "cohort.time" & length(obj$elements$a) >
1) {
print(obj$tab)
cat("\nPoint estimates and", obj$conf.level * 100, "%",
"CIs:")
cat("\n---------------------------------------------------------")
with(obj$res, {
cat(sprintf("\nInc rate ratio (crude) %.2f (%.2f, %.2f)",
IRR.crude$est, IRR.crude$lower, IRR.crude$upper))
cat(sprintf("\nInc rate ratio (M-H) %.2f (%.2f, %.2f)",
IRR.mh$est, IRR.mh$lower, IRR.mh$upper))
cat(sprintf("\nInc rate ratio (crude:M-H) %.2f",
IRR.conf$est))
cat(sprintf("\nAttrib rate (crude) * %.2f (%.2f, %.2f)",
ARate.crude$est, ARate.crude$lower, ARate.crude$upper))
cat(sprintf("\nAttrib rate (M-H) * %.2f (%.2f, %.2f)",
ARate.mh$est, ARate.mh$lower, ARate.mh$upper))
cat(sprintf("\nAttrib rate (crude:M-H) %.2f",
ARate.conf$est))
})
cat("\n---------------------------------------------------------")
cat("\n", "*", obj$res$time.units, "\n")
}
if (obj$method == "case.control" & length(obj$elements$a) ==
1) {
print(obj$tab)
cat("\nPoint estimates and", obj$conf.level * 100, "%",
"CIs:")
cat("\n---------------------------------------------------------")
with(obj$res, {
cat(sprintf("\nOdds ratio %.2f (%.2f, %.2f)",
OR.crude$est, OR.crude$lower, OR.crude$upper))
cat(sprintf("\nAttrib prevalence * %.2f (%.2f, %.2f)",
ARisk.crude$est, ARisk.crude$lower, ARisk.crude$upper))
cat(sprintf("\nAttrib prevalence in population * %.2f (%.2f, %.2f)",
PARisk.crude$est, PARisk.crude$lower, PARisk.crude$upper))
cat(sprintf("\nAttrib fraction (est) in exposed (%%) %.2f (%.2f, %.2f)",
AFest.crude$est * 100, AFest.crude$lower * 100,
AFest.crude$upper * 100))
cat(sprintf("\nAttrib fraction (est) in population (%%) %.2f (%.2f, %.2f)",
PAFest.strata$est * 100, PAFest.strata$lower *
100, PAFest.strata$upper * 100))
})
cat("\n---------------------------------------------------------")
cat("\n", "*", obj$res$count.units, "\n")
}
if (obj$method == "case.control" & length(obj$elements$a) >
1) {
print(obj$tab)
cat("\nPoint estimates and", obj$conf.level * 100, "%",
"CIs:")
cat("\n---------------------------------------------------------")
with(obj$res, {
cat(sprintf("\nOdds ratio (crude) %.2f (%.2f, %.2f)",
OR.crude$est, OR.crude$lower, OR.crude$upper))
cat(sprintf("\nOdds ratio (M-H) %.2f (%.2f, %.2f)",
OR.mh$est, OR.mh$lower, OR.mh$upper))
cat(sprintf("\nOdds ratio (crude:M-H) %.2f",
OR.conf$est))
cat(sprintf("\nAttrib prevalence (crude) * %.2f (%.2f, %.2f)",
ARisk.crude$est, ARisk.crude$lower, ARisk.crude$upper))
cat(sprintf("\nAttrib prevalence (M-H) * %.2f (%.2f, %.2f)",
ARisk.mh$est, ARisk.mh$lower, ARisk.mh$upper))
cat(sprintf("\nAttrib prevalence (crude:M-H) %.2f",
ARisk.conf$est))
})
cat("\n---------------------------------------------------------")
cat("\n", "*", obj$res$count.units, "\n")
}
if (obj$method == "cross.sectional" & length(obj$elements$a) ==
1) {
print(obj$tab)
cat("\nPoint estimates and", obj$conf.level * 100, "%",
"CIs:")
cat("\n---------------------------------------------------------")
with(obj$res, {
cat(sprintf("\nPrevalence ratio %.2f (%.2f, %.2f)",
RR.crude$est, RR.crude$lower, RR.crude$upper))
cat(sprintf("\nOdds ratio %.2f (%.2f, %.2f)",
OR.crude$est, OR.crude$lower, OR.crude$upper))
cat(sprintf("\nAttrib prevalence * %.2f (%.2f, %.2f)",
ARisk.crude$est, ARisk.crude$lower, ARisk.crude$upper))
cat(sprintf("\nAttrib prevalence in population * %.2f (%.2f, %.2f)",
PARisk.crude$est, PARisk.crude$lower, PARisk.crude$upper))
cat(sprintf("\nAttrib fraction in exposed (%%) %.2f (%.2f, %.2f)",
AFRisk.crude$est * 100, AFRisk.crude$lower *
100, AFRisk.crude$upper * 100))
cat(sprintf("\nAttrib fraction in population (%%) %.2f (%.2f, %.2f)",
PAFRisk.strata$est * 100, PAFRisk.strata$lower *
100, PAFRisk.strata$upper * 100))
})
cat("\n---------------------------------------------------------")
cat("\n", "*", obj$res$count.units, "\n")
}
if (obj$method == "cross.sectional" & length(obj$elements$a) >
1) {
print(obj$tab)
cat("\nPoint estimates and", obj$conf.level * 100, "%",
"CIs:")
cat("\n---------------------------------------------------------")
with(obj$res, {
cat(sprintf("\nPrevalence ratio (crude) %.2f (%.2f, %.2f)",
RR.crude$est, RR.crude$lower, RR.crude$upper))
cat(sprintf("\nPrevalence ratio (M-H) %.2f (%.2f, %.2f)",
RR.mh$est, RR.mh$lower, RR.mh$upper))
cat(sprintf("\nPrevalence ratio (crude:M-H) %.2f",
RR.conf$est))
cat(sprintf("\nOdds ratio (crude) %.2f (%.2f, %.2f)",
OR.crude$est, OR.crude$lower, OR.crude$upper))
cat(sprintf("\nOdds ratio (M-H) %.2f (%.2f, %.2f)",
OR.mh$est, OR.mh$lower, OR.mh$upper))
cat(sprintf("\nOdds ratio (crude:M-H) %.2f",
OR.conf$est))
cat(sprintf("\nAtributable prevalence (crude) * %.2f (%.2f, %.2f)",
ARisk.crude$est, ARisk.crude$lower, ARisk.crude$upper))
cat(sprintf("\nAtributable prevalence (M-H) * %.2f (%.2f, %.2f)",
ARisk.mh$est, ARisk.mh$lower, ARisk.mh$upper))
cat(sprintf("\nAtributable prevalence (crude:M-H) %.2f",
ARisk.conf$est))
})
cat("\n---------------------------------------------------------")
cat("\n", "*", obj$res$count.units, "\n")
}
}
## Summary method
summary.epi.2by2
function (obj)
{
return(obj$rval)
}