# Is SPSS giving the wrong residuals for a repeated measures design?

When running a repeated measures ANOVA in SPSS, it's possible to 'Save' the residuals as new variables in the data editor.

But the values output do not match the residuals given in R, and seem to be residuals for a between-subjects model. Unless I am missing something? Is SPSS giving the wrong residuals?

Example in R:

set.seed(1)  # hopefully this keeps things the same every time!

# create a data frame with each line representing one subject,
# and create first and second observations for some experiment

DF <- data.frame(participant=factor(1:5), first=rnorm(5, 10, 5), second=rnorm(5, 20, 5))

DF


-

  participant     first   second
1           1  6.867731 15.89766
2           2 10.918217 22.43715
3           3  5.821857 23.69162
4           4 17.976404 22.87891
5           5 11.647539 18.47306


-

  # reshape it for an ANOVA in R
DFlong <- reshape(DF, direction="long", varying=c("first", "second"), v.names="value", idvar="participant", times=c(1, 2), timevar="group")

DFlong


-

    participant group     value
1.1           1     1  6.867731
2.1           2     1 10.918217
3.1           3     1  5.821857
4.1           4     1 17.976404
5.1           5     1 11.647539
1.2           1     2 15.897658
2.2           2     2 22.437145
3.2           3     2 23.691624
4.2           4     2 22.878907
5.2           5     2 18.473058


-

my.aov <- aov(value ~ group + Error( participant / group ), DFlong)
summary(my.aov)


-

Error: participant
Df Sum Sq Mean Sq F value Pr(>F)
Residuals  4 86.474  21.619

Error: participant:group
Df  Sum Sq Mean Sq F value  Pr(>F)
group      1 251.469 251.469  19.871 0.01118 *
Residuals  4  50.619  12.655
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


-

my.aov$"participant:group"$residuals


-

        6          7          8          9         10
0.7066837 -1.0533061 -5.5440267  3.6252135 -2.2654355


-

# import into SPSS:
write.table(DF, "C:/test.txt", row.names=FALSE)


GET DATA  /TYPE = TXT
/FILE = 'C:\test.txt'
/DELCASE = LINE
/DELIMITERS = " "
/QUALIFIER = '"'
/ARRANGEMENT = DELIMITED
/FIRSTCASE = 2
/IMPORTCASE = ALL
/VARIABLES =
participant F1.0
first F16.14
second F16.13
.
CACHE.
EXECUTE.
DATASET NAME DataSet1 WINDOW=FRONT.


