1
$\begingroup$

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
$\endgroup$

closed as off-topic by mdewey, Michael Chernick, John, Peter Flom Feb 16 '17 at 13:59

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – mdewey, Michael Chernick, Peter Flom
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ I suspect you would get better responses on StackOverflow or R-help $\endgroup$ – mdewey Feb 14 '17 at 16:37
0
$\begingroup$

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)

Then you have the results that match up your SPSS output:

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
$\endgroup$
  • $\begingroup$ Many thanks for your answer @bluepole. That's exactly what I needed. Do you know if I can also use this to compute the polynomial contrasts for the time * cond interaction. Would be nice if I can compute the polynomial contrast for time and time * cond with the same function. $\endgroup$ – Steve Feb 16 '17 at 9:36

Not the answer you're looking for? Browse other questions tagged or ask your own question.