[R] puzzle about contrasts
Peter Dalgaard
P.Dalgaard at biostat.ku.dk
Tue Sep 9 18:35:05 CEST 2008
Prof Brian Ripley skrev:
> -0.5*(A+B) is not a contrast, which is the seat of your puzzlement.
>
> All you can get from y ~ x is an intercept (a column of ones) and a
> single 'contrast' column for 'x'.
>
> If you use y ~ 0+x you can get two columns for 'x', but R does not
> give you an option of what columns in the case: see the source of
> contrasts(). So you would need to replace contrasts(), which I think
> will be hard as model.matrix.default will look in the 'stats'
> namespace. It would probably be easier to create the model matrix
> yourself.
>
Or accept the default and do the parameter transformations yourself.
l <- lm(y~x)
T <- rbind(
c(-1,-.5),
c(0,1))
c2 <- T%*%coef(l)
V2 <- T%*%vcov(l) %*% t(T)
cbind(coef=c(c2), s.e.=sqrt(diag(V2)))
> On Tue, 9 Sep 2008, Kenneth Knoblauch wrote:
>
>> 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
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
> ------------------------------------------------------------------------
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
--
O__ ---- Peter Dalgaard Øster Farimagsgade 5, Entr.B
c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
More information about the R-help
mailing list