R help - Apr 2008 - Post hoc tests with lme

Using the "ergoStool" data cited in Mixed-Effects Models in S and
S-PLUS by Pinheiro and Bates as an example, we have

=======> library(nlme)
> fm <- lme(effort~Type-1, data=ergoStool, random=~1|Subject)
> summary(fm)
Linear mixed-effects model fit by REML
  Data: ergoStool
       AIC      BIC   logLik
  133.1308 141.9252 -60.5654

Random effects:
 Formula: ~1 | Subject
        (Intercept) Residual
StdDev:    1.332465 1.100295

Fixed effects: effort ~ Type - 1
           Value Std.Error DF  t-value p-value
TypeT1  8.555556 0.5760123 24 14.85308       0
TypeT2 12.444444 0.5760123 24 21.60448       0
TypeT3 10.777778 0.5760123 24 18.71102       0
TypeT4  9.222222 0.5760123 24 16.01046       0
 Correlation:
       TypeT1 TypeT2 TypeT3
TypeT2 0.595
TypeT3 0.595  0.595
TypeT4 0.595  0.595  0.595

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max
-1.80200345 -0.64316591  0.05783115  0.70099706  1.63142054

Number of Observations: 36
Number of Groups: 9

=======
Now suppose I want to test the following hypothesis

H0: TypeT1 =0 and TypeT2 = 0

I've tried estimable() and glh.test() in package gmodels, esticon() in
package boBy, and linear.hypothesis() in package car, but it seems
none of them would work with objects from lme:

=======> library(gmodels)
> estimable(fm, rbind(c("TypeT1"=1), c("TypeT2"=1)))
Error in FUN(newX[, i], ...) :
  `param' has no names and does not match number of coefficients of
model. Unable to construct coefficient vector
> glh.test(fm, rbind(c("TypeT1"=1), c("TypeT2"=1)))
Error in glh.test(fm, rbind(c(TypeT1 = 1), c(TypeT2 = 1))) :
  Only defined for lm,glm objects
> library(doBy)
> esticon(fm, rbind(c("TypeT1"=1), c("TypeT2"=1)))
Error in t(abs(t(tmp) * obj$fixDF$X)) :
  dims [product 2] do not match the length of object [4]
In addition: Warning message:
In esticon.lme(fm, rbind(c(TypeT1 = 1), c(TypeT2 = 1))) :
  The esticon function has not been thoroughly teste on 'lme' objects
> library(car)
> linear.hypothesis(fm, rbind(c("TypeT1"=1),
c("TypeT2"=1)))
Error in L %*% b : requires numeric matrix/vector arguments
=======
So is there any other package with which I can run this kind of tests?

Thanks,
Gang

At the risk of annoying Doug Bates, I'll point out that "glht" in
the multcomp package works with lmer objects.  In fact, you can
supply your own degrees of freedom value via the "df" argument
(a fact which is not immediately obvious in the glht help page).
If you don't supply a df value, it will use "Inf" (i.e., normal
instead of t).

Regards,   Rob Kushler


Gang Chen wrote:
> Using the "ergoStool" data cited in Mixed-Effects Models in S and
> S-PLUS by Pinheiro and Bates as an example, we have
> 
> =======>> library(nlme)
>> fm <- lme(effort~Type-1, data=ergoStool, random=~1|Subject)
>> summary(fm)
> 
> Linear mixed-effects model fit by REML
>   Data: ergoStool
>        AIC      BIC   logLik
>   133.1308 141.9252 -60.5654
> 
> Random effects:
>  Formula: ~1 | Subject
>         (Intercept) Residual
> StdDev:    1.332465 1.100295
> 
> Fixed effects: effort ~ Type - 1
>            Value Std.Error DF  t-value p-value
> TypeT1  8.555556 0.5760123 24 14.85308       0
> TypeT2 12.444444 0.5760123 24 21.60448       0
> TypeT3 10.777778 0.5760123 24 18.71102       0
> TypeT4  9.222222 0.5760123 24 16.01046       0
>  Correlation:
>        TypeT1 TypeT2 TypeT3
> TypeT2 0.595
> TypeT3 0.595  0.595
> TypeT4 0.595  0.595  0.595
> 
> Standardized Within-Group Residuals:
>         Min          Q1         Med          Q3         Max
> -1.80200345 -0.64316591  0.05783115  0.70099706  1.63142054
> 
> Number of Observations: 36
> Number of Groups: 9
> 
> =======> 
> Now suppose I want to test the following hypothesis
> 
> H0: TypeT1 =0 and TypeT2 = 0
> 
> I've tried estimable() and glh.test() in package gmodels, esticon() in
> package boBy, and linear.hypothesis() in package car, but it seems
> none of them would work with objects from lme:
> 
> =======>> library(gmodels)
>> estimable(fm, rbind(c("TypeT1"=1), c("TypeT2"=1)))
> Error in FUN(newX[, i], ...) :
>   `param' has no names and does not match number of coefficients of
> model. Unable to construct coefficient vector
>> glh.test(fm, rbind(c("TypeT1"=1), c("TypeT2"=1)))
> Error in glh.test(fm, rbind(c(TypeT1 = 1), c(TypeT2 = 1))) :
>   Only defined for lm,glm objects
> 
>> library(doBy)
>> esticon(fm, rbind(c("TypeT1"=1), c("TypeT2"=1)))
> Error in t(abs(t(tmp) * obj$fixDF$X)) :
>   dims [product 2] do not match the length of object [4]
> In addition: Warning message:
> In esticon.lme(fm, rbind(c(TypeT1 = 1), c(TypeT2 = 1))) :
>   The esticon function has not been thoroughly teste on 'lme'
objects
> 
>> library(car)
>> linear.hypothesis(fm, rbind(c("TypeT1"=1),
c("TypeT2"=1)))
> Error in L %*% b : requires numeric matrix/vector arguments
> =======> 
> So is there any other package with which I can run this kind of tests?
> 
> Thanks,
> Gang
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
>

Try glht in package multcomp.

Simon.

On Wed, 2008-04-16 at 12:00 -0400, Gang Chen wrote:
> = 1), c(TypeT2 = 1))) :
>   Only defined for lm,glm objec
-- 
Simon Blomberg, BSc (Hons), PhD, MAppStat. 
Lecturer and Consultant Statistician 
Faculty of Biological and Chemical Sciences 
The University of Queensland 
St. Lucia Queensland 4072 
Australia
Room 320 Goddard Building (8)
T: +61 7 3365 2506
http://www.uq.edu.au/~uqsblomb
email: S.Blomberg1_at_uq.edu.au

Policies:
1.  I will NOT analyse your data for you.
2.  Your deadline is your problem.

The combination of some data and an aching desire for 
an answer does not ensure that a reasonable answer can 
be extracted from a given body of data. - John Tukey.