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")
}