Skip to content

Problems with glht function for lme object

5 messages · anord, Ben Bolker, David Winsemius +1 more

anord
# 
Dear all, 

I'm trying to make multiple comparisons for an lme-object. The data is for
an experiment on parental work load in birds, in which adults at different
sites were induced to work at one of three levels ('treat'; H, M, L). The
response is 'feedings', which is a quantitative measure of nest provisioning
per parent per chick per hour. Site is included as a random effect (one
male/female pair per site). 

My final model takes the following form:
feedings ~ treat + year + data^2, random = ~1|site,data=feed.df

For this model, I would like to do multiple comparisons on 'treat', using
the multcomp package:
summary(glht(m4.feed,linfct=mcp(treat="Tukey")))

However, this does not work, and I get the below error message.

Error in if (is.null(pkg) | pkg == "nlme") terms(formula(x)) else slot(x,  : 
  argument is of length zero
Error in factor_contrasts(model) : 
  no ?model.matrix? method for ?model? found!
 
I suspect this might have quite a straightforward solution, but I'm stuck at
this point.
Any help would be most appreciated. Sample data below.

Kind regards, 
Andreas Nord
Sweden

==============
feedings     sex site  treat year  date^2
1.8877888   M  838     H 2009      81
1.9102787   M  247     H 2009      81
1.4647229   M  674     H 2010      121
1.4160590   M 7009     M 2009      144
1.3106749   M  863     M 2010      196
1.2718121   M   61     M 2009      225
1.2799263   M  729     L 2009      256
1.5829564   M  629     L 2009      256
1.4847251   M  299     L 2010      324
1.2463151   M  569     L 2010      324
2.1694169   F  838     H 2009      81
1.5966899   F  247     H 2009      81
2.4136983   F  674     H 2010      121
1.7784873   F 7009     M 2009      144
1.6681317   F  863     M 2010      196
2.3691275   F   61     M 2009      225
2.0672192   F  729     L 2009      256
1.6389902   F  629     L 2009      256
0.9307536   F  299     L 2010      324
1.6786767   F  569     L 2010      324
==============
Ben Bolker
# 
anord <andreas.nord <at> zooekol.lu.se> writes:
(For more complicated mixed model questions you might want
to try the r-sig-mixed-models list instead.)

It works for me:

feed.df <- read.table(textConnection("
feedings     sex site  treat year  date
1.8877888   M  838     H 2009      81
1.9102787   M  247     H 2009      81
1.4647229   M  674     H 2010      121
1.4160590   M 7009     M 2009      144
1.3106749   M  863     M 2010      196
1.2718121   M   61     M 2009      225
1.2799263   M  729     L 2009      256
1.5829564   M  629     L 2009      256
1.4847251   M  299     L 2010      324
1.2463151   M  569     L 2010      324
2.1694169   F  838     H 2009      81
1.5966899   F  247     H 2009      81
2.4136983   F  674     H 2010      121
1.7784873   F 7009     M 2009      144
1.6681317   F  863     M 2010      196
2.3691275   F   61     M 2009      225
2.0672192   F  729     L 2009      256
1.6389902   F  629     L 2009      256
0.9307536   F  299     L 2010      324
1.6786767   F  569     L 2010      324"),
           header=TRUE)

library(nlme)
m4.feed <- lme(feedings ~ treat + year + date, random = ~1|site,data=feed.df)

library(multcomp)
gg <- glht(m4.feed,linfct=mcp(treat="Tukey"))
plot(gg)
summary(gg)

It works for me (although it doesn't look like there's anything
going on there ...)

=========

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lme.formula(fixed = feedings ~ treat + year + date, data = feed.df, 
    random = ~1 | site)

Linear Hypotheses:
           Estimate Std. Error z value Pr(>|z|)
L - H == 0  -0.6452     0.7667  -0.842    0.592
M - H == 0  -0.3996     0.4276  -0.934    0.529
M - L == 0   0.2456     0.4229   0.581    0.771
(Adjusted p values reported -- single-step method)
sessionInfo()
R version 2.12.1 (2010-12-16)
Platform: i486-pc-linux-gnu (32-bit)

locale:
 [1] LC_CTYPE=en_CA.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_CA.UTF-8        LC_COLLATE=en_CA.UTF-8    
 [5] LC_MONETARY=C              LC_MESSAGES=en_CA.UTF-8   
 [7] LC_PAPER=en_CA.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] splines   stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
[1] multcomp_1.2-4  survival_2.36-3 mvtnorm_0.9-95  nlme_3.1-97    

loaded via a namespace (and not attached):
[1] grid_2.12.1     lattice_0.19-18
David Winsemius
#
On Jan 7, 2011, at 8:28 AM, anord wrote:

            

      
      
      
    

        
I am guessing that problems could easily arise as a result of your  
variable name containing "^". That is an invalid name in an un-back- 
quoted state. You have clearly not given us a copy of a console  
session since your variable is date^2  below and data^2 above.

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Dieter Menne
    

    

      
      #
    

  


  

      
anord wrote:

            

      
      
      
    

        
1) Try I(date^2)
2) Check the version you are using. For some older versions, multcomp and
lme were not close friends (even if that is some time ago now)

Dieter

  
    

      
      
    

  


  


      

    

    
1 day later

    

      
      

        


  

        

      anord
    

    

      
      #
    

  


  

      
Dear all, 

Thanks for your input. It works fine now. All it took was for me to tidy up
the workspace. Simple solution, that I should have considered earlier.

Regards, 
Andreas