# Solved – Within-subjects contrasts in repeated measures anova with unbalanced design

I would like to reproduce the standard output of a repeated measures anova in SPSS with R. For a balanced design (i.e. equal group sizes) I can exactly reproduce the SPSS output of the within-subjects contrasts. But when using an unbalanced (i.e. different group sizes) the results of the within-subjects contrasts are not exactly the same, especially for the main effect of the within-subject factor.

Create a test dataset: (un)comment the line `cond` to create (un)balanced dataset

``set.seed(1234) t <- 4 idmax <- 40  myData <- data.frame(PID = rep(seq(from = 1, to = idmax, by = 1), t),                  #cond = rep(c(rep(c("A"), 20), rep(c("B"), 20)), t), #BALANCED                  cond = rep(c(rep(c("A"), 22), rep(c("B"), 18)), t), #UNBALANCED                  time = rep(x = 1:t, each = idmax),                  acc = sample(x = 1:100, size = idmax*t, replace = TRUE) )  myData <- within(myData, {                  PID   <- factor(PID)                  cond <- factor(cond)                  time <- factor(time)                  acc <- acc/100 }) ``

Using the `ezAnova` function of the `ez` package I can exactly reproduce the standard SPSS output of the Between-Subject Effects and the Within-Subject Effects for both balanced and unbalanced designs.

``library(ez) ezANOVA( data = myData, dv = acc, wid = PID, within = c("time"), between = cond, type = 3) ``

But I cannot reproduce the Within-subjects Contrasts table of the SPSS output for unbalanced designs.

For a balanced design I can reproduce the results using the `aov` function.

``time.L <- C(myData\$time, poly, 1) time.Q <- C(myData\$time, poly, 2) time.C <- C(myData\$time, poly, 3) mdl_balanced <- aov(acc ~ cond*time + Error(factor(PID)/(time.L + time.Q + time.C)), data = myData) summary(mdl_balanced) ``

R Output:

``Error: factor(PID)           Df Sum Sq Mean Sq F value Pr(>F) cond       1  0.128 0.12826   1.524  0.225 Residuals 38  3.198 0.08416                 Error: factor(PID):time.L           Df Sum Sq Mean Sq F value Pr(>F)   time       1 0.2999 0.29993   4.213 0.0471 *                   <--1 cond:time  1 0.0000 0.00003   0.000 0.9842   Residuals 38 2.7055 0.07120                  --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Error: factor(PID):time.Q           Df Sum Sq Mean Sq F value Pr(>F)   time       1  0.553  0.5534   6.210 0.0172 *                   <--2 cond:time  1  0.014  0.0135   0.152 0.6992   Residuals 38  3.386  0.0891                  --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Error: factor(PID):time.C           Df Sum Sq Mean Sq F value Pr(>F) time       1 0.1208 0.12079   1.507  0.227                     <--3 cond:time  1 0.0259 0.02588   0.323  0.573 Residuals 38 3.0454 0.08014 ``

SPSS output

``Tests of Within-Subjects Contrasts                       Measure:   MEASURE_1  Source      time    Type III SS df  Mean Square    F    Sig. time        Linear      ,300    1   ,300       4,213    ,047   <--1             Quadratic   ,553    1   ,553       6,210    ,017   <--2             Cubic       ,121    1   ,121       1,507    ,227   <--3 time * cond Linear  2,812E-005  1   2,812E-005  ,000    ,984             Quadratic   ,014    1   ,014        ,152    ,699             Cubic       ,026    1   ,026        ,323    ,573 Error(time) Linear     2,706    38  ,071                     Quadratic  3,386    38  ,089                     Cubic       3,045   38  ,080 ``

But if I run the same `aov` code for the unbalanced design. The results in the R output for the within-subject factor `time` clearly differ from the SPSS output (indicated with `<--`). The R output for `cond:time` match the SPSS output.

``Error: factor(PID)           Df Sum Sq Mean Sq F value Pr(>F) cond       1  0.114 0.11398   1.348  0.253 Residuals 38  3.213 0.08454                 Error: factor(PID):time.L           Df Sum Sq Mean Sq F value Pr(>F)   time       1 0.2999 0.29993   4.229 0.0466 *                   <--1 cond:time  1 0.0106 0.01063   0.150 0.7009   Residuals 38 2.6949 0.07092                  --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Error: factor(PID):time.Q           Df Sum Sq Mean Sq F value Pr(>F)   time       1  0.553  0.5534   6.202 0.0172 *                   <--2 cond:time  1  0.009  0.0093   0.104 0.7491   Residuals 38  3.391  0.0892                  --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Error: factor(PID):time.C           Df Sum Sq Mean Sq F value Pr(>F) time       1 0.1208 0.12079   1.527  0.224                     <--3 cond:time  1 0.0649 0.06486   0.820  0.371 Residuals 38 3.0064 0.07912                ``

SPSS output

``Tests of Within-Subjects Contrasts                       Measure:   MEASURE_1  Source      time    Type III SS df  Mean Square    F    Sig. time        Linear      ,286    1   ,286       4,030    ,052   <--1             Quadratic   ,562    1   ,562       6,301    ,016   <--2             Cubic       ,103    1   ,103       1,297    ,262   <--3 time * cond Linear      ,011    1   ,011        ,150    ,701             Quadratic   ,009    1   ,009        ,104    ,749             Cubic       ,065    1   ,065        ,820    ,371 Error(time) Linear     2,695    38  ,071                     Quadratic  3,391    38  ,089                     Cubic      3,006    38  ,079 ``

So, my question is: How can I compute the correct F-values for polynomial contrasts for the factor `time` in an unbalanced repeated measures anova design using R?

PS. I can use a linear-mixed effects model for the unbalanced design in combination with `fit.contrast`, but this produces t-values for the contrasts instead of F-values.

``options(contrasts = c("contr.sum", "contr.poly")) library(nlme) model.lme <- lme(acc ~ cond*time, random = ~ 1 | PID, data = myData) anova(model.lme, type="marginal")  library(gmodels) fit.contrast(model.lme, "time", t(contr.poly(4, c(1, 2, 3, 4))))           Estimate Std. Error  t-value    Pr(>|t|) time.L 0.08495364  0.0448776 1.893008 0.060892471 time.Q 0.11915404  0.0448776 2.655089 0.009063148 time.C 0.05090443  0.0448776 1.134295 0.259050482 ``
Contents

I'm not familiar with the R package `ez`, but I can offer an alternative approach with the R package `afex`:
``require(afex) fm <- aov_car(acc ~ cond+Error(PID/time), data = myData) ``
``require(phia) testInteractions(fm\$lm, custom=list(time=contr.poly(4, c(1, 2, 3, 4))[,1]), idata = fm\$data[["idata"]])           Value Df test stat approx F num Df den Df  Pr(>F)   time1 0.084954  1  0.095882   4.0299      1     38 0.05185 .  testInteractions(fm\$lm, custom=list(time=contr.poly(4, c(1, 2, 3, 4))[,2]), idata = fm\$data[["idata"]])          Value Df test stat approx F num Df den Df  Pr(>F)   time1 0.11915  1   0.14223   6.3011      1     38 0.01644 *  testInteractions(fm\$lm, custom=list(time=contr.poly(4, c(1, 2, 3, 4))[,3]), idata = fm\$data[["idata"]])           Value Df test stat approx F num Df den Df Pr(>F) time1 0.050904  1  0.033005    1.297      1     38 0.2619 ``