[R-sig-ME] contrast of predicted slopes from MCMCglmm model

[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

> On Apr 10, 2015, at 5:39 PM, Daniel Fulop <dfulop.ucd at gmail.com> wrote:
> 
> 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:
>> 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, Ph.D.
> Postdoctoral Scholar
> Dept. Plant Biology, UC Davis
> Maloof Lab, Rm. 2220
> Life Sciences Addition, One Shields Ave.
> Davis, CA 95616
> 
> 510-253-7462
> dfulop at ucdavis.edu
>