Now change the variable types to scale (in the 'variables' tab - I don't know the syntax for this). Then run:

GLM
first second
/WSFACTOR = factor1 2 Polynomial
/METHOD = SSTYPE(3)
/SAVE = RESID
/CRITERIA = ALPHA(.05)
/WSDESIGN = factor1 .


Or, do the above SPSS commands using the GUI: File->Read text data... find C:\test.txt, import it, remember to specify that the file has variable names as the first case, and run:

1. Analyze->General Linear Model->Repeated Measures...

2. Set number of levels to 2

3. Put variables into analysis, 'first' and 'second'.

4. Open 'Save...' dialog box, check 'Residuals->Unstandardized'

5. Run analysis, SPSS creates two variables of residuals:

RES_1    RES_2
-3.78    -4.78
.27     1.76
-4.82     3.02
7.33     2.20
1.00    -2.20


Note these values are different to R. So has SPSS got it wrong?

• You might as well paste your SPSS syntax to make it completely self contained. – Andy W Aug 30 '11 at 15:45

SPSS is giving the residuals from the group means without correcting for individual error.

> my.lm <- lm(value ~ group, DFlong)
> round(matrix(residuals(my.lm),ncol=2),2)
[,1]  [,2]
[1,] -3.78 -4.78
[2,]  0.27  1.76
[3,] -4.82  3.02
[4,]  7.33  2.20
[5,]  1.00 -2.20


The residuals after correcting for individual error are as follows; they're not what you're finding in the aov fit either.

> my.lm <- lm(value ~ group + participant, DFlong)
> round(matrix(residuals(my.lm),ncol=2),2)
[,1]  [,2]
[1,]  0.50 -0.50
[2,] -0.74  0.74
[3,] -3.92  3.92
[4,]  2.56 -2.56
[5,]  1.60 -1.60


The aov fit doesn't have a residuals method, which is a big hint that the residuals it's calculating are probably not what most end users want.

• Yes, you're right about how the residuals are calculated in SPSS, and how to get them in R, thanks a lot for your help! Now, since as you say aov doesn't have a residuals method, I would like to know how to get the residuals (in R) in a more complex design... I think I'll post this as a separate question! – trev Aug 31 '11 at 12:52
• Yes, better as a separate question. I just fit the model that I want the residuals for, as above; we'll see if others have other ideas. – Aaron Aug 31 '11 at 14:12

I've read more now on the repeated measures general linear model, and can explain how the residuals are calculated that R is reporting in my.aov.

A repeated measures ANOVA is applied to a transformation of the within-subject variables, and the residuals of this model are the ones returned from the aov fit.

The transformation multiplies the original variables by a transformation matrix, which can be output in SPSS by adding /PRINT = TEST(MMATRIX) to the GLM command. SPSS doesn't seem to allow the actual transformed variables to be saved, but we can calculate them ourselves, in R, see below.

The way the transformation matrix is generated follows certain rules, outlined here: http://www.uccs.edu/~faculty/lbecker/SPSS/glm_1withn.htm

i.e. the sum of coefficients should be zero, the sum of squares of coefficients should be 1. But for now I don't deal with generating the transformation matrix, I simply apply the matrix output from SPSS to our original variables to create the transformed variable, and then check it's the model being used by aov by comparing the residuals from my.aov to an lm call on our own model.

# As before, continuing from the code in the original question:
# my.aov <- aov(value ~ group + Error( participant / group ), DFlong)
summary(my.aov)

# create our transformed variable
# (there will always be one less new variable
# than the number of within subject variables)
# We are applying here the transformation matrix that is output from SPSS,
DF$DV <- sqrt(0.5) * DF$first + -sqrt(0.5) * DF$second residuals(lm(DV ~ 1, DF)) # Our self-calculated residuals should be the same as the ones returned by the aov fit. # They are the same except for the last one, # which is the negative of the residual calculated my.aov$"participant:group"$residuals # Again fitted values are the same except for the *-1 at the end fitted.values(lm(DV ~ 1, DF)) my.aov$"participant:group"\$fitted.values


I am not sure how aov calculates the participant residuals, but presumably that is why they have the *-1, because I have not included individual differences in the model above.

In summary the residuals saved by SPSS are not from the repeated-measures GLM.

• Nice explanation. Did you find anything about why someone might be interested in these residuals rather than the ones from the fit for each strata? – Aaron Sep 2 '11 at 12:54
• As far as I know, these are the ones people should be interested in for a repeated-measures design, because they are the residuals from the model that we are applying to our data. The df from the anova = 4, so we would expect 5 residuals (only one variable is being modelled in this case, that represents the difference between the two groups, i.e. the within-subjects effect). I am sure the values are somehow mathematically related to the ones you calculated above. – trev Sep 2 '11 at 13:41
• With regards to why we apply the model in this way to our data, I think it's because we can specify our predictions by adjusting the transformation matrix, e.g. to predict for a quadratic or cubic effect of time. But I've certainly never done that in an analysis before. – trev Sep 2 '11 at 13:44
• I guess it depends what you're doing with the residuals. If you're doing diagnostics, I would think you'd want the ones from each strata so you could easily identify which observation might be troublesome. – Aaron Sep 2 '11 at 13:49
• You mean after correcting for individual error? These are symmetrical around zero, so you could not say either 'first' or 'second' measurements were troublesome. Essentially with two time points there is only one within-subject error per participant. Or do you mean finding which participant might be troublesome? In which case you can do that from the 5 residuals from the model: each one corresponds to the within-subject error for each participant. If one is particularly high it could be an outlier. – trev Sep 2 '11 at 13:56