R code for GLMMs and randomized permutations for a BW - one between and one within factor interaction
The Akaike Information Criterion may be compared for the observed data and bootstrap samples using the Generalized linear mixed Poisson model below if we can assume that the repeated measures factor can be regarded as an ordinal variable (e.g. lying < sitting < standing). . In this example (diagnostic) group is between subjects and position within subjects (with three levels: lying, sitting and standing). The t-value for the interaction can then be read off.
library(foreign)
x <- read.spss("U://My Documents//WHIMINT.sav")
x1 <- attach(x)
library(foreign)
library(lme4)
int <-glmer( resp ~ position + grp + position*grp + (1|sub), data = x, family = poisson)
noint <-glmer( resp ~ position + grp + (1|sub), data = x, family = poisson)One could assess the importance of the interaction treating position as a factor by comparing the GLMMs with and without the interaction present to see how much, if at all, the AIC decreases with the addition of the interaction.
library(foreign)
x <- read.spss("U://My Documents//WHIMINT.sav")
x1 <- attach(x)
library(lme4)
int <-lmer( log(resp) ~ as.factor(position) + as.factor(grp) + as.factor(position)*as.factor(grp) + (1|sub), data = x)
int2 <-lmer( log(resp) ~ as.factor(position) + as.factor(grp) + (1|sub), data = x)
diff <- AIC(int) - AIC(int2)Randomized permutation tests
Edgington (1995) presents a method of obtaining reliable p-values for small samples sizes using all possible permutations of responses. In addition on pages 50-51 he shows that the p-values from using all possible permutations are close to those generated using a large number (N=10000) of random bootstrap permutations which corresponds to randomly sampling the observed responses (with replacement). R is a natural choice to perform randomized permutations as it has a straightforward concise procedure, sample(), to produce bootstrap samples. A variety of fit diagnostics may be used to assess the importance of the interaction.
For repeated measures ANOVA treating the response as continuous the syntax below uses a logged response to handle any skew in the response with the difference in the Akaike Information Criterion (AIC) comparing model fits of the ANOVA (fitted as a mixed model) with and without the interaction term as the fit statistic of interest.
library(foreign)
x <- read.spss("U://My Documents//WHIMINT.sav")
x1 <- attach(x)
library(lme4)
int <-lmer( log(resp) ~ as.factor(position) + as.factor(grp) + as.factor(position)*as.factor(grp) + (1|sub), data = x)
int2 <-lmer( log(resp) ~ as.factor(position) + as.factor(grp) + (1|sub), data = x)
diff <- AIC(int) - AIC(int2)
nb<-10000
n <- length(resp)
boot<-matrix(NA,nb,6)
attributes(boot)
for (i in 1:nb)
{
bs<-sample(resp,n,replace =F)
int <-lmer( log(bs) ~ as.factor(position) + as.factor(grp) + as.factor(position)*as.factor(grp) + (1|sub), data = x)
boot[i,1] <- AIC(int)
int2 <-lmer( log(bs) ~ as.factor(position) + as.factor(grp) + (1|sub), data = x)
boot[i,2] <- AIC(int2)
boot[i,3] <- AIC(int2) - AIC(int)
}
sum(boot[,3]<=diff)/nb # GIVES randomization test p-valueOne could alternatively compute a randomized permutation p-value for the differences in AICs with and without the interaction from fitting a Poisson GLMM although this takes much longer to run than using the above lmer() procedure which assumes a continuous response.
library(foreign)
x <- read.spss("U://My Documents//WHIMINT.sav") #move the interaction SPSS file to the My Documents folder
x1 <- attach(x)
library(lme4)
int <-glmer( resp ~ as.factor(position) + as.factor(grp) + as.factor(position)*as.factor(grp) + (1|sub), data = x, family = poisson)
int2 <-glmer( resp ~ as.factor(position) + as.factor(grp) + (1|sub), data = x, family = poisson)
diff <- AIC(int) - AIC(int2)
nb<-10000
n <- length(resp)
boot<-matrix(NA,nb,6)
attributes(boot)
for (i in 1:nb)
{
bs<-sample(resp,n,replace =F)
int <-glmer( bs ~ as.factor(position) + as.factor(grp) + as.factor(position)*as.factor(grp) + (1|sub), data = x, family = poisson)
boot[i,1] <- AIC(int)
int2 <-glmer( bs ~ as.factor(position) + as.factor(grp) + (1|sub), data = x, family = poisson)
boot[i,2] <- AIC(int2)
boot[i,3] <- AIC(int2) - AIC(int)
}
sum(boot[,3]<=diff)/nbAlternatively the glm() procedure, for the general linear model, may be used comparing models with and without the interaction term on changes in either log likelihoods or in the interaction sum of squares defined as the difference in the residual sum of squares from fitting a model with and without the interaction term (see R syntax below). In practice the choice of fit diagnostic should not make a lot of difference and the different diagnostics should give very similar p-values.
library(multcomp)
library(MASS)
library(nlme)
library(foreign)
library(car)
library(stats)
x <- read.spss("U://My Documents//WHIMINT.sav")
out <- glm(x$resp ~ as.factor(x$grp)*as.factor(x$position))
logLik(out)
out2 <- glm(x$resp ~ as.factor(x$grp)+as.factor(x$position))
logLik(out2)
rssint <-sum(out2$residuals*out2$residuals)- sum(out$residuals*out$residuals)
diff <- mean(logLik(out))-mean(logLik(out2)) # change in model fit when adding in the interaction term
nb<-10000
n <- length(x$resp)
boot<-matrix(NA,nb,6)
attributes(boot)
for (i in 1:nb)
{
bs<-sample(x$resp,n,replace =F)
int <- glm(bs ~ as.factor(x$grp)*as.factor(x$position))
boot[i,1] <- mean(logLik(int))
int2 <- glm(bs ~ as.factor(x$grp)+as.factor(x$position))
boot[i,2] <- mean(logLik(int2))
boot[i,3] <- mean(logLik(int))-mean(logLik(int2))
boot[i,4] <-sum(int2$residuals*int2$residuals)- sum(int$residuals*int$residuals)
}
sum(boot[,3]>=diff)/nb # Gives p-value using changes in log likelihood
sum(boot[,4]>=rssint)/nb # Gives p-value for residual SS > than the observed RSS (=rssint). Reference
Edgington, ES (1995). Randomization tests. Third Edition. Marcel Dekker Inc:New York.