[R] puzzle about contrasts

Kenneth Knoblauch ken.knoblauch at inserm.fr
Tue Sep 9 17:27:17 CEST 2008


Hi,

I'm trying to redefine the contrasts for a linear model.
With a 2 level factor, x, with levels A and B, a two level
factor outputs A and B - A from an lm fit, say
lm(y ~ x). I would like to set the contrasts so that
the coefficients output are -0.5 (A + B) and B - A,
but I can't get the sign correct for the first coefficient
(Intercept).

Here is a toy example,

set.seed(12161952)
y <- rnorm(10)
x <- factor(rep(letters[1:2], each = 5))
##  so  A and B =
tapply(y, x, mean)

         a          b
-0.7198888  0.8323837

## and with treatment contrasts
coef(lm(y ~ x))  ## A and B - A

(Intercept)          xb
  -0.7198888   1.5522724

Then, I try to redefine the contrasts

### would like contrasts: -0.5 (A + B) and B - A
D1 <- matrix( c(-0.5, -0.5,
		-1, 	1),
		2, 2, byrow = TRUE)
C1 <- solve(D1)
Cnt <- C1[, -1]
contrasts(x) <- Cnt
coef(lm(y ~ x))

(Intercept)          x1
  0.05624745  1.55227241

but note that the desired value is
-0.5 * sum(tapply(y, x, mean))

[1] -0.05624745

I note that the first column of C1 is -1's not +1's
and that working by hand, if I tamper with the model matrix

mm <- model.matrix(y ~ x)
mm[, 1] <- -1

mm
    (Intercept)   x1
1           -1 -0.5
2           -1 -0.5
3           -1 -0.5
4           -1 -0.5
5           -1 -0.5
6           -1  0.5
7           -1  0.5
8           -1  0.5
9           -1  0.5
10          -1  0.5
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$x
   [,1]
a -0.5
b  0.5

solve(t(mm) %*% mm) %*% t(mm) %*% y  ##Yes, I know. Use QR
                    [,1]
(Intercept) -0.05624745
x1           1.55227241

gives the correct sign.

So, I guess my question reduces to how one would set the
contrasts for the model.matrix to be correct
for this to work out correctly?

Thank you.

Ken


-- 
Ken Knoblauch
Inserm U846
Institut Cellule Souche et Cerveau
Département Neurosciences Intégratives
18 avenue du Doyen Lépine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: +33 (0)6 84 10 64 10
http://www.sbri.fr



More information about the R-help mailing list