Online Trick: permutation test

Showing posts with label permutation test. Show all posts

Sunday, September 20, 2015

Permutation Test with Stratified Data and Repeated Measurements

This is an example for a permutation test on stratified samples with repeated measurements. Samples are interdependent firstly because they come from several sites and secondly because the sampling was repeated a second time. That is samples of the same sites are dependent and sample t1 and sample t2, taken from the very same places are dependent.

What I want to test is whether there is a difference between timepoint one (t1) and two (t2) or not. A hypothesis could be: the average difference t1-t2 is sign. larger than zero (a one-sided test). Another hypothesis could be: the average difference is sign. different from zero, either larger or smaller (a two-sided test).

If you deal with a distribution of your data that ordinary Linear Mixed Models (LMMs) or Generalized LMMs (GLMMs) can handle you should vote for this option - but sometimes you deal with awkard data and permutation tests may the only thing to bail you out...

library(vegan)

### the data:

sites <- sort(rep(letters[1:8],6))
tm <- rep(1:2,24)
pairs <- gl(24, 2)
n <- length(sites)
id <- 1:n

print(data <- data.frame(id, sites, tm, pairs))

### set up a Control object:

ctrl <- permControl(strata = sites, type = "series")

### check permutation scheme -
### by permuting the ids of the pair vector with the above
### control we assure that time points are either changed
### or unchanged within each stratum - this allows for the
### dependency within strata - the nr. of possible permutations
### is 2^n(strata) which is 2^8 = 256, meaning that the smallest
### p-value that we will be able to obtain is 1/256 = 0.0039

### in this frame you can retry and see that what i described above
### really happens:

data.frame(id, sites, pairs, tm, tm_p = tm[permuted.index2(n, control = ctrl)])

### now what i want to test is pair differences - i will
### yield differences by aggregating pairs, i will use the
### aggregate function with FUN = sum, to yield differences
### i set up a "sign" variable, with the negative sign aligned
### to the second time point, that is a positive difference
### means a decrease with time and vice versa:

sign <- rep(c(1,-1), 24)

### this will be a simulated dependent variable with a time-effect
### at the 2nd time point var is decreased

var <- rnorm(48, 10, 2)
var <- ifelse(tm == 2, var * 0.7, var * 1)

### in a frame:
data.frame(id, sites, pairs, tm, var)

### x will be the differences:
### you see, the added effect simulates a decrease

print(true <- aggregate(x = var*sign, by = list(sites = sites, pairs = pairs), FUN = sum))

### and now the actual permutation comes in -
### we permute the sign and than aggregate which
### gives the differences of permuted pairs

aggregate(x = var*sign[permuted.index2(n, control = ctrl)],
by = list(pairs = pairs), FUN = sum)

### and now we repeat the above step and record the mean differences
### which will give us a null-population - this we will compare to
### the mean of the observed differences

m.true <- mean(true$x)

### set number of perms

B <- 1999

### setting up frame for the null-population of differences:

pop <- rep(NA, B + 1)
pop[1] <- m.true

for(i in 2:(B+1)){
perm <- aggregate(x = var*sign[permuted.index2(n, control = ctrl)],
by = list(pairs = pairs), FUN = sum)
pop[i] <- mean(perm$x)
}

### show the null-population and the observed difference:

hist(pop, xlab = "Differences")
abline(v = pop[1], col = "red", lty = 3)
text(x = pop[1]*0.9, y = 150, srt = 90, "the true Difference")

### depending on your hypothesis you will take the one-sided p-value:

print(pval_1 <- sum(pop >= pop[1]) / (B + 1))

### or the two-sided, which simply is:

print(pval_2 <- sum(abs(pop) >= abs(pop[1])) / (B + 1))

### the one-sided test would be appropiate if your hypothesis
### was taht the observed difference is sign. larger than the null.
### the two-sided is appropiate if your hypothesis was that
### the obs. diff. is different from the null - that is, it could be larger
### or smaller than the null.

### You can check the result without time effect by resetting
### the var like beneath, then re-run the test:
var <- rnorm(48, 10, 2)

To cite package ‘vegan’ in publications use:

Jari Oksanen, F. Guillaume Blanchet, Roeland Kindt, Pierre Legendre, R.

B. O'Hara, Gavin L. Simpson, Peter Solymos, M. Henry H. Stevens and

Helene Wagner (2011). vegan: Community Ecology Package. R package

version 1.17-9. http://CRAN.R-project.org/package=vegan

Legendre & Borcard: Nested ANOVA by Permutation

..here's a very illustrative R-Script example by Legendre & Borcard showing how a Nested ANOVA can be done by permutation.

