Barnard.R

Barnard <- function(data, Tbx = 100, to.plot = TRUE, 
    to.print = TRUE) {
    #Examples:
    #Convictions <-matrix(c(2, 10, 15, 3), nrow = 2, dimnames =
    #   list(c('Dizygotic', 'Monozygotic'), c('Convicted', 'Not
    #   convicted')))
    #Barnard(Convictions)
    #Compare this to:
    #fisher.test(Convictions, alternative = 'less')
    #This table did not work on the previous code due to
    #   computational difficulty:
    #Barnard(matrix(c(10,1,200,5),2,2))
    if (sum(dim(data) == c(2, 2)) != 2) {
        stop("Input matrix must be a 2x2 matrix")
    }
    if (sum(is.finite(data)) != 4 || !is.numeric(data)) {
        stop("All X values must be numeric and finite")
    }
    if (sum(data[1, ]) <= 0 || sum(data[2, ]) <= 0 || sum(data[, 
        1]) <= 0 || sum(data[, 2]) <= 0) {
        stop("Need at least one observation in each row or column")
    }
    Cs <- c(sum(data[, 1]), sum(data[, 2]))
    N <- sum(Cs)
    I <- matrix(0, Cs[1] + 1, Cs[2] + 1)
    for (i in 1:(Cs[1] + 1)) {
        I[i, ] <- rep(i, 1, Cs[2] + 1) - 1
    }
    J <- matrix(0, Cs[1] + 1, Cs[2] + 1)
    for (j in 1:(Cs[2] + 1)) {
        J[, j] <- t(rep(j, 1, Cs[1] + 1) - 1)
    }
    TX <- (I/Cs[1] - J/Cs[2])/sqrt(((I + J)/N) * (1 - ((I + J)/N)) * 
        sum(1/Cs))
    TX[which(is.na(TX))] = 0
    #Wald Statistic:
    TXO <- abs(TX[data[1] + 1, data[3] + 1])
    idx <- matrix(0, Cs[1] + 1, Cs[2] + 1)
    idx[which(TX >= TXO)] <- 1
    B <- Cs + 1
    npa <- seq(1e-04, 0.9999, length = Tbx)
    LP <- log(npa)
    ALP <- log(1 - npa)
    E <- list(I + J)
    for (i in 2:Tbx) {
        E[[i]] <- E[[1]]
    }
    F <- list(N - E[[1]])
    for (i in 2:Tbx) {
        F[[i]] <- F[[1]]
    }
    CF <- list(sum(lgamma(B)) - (lgamma(I + 1) + lgamma(J + 1) + 
        lgamma(B[1] - I) + lgamma(B[2] - J)))
    for (i in 2:Tbx) {
        CF[[i]] <- CF[[1]]
    }
    replaced1 <- list(0)
    for (i in 1:Tbx) {
        replaced1[[i]] <- matrix(LP[i], Cs[1] + 1, Cs[2] + 1)
    }
    replaced2 <- list(0)
    for (i in 1:Tbx) {
        replaced2[[i]] <- matrix(ALP[i], Cs[1] + 1, Cs[2] + 1)
    }
    S <- list(0)
    for (i in 1:Tbx) {
        S[[i]] <- exp(CF[[i]] + E[[i]] * replaced1[[i]] + F[[i]] * 
            replaced2[[i]])
    }
    replaced3 <- list(0)
    for (i in 1:Tbx) {
        replaced3[[i]] <- idx
    }
    Snew <- {
    }
    for (i in 1:Tbx) {
        Snew <- c(Snew, S[[i]])
    }
    nidx <- (Cs[1] + 1) * (Cs[2] + 1)
    dummy1 <- {
    }
    for (i in 1:Tbx) {
        dummy1 <- c(dummy1, Snew[which(replaced3[[i]] > 0) + 
            nidx * (i - 1)])
    }
    cols <- sum(idx[idx == 1])
    dummy2 <- matrix(0, cols, Tbx)
    for (i in 1:Tbx) {
        dummy2[, i] <- dummy1[(cols * i):(cols * (i + 1) - 1) - 
            (cols - 1)]
    }
    P <- {
    }
    for (i in 1:Tbx) {
        P <- c(P, sum(dummy2[, i]))
    }
    #1-tailed p-value:
    PV1 <- max(P)
    #2-tailed p-value:
    PV2 <- min(2 * PV1, 1)
    #Nuisance parameter:
    np <- npa[P == PV1]
    if (to.print) {
        cat("\n", noquote(paste(paste(paste("2x2 matrix Barnard's exact test:", 
            Tbx, seq = ""), paste(paste(Cs[1] + 1, "x", sep = ""), 
            Cs[2] + 1, sep = ""), sep = ""), "tables were evaluated")))
        cat("\n", noquote("-----------------------------------------------------------"))
        cat("\n", noquote(c("Wald statistic = ", format(TXO, 
            digits = 5, scientific = FALSE))))
        cat("\n", noquote(c("Nuisance parameter = ", format(np, 
            digits = 5, scientific = FALSE))))
        cat("\n", noquote(c("p-values: ", "1-tailed = ", format(PV1, 
            digits = 5, scientific = FALSE), "2-tailed = ", format(PV2, 
            digits = 5, scientific = FALSE))))
        cat("\n", noquote("-----------------------------------------------------------"), 
            "\n", "\n")
    }
    if (to.plot) {
        plot(npa, P, type = "l", main = "Barnard's exact P-value", 
            xlab = "Nuisance parameter", ylab = "P-value")
        points(np, PV1, col = 2)
    }
}


barnard.test <- Barnard # just so the option will be available...


#Examples:
# Convictions <-matrix(c(2, 10, 15, 3), nrow = 2, dimnames
#   = list(c('Dizygotic', 'Monozygotic'), c('Convicted',
#   'Not convicted')))
# Barnard(Convictions)
# fisher.test(Convictions, alternative = 'less')