Repeated Measure Randomization Test

Setting: two within-subject repeated measures and one between-subject measure

Loading the data and “fixing” the data into R-useable format

Original Version from Jeff Witmer

Here’s Jeff Witmer’s code for this.

HowellData.orig <- 
  read.table("http://www.uvm.edu/~dhowell/methods8/DataFiles/Tab14-11.dat", header = TRUE)
dv <- numeric()
for (i in 1:24) {
  w <- with( HowellData.orig, c(C1P1[i], C1P2[i], C2P1[i], C2P2[i], C3P1[i], C3P2[i], C4P1[i], C4P2[i]))
  dv <- c(dv, w)
}

Phase <- factor(rep(1:2, times = 96))
Cycle <- factor(rep(1:4, each = 2, times = 24))
Group = factor(rep(1:3, each = 64))
Subj <- factor(rep(1:24, each = 8))

HowellDataFramed <- data.frame(dv, Phase, Cycle, Group, Subj)
HowellDataFramed %>% filter(Subj == "1")
##   dv Phase Cycle Group Subj
## 1  1     1     1     1    1
## 2 28     2     1     1    1
## 3 22     1     2     1    1
## 4 48     2     2     1    1
## 5 22     1     3     1    1
## 6 50     2     3     1    1
## 7 14     1     4     1    1
## 8 48     2     4     1    1

Using the Hadleyverse

This method (suggested by R Pruim) uses tidyr and dplyr to make the data transformations clearer.

HowellData.orig <- 
  read.table("http://www.uvm.edu/~dhowell/methods8/DataFiles/Tab14-11.dat", header = TRUE)
HowellData <-
  HowellData.orig %>% 
  mutate(Subj = factor(1:nrow(HowellData.orig))) %>%
  gather(key = "key", value = "dv", C1P1:C4P2) %>%
  separate(key, into = c("Cycle", "Phase"), sep = 2) %>%
  mutate(Group = factor(Group),
         Cycle = factor(readr::parse_number(Cycle)), # extract number
         Phase = factor(readr::parse_number(Phase))) # extract number
HowellData %>% filter(Subj == "1")
##   Group Subj Cycle Phase dv
## 1     1    1     1     1  1
## 2     1    1     1     2 28
## 3     1    1     2     1 22
## 4     1    1     2     2 48
## 5     1    1     3     1 22
## 6     1    1     3     2 50
## 7     1    1     4     1 14
## 8     1    1     4     2 48

What we do with one model

# pull the 7 F statistics out of a model object
extractStats <- function(model) {
  smodel <- summary(model)
  c(
    group       = smodel[[1]][[1]][[4]][[1]],
    cycle     = smodel[[2]][[1]][[4]][[1]],
    grpxcycle     = smodel[[2]][[1]][[4]][[2]],
    phase     = smodel[[3]][[1]][[4]][[1]],
    grpxphase     = smodel[[3]][[1]][[4]][[2]],
    cycxphs     = smodel[[4]][[1]][[4]][[1]],
    grpxcycxphs = smodel[[4]][[1]][[4]][[2]]
  )
}

# This can be omitted by setting contrasts in each call to aov()
# options.before <- options(contrasts = c("contr.sum","contr.poly"))

model1 <- 
  aov(dv ~ (Group*Cycle*Phase) + Error(Subj/(Cycle*Phase)), data = HowellData,  
      contrasts = contr.sum)
print(summary(model1))   # Standard repeated measures anova  
## 
## Error: Subj
##           Df Sum Sq Mean Sq F value Pr(>F)  
## Group      2   4617  2308.4   3.083  0.067 .
## Residuals 21  15723   748.7                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Subj:Cycle
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Cycle        3   2727   909.0  12.027 2.53e-06 ***
## Group:Cycle  6   1047   174.5   2.309   0.0445 *  
## Residuals   63   4761    75.6                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Subj:Phase
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Phase        1  11703   11703  129.85 1.88e-10 ***
## Group:Phase  2   4054    2027   22.49 6.01e-06 ***
## Residuals   21   1893      90                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Subj:Cycle:Phase
##                   Df Sum Sq Mean Sq F value  Pr(>F)   
## Cycle:Phase        3    742  247.17   4.035 0.01090 * 
## Group:Cycle:Phase  6   1274  212.30   3.466 0.00505 **
## Residuals         63   3859   61.26                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
observed <- extractStats(model1)

Generation of random data

We use a helper function that samples “block-wise”.

# x is a vector
# block is a vector of equal length
# x[i] and x[j] should match if block[i] == block[j]
# the result is a vector where x values have been permuted "block-wise"

blockSample <- function(x, blocks) {
  id <- as.numeric(factor(blocks))
  blocks <- aggregate(x, by = list(blocks), FUN = function(x) head(x, 1))$x
  sample(blocks)[id]
}

(This is still brittle with respect to the variable names used, and one could imagine fancier versions that avoid this.)

randomData <- function(data) {
  data %>% 
    mutate(
      dv.orig = dv,
      dv = mosaic::shuffle(dv, groups = Subj),  # this is using mosaic::shuffle()
      Group = blockSample(Group, blocks = Subj)
    )
}
randomModel <- function() {
  aov(dv ~ (Group*Cycle*Phase) + Error(Subj/(Cycle*Phase)), 
      data = randomData(HowellDataFramed), 
      contrasts = contr.sum)
}

Doing it lots of times

SimStats <- do(1000) * extractStats(randomModel())

# Summarise permuted data
SimStats %>%
  summarise(
    probgroup   = prop(group >= observed["group"]),
    probcycle = prop(cycle >= observed["cycle"]),
    probgrpxcycle   = prop(grpxcycle >= observed["grpxcycle"]),
    probphase = prop(phase >= observed["phase"]),
    probgrpxphase  = prop(grpxphase >= observed["grpxphase"]),
    probcycxphs = prop(cycxphs >= observed["cycxphs"]),
    probgrpxcycxphs  = prop(grpxcycxphs >= observed["grpxcycxphs"])
  )
##   probgroup probcycle probgrpxcycle probphase probgrpxphase probcycxphs
## 1     0.064         0         0.044         0             0       0.008
##   probgrpxcycxphs
## 1           0.005