[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