[R] How to test combined effects?
Gregory Warnes
gregory.warnes at mac.com
Fri Nov 2 18:35:15 CET 2007
On Nov 2, 2007, at 12:20PM , Gang Chen wrote:
> 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?
I'll have to defer to others on the validity of the LR test in this
context, as I've got my programmer hat on today and not my stat
theory hat ... ;-)
>
>
>> 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?
Yes, this provides a 1-degree of freedom contrast between an
individual with ('IQ:age'=0, 'IQ:I(age^2)'=0, 'IQ:I(age^3)'=0) and
another with ('IQ:age'=1, 'IQ:I(age^2)'=1, 'IQ:I(age^3)'=1).
Although we do provide a joint Wald test for 'lm', 'glm', 'aov',
'gee' or 'geese' objects, we have not done so for 'lme' objects.
The code to do so is straightforward, and I'll send you a copy
privately.
Perhaps someone who has the stat theory hat on today can comment on
the validity of the Wald X^2 here, and on the relative merits of the
LR test vs the Wald test for fixed effects of an LME.
In any case here is how you use the enhanced code.
> library(nlme)
> library(gmodels)
> Orthodont$Rand <- rnorm(nrow(Orthodont)) # add a noise variable
> fm2 <- lme(distance ~ age * Rand + Sex * Rand, data = Orthodont,
random = ~ 1)
> cmat <- cbind( "Rand"=c(1,0,0), "age:Rand"=c(0,1,0),
"Rand:SexFemale"=c(0,0,1) )
> cmat
Rand age:Rand Rand:SexFemale
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
> estimable( fm2, cmat, joint.test=FALSE ) # individual contrasts
Estimate Std. Error t value DF Pr(>|t|)
(0 0 1 0 0 0) 0.39919525 0.9263579 0.4309298 77 0.6677233
(0 0 0 0 1 0) -0.07152143 0.0783120 -0.9132883 77 0.3639418
(0 0 0 0 0 1) 0.47758669 0.2993967 1.5951635 77 0.1147722
> attach(environment(estimable)) # make the estimable functions
available to local code
The following object(s) are masked from package:gmodels :
CrossTable ci ci.binom coefFrame estimable fast.prcomp fast.svd
fit.contrast glh.test make.contrasts print.glh.test summary.glh.test
> source("~/src/r-gregmisc/gmodels/R/estimable.R") # load the
modified code
> estimable( fm2, cmat, joint.test=TRUE ) # joint Wald X^2
X2.stat DF Pr(>|X^2|)
1 4.876038 3 0.1811025
> q()
-Greg
> 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