Permutation test explained, briefly

Typically, if we are comparing outcomes for two treatment arms, we calculate a statistic that quantifies the comparison. For example, this could be a difference in group means, a risk ratio, or a log-odds ratio. For whatever statistic we use, there would be an underlying sampling distribution of that statistic under the assumption that there is no difference between the two groups. Typically, the sampling distribution would be estimated analytically using additional assumptions about the underlying distributions of the observed data, such as normal or Poisson. We then use the sampling distribution to calculate a p-value for the observed value of the statistic.

The permutation approach is an alternative way to generate the sampling distribution of the statistic under an assumption of no group differences without making any assumptions about the distributions of the data. If group membership does not influence the outcome, it wouldn’t matter if we re-arranged all the treatment assignments in the data. We could do that and estimate the statistic. In fact, we could do that for all the possible arrangements, and that would give us a distribution of the statistic for that sample under the assumption of no group effect. (If the number of possible arrangements is excessive, we could just take a large sample of those possible arrangements, which is what I do below.) To get a p-value, we compare the observed statistic to this manufactured distribution.

Now to the simulations.

The data generation process

In this proposed study, we are interested in measuring the rate of Covid-19 infection in 8 nursing homes. Given the nature of the spread of disease and the inter-relationship of the infections between residents, the nursing home is the logical unit of analysis. So, we will only have 8 observations – hardly anything to hang your hat on. But, can the permutation provide any useful information?

The data generation process starts with generating a pool of residents at each site, about 15 per home. The study will run for followed for 4 months (or 120 days), and residents will come and go during that period. We are going to assume that the average residence time is 50 days at each home, but there will be some variability. Based on the number of patients and average length of stay, we can calculate the number of patient-days per site. The number of infected patients \(y\) at a site is a function of the intervention and the time of exposure (patient-days). We will be comparing the average rates (\(y/patDays\)) for the two groups.

In the first simulation, I am assuming no treatment effect, because I want to assess the Type 1 error (the probability of concluding there is an effect given we know there is no effect).

Here is a function to generate the data definitions and a second function to go through the simple data generation process:

library(simstudy)
library(parallel)
defs <- function() {
  
  def <- defDataAdd(varname = "nRes", formula = 15, dist = "poisson")
  def <- defDataAdd(def, varname = "nDays", formula = 50, dist = "poisson")
  def <- defDataAdd(def, varname = "patDays", 
                    formula = "nRes * pmin(120, nDays)",
                    dist = "nonrandom")
  def <- defDataAdd(def, varname = "y", 
                    formula = "-4 - 0.0 * rx + log(patDays)", 
                    variance = 1,
                    dist = "negBinomial", link = "log")
  def <- defDataAdd(def, varname = "rate", 
                    formula = "y/patDays", 
                    dist = "nonrandom")
  
  return(def[])

}
gData <- function(n, def) {
  
  dx <- genData(n)
  dx <- trtAssign(dx, grpName = "rx")
  dx <- addColumns(def, dx)
  
  dx[]
  
}

And here we actually generate a single data set:

RNGkind(kind = "L'Ecuyer-CMRG")
#set.seed(72456)
set.seed(82456)

def <- defs()
dx <- gData(8, def)
dx
##    id rx nRes nDays patDays  y        rate
## 1:  1  1   16    40     640  0 0.000000000
## 2:  2  0    9    59     531  8 0.015065913
## 3:  3  1   20    57    1140 58 0.050877193
## 4:  4  1   14    51     714 14 0.019607843
## 5:  5  1    7    59     413  4 0.009685230
## 6:  6  0   14    38     532  4 0.007518797
## 7:  7  0   10    40     400 19 0.047500000
## 8:  8  0   11    56     616 22 0.035714286

The observed difference in rates is quite close to 0:

dx[, mean(rate), keyby = rx]
##    rx         V1
## 1:  0 0.02644975
## 2:  1 0.02004257
obs.diff <- dx[, mean(rate), keyby = rx][, diff(V1)]
obs.diff
## [1] -0.006407182

The permutation test

With 8 sites, there are \(8!\) possible permutations, or a lot ways to scramble the treatment assignments.

factorial(8)
## [1] 40320

I decided to implement this myself in a pretty rudimentary way, though there are R packages out there that can certainly do this better. Since I am comparing averages, I am creating a vector that represents the contrast.

rx <- dx$rx/sum(dx$rx)
rx[rx==0] <- -1/(length(dx$rx) - sum(dx$rx))

rx
## [1]  0.25 -0.25  0.25  0.25  0.25 -0.25 -0.25 -0.25

I’m taking a random sample of 5000 permutations of the contrast vector, storing the results in a matrix:

perm <- t(sapply(1:5000, function(x) sample(rx, 8, replace = FALSE)))
head(perm)
##       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]
## [1,] -0.25 -0.25  0.25 -0.25  0.25  0.25 -0.25  0.25
## [2,] -0.25 -0.25  0.25 -0.25  0.25 -0.25  0.25  0.25
## [3,]  0.25 -0.25 -0.25  0.25 -0.25 -0.25  0.25  0.25
## [4,]  0.25 -0.25  0.25  0.25 -0.25  0.25 -0.25 -0.25
## [5,] -0.25 -0.25 -0.25  0.25  0.25  0.25  0.25 -0.25
## [6,]  0.25 -0.25  0.25 -0.25  0.25 -0.25  0.25 -0.25

Using a simple operation of matrix multiplication, I’m calculating a rate difference for each of the sampled permutations:

perm.diffs <- perm %*% dx$rate
head(perm.diffs)
##              [,1]
## [1,]  0.005405437
## [2,]  0.025396039
## [3,]  0.004918749
## [4,] -0.007490399
## [5,] -0.004336380
## [6,]  0.007538896

Here is an estimate of the 2-sided p-value:

mean(abs(perm.diffs) > abs(obs.diff))
## [1] 0.7166

And finally, here is a histogram of the permuted rate differences, with the observed rate different overlaid as a red line. The observed value lies pretty much right in the middle of the distribution, which is what the p-value suggests:

ggplot(data = data.frame(perm.diffs), aes(x = perm.diffs)) +
  geom_histogram(binwidth = .005, color = "white") +
  geom_vline(xintercept = obs.diff, color = "red") +
  scale_y_continuous(expand = c(0,0), breaks = seq(2500, 10000, 2500)) +
  theme(panel.grid = element_blank())

Operating characteristics

As in a power analysis by simulation, we can estimate the Type 1 error rate by generating many data sets, and for each one calculate a p-value using the permutation test. The proportion of p-values less than 0.05 would represent the Type 1 error rate, which should be close to 0.05.

iter <- function(n) {
  
  dx<- gData(n, defs())
  obs.diff <- dx[, mean(rate), keyby = rx][, diff(V1)]

  rx <- dx$rx/sum(dx$rx)
  rx[rx==0] <- -1/(n - sum(dx$rx))
  
  perm <- t(sapply(1:20000, function(x) sample(rx, n, replace = FALSE)))
  
  perm.diffs <- perm %*% dx$rate
  mean(abs(perm.diffs) > abs(obs.diff))
}

Here we use 5000 data sets to estimate the Type 1 error rate under the data generating process we’ve been using all along, and for each of those data sets we use 5000 permutations to estimate the p-value.

res <- unlist(mclapply(1:5000, function(x) iter(8), mc.cores = 4))
mean(res < .05)
## [1] 0.0542

The risks of using a model (with assumptions)

If we go ahead and try to find a signal using a parametric model, there’s a chance we’ll be led astray. These data are count data, so it would not be strange to consider Poisson regression model to estimate the treatment effect (in this case, the effect would be a rate ratio rather than a rate difference). Given that the data are quite limited, we may not really be in a position to verify whether the Poisson distribution is appropriate; as a result, it might be hard to actually select the right model. (In reality, I know that this model will lead us astray, because we used a negative binomial distribution, a distribution with more variance than the Poisson, to generate the count data.)

Just as before, we generate 5000 data sets. For each one we fit a generalized linear model with a Poisson distribution and a log link, and store the effect estimate along with the p-value.

chkglm <- function(n) {
  
  dx <- gData(n, defs())
  glmfit <- glm( y ~ rx + offset(log(patDays)), family = poisson, data = dx)
  
  data.table(t(coef(summary(glmfit))["rx",]))
  
}

glm.res <- rbindlist(mclapply(1:5000, function(x) chkglm(8)))

The estimated Type 1 error is far greater than 0.05; there would be a pretty good chance that we will be over-enthusiastic about the potential success of our new nursing home strategy if it was not actually effective.

glm.res[, .(mean(`Pr(>|z|)` < 0.05))]
##        V1
## 1: 0.6152

When there is a treatment effect

In the case where there is actually a treatment effect, the observed effect size is more likely to fall closer to one of the distribution’s tails, depending on the direction of the effect. If the treatment reduces the number infections, we would expect the rate difference to be \(< 0\), as it is in this particular case:

def <- updateDefAdd(def, changevar = "y", 
          newformula = "-4 - 1.2 * rx + log(patDays)" )

At the end of the day, if you feel like you must estimate the treatment effect in a pilot study before moving on to the larger trial, one option is to use a non-parametric approach like a permutation test that requires fewer assumptions to lead you astray.

In the end, though, we opted for a different model. If we do get the go ahead to conduct this study, we will fit a Bayesian model instead. We hope this will be flexible enough to accommodate a range of assumptions and give us a potentially more informative posterior probability of a treatment effect. If we actually get the opportunity to do this, I’ll consider describing that model here.

Support: