Nested ANOVA permutation test in R

A R function for the permutation test for onw-way nested design ANOVA.

This work is inspired by Dr. David C. Howell, University of Vermont. In Dr. Howell’s webpage “Permutation Tests for Factorial ANOVA Designs” (http://www.uvm.edu/~dhowell/StatPages/More_Stuff/Permutation%20Anova/PermTestsAnova.html [access on May 28, 2012]), he described the algorithm of permutation test for a one-way nested design ANOVA by using R language. I followed this algorithm and made a R function to do this test.

	# Copyright 2012 Chen-Pan Liao
	# This program is free software: you can redistribute it and/or modify
	# it under the terms of the GNU General Public License as published by
	# the Free Software Foundation, either version 2 of the License, or
	# any later version.
	#
	# This program is distributed in the hope that it will be useful,
	# but WITHOUT ANY WARRANTY; without even the implied warranty of
	# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	# GNU General Public License for more details.
	#
	# You should have received a copy of the GNU General Public License
	# along with this program. If not, see <http://www.gnu.org/licenses/>.
	# Author: Chen-Pan Liao
	# Date: May 28, 2012
	#
	# This work is inspirated by Dr. David C. Howell, University of Vermont.
	# See Dr. Howell's webpage to learn more:
	# http://www.uvm.edu/~dhowell/StatPages/More_Stuff/
	# Permutation%20Anova/PermTestsAnova.html
	nestedPermutationAnova <- function(
	dv, # Dependent variable
	treatment, # Treatment
	subject, # Experimental unit
	nreps = 1999 # Number of permutation
	) {
	# initiation
	dv <- as.numeric(dv)
	treatment <- as.factor(treatment)
	subject <- as.factor(subject)
	nreps <- as.integer(nreps)
	subject <- as.factor(paste(subject,treatment,sep="_in_"))
	# rename subject as "subject_in_treatment"
	treatment.name <- levels(treatment)
	subject.name <- levels(subject)
	treatment.n <- nlevels(treatment)
	subject.n <- nlevels(subject)
	f.treatment <- numeric(nreps)
	f.subject <- numeric(nreps);
	# original model
	model.0 <- summary(aov(dv ~ treatment / subject))
	f.subject.0 <- model.0[[1]][1,4]
	f.treatment.0 <- model.0[[1]][1,3] / model.0[[1]][2,3]
	# permutation
	for (i in 1:nreps) {
	# permutation for effect of subject
	dv.subject <- NULL
	for (j in 1:treatment.n){
	dv.subject <- c(
	dv.subject,
	sample(dv[ treatment == treatment.name[j] ])
	)
	}
	model <- summary(aov(dv.subject ~ treatment / subject))
	f.subject[i] <- model[[1]][1,4]
	# permutation for effect of treatment
	dv.treatment <- NULL
	subject.name.reorder <- sample(subject.name)
	for (j in 1:subject.n){
	dv.treatment <- c(
	dv.treatment,
	dv[ subject == subject.name.reorder[j] ]
	)
	}
	model <- summary(aov(dv.treatment ~ treatment / subject))
	f.treatment[i] <- model[[1]][1,3] / model[[1]][2,3]
	}
	# p-value
	p.treatment <- (sum(f.treatment >= f.treatment.0)+1) / (nreps+1)
	p.subject <- (sum(f.subject >= f.subject.0)+1) / (nreps+1)
	# plot
	par(mfrow = c(1,2))
	hist(f.treatment, breaks=40, xlab="Sudo f-value for treatment", main="")
	abline(v=f.treatment.0, lty=2)
	hist(f.subject, breaks=40, xlab="Sudo f-value for subject", main="")
	abline(v=f.subject.0, lty=2)
	# return
	return (
	list(
	p.treatment = p.treatment,
	p.subject = p.subject,
	nreps = nreps
	)
	)
	}

view raw nestedPermutationAnova.R hosted with ❤ by GitHub

Example

dat <- read.csv(textConnection("
 trt , unit , obs
 1 , 1 , 11
 1 , 1 ,  9
 1 , 1 ,  9
 1 , 2 ,  8
 1 , 2 ,  7
 1 , 2 ,  6
 1 , 3 ,  8
 1 , 3 , 10
 1 , 3 , 11
 2 , 4 , 11
 2 , 4 ,  8
 2 , 4 ,  7
 2 , 5 , 10
 2 , 5 , 14
 2 , 5 , 12
 2 , 6 ,  9
 2 , 6 , 10
 2 , 6 ,  8
"))

## traditional nested ANOVA
mod.2 <- aov(
 obs ~ factor(trt) + Error(factor(unit)),
 data = dat
)
summary(mod.2) # alternative

## permutation nested ANOVA
nestedPermutationAnova(dat$obs, dat$trt, dat$unit, 499)