B Código da função paired_proportions
Abaixo é exibido o código da função paired_proportions, utilizada no capítulo 17, seção 17.3.3.
library(dplyr)
paired_proportions <- function(data, id, row, col, row_ref, row_trt, col_ref,
col_out, case_control = FALSE, alpha = 0.05, ...) {
id_subj99 = id
summary = mutate(data,
cell = case_when(
data[, row] == row_ref & data[, col] == col_ref ~ 1,
data[, row] == row_trt & data[, col] == col_ref ~ 2,
data[, row] == row_ref & data[, col] == col_out ~ 4,
data[, row] == row_trt & data[, col] == col_out ~ 8,
TRUE ~ 0)) %>% group_by(data[, id_subj99]) %>%
summarize(cell_class = sum(cell))
tab = table(summary$cell_class)
if (case_control == TRUE) {
#
# tab["5"] -> (row_ref in both),
# tab["6"] -> (row_ref in col_out and row_trt in col_ref),
# tab["9"] -> (row_trt in col_out and row_ref in col_ref),
# tab["10"] -> (row_trt in both),
#
if (is.na(tab["5"])) tab["5"] = 0
if (is.na(tab["6"])) tab["6"] = 0
if (is.na(tab["9"])) tab["9"] = 0
if (is.na(tab["10"])) tab["10"] = 0
mat = matrix(c(tab["10"], tab["6"], tab["9"], tab["5"]), nrow = 2)
colnames(mat) = c(row_trt, row_ref)
rownames(mat) = c(row_trt, row_ref)
names(dimnames(mat)) = c(col_out,col_ref)
} else {
#
# tab["3"] -> (col_ref in both),
# tab["6"] -> (col_ref in row_trt and col_out in row_ref),
# tab["9"] -> (col_out in row_trt and col_ref in row_ref),
# tab["12"] -> (col_out in both),
#
if (is.na(tab["3"])) tab["3"] = 0
if (is.na(tab["6"])) tab["6"] = 0
if (is.na(tab["9"])) tab["9"] = 0
if (is.na(tab["12"])) tab["12"] = 0
mat = matrix(c(tab["12"], tab["6"], tab["9"], tab["3"]), nrow = 2)
colnames(mat) = c(col_out, col_ref)
rownames(mat) = c(col_out, col_ref)
names(dimnames(mat)) = c(row_trt, row_ref)
}
#
# IC proportion differences - Wald with Bonett–Price Laplace adjustment
#
n = sum(mat)
r = mat[1,1]
t = mat[2,1]
s = mat[1,2]
n1. = r+s
n.1 = r+t
p2 = n1. / n
p1 = n.1 / n
p2_ = (s + 1) / (n + 2)
p1_ = (t + 1) / (n + 2)
z = qnorm(1-alpha/2)
D = p2 - p1
Dinf = max(-1, (p2_-p1_) - z * sqrt((p2_+p1_-(p2_-p1_)^2)/(n+2)))
Dsup = min((p2_-p1_) + z * sqrt((p2_+p1_-(p2_-p1_)^2)/(n+2)), 1)
prop_diff = c(D, Dinf, Dsup)
#
# IC odds ratio - Transformed Wilson score
#
or_est = s/t
nd = t+s
div = 2*(nd+z^2)
num1 = (2*s + z^2)
num2 = z*sqrt(z^2 + 4*s*(1-s/nd))
pu = (num1 + num2) / div
pi = (num1 - num2) / div
or = c(or_est, pi/(1-pi), pu/(1-pu))
#
# IC relative risk - Bonett–Price hybrid Wilson score
#
rr_est = p2/p1
nn = r + t + s
A = sqrt((nd + 2) / ((n1.+1)*(n.1+1)))
B = sqrt((1-(n1.+1)/(nn+2))/(n1.+1))
C = sqrt((1-(n.1+1)/(nn+2))/(n.1+1))
z = A/(B+C) * z
den =2*(nn+z^2)
l1 = 2*n1.+z^2 - z*sqrt(z^2+4*n1.*(1-n1./nn))
u1 = 2*n1.+z^2 + z*sqrt(z^2+4*n1.*(1-n1./nn))
l2 = 2*n.1+z^2 - z*sqrt(z^2+4*n.1*(1-n.1/nn))
u2 = 2*n.1+z^2 + z*sqrt(z^2+4*n.1*(1-n.1/nn))
rr = c(rr_est, l1/u2, u1/l2)
print(mat)
print(mcnemar.test(mat, ...))
text_low = paste("Lower ", sprintf("%.0f", (1-alpha)*100),"% CI",
sep = '')
text_upper = paste("Upper ", sprintf("%.0f", (1-alpha)*100),"% CI",
sep = '')
if (case_control == FALSE) {
names(prop_diff) = c("proportion differences", text_low, text_upper)
print(prop_diff)
cat("", sep = "\n")
names(rr) = c("relative risk", text_low, text_upper)
print(rr)
cat("", sep = "\n")
}
names(or) = c("odds ratio", text_low,text_upper)
print(or)
cat("", sep="\n")
}