nest.anova.perm <- function(Y, Xmain, Xnes, nperm=0)
#
##########################################################################
#
# Nested anova with a main factor and a nested factor.
# Balanced design, i.e., equal replication in each level of the nested factor.
#
# Arguments
#
#   Y     : Vector containing the response variable.
#
#   Xmain : Vector or data.frame containing the values of the main factor for
#           each observation.
#
#   Xnes : Vector or data.frame containing the values of the nested factor for
#           each observation.
#
# Note    : Xmain and Xres are transformed to factors inside the function.
#
#                             Daniel Borcard and Pierre Legendre, September 2007
#
##########################################################################
#
## Example from Sokal & Rohlf (1981, 1995), Box 10.1
#
# Var   = c(58.5,59.5,77.8,80.9,84.0,83.6,70.1,68.3,69.8,69.8,56.0,54.5, 50.7,49.3,63.8,65.8,56.6,57.5,77.8,79.2,69.9,69.2,62.1,64.5)
# Xmain = c(1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3)
# Xnes = c(1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,11,11,12,12)
#
# nest.out = nest.anova.perm(Var, Xmain, Xnes, nperm=999)
#
##########################################################################
{
a.fac = as.factor(Xmain)
b.fac = as.factor(Xnes)

dat = as.data.frame(cbind(Y, a.fac, b.fac))

# Check model balance (equal replication in each level of the nested factor)
balance = var(as.vector(table(dat[,2])))
if(balance > 0)
stop("Design unbalanced. This function can only handle balanced designs.")

# The nesting is done by %in% in 'lm'
anova.res = anova(lm(Y ~ a.fac + b.fac %in% a.fac))
if(nrow(anova.res)<3) cat("Problem with the data, probably no replicates",'\n',"Permutation test on the main factor only",'\n')

# Recompute parametric F-statistic and associated probability for the main factor
F.main = anova.res[1,3] / anova.res[2,3]
P.main = pf(F.main, anova.res[1,1], anova.res[2,1], lower.tail=FALSE)
anova.res[1,4] = F.main
anova.res[1,5] = P.main

# Permutation test (if requested)

if(nperm>0){

n = nrow(dat)
vec = seq(1:n)
nblock = length(levels(a.fac))
nobs.block = nrow(dat)/nblock
k = nrow(anova.res) - 1
Pperm = c(rep(0,k), NA)

# Permutation test of the nested factor. The data are permuted within the levels
# of the main factor. Done only if replicates are available in the nested factor.
GEn = NA
if(nrow(anova.res)>=3){

# First, make sure that the data are ordered on the main factor.
data2 = dat[order(dat[,2]),]

# Restricted permutation
GEn = 1
for(i in 1:nperm){
vecperm = restrictedPerm(nobs.block, nblock, n, restPerm=TRUE, vec)
Y.perm = data2[order(vecperm),1]
anova.perm = anova(lm(Y.perm ~ a.fac + b.fac %in% a.fac))
if(anova.perm[2,4] >= anova.res[2,4]) GEn = GEn + 1
                 }
    Pperm[2] = GEn/(nperm+1)
                      }

# Permutation test of the main factor
GEm = 1
for(j in 1:nperm){
Y.perm = sample(Y, n)
anova.perm = anova(lm(Y.perm ~ a.fac + b.fac %in% a.fac))
F.main.perm = anova.perm[1,3] / anova.perm[2,3]
if(F.main.perm >= F.main) GEm = GEm + 1
                 }
    Pperm[1] = GEm/(nperm+1)

anova.res = data.frame(anova.res, Pperm)
colnames(anova.res) = c("Df", "Sum Sq", "Mean Sq", "F value", "Prob(param)", "Prob(perm)")
note = "Nested anova, parametric and permutation tests"

          }

# return(anova.res)

if(nperm <= 0) {
   note = "Nested anova, parametric tests"
   return(list(anova.type=note,anova.table=anova.res))
   } else {
   return(list(anova.type=note, nperm=nperm, anova.table=anova.res))
   }

}

restrictedPerm <- function(nobs.block, nblock, n, restPerm=TRUE, vec)
#
# restPerm == F: Unrestricted permutation.
#
# restPerm == T: Restricted permutation of the observations within each block.
# The data are arranged by blocks, with all observations forming a block
# placed in adjacent positions in the file. Example:
# BLOCK-1: obs1, obs2, ... obs-nobs.block; BLOCK-2: obs1, obs2, ... obs-nobs.block; etc.
#
# Vector 'vec' contains the initial order of the objects, e.g. vec=c(1:n).
# At the end of the function, it gives the permuted order of the objects.
#
# Examples: toto0 <- restrictedPerm(6,4,24,FALSE,c(1:24))
#            toto1 <- restrictedPerm(6,4,24,TRUE, c(1:24))
#
#                                       Pierre Legendre, January 2006
{
if(restPerm == FALSE) { vec <- sample(vec[1:n],n)

   } else {

   for(j in 1:nblock) {
      i1 <- nobs.block*(j-1)+1
      i2 <- nobs.block*j
      vec[i1:i2] <- sample(vec[i1:i2],nobs.block)
      }
   }

return(vec)
}

See this and other examples at:

http://www.bio.umontreal.ca/legendre/indexEn.html#RFunctions and
http://www.bio.umontreal.ca/numecolR/