R/mhodds.R

In myominnoo/mStats_beta: A tool for data management and statistical analysis

#' @title Tabulation to display Odds Ratios
#' @description
#' \code{mhodds} generates cross-tabulation between two variables and
#' display odds as well as Odds Ratios of failure \code{var_case}
#' against a categorical explanatory variable \code{var_exp}.
#'
#' It is used with cross-sectional as well as case-control data.
#' @param var_case Case or outcome variable should be binary vector.
#' @param var_exp Exposure variable.
#' @param data a data frame object (Optional)
#' @param exposed a number or character representing reference value
#' @param strata variable to stratify by Mantel-Haenszel method
#' @param na.rm A logical value to specify missing values, <NA> in the table
#' @param rnd specify rounding of numbers. See \code{\link{round}}.
#' @param print.table logical value to display formatted outputs
#' @details
#' The variable \code{var_case} should coded 1 for case and 0 for non-case.
#' The arugment \code{exposed} can be used to set the reference value for
#' comparison of odds ratios.
#'
#' It produces a simple table with odds, Odds Ratio,95% CI as well as
#' p-value. If \code{strata} is specified, \code{Mantel-Haenzsel\'s Odds
#' Ratio} is generated along with chi-sqaured test for heterogeneity.
#'
#' \strong{Odds Ratio}
#'
#' \deqn{p(x) = \frac{\lambda\^{x} e^{-\lambda}}{x!}}{%
#' p(x) = \lambda^{x} exp(-\lambda)/x!}
#' for \eqn{x = 0, 1, 2, \ldots}.
#'
#' \deqn{\square_{i=1}^n X_i$$}
#'
#' \strong{Error Factor}
#'
#' is calcluated using the following formula:
#' \deqn{EF = exp(log(Odds Ratio) \pm{} 1.96 * SE}
#'
#' \deqn{SE = \sqrt{1/D1 + 1/D0 + 1/H1 + 1/H0}}
#'
#' \strong{References:}
#' \enumerate{
#'   \item Essential Medical Statistics, Betty R. Kirkwood & Jonathan A.C. Sterne,
#'   Second Edition. Chapter 3
#' }
#' @seealso \code{\link{tabodds}}, \code{\link{expandTables}}
#' @keywords odds, odds ratio, frequency table, statistics, descriptive
#' @author Myo Minn Oo (Email: \email{dr.myominnoo@@gmail.com} |
#' Website: \url{https://myominnoo.github.io/})
#' @examples
#' \dontrun{
#' ## Create stratified 2x2 tables from page 178 (Essential Medical
#' Statistics)
#' ## See ?expandTables
#' ab <- expandTables(male = c(36, 14, 50, 50),
#'     female = c(24, 126, 10, 90),
#'     exp_name = "areaType", exp_lvl = c("Rural", "Urban"),
#'     case_name = "Antibodies", case_lvl = c("Yes", "No"),
#'     strata_name = "gender")
#'
#' ## Labelling data and variables
#' ab <- labelData(ab, "Prevalence survey of leptospirosis in West Indies")
#' ab <- labelVars(ab, areaType, "Type of area in west indies")
#' ab <- labelVars(ab, Antibodies, "Leptospirosis antibodies")
#' ab <- labelVars(ab, gender, "Male or Female")
#'
#' ## Claculating Odds Ratios
#' mhodds(areaType, Antibodies, ab)
#' mhodds(areaType, Antibodies, ab, exposed = "Urban")
#'
#' ## Claculating Mantel-Haenzsel Odds Ratio
#' mhodds(areaType, Antibodies, ab, exposed = "Urban",
#'     strata = gender)
#'
#' ## Breast Cancer example (page 162, Essential Medical Statistics)
#' breast <- expandTables(c(88, 12, 71, 29),
#'                        exp_name = "drug",
#'                        exp_lvl = c("new", "standard"),
#'                        case_name = "nausea",
#'                        case_lvl = c("severe", "none"))
#'
#' tabodds(drug, nausea, breast)
#' mhodds(drug, nausea, breast, exposed = "standard")
#'
#'
#' ## Asthma Example (page 160, Essential Medical Statistics)
#' asthma <- expandTables(c(81, 995, 57, 867),
#'              exp_name = "Sex",
#'              exp_lvl = c("Women", "Man"),
#'              case_name = "Asthma",
#'              case_lvl = c("Yes", "No"))
#'
#' ## labelling
#' asthma <- labelData(asthma, "Hypothetical Data of Asthma Prevalence")
#' asthma <- labelVars(asthma,
#'     c(Sex, Asthma), c("Man or Woman", "Asthma or No Asthma"))
#'
#' ## Checking codebook
#' codebook(asthma)
#'
#' ## Odds
#' tabodds(Sex, Asthma, asthma)
#'
#' ## Odds ratios
#' mhodds(Sex, Asthma, asthma)
#' mhodds(Sex, Asthma, asthma, "Women")
#'
#' ## The odds ratio, J Martin Bland, Douglas G Altman, BMJ 2000;320:1468
#' hay <- expandTables(c(141, 420, 928, 13525),
#'                     exp_name = "eczema",
#'                     exp_lvl = c("Yes", "No"),
#'                     case_name = "hayFever",
#'                     case_lvl = c("Yes", "No"))
#' hay <- labelData(hay, "hay fever and eczema in 11 year old children")
#' hay <- labelVars(hay,
#'                  c(eczema, hayFever),
#'                 c("prevalence of eczema", "prevalence of hay fever"))
#'
#'
#' tabodds(eczema, hayFever, hay)
#' mhodds(eczema, hayFever, hay, "Yes")
#' }

#' @export
mhodds <- function(var_case, var_exp, data = NULL,
                   exposed = NULL, strata = NULL,
                   na.rm = FALSE, rnd = 3,
                   print.table = TRUE)
{
    arguments <- as.list(match.call())
    if (!is.null(data)) {
        var_case <- eval(substitute(var_case), data)
        var_exp <- eval(substitute(var_exp), data)
        strata <- eval(substitute(strata), data)
    }

    if (is.null(strata)) {
        x <- as.character()
    } else {
        x <- list()
    }

    UseMethod("mhodds", x)
}

#' @rdname mhodds
#' @export
mhodds.character <- function(var_case, var_exp, data = NULL,
                             exposed = NULL, strata = NULL,
                             na.rm = FALSE, rnd = 3,
                             print.table = TRUE)
{
    arguments <- as.list(match.call())
    exp_name <- (deparse(substitute(var_exp)))
    case_name <- (deparse(substitute(var_case)))

    if (!is.null(data)) {
        var_exp <- eval(substitute(var_exp), data)
        var_case <- eval(substitute(var_case), data)
    }

    na.rm <- ifelse(na.rm, "no", "ifany")

    t <- tabodds(var_case, var_exp, data = NULL,
                 plot = FALSE, print.table = FALSE)
    t <- t[, 1:3]
    lvl <- row.names(t)

    if (is.null(exposed)) {
        exposed <- lvl[1]
        ref_lvl <- lvl %in% exposed
    } else {
        ref_lvl <- lvl %in% exposed
        if (!any(ref_lvl)) {
            stop(" >>> Wrong Reference Value <<< ")
        }
    }

    if (length(lvl) == 2) {
        t <- calc_OR(t, ref_lvl, rnd)
    } else {
        or <- do.call(
            rbind,
            lapply(1:nrow(t), function(z) {
                ref_value <- row.names(t)[z]
                if (ref_value != exposed) {
                    t <- rbind(t[row.names(t) == ref_value,],
                               t[row.names(t) == exposed, ])
                    t <- calc_OR(t, c(TRUE, FALSE), rnd)[-1, ]
                }
            })
        )
        v <- cbind(t[ref_lvl, ], "Reference", "", "", "")
        names(v) <- names(or)
        t <- rbind(v, or)
    }

    if (print.table) {
        printText(
            t,
            paste0("Maximum likelihood estimates of Odds Ratios\n",
                   "Comparing odds of '", case_name, "' among '",
                   exp_name, "'"),
            split = "Comparing"
        )

        if (!is.null(attr(var_case, "label")) |
            !is.null(attr(var_exp, "label"))) {
            printMsg("Labels:")
            printMsg(paste0(exp_name, ": ",
                            attr(var_exp, "label"), collapse = ""))
            printMsg(paste0(case_name, ": ",
                            attr(var_case, "label"), collapse = ""))
        }
    }

    invisible(t)
}

