Use a baseline including all observations in lm()

When using factors in a linear model, I would like to retrieve the 'average' effect as baseline (intercept), rather than level 1 of all factors. Is this possible?

tab <- data.frame(gl(2,5),rnorm(10)+as.numeric(gl(2,5)))
colnames(tab) <- c('pred','dat')
mylm <- lm(dat~pred,data=tab)
coefficients(mylm)


returns something like

(Intercept)       pred2
0.5879712   1.8542153


coefficients are on average 1 for both intercept and pred2. I would like to rather have intercept = 1.5, pred1 = -0.5 and pred2 = +0.5

model.matrix(mylm)


returns

         (Intercept) pred2
1            1     0
2            1     0
3            1     0
4            1     0
5            1     0
6            1     1
7            1     1
8            1     1
9            1     1
10           1     1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$pred [1] "contr.treatment"  I would like a matrix like this  (Intercept) pred1 pred2 1 1 1 0 2 1 1 0 3 1 1 0 4 1 1 0 5 1 1 0 6 1 0 1 7 1 0 1 8 1 0 1 9 1 0 1 10 1 0 1  Does that make sense? if not, why? Is there some command to generate such a design matrix? I've been looking for a way to do this without success... - Why not just subtract the mean effect from the coefficients? – whuber Feb 6 '14 at 15:45 Because it seems to me that this would work only if there is the same amount of observations for each predictor level. Should I then weigh this mean effect by the number of observations? How about if I have also covariates? would that change something? I don't feel very comfortable doing this 'by hand'... – Max Feb 6 '14 at 20:36 What, exactly, do you mean by "average effect"? Is it the average over the effect coefficients or the average over the data? Take your pick. When there are covariates, there is no problem: the effects are still the effects. If, however, there are interactions between them and your effects, then the effects depend on levels of the covariates and you will have to average them over all the data to obtain anything meaningful. – whuber Feb 6 '14 at 21:01 Well, I mean the intercept, actually. The average of dat across observations, independently of pred... – Max Feb 7 '14 at 13:05 If pred were a continuous variable, I would get an intercept that would be the modeled level of dat when pred == 0. Why can't I here retrieve the modeled level of dat for all (i.e. the average of all) levels of pred? – Max Feb 7 '14 at 13:11 2 Answers It does not make sense to have such a design matrix because the columns are linearly dependent (specifically, column 1 = column 2 + column 3) so you cannot compute the OLS estimator, which requires inversion of$X'X$, where$X$is said design matrix. Your proposed design matrix falls under what is sometimes called the "dummy variable trap". What you can do, if you want to include both 'pred1' and 'pred2', is to drop the intercept instead, but you cannot have exhaustive dummies and an intercept. Addition To do this in R, and with your example code, run the following: lm(dat ~ -1 + pred, data=tab)  The '-1' removes the intercept and then you interpret the coefficients as the mean in the respective group, that is, there is no baseline group. - Ok, I understand. Now if I want the effect of pred1, how can I do that? – Max Feb 5 '14 at 14:33 @Max I've updated with a command. – Student001 Feb 6 '14 at 9:34 Your coding does not make sense. It sound like you want deviation coding (or contrast) coding rather than traditional dummy coding. The UCLA website has code to create this: http://www.ats.ucla.edu/stat/r/library/contrast_coding.htm Or assuming you have a four level variable (unlike your example): #the contrast matrix for categorical variable with four levels > contr.sum(4) [,1] [,2] [,3] 1 1 0 0 2 0 1 0 3 0 0 1 4 -1 -1 -1 #assigning the deviation contrasts to pred.f contrasts(tab$pred.f) = contr.sum(4)
#the regression
model1=lm(dat ~ pred.f, tab)
summary(model1)


The annoying thing about R is that it will only show you the coefficients for pred2. To get all the coefficients you then have to either use the fact that the coefficients for deviation coding will add up to zero or try.

dummy.coef(model1)

-
I don't fully understand this contr.sum(). I see it allows to extract the effect specific to each of the factor levels, but why isn't there a way to retrieve the average of all levels (i.e. the 'real' intercept)? – Max Feb 7 '14 at 13:16
I might be misunderstanding you. Sorry. But my understanding is the 'average' effect is the intercept which you should get after summary(model1). – charles Feb 7 '14 at 20:55