[R] How to test combined effects?

Gang Chen gangchen at mail.nih.gov
Fri Nov 2 17:20:15 CET 2007


Thanks a lot for the help!


> First, if you would like to performa an overall test of whether the  
> IQ interactions are necessary, you may find it most useful to use  
> anova to compare a full and reduced model.  Something like:
>
> 	ModelFit.full <-lme(mct~ IQ*age+IQ*I(age^2)+IQ*I(age^3), MyData,   
> random=~1|ID)
> 	ModelFit.reduced <-lme(mct~ IQ + age+I(age^2)+I(age^3), MyData,   
> random=~1|ID)
> 	anova(ModelFit.full, ModelFit.reduced, test="F")


I had done this before, but it seems that I would only get a  
likelihood ratio test, not a partial F test, out of the anova. The  
'test' option with logical value in anova seems to take TRUE or  
FALSE, thus either I get a likelihood ration test or not. I guess a  
likelihood ratio test with ML method is legitimate in this context?


> Second, you don't have the syntax right for estimable().  As  
> described and shown by example in the manual page.  The correct  
> syntax is:
>
> 	library(gmodels)
> 	estimable(ModelFit, c('IQ:age'=1, 'IQ:I(age^2)'= 1, 'IQ:I(age^3)'  
> = 1))
>
> Note the pattern of quoted name, followed by '=', and then the  
> value 1 (not zero).  This will perform a single joint test whether  
> these three coefficients are zero.


Thanks for catching the errors! Yes this works. However it gives a t  
test with 1 degree of freedom, not exactly a partial F test. So does  
it mean that this only tests the average effect of those 3 terms? If  
so, that would be slightly different from the partial F test I was  
looking for, no?

Gang


> -G
>
>
>
> On Oct 30, 2007, at 5:26PM , Gang Chen wrote:
>
>> 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
>>>
>>
>> ______________________________________________
>> 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.
>



More information about the R-help mailing list