Skip to content

contrast of predicted slopes from MCMCglmm model

4 messages · Daniel Fulop, Lenth, Russell V

Original
Lenth, Russell V
# 
Daniel,

As long as you can get the regression coefficients, covariance of the coefficients, and a suitable grid of linear functions, you can obtain an lsmobj object that can be summarized and further analyzed in the 'lsmeans' package. For your question, I think something like this will work (I'm assuming your fitted mcmcGLMM object is named 'mod'):

# Get coefs and covariance matrix
bhat <- apply(mod$Sol, 2, mean)
V <- cov(mod$Sol)

# Set up desired factor levels (you likely want to change this)
levels <- list(day = 1:4, treatment = c("control", "cold"), species = factor(1:10))

# Obtain a differece quotient of model matrices 
lv1 <- lv2 <- levels
lv1$day = lv1$day - .0001
lv2$day = lv2$day + .0001
grid1 <- do.call(expand.grid, lv1)
grid2 <- do.call(expand.grid, lv2)
lf1 <- model.matrix(fixedForm, data = grid1)
lf2 <- model.matrix(fixedForm, data = grid2)
linfct <- (lf2 - lf1) / .0002

# Create the needed lsmobj
library(lsmeans)
mod.lsm <- lsmobj(bhat = bhat, V = V, levels = levels, 
    linfct = linfct, by = "species")

# Can now use summary, contrast, pairs, etc. on this result


I hope this helps

Russ

Russell V. Lenth  -  Professor Emeritus
Department of Statistics and Actuarial Science   
The University of Iowa  -  Iowa City, IA 52242  USA   
Voice (319)335-0712 (Dept. office)  -  FAX (319)335-3017


-----Original Message-----
Date: Thu, 09 Apr 2015 18:32:11 -0700
From: Daniel Fulop <dfulop.ucd at gmail.com>
To: R_MixedModels <r-sig-mixed-models at r-project.org>
Subject: [R-sig-ME] contrast of predicted slopes from MCMCglmm model
Message-ID: <5527281B.3080607 at gmail.com>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed

Hi mixed modelers,

I would like to contrast average slopes using predictions from an LMM model fit with MCMCglmm.  In essence I want to do something analogous to lsmeans' lstrends() function.  The predicted slopes are growth rates ...more below.

I have 12 days of plant growth data for 10 closely related species at 2 temperatures.  There are several individuals per species at each temperature, and they're fully randomized.  I want to asses whether each species' growth rate differs between temperatures.  What I mean by that is that I would like to contrast the average slopes for each species in control vs. cold temperature; these slopes are the growth rates.

I am modeling the data in MCMCglmm so I can account for phylogeny and also to be able to try multi-response models (to jointly model the stem, plastochrons, root, etc).  I am starting out with the stem length data and after trying different time dependencies settled on this 3rd degree polynomial model:

stemLen ~ poly(day, 3, raw=TRUE) * treatment * species + (1 + day | indiv)

where day is a numeric time variable, treatment a 2 level factor (control and cold), and species a 10 level factor.

I have a slight idea of how I would go about contrasting the average slopes for each species and that I would use the predict.MCMCglmm() function, but I would really appreciate some guidance before I go down the wrong the wrong path.

Thanks for your help!
Dan.

--
Daniel Fulop, Ph.D.
Postdoctoral Scholar
Dept. Plant Biology, UC Davis
Maloof Lab, Rm. 2220
Life Sciences Addition, One Shields Ave.
Davis, CA 95616
Daniel Fulop
# 
Hi Russ,

Thanks a ton!  I did notice that lsmeans could be extended, but wasn't 
sure of how to get all the pieces from my MCMCglmm model object.  I was 
going to start getting the slopes from predictions on a grid of 'day' 
values to then average and finally subtract averages to contrast, but 
you've spared me that trouble.  Plus, I already have code to summarize 
and plot growth rate results from lsmeans, which I can now reuse.

I have one question about your example code, though.  Since I have 
posterior estimates, wouldn't it be better to get the bhat vector by 
using posterior.mode() instead of mean()?

Thanks,
Dan.
Lenth, Russell V wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Lenth, Russell V
    

    

      
      #
    

  


  

      
On your question, maybe, but I'm not so sure the covariance matrix is right for that. Perhaps others can comment. It'd be good also to do some plots and perhaps other diagnostics to see if the sampler results are reasonably modeled as multivariate normal, as that is an underlying assumption of all the frequentist-style tests and CIs that you obtain.

Russ

Sent from my iPad

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Daniel Fulop
    

    

      
      #
    

  


  

      
Hi Russ,

Thanks again for the feedback and quick response!  My primary goal is 
simply to compare slopes within species and slope differences among 
species, the tests and CIs are a secondary concern.  lstrends() provides 
an easy means of estimating average slopes andslope differences .  I 
take your point about diagnostics and plotting, though, and will do some 
checking such as comparing the lstrends estimates to analogous estimates 
calculated using predict.MCMCglmm output.

On the point about the covariance matrix, the MCMCglmm output includes 
VCV samples.  So, I suppose I could calculate posterior modes of each 
cell in the vcv and use that instead?  Otherwise using the means should 
be fine.  From visual inspection of the posterior samples they're quite 
normal in appearance, so the means and modes should be quite close to 
each other.

Cheers,
Dan.

        
Lenth, Russell V wrote: