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

[Daniel Fulop]

Hi Russ,

Okay, I'll stick to your code.

As to the alternative, yes that would be better, but I'm afraid I have 
to sit down with pen and paper to figure out how to calculate the slope 
from each sample.  That's why I thought of resorting to using 
predict.MCMCglmm instead.  I guess I can use predict() and lstrends() to 
check that my slope samples are well calculated.

Thanks,
Dan.


Lenth, Russell V wrote:
> Dan,
>
> If you're talking about using the posterior mode of the variances and covariances, I think a whole big matrix of those would likely be very flaky. I think you are better off with my original suggestion.
>
> Alternatively, the really right thing to do is compute the slopes you need for *each* MCMC sample. Then you'll have a sample of each slope, from which you can obtain modes, HDRs, etc.
>
> Russ
>
> Sent from my iPhone
>
>> On Apr 10, 2015, at 6:51 PM, Daniel Fulop<dfulop.ucd at gmail.com>  wrote:
>>
>> 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:
>>> 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
>> -- 
>> 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
>>

-- 
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