#' @rdname mhodds
#' @export
mhodds.list <- function(var_case, var_exp, data = NULL,
                        exposed = NULL, strata = NULL,
                        na.rm = FALSE, rnd = 3,
                        print.table = TRUE)
{
    arguments <- as.list(match.call())
    case_name <- (deparse(substitute(var_case)))
    exp_name <- (deparse(substitute(var_exp)))
    strata_name <- (deparse(substitute(strata)))

    if (!is.null(data)) {
        var_exp <- eval(substitute(var_exp), data)
        var_case <- eval(substitute(var_case), data)
        strata <- eval(substitute(strata), data)
    }

    t <- mhodds(var_exp = var_exp, var_case = var_case,
                data = NULL, exposed = exposed,
                strata = NULL, na.rm = na.rm, rnd = rnd,
                print.table = FALSE)

    lvl <- row.names(t)

    if (length(lvl) > 2) {
        stop(" >>> Only Stratified 2x2 tables are supported. <<< ")
    }

    if (is.null(exposed)) {
        exposed <- lvl[1]
        ref_lvl <- lvl %in% exposed
    } else {
        ref_lvl <- lvl %in% exposed
        if (!any(ref_lvl)) {
            stop(" >>> Wrong Reference Value <<< ")
        }
    }

    strata.lvl <- sort(unique(strata))


    f <- do.call(
        rbind,
        lapply(strata.lvl, function(z) {
            t.var_case <- var_case[strata == z]
            t.var_exp <- var_exp[strata == z]
            v <- mhodds(t.var_case, t.var_exp, NULL, exposed, NULL,
                        na.rm, rnd, print.table = FALSE)
            d1i <- as.numeric(as.character(v[2, 2]))
            h1i <- as.numeric(as.character(v[2, 1]))
            d0i <- as.numeric(as.character(v[1, 2]))
            h0i <- as.numeric(as.character(v[1, 1]))
            ni <- sum(d1i, h1i, d0i, h0i)
            di <- sum(d1i, d0i)
            hi <- sum(h1i, h0i)
            n1i <- sum(d0i, h0i)
            n0i <- sum(d1i, h1i)

            Q <- d1i * h0i / ni
            R <- d0i * h1i / ni
            V <- (di * hi * n0i * n1i) / ((ni ^ 2) * (ni - 1))
            E <- (di * n1i) / ni
            O <- d1i
            v <- cbind(v[ref_lvl, ], Q, R, V, E, O,
                       d1i, h0i, d0i, h1i, ni)
            row.names(v) <- z
            v
        })
    )

    Q <- sum(f[, "Q"])
    R <- sum(f[, "R"])
    V <- sum(f[, "V"])
    E <- sum(f[, "E"])
    O <- sum(f[, "O"])
    mhor <- Q / R
    mhse <- sqrt(V / (Q * R))
    mhef <- exp(1.96 * mhse)
    chi <- (O - E) ^ 2 / V
    chi <- pchisq(chi, df = 1, lower.tail = FALSE)

    ## Calculations for the Chi-squared test of heterogeneity
    hetero <- ((f[, "d1i"] * f[, "h0i"]) -
                   (mhor * f[, "d0i"] * f[, "h1i"]))^2 /
        (mhor * f[, "V"] * (f[, "ni"]) ^ 2)
    hetero <- sum(hetero)
    hetero <- pchisq(hetero, df = 1, lower.tail = FALSE)

    mhor <- data.frame(cbind("Haenzsel's", "OR", ">>>",
                             sprintf(mhor, fmt = paste0('%#.', rnd, 'f')),
                             sprintf(mhor / mhef, fmt = paste0('%#.', rnd, 'f')),
                             sprintf(mhor * mhef, fmt = paste0('%#.', rnd, 'f')),
                             sprintf(chi, fmt = '%#.5f')
    ))
    f <- f[, 1:7]
    names(mhor) <- names(f)
    row.names(mhor) <- "Mantel"

    f <- rbind(f, mhor)


    if (print.table) {
        printText(
            f,
            paste0("Mantel-Haenzsel estimate of Odds Ratios of ",
                   exp_name, " ~ ", case_name,
                   "\n Stratified by ", strata_name,
                   "\nReference: ", exp_name, " == ", lvl[ref_lvl]),
            split = "Comparing")

        printMsg(paste0("Chi-squared test of heterogeneity [P>Chi2]: ",
                        sprintf(hetero, fmt = '%#.5f')))

        if (!is.null(attr(var_case, "label")) |
            !is.null(attr(var_exp, "label")) |
            !is.null(attr(strata, "label"))) {
            printMsg("Labels:")
            printMsg(paste0(exp_name, ": ",
                            attr(var_exp, "label"), collapse = ""))
            printMsg(paste0(case_name, ": ",
                            attr(var_case, "label"), collapse = ""))
            printMsg(paste0(strata_name, ": ",
                            attr(strata, "label"), collapse = ""))
        }
    }

    invisible(f)
}


#' @param t 2x2 table input to calaculate odds ratio of non-reference
#' @param ref_lvl logical vector to indicate where reference category is
#' @rdname mhodds
#' @export
calc_OR <- function(t, ref_lvl, rnd)
{
    d1 <- as.numeric(as.character(t[ref_lvl, 2]))
    d0 <- as.numeric(as.character(t[!ref_lvl, 2]))
    h1 <- as.numeric(as.character(t[ref_lvl, 1]))
    h0 <- as.numeric(as.character(t[!ref_lvl, 1]))
    SE <- sqrt((1/d1) + (1/h1) + (1/d0) + (1/h0))
    EF <- exp(1.96 * SE)
    OR <- ((d1 * h0) / (d0 * h1))
    lower <- OR / EF
    upper <- OR * EF
    v <- as.table(cbind(as.numeric(as.character(t[, 1])),
                        as.numeric(as.character(t[, 2]))))
    chi <- tryCatch({
        suppressWarnings(chisq.test(v, correct = FALSE)$p.value)
    }, error = function(err) {
        return(NA)
    })
    f <- as.data.frame(cbind(
        sprintf(OR, fmt = paste0('%#.', rnd, 'f')),
        sprintf(lower, fmt = paste0('%#.', rnd, 'f')),
        sprintf(upper, fmt = paste0('%#.', rnd, 'f')),
        sprintf(chi, fmt = '%#.5f')
    ))
    f <- cbind(t[ref_lvl, ], f)
    v <- cbind(t[!ref_lvl, ], "Reference", "", "", "")
    names(f) <- c(names(t), "Odds Ratio", "[95% Conf.",
                  "Interval]", "P>chi2")
    names(v) <- names(f)
    t <- rbind(v, f)
    return(t)
}