## Translating SPSS to R: Non-orthogonal (or unbalanced) Factorial Between-Subjects ANOVA

First, hello, yes, the blog series is back. I got a little sidetracked actually teaching the course and then didn’t pick it up last summer, and then the cycle repeated itself. Part of the reason is also that I hit something of a snag along the way and was unable to get SPSS and R to produce the same answers for certain analyses, and it took me a while to figure out what was going on. Now that I believe I understand the problem (and the solution to it), it’s time to plunge back in.

So, now, the situation: non-orthogonal (or unbalanced) factorial between-subjects ANOVA. That alone is quite a mouthful. What that means is a two (or more) way ANOVA with one observation per subject where the number of subjects in each cell is not equal.

For most people who might read this blog (all three of you), the implications of this are generally pretty well-known: it means the different ANOVA effects are not statistically independent of each other. This means that something generally has to be done to correct. This is a bit of a controversial point in the statistics universe. Many statistical purists (and, by my reading of it, most of the R community) think the right approach to this problem is to use Type II sums of squares. There is another contingent, however—and most psychologists, as far as I can tell—that argues that Type III sums of squares are the way to go. This is the kind of thing that generates name-calling on statistical Web sites (particularly from the R side).

I am not interested in having this argument here. I just care about getting the statistics software to produce whichever analysis is preferred by the analyst. It seems to me this ought to be the focus of discussions on the software side of things, but that doesn’t seem to be how it usually plays out. Regardless, I’m going to present both without discussion of the pros and cons of the different approaches, just the “how to” part.

And the good news is that both packages can do this, although some of the default behaviors by the two packages are somewhat odd. I’ll get to the details soon.

The data for these examples is in nonorth.xls, which has three variables: “Dep,” the dependent variable, and “A” and “B,” the two independent variables. Noting fancy, no cover story, just the data.

**Basic ANOVA and Sums of Squares**

First, in SPSS. The default behavior in SPSS is to produce Type III sums of squares. This is entirely reasonable given that this is still the majority viewpoint in the social sciences—or, at least, it’s what I was taught in grad school in Psycholgy in the early 1990s, anyway. So, the simplest command:

`UNIANOVA Dep BY A B.`

produces this ANOVA table:

And it even nicely notes that it’s using Type III. Great so far. If you’d prefer Type II sums of squares, you just tell SPSS to override the default:

`UNIANOVA Dep BY A B`

/METHOD = SSTYPE(2).

Note that it is possible to request Type I sums of squares, but nobody uses that. Anyway, the output:

Notice that the sums of squares for both A and B have increased, as you’d expect with Type II sums of squares. Very straightforward.

Now, R. Despite the fact that nobody I can find anywhere actually recommends the use of Type I sums of squares, that’s the default in R. Yes, it defaults to a behavior that nobody wants—in what universe that makes sense, I’m sure I don’t know. So, assuming the data frame is already attached, first create factors and then run the ANOVA with this code:

`AF = factor(A)`

BF = factor(B)

summary(aov(Dep ~ AF * BF))

The following output is produced:

Df Sum Sq Mean Sq F value Pr(>F) AF 1 11438 11438 42.846 1.20e-09 *** BF 3 9791 3264 12.226 4.14e-07 *** AF:BF 3 496 165 0.619 0.604 Residuals 132 35237 267 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Nowhere in the output does it identify that these are Type I sums of squares, but if you reverse the order of the factors in the “aov()” command, you do indeed get different output. That is, this code:

`summary(aov(Dep ~ BF * AF))`

yields this ANOVA table:

Df Sum Sq Mean Sq F value Pr(>F) BF 3 17521 5840 21.879 1.45e-11 *** AF 1 3707 3707 13.888 0.000287 *** BF:AF 3 496 165 0.619 0.603611 Residuals 132 35237 267 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The sums of squares for A and B change dramatically, as expected. All of this is well and good, but not very helpful, since (almost) nobody actually uses Type I sums of squares. In order to get Type II or III sums of squares, the “Anova()” command—note the capital “A”—from the “car” package is required. In addition, certain options in R must be set in order for the internal contrasts used by the system to generate the correct output, so again, the defaults produce things other than what is desired. (This one is a little more defensible, since internally R defaults to dummy coding of categorical variables, which is not a wholly unreasonable choice, but it does mess up SS computations in unbalanced ANOVAs when done that way.)

So, first, load the “car” package, set the relevant options, create an aov object, and then run that through “Anova()” with the options set:

`library(car)`

options(contrasts = c("contr.sum","contr.poly"))

Anv = aov(Dep ~ AF * BF)

Anova(Anv, type=2)

This yields output that matches SPSS (whew):

Anova Table (Type II tests) Response: Dep Sum Sq Df F value Pr(>F) AF 3707 1 13.8883 0.0002866 *** BF 9791 3 12.2259 4.137e-07 *** AF:BF 496 3 0.6194 0.6036113 Residuals 35237 132 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

It does also note that it’s using Type II sums of squares. To get Type III, the “type” parameter in “Anova()” is just changed appropriately:

`Anova(Anv, type=3)`

And that produces output that again matches SPSS:

Anova Table (Type III tests) Response: Dep Sum Sq Df F value Pr(>F) (Intercept) 31112 1 116.5482 < 2e-16 *** AF 1192 1 4.4656 0.03647 * BF 2892 3 3.6113 0.01509 * AF:BF 496 3 0.6194 0.60361 Residuals 35237 132 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

So while it is mildly irritating that R has a useless default, it’s at least not actually difficult to get it to produce something reasonable.

**Contrasts**

This is where things get a little more interesting. Let’s say that independent variable B is an interval variable so that polynomial contrast analysis is a reasonable approach. SPSS can generate this pretty easily:

