Hi all,
I was trying to use glht() from multcomp package to construct a contrast on
interaction term
in a linear model to do some comparisons. I am little uncertain on how to
construct contrasts on a 3-way interaction containing a continuous variable, and
hope someone can confirm what I did is correct or wrong:
The linear model has a continuous dependent
variable “y”, with treatment factor “Trt” with value 0 and 1, a factor variable
“A” with value 0 and 1, a continuous variable “x”.
A simpler model is:
set.seed(10)
dat <- cbind(y=c(rnorm(10,3),rnorm(10,4),rnorm(10,3.1),rnorm(10,6)),
x=runif(40,5,15), expand.grid(A=rep(factor(0:1),each=10),Trt=factor(0:1)))
fit <- lm(y
~ x + Trt * A,dat)
My purpose is to test whether treatment effect is significant given
each level of factor A, so I used contrasts:
library(multcomp)
K <- rbind(c(0,0,1,0,0), c(0,0,1,0,1))
rownames(K) <- c('Trt 1-0|A=0','Trt 1-0|A=1')
colnames(K) <- names(coef(fit))
K
(Intercept) x Trt1 A1 Trt1:A1
Trt 1-0|A=0 0 0 1 0 0
Trt 1-0|A=1 0 0 1 0 1
(glht.fit <- summary(glht(fit, linfct K),
test=adjusted(type='none')))
Linear
Hypotheses:
Estimate Std. Error t value
Pr(>|t|)
Trt 1-0|A=0 =0 -0.2720 0.3616 -0.752 0.45701
Trt 1-0|A=1 =0 1.0690 0.3564 2.999 0.00496 **
Now I suspect independent variable “x” may play a role in the treatment
effect at each level of A, so I would like to add in a 3-way interaction
between Trt, A and x:
fit <- lm(y
~ x * Trt * A,dat)
If my purpose is to test whether treatment is significant at each level
of factor A and certain value of covariate “x”, for example, when x=10, would
following code give me what I wanted?
K <- rbind(c(0,0,1,0,10,0,0,0), c(0,0,1,0,10,0,1,10))
rownames(K) <- c('Trt 1-0|A=0 x=10','Trt
1-0|A=1 x=10')
colnames(K) <- names(coef(fit))
K
(Intercept) x Trt1 A1 x:Trt1 x:A1 Trt1:A1 x:Trt1:A1
Trt 1-0|A=0 x=10 0 0 1 0 10 0 0 0
Trt 1-0|A=1 x=10 0 0 1 0 10 0 1 10
(glht.fit <- summary(glht(fit, linfct K),
test=adjusted(type='none')))
Linear
Hypotheses:
Estimate Std. Error t
value Pr(>|t|)
Trt 1-0|A=0
x=10 == 0 -0.3526 0.3254 -1.083 0.286731
Trt 1-0|A=1
x=10 == 0 1.4621 0.3328 4.394 0.000115 ***
So the above test was testing whether treatment effect is significant at each
level of factor A when x=10, am I correct?
Appreciate if someone would confirm this?
Thanks
John
[[alternative HTML version deleted]]
