Setting up contrasts for means in R can be very simple, or it can me much more complex than you would expect. It is quite simple to compare each mean with the one before it or compare all means to a designated base mean. It can be almost as easy to create and test a set of orthogonal mean contrasts. But it gets quite a bit messier if you have a set of contrasts that are not orthogonal.
If you consider what you do in SPSS or SAS (a comparison that usually upsets true R devotees), you do some of the same things for some contrasts. For example, in SPSS you can have a statement that reads "/CONTRAST(Group)=Helmert", or "/CONTRAST(Group)=Simple", etc. Helmert contrasts compare the 1st group with the all subsequent levels, the 2rd with the average of the levels 3, 4, 5, etc., and so on. (R reverses the ordering of groups (1st with 2nd, 1st and 2nd with 3rd, etc.), "/CONTRAST(Group) = Simple" compares each group with a reference group. R allows basically the same thing. You can simply state " contrasts(dat$Group) = contr.helmert" to get helmert contrasts. Similar statements will give you other standard contrasts. So that part is easy. You just embed such statements in
contrasts(dat$Group) <- contr.helmert model1 <- lm(Time~Group, data = dat) summary(model1)
and get the result you seek.
Perhaps I just work in the wrong field, or perhaps I am just a crab, but I find that whole set of possible contrasts, such as Helmert, to be less than helpful. In my experience, the ordering of treatments is often quite arbitrary, and I can't say that the second should be compared to the first or the third compared to the average of the first two. That might happen, but I haven't yet seen a case. The same applies to similar kinds of contrasts. In my experience a set of five means might have two, three, or four useful specific contrasts, but I want to be able to specify what they are and not let the program do that for me or lock me in to a specific set.
There are a couple of choices. In SPSS, SAS, and R you can specify a set of orthogonal contrasts. That is certainly getting closer to my ideal.
In R you can do so with something like
c1 <- c(1/3, 1/3, 1/3, -.5, -.5)
c2 <- c(1,1,-2,0,0)
c3 <- c(1,-1,0,0,0)
c4 <- c(0,0,0,1,-1)
mat <- cbind(c1,c2,c3,c4)
contrasts(dat$Group) <- mat
model1 <- lm(Time ~ Group, data = dat)
summary(model1, split=list(Group=list("First 3 vs last 2"=1, "First and Second vs Third" = 2,
"One vs two"=3, "Four vs Five" = 4)))
The first four lines specify the groups to be compared, the fifth combines them into a matrix, the 6th tells R that the matrix gives the contrasts you want (sort of), and the next two lines give you your result. That's pretty easy and makes far more sense to me than to select one of the previous versions. You do something very similar to this in SPSS or SAS. But in R there is a catch. This system only works if the contrasts are orthogonal--i.e. independent. If they are not orthogonal you will get a result, but it is not the result you thought that you were getting.
So what is the problem? Well in SPSS or SAS, the program takes your matrix and automatically adjust it to be what you need. You don't see this happening, but it happens anyway. That's good. But R takes itself more literally. It uses the matrix you give it, but it doesn't make a necessary adjustment because you didn't tell it to and it only does what you tell it. Now this doesn't make any real difference if your contrasts are orthogonal, but it does if they aren't.
What SPSS really does is to take the matrix of contrasts that you give it, augment it temporaily, takes the inverse of that matrix, deaugments that (if that is a word) and works with that inverse. It just happens by magic. If the original matrix was orthogonal, the inverse version is just a linear function of that and it won't make any difference to the outcome which matrix you work with. So you specify the code that I gave above and you will get the right answer. But if the matrix is not a set of orthogonal coefficients, working directly with that matrix will give you the wrong answers for some of the contrasts. And they will look just fine, leaving you to think that you got what you wanted. So what do we do? We get the inverse ourselves.
Let's start with a matrix of contrasts that are not orthogonal. Suppose that we have the same five means that we used above, but we set up a different matrix of contrasts. We can do this as follows (or we could use the pattern above if you prefer).
You may wonder about the first row of that command, but it sets up a column that represents the grand mean. (you could use 1,1,1,1,1 in place of the 1/5s if yu prefer) It is there largely to make the matrix square so that we can take its inverse. We do that by So what is all that. The first command takes the inverse (R likes the word "solve") of the augmented matrix. The next line removes the first column from that matrix--what I called "deaugmenting". The third line sets that new matrix up as the specified contrast matrix. This is what SPSS does behind the scenes where you don't see it. Now we are ready to find out what differences are significant. But we will do so with an actual example. From page 381 in my book Statistical Methods for Psychology, 8th ed. we can read in the data for five groups. Since writing this page I have found a very nice discussion of the issue at http://sites.stat.psu.edu/~jls/stat512/lectures/ContrastsInR.pdf. I believe that the author is Joe Shafer at Penn State, and it is quite a good document. That's all there is to it. And don't think that you need four contrasts just because you have four df for Groups. Don't run things just because you can--run them because they ask relevant questions. The more tests you run, the more you increase the probability of a type I error. Why do that?
LInv <- solve(t(L))
Nonorth.Contrasts <- LInv[,2:5]
contrasts(dat$Group) = Nonorth.Contrasts
###______________________________________________________
# The following gives correct solution for nonorthogonal design
# Notice that it uses the Linverse
# I need to use the Linverse with nonorthogonal designs
setwd ("www.uvm.edu/~dhowell/StatPages/methods8/DataFiles")
dat <- read.table("Tab12-1.dat", header = T)
dat$Group <- factor(dat$Group)
# I am specifying the contrast matrix as L and its inverse as LInv
L <- matrix(c(1,1,1,1,1,-.333, -.333, -.333, .5, .5, 0, -1, 0,0,1, -1,0,1,0,
0,0,1,-1,0,0), ncol = 5)
Linv <- solve(t(L))
new.contrasts <- Linv[,2:5]
contrasts(dat$Group) = new.contrasts
model5 <- lm(Time~Group, data = dat )
summary(model5,split=list(Group=list("First 3 vs last 2"=1, "Second vs Fifth" = 2,
"First vs Third"=3, "Second vs Third" = 4)))
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.6000 0.8944 17.441 < 2e-16 ***
Group1 18.1594 1.8250 9.950 9.66e-12 ***
Group2 19.0000 2.8284 6.718 8.87e-08 ***
Group3 7.0000 2.8284 2.475 0.0183 *
Group4 -1.0000 2.8284 -0.354 0.7258
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 5.657 on 35 degrees of freedom
Multiple R-squared: 0.7574, Adjusted R-squared: 0.7297
F-statistic: 27.32 on 4 and 35 DF, p-value: 2.443e-10
dch:
David C. Howell
University of Vermont
David.Howell@uvm.edu