# Within-subjects contrasts in repeated measures anova with unbalanced design [closed]

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  ## 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. • I suspect you would get better responses on StackOverflow or R-help – mdewey Feb 14 '17 at 16:37 ## 1 Answer 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

• 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. – Steve Feb 16 '17 at 9:36