R/CTA.R

In evidence: Analysis of Scientific Evidence Using Bayesian and Likelihood Methods

CTA <- function(X, extBayes = FALSE) {
    if (length(dim(X)) != 2) 
        stop("X must be a matrix")
    n <- nrow(X)
    m <- ncol(X)
    N <- n * m
    cat("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n")
    cat("               CTA: contingency table analysis\n")
    cat("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n")
    Xname <- deparse(substitute(X))
    cat("Data analyzed:", Xname, ":\n")
    print(X)
    dmnmsX <- dimnames(X)
    par(mfrow = c(2, 1))
    a <- matrix(1, nrow = n, ncol = m)
    y <- X
    ni. <- rowSums(X)
    nj. <- colSums(X)
    bf <- LearnBayes::ctable(X, a)
    cat("\nBayesian analysis (J. Albert):\n")
    cat("Bayes-factor against independence:", round(bf, 
        2), "\n")
    if(extBayes) {
    cat("analyzing...\n")
    logK <- seq(0, 10, 0.2)
    logBF <- 0 * logK
    for (j in 1:length(logK)) {
        q <- LearnBayes::bfindep(X, exp(logK[j]), 1e+05)
        logBF[j] <- log(q$bf)
    }
    logK <- as.numeric(logK)
    mBF <- exp(max(logBF))
    cat("Bayes Factor against independence (near-independence prior):", 
        round(mBF, 2), "\n")
    }
    
    cat("\nModel comparison: \n")
    if (is.null(rownames(X))) 
        rnames <- paste("R", as.character(1:n), sep = "") else rnames <- rownames(X)
    if (is.null(colnames(X))) 
        cnames <- paste("C", as.character(1:m), sep = "") else cnames <- colnames(X)
    cn <- character(N)
    R <- factor(rep(rnames, m), labels = rnames)
    C <- factor(sort(rep(cnames, n)), labels = cnames)
    Y <- as.numeric(X)
    fit0 <- glm(Y ~ R + C, family = poisson)
    fit1 <- glm(Y ~ R * C, family = poisson)
    cat("Null model (independence):\n")
    print(logLik(fit0))
    cat("Alt. model (with dependence):\n")
    print(logLik(fit1))
    cat("logLikelihood ratio against independence:")
    lL <- 2 * (logLik(fit1) - logLik(fit0))
    df.lL <- (n - 1) * (m - 1)
    cat(round(lL, 4), "\n")
    cat("AIC of null (independence) and alt. (dependence) models:\n")
    print(AIC(fit0, fit1))
    Lg0 <- exp(-0.5 * (AIC(fit0) - AIC(fit1)))
    Lg1 <- exp(-0.5 * (AIC(fit1) - AIC(fit0)))
    w0 <- Lg0/(Lg0 + Lg1)
    w1 <- Lg1/(Lg0 + Lg1)
    cat("Evidence-ratio against independence:", round(Lg1/Lg0, 
        3), "\n")
    cat("Predicted counts for null model:\n")
    prd0 <- predict(fit0, type = "response")
    prd0mtrx <- matrix(prd0, n, m)
    dimnames(prd0mtrx) <- dmnmsX
    print(prd0mtrx)
    cat("Predicted counts for alt. model:\n")
    prd1 <- predict(fit1, type = "response")
    prd1mtrx <- matrix(prd1, n, m)
    dimnames(prd1mtrx) <- dmnmsX
    print(prd1mtrx)
    cat("\nFrequentist Chi-squared test:\n")
    print(X2t <- chisq.test(X))
    assocplot(X)
    mosaicplot(X, color = TRUE, type = "dev", main = "")
    if (sum(X) >= 10 * N) 
        cat("Very large table, asymptotic tests reasonable.\n") else cat("Table not very large: asymtotic tests suspect!\n")
    cat("\nVerdict:\n")
    if (bf < 4) 
        cat("Bayes: inconclusive\n")
    if (bf >= 4) 
        cat("Bayes: independence rejected\n")
    if (extBayes && mBF < 4) 
        cat("Bayes with near-independence prior: inconclusive\n")
    if (extBayes && mBF >= 4) 
        cat("Bayes with near-independence prior: independence rejected\n")
    if (AIC(fit1) - AIC(fit0) > 3) 
        cat("AIC: independence established\n")
    if (AIC(fit1) - AIC(fit0) <= 3) 
        cat("AIC: independence rejected\n")
    if (X2t$p.value < 0.05) 
        cat("Frequentist test: independence rejected\n") else cat("Frequentist test: insufficient evidence to reject independence\n")
    
    if (nrow(X) == 2 && ncol(X) == 2) {
        cat("Because the matrix you analysed is a 2 x 2 one you should also\n")
        cat("consider using function Bft2x2.\n")
    }
    par(ask = FALSE)
    par(mfrow = c(1, 1))
}