mcnemar.R

# k-by-k version of a mcnemar test

mcnemar <-
    function(tab, n_perm=1000)
{
    if(diff(dim(tab)) != 0 || any(tab < 0) || any(round(tab) != tab)) {
        stop("tab should be a square matrix with non-negative integers")
    }

    calc_stat <- function(tab) {
        marg <- rbind(rowSums(tab), colSums(tab))
        expected <- rbind(colMeans(marg), colMeans(marg))
        sum( (marg - expected)^2/expected )
    }

    obs <- calc_stat(tab)
    if(n_perm==0) return(obs)

    perm_stat <- function(tab) {
        # constrain sums of paired off-diagonal elements
        n <- (tab + t(tab))[lower.tri(tab)]

        # sample a new table
        x <- rbinom(length(n), n, 0.5)

        # a bit of contortion to get the values in the right places
        tab[lower.tri(tab)] <- x
        tab <- t(tab)
        tab[lower.tri(tab)] <- n-x

        # make sure the sums of the corresponding off-diagonal elements haven't changed
        if( any((tab + t(tab))[lower.tri(tab)] != n) ) {
            stop("Something got messed up; sums are not identical")
        }

        calc_stat(tab)
    }

    perms <- replicate(n_perm, perm_stat(tab))

    list(observed=obs,
         pvalue=mean(perms >= obs),
         perms=perms)
}