`UNIANOVA Dep BY A B`

/CONTRAST(B) = POLYNOMIAL.

What’s unfortunate is that the output is so dumb. Here’s what the relevant part of the output looks like:

Hey, SPSS, how about the actual *t*-values? It doesn’t print the relevant *F*-values, either, just the *p*-values. Grr. The corresponding t-values are -1.27 for the linear contrast, -4.55 for the quadratic, and 3.01 for the cubic. (Fs would be the squares of those: 1.61, 20.71, and 9.06.)

What’s much more interesting about this is that SPSS doesn’t care what you set for the sum of squares method. If you set the sum of squares method for Type I, the *contrast output is not affected*. So it’s not sensitive to what’s being done in terms of variable entry.

Now, let’s take a look at what R does with the same basic idea. This R code first makes sure that the B factor uses polynomial contrasts, then sets up labels for those contrasts, runs the ANOVA, and prints out a summary table with the B factor split and labelled:

`contrasts(BF) = contr.poly(4)`

CondLab = list("Linear"=1, "Quad"=2, "Cubic"=3)

AnvL1 = aov(Dep ~ BF*AF)

summary(AnvL1, split=list(BF=CondLab))

Here’s the output:

Df Sum Sq Mean Sq F value Pr(>F) BF 3 17521 5840 21.879 1.45e-11 *** BF: Linear 1 4668 4668 17.486 5.24e-05 *** BF: Quad 1 7797 7797 29.206 2.94e-07 *** BF: Cubic 1 5057 5057 18.943 2.68e-05 *** AF 1 3707 3707 13.888 0.000287 *** BF:AF 3 496 165 0.619 0.603611 BF:AF: Linear 1 68 68 0.254 0.615356 BF:AF: Quad 1 75 75 0.282 0.596155 BF:AF: Cubic 1 353 353 1.322 0.252225 Residuals 132 35237 267 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Hmm, none of those match SPSS at all. If you look carefully you’ll note that the sums of squares for the three B contrasts do in fact add up to the sum of squares for B, but this is the Type I sum of squares for B, with B added first. These numbers are inflated! This is probably pretty dangerous.

Note that, because of this Type I computation, R is sensitive to the order in which variables are listed in the model specification. So, this code:

`AnvL2 = aov(Dep ~ AF*BF)`

summary(AnvL2, split=list(BF=CondLab))

produces very different results:

Df Sum Sq Mean Sq F value Pr(>F) AF 1 11438 11438 42.846 1.20e-09 *** BF 3 9791 3264 12.226 4.14e-07 *** BF: Linear 1 1 1 0.004 0.950627 BF: Quad 1 6330 6330 23.711 3.15e-06 *** BF: Cubic 1 3461 3461 12.963 0.000448 *** AF:BF 3 496 165 0.619 0.603611 AF:BF: Linear 1 68 68 0.254 0.615356 AF:BF: Quad 1 75 75 0.282 0.596155 AF:BF: Cubic 1 353 353 1.322 0.252225 Residuals 132 35237 267 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

These aren’t wrong in the sense that R is making some kind of computational error, but again, these are probably not the results you want to report in a journal article.

So, how do you get R to produce something that makes more sense? And what exactly is SPSS doing internally to get those numbers? It took me a while to figure this out. I recoded the thing as a straightforward multiple regression in SPSS and tried doing it a whole bunch of different ways in terms of order of variable entry—I assumed it was doing some kind of change in F if the term is added last kind of thing—and then hit it completely by accident when I just ran the regression with all variables in it. It’s nice and clean in R as a simple regression:

`AnvL3 = lm(Dep ~ AF * BF)`

summary(AnvL3)

Note the use of “lm()” there instead of an ANOVA command. Since all our independent variables now have a single degree of freedom, it’s just a regression, right? And that’s when I realized what SPSS is doing internally, as here’s the result:

Call: lm(formula = Dep ~ AF * BF) Residuals: Min 1Q Median 3Q Max -39.120 -11.779 2.146 11.190 45.396 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 94.779 1.872 50.643 < 2e-16 *** AF1 -6.397 1.872 -3.418 0.000838 *** BF.L -4.918 3.874 -1.270 0.206471 BF.Q -17.023 3.743 -4.548 1.21e-05 *** BF.C 10.867 3.608 3.012 0.003108 ** AF1:BF.L 2.309 3.874 0.596 0.552109 AF1:BF.Q 1.242 3.743 0.332 0.740653 AF1:BF.C 4.149 3.608 1.150 0.252225 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 16.34 on 132 degrees of freedom Multiple R-squared: 0.3814, Adjusted R-squared: 0.3486 F-statistic: 11.63 on 7 and 132 DF, p-value: 1.857e-11

If you’re not sure what this means, look again at the t-values for B, noting that “.L” denotes linear, “.Q” quadratic, and so on. These match what SPSS produces (though of course SPSS doesn’t print out the actual t-values). So, what SPSS is doing internally is simultaneous entry of all terms into a regression, and tests of the coefficients for the contrast variables. This is actually a pretty reasonable thing to do, but is completely opaque to the SPSS user and is not explained anywhere that I could find in the SPSS documentation.

Note that in an unbalanced ANVOA, this is almost certainly the right thing to do, too, since the Type I approach is potentially more than a little permissive depending on the order of entry of the variables. So, if you’re going to use R (perhaps migrating from SPSS like me), then this is probably how contrasts should be handled if the ANOVA is unbalanced.

Note that neither package produces effect sizes (e.g., eta squared, Cohen’s ƒ) for contrasts. Bummer.

Next time: one-way repeated measures. There are many ways to do this in R, but some are better than others…