[R] a question about glht function

[Achim Zeileis]

On Tue, 26 Jul 2011, Peter Ehlers wrote:

> On 2011-07-26 01:13, Achim Zeileis wrote:
>> 
>> 
>> On Tue, 26 Jul 2011, Lao Meng wrote:
>> 
>>> Hi all:
>>> There's a question about glht function.
>>> 
>>> My data:data_ori,which inclue CD4, GROUP,time.
>>> 
>>> 
>>> f_GROUP<-factor(data_ori$GROUP)
>>> 
>>> f_GROUP is a factor of 3 levels(0,1,2,3)
>>> 
>>> 
>>> result<- lme(sqrt(CD4) ~ f_GROUP*time ,random = ~time|ID,data=data_ori)
>>> 
>>> glht(result, linfct = mcp(f_GROUP="Tukey") )
>>> Error in `[.data.frame`(mf, nhypo[checknm]) : undefined columns selected
>>> 
>>> 
>>> I can't find out the reason for Error.
>> 
>> We can't either based on the information you provide. If I set up random
>> data based on your specification, then everything works fine for me in R
>> 2.13.1 with nlme 3.1-102 and multcomp 1.2-7.
>> 
>> ## random data
>> set.seed(0)
>> data_ori<- data.frame(
>>     CD4 = rnorm(100)^2,
>>     time = runif(100),
>>     f_GROUP = factor(sample(0:3, 100, replace = TRUE)),
>>     ID = factor(rep(1:10, each = 10))
>> )
>> 
>> ## model fitting
>> library("nlme")
>> result<- lme(sqrt(CD4) ~ f_GROUP * time, random = ~ time | ID,
>>     data = data_ori)
>> 
>> ## multiple comparisons
>> library("multcomp")
>> glht(result, linfct = mcp(f_GROUP = "Tukey"))
>> 
>> Everything works ok. glht() just issues an appropriate warning that Tukey
>> contrasts of f_GROUP might not be what you want in the presence of
>> interactions.
>
> Yes, but try this:
>
> set.seed(0)
> data_ori <- data.frame(
>    CD4 = rnorm(100)^2,
>    time = runif(100),
>    GROUP = sample(0:3, 100, replace = TRUE),
>    ID = factor(rep(1:10, each = 10))
> )
>
> f_GROUP<-factor(data_ori$GROUP)
> result <- lme(sqrt(CD4) ~ f_GROUP * time, random = ~ time | ID,
>    data = data_ori)
> glht(result, linfct = mcp(f_GROUP = "Tukey"))
> #Error in `[.data.frame`(mf, nhypo[checknm]) : undefined columns
> # selected
>
> (Note that I changed data_ori$f_GROUP to data_ori$GROUP.)

Thanks for pointing this out. So the workaround seems to be to simply put 
f_GROUP into data_ori.

The source of the problem is that because "nlme" does not provide a 
model.frame() method for "lme" objects, "multcomp" brings its own 
unexported multcomp:::model.frame.lme() which simply extracts $data. 
However, this is not always the model.frame - as in this example. It would 
probably be safer to modify the $call and reevaluate the model frame as 
lme() itself does.

> Peter Ehlers
>
>> 
>>> 
>>> Any help from you are welcome.
>>> 
>>> 
>>> Many thanks!
>>>
>>> 	[[alternative HTML version deleted]]
>>> 
>>> ______________________________________________
>>> 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.
>>> 
>> 
>> ______________________________________________
>> 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.
>
>