[R] How to test combined effects?

Gang Chen gangchen at mail.nih.gov
Tue Oct 30 22:26:36 CET 2007


Dieter,

Thank you very much for the help!

I tried both glht() in multcomp and estimable() in gmodels, but  
couldn't get them work as shown below. Basically I have trouble  
specifying those continuous variables. Any suggestions?

Also it seems both glht() and estimable() would give multiple t  
tests. Is there a way to obtain sort of partial F test?


 > ModelFit<-lme(mct~ IQ*age+IQ*I(age^2)+IQ*I(age^3), MyData,  
random=~1|ID)
 > anova(ModelFit)

                mDF denDF  F-value p-value
(Intercept)     1   257 54393.04  <.0001
IQ              1   215     3.02  0.0839
age             1   257    46.06  <.0001
I(age^2)        1   257     8.80  0.0033
I(age^3)        1   257    21.30  <.0001
IQ:age          1   257     1.18  0.2776
IQ:I(age^2)     1   257     0.50  0.4798
IQ:I(age^3)     1   257     0.23  0.6284

 > library(multcomp)
 > glht(ModelFit, linfct = c("IQ:age = 0", "IQ:I(age^2) = 0", "IQ:I 
(age^3) = 0"))
Error in coefs(ex[[3]]) :
   cannot interpret expression 'I''age^2' as linear function

 > library(gmodels)
 > estimable(ModelFit, rbind('IQ:age'=0, 'IQ:I(age^2) = 0', 'IQ:I 
(age^3) = 0'))
Error in FUN(newX[, i], ...) :
   `param' has no names and does not match number of coefficients of  
model. Unable to construct coefficient vector

Thanks,
Gang


On Oct 30, 2007, at 9:08 AM, Dieter Menne wrote:


> Gang Chen <gangchen <at> mail.nih.gov> writes:
>
>
>>
>> Suppose I have a mixed-effects model where yij is the jth sample for
>> the ith subject:
>>
>> yij= beta0 + beta1(age) + beta2(age^2) + beta3(age^3) + beta4(IQ) +
>> beta5(IQ^2) + beta6(age*IQ)  + beta7(age^2*IQ)  +  beta8(age^3 *IQ)
>> +random intercepti + eij
>>
>> In R how can I get an F test against the null hypothesis of
>> beta6=beta7=beta8=0? In SAS I can run something like contrast  age*IQ
>> 1, age^2*IQ 1, age^3 *IQ 1, but is there anything similar in R?
>>
>
> Check packages multcomp and gmodels for contrast tests that work  
> with lme.
>
> Dieter
>



More information about the R-help mailing list