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

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:

wpid-pastedgraphic-2014-06-13-14-491.png

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:

wpid-pastedgraphic1-2014-06-13-14-49.png

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:

wpid-pastedgraphic-2014-06-13-14-49.png

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…

Advertisements
Categories: Uncategorized Tags: , ,
  1. No comments yet.
  1. No trackbacks yet.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: