Trying to answer this question brought up my own question on the workings
of the C function for defining contrasts. Others in the group, please
skim to the bottom of my response to look at my question. Thanks.
On Fri, 15 Feb 2002, Rene Holst wrote:
> Hi'
>
> S+ 4.5 R.2, Win2K
>
> I want to fit a model y~-1+a/x, where a is a factor with three levels and x
> is a continuous variable. The model corresponds to fitting a line for each
> level of a:
>
> a1: b0.1+b1.1*x
> a2: b0.2+b1.2*x
> a3: b0.3+b1.3*x
>
> How can I fit the model subject to the following constraints
> 1) b1.1=b1.3
> 2) b1.1=b1.2=b1.3
> 3) b0.1=b0.2=0
To put these types of constraints on the coefficients you need to think
about how your final x-matrix could change to impose this. For example in
the first constraint you want 2 of the columns to have the same
coefficient, this is the same as multiplying that coefficient by the sum
of the 2 columns, so if you can make S-PLUS create an x-matrix with 1
column equaling the sum of the 2 original columns, then you have your
constraint.
The following commands illustrate a way to do this:
a <- factor(rep(c("a","b","c"), each=3))
x <- rep(1:3, 3)
y <- c(1,2,3, 1.2,2.3,3.4, 1,1.9,2.8)
plot(x,y, type="n")
text(x,y, as.character(a))
abline(0,1)
abline(.1,1.1)
abline(.1,.9)
options(contrasts=c("contr.treatment","contr.poly"))
temp1 <- lm(y~-1+a/x,x=T)
temp1
temp1$x
temp2 <- lm(y ~ -1 + a + a:x, x=T)
temp2
temp2$x
# constraint 1, b1.1=b1.3
b <- cbind( c(1,0,1), c(0,1,0) )[codes(a),]
b
temp3 <- lm( y~ -1 + a + b:x, x=T)
temp3
temp3$x
tempc3 <- coef(temp3)
plot(x,y, type="n")
text(x,y, as.character(a))
abline(tempc3[1], tempc3[4])
abline(tempc3[2], tempc3[5])
abline(tempc3[3], tempc3[4])
# constraint 2, b1.1=b1.2=b1.3
temp4 <- lm( y~ -1 + a + x, x=T)
temp4
temp4$x
tempc4 <- coef(temp4)
plot(x,y, type="n")
text(x,y, as.character(a))
abline(tempc4[1], tempc4[4])
abline(tempc4[2], tempc4[4])
abline(tempc4[3], tempc4[4])
# constraint 3, b0.1=b0.2=0
b <- cbind( c(0,0,1) )[codes(a),]
temp5 <- lm( y~ -1 + b + a:x, x=T)
temp5
temp5$x
tempc5 <- coef(temp5)
plot(x,y, type="n")
text(x,y, as.character(a))
abline(0, tempc5[2])
abline(0, tempc5[3])
abline(tempc5[1], tempc5[4])
Now my question for the group:
I originally thought that this could also be done using the contrasts or C
functions, but if I try the following code (S-PLUS 6.0.3, R 1.3.1)
b <- a
contrasts(b,2) <- cbind( c(1,0,1), c(0,1,0) )
attributes(b)
temp6 <- lm( y ~ -1 + a + b:x, x=T)
temp6$x
then there is no evidence that the contrasts were ever used, I expected
the same x-matrix as temp3$x above. Why are the contrasts not used in
creating the interaction terms?
I thought that maybe this was just something to do with the code that
generates the interactions so I tried:
temp7 <- lm( y ~ -1 + b + b:x, x=T, singular.ok=T)
temp7$x
Now in S-PLUS 6 the intercept has come back despite the -1 (needing the
singular.ok=T to compute and not giving the desired output, though
dummy.coef does give the answers), and in R 1.3.1 the contrasts are again
ignored.
Is this a bug, or am I missing a reason why the contrasts should be
working the way they are?
Thanks,
--
Greg Snow, PhD Office: 223A TMCB
Department of Statistics Phone: (801) 378-7049
Brigham Young University Dept.: (801) 378-4505
Provo, UT 84602 email: gls@byu.edu
|