R/permutationTest.R

In surveillance: Temporal and Spatio-Temporal Modeling and Monitoring of Epidemic Phenomena

################################################################################
### Permutation test to compare the means of paired samples
###
### Copyright (C) 2011-2012 Michaela Paul, 2013-2015 Sebastian Meyer
###
### This file is part of the R package "surveillance",
### free software under the terms of the GNU General Public License, version 2,
### a copy of which is available at https://www.R-project.org/Licenses/.
################################################################################

permutationTest <- function(score1, score2, nPermutation = 9999,
                            plot = FALSE, verbose = FALSE)
{
    stopifnot((nTime <- length(score1)) == length(score2),
              !is.na(score1), !is.na(score2))

    meanScore1 <- mean(score1)
    meanScore2 <- mean(score2)
    diffObserved <- meanScore1 - meanScore2

    diffMean <- replicate(nPermutation, {
        sel <- rbinom(nTime, size=1, prob=0.5)
        g1 <- (sum(score1[sel==0]) + sum(score2[sel==1]))/nTime
        g2 <- (sum(score1[sel==1]) + sum(score2[sel==0]))/nTime
        g1 - g2
    })

    if (isTRUE(plot)) plot <- list()
    if (is.list(plot)) {
        do.call("permtestplot", args = modifyList(
            list(permstats = diffMean, xmarks = c("observed" = diffObserved),
                 xlab = "Difference between means", ylab = "Density", main = ""),
            plot))
    }

    pVal <- (1+sum(abs(diffMean)>=abs(diffObserved))) / (nPermutation+1)

    pTtest <- t.test(score1, score2, paired=TRUE)$p.value

    if (verbose)
        cat("mean difference =", diffObserved,
            "\tp(permutation) =", pVal,
            "\tp(paired t-test) =", pTtest, "\n")

    list(diffObs=diffObserved, pVal.permut=pVal, pVal.t=pTtest)
}