[R-sig-ME] Syntax for a multivariate & multilevel MCMCglmm model

Fri Mar 4 20:54:06 CET 2016

Hi Allan,

Maybe take a look at the residuals, if they show independence and a lovely  
distribution for each trait (or even across traits after standardisation)  
when plotted the "us(trait):units" works for you. You assume each trait  
has a different residual variance and residuals are correlated between  
trait pairs within individuals. You can, of course, relax this or make  
this more complicated - only model diagnostics and comparisons can tell  
you.

I don't dare giving advice on priors (being completely clueless myself  
once correlations are involved), but the "course notes" by Jarrod Hadfield  
are extremely informative on many modeling matters and also include  
information about covariance structures and priors. You can find them here:

https://cran.r-project.org/web/packages/MCMCglmm/vignettes/CourseNotes.pdf

Best,
Paul

On Fri, 04 Mar 2016 19:00:21 +0200, Allan Debelle <a.debelle at sussex.ac.uk>  
wrote:

> Hi Paul,
>
> Thank you very much for taking the time to have a look at this.
> No need to modify the residual variance structure then? Just "rcov = ~  
> us(trait):units" ? And I don't need to change the syntax of the prior?
>
> Thank you.
>
> Al
>>> 	>>Paul Debes	>>4 March 2016 14:59

>> Hi Allan,
>> This model (explicitly nesting random replicates within treatments):
>> fixed = var ~ 1 + treatment + cov1 + cov2,random = ~  
>> treatment:replicate,rcov = ~ units
>> expanded to a multivariate model with unstructured covariances in G and  
>> R may be (you missed the "trait" interaction in the fixed >>statement  
>> of your first model):
>> fixed = cbind(varA, varB, varC, varD) ~ trait + trait:treatment +  
>> trait:cov1 + trait:cov2,random = ~ us(trait):treatment:replicate,rcov =  
>> ~ us(trait):units
>> Adding the trait interactions probably gets you back the treatment  
>> effects you did see in the univariate models.
>> How you add an additional random term, nested within replicates, will  
>> depend on how your random treatment:replicate effects >>are correlated  
>> between traits (e.g., "us" above may be reduced to "diag" if not  
>> correlated) and how the sub-replicates are >>correlated with the  
>> replicates and the traits. Plotting random effects and nested model  
>> testing may be required to find a useful >>solution.
>> Paul
>>
>> On Fri, 04 Mar 2016 15:30:19 +0200, Allan Debelle  
>> <a.debelle at sussex.ac.uk> wrote:
>>
>>>> 	>>Allan Debelle	>>4 March 2016 13:30

>> Hello everyone,
>>
>> I am trying to analyse a dataset using a multivariate MCMCglmm model, 
>> but I cannot solve a syntax problem related to the structure of my  
>> dataset.
>>
>> Basically, 50 individuals were exposed to either an A or a B treatment, 
>> and then 4 traits were measured on each individual, as well as 2 
>> covariates. This experiment was replicated 4 times, meaning that I have 
>> 4 replicates (rep1, rep2, rep3 and rep4) for each of the two treatments 
>> (A and B).
>>
>> The dataset looks something like this, with 'treat_rep' just being a 
>> unique identifier for each combination of treatment x replicate (I 
>> simplified it and generated random numbers, just FYI):
>>
>> varA varB varC varD cov1 cov2 treatment replicatetreat_rep
>>
>> 0.109488619 0.675591081 0.940580782 0.3669807360.911079042 0.468285642  
>> A rep1 A1
>> 0.400887809 0.43162503 0.823351715 0.76152362 0.0032215530.622398329 A  
>> rep2 A2
>> 0.176948608 0.074676865 0.91156597 0.7476630210.028740661 0.856395358 A  
>> rep3 A3
>> 0.16468968 0.776755434 0.367321135 0.8863521010.083174005 0.246884218 A  
>> rep4 A4
>> 0.090596435 0.785823554 0.061103075 0.8944361440.989249957 0.50533554 B  
>> rep1 B1
>> 0.839349822 0.639784216 0.293986243 0.558128180.307394708 0.036590982 B  
>> rep2 B2
>> 0.711702334 0.174145468 0.634296876 0.4421127730.480754889 0.863199351  
>> B rep3 B3
>> 0.052352675 0.492622311 0.668647151 0.6832434550.958397833 0.857768988  
>> B rep4 B4
>> ...
>>
>> ##############################
>>
>> What I want to do is to test for the effect of the treatment on the 
>> traits I measured, while taking into account the effects of my 
>> covariates, so something like this:
>>
>> varA + varB +varC +varD ~ treatment + cov1 + cov2
>>
>> The problem is that I have to deal with a nested data structure. Each 
>> individual belongs to a replicate of a treatment (e.g. rep1 of treatment 
>> A, or rep3 of treatment B, etc.), and I am not sure how to specify this 
>> correctly in MCMCglmm. I tried a few things, but something is wrong 
>> either in the prior specification or/and in my specification of the 
>> variance structure of the random effects (and potentially of the 
>> residuals too).
>>
>> ##############################
>>
>> Among other things, I tried to use the treat_rep variable as a grouping 
>> factor:
>>
>> prior<- list(R=list(V=diag(4),nu=4),G=list(G1=list(V=diag(4),nu=4)))
>>
>> model <- MCMCglmm(
>>
>> fixed = cbind(varA, varB, varC, varD) ~ -1 + treatment + cov1 + cov2,
>>
>> random = ~ us(trait):treat_rep,
>>
>> rcov = ~ us(trait):units,
>>
>> prior = prior,
>>
>> family = c("gaussian", "gaussian", "gaussian", "gaussian"), nitt = 
>> 100000, burnin = 5000,
>>
>> thin = 25, data = dataset)
>>
>> However, and unless I am wrong, it is not nesting 'replicate' within 
>> 'treatment'. Which means I end up with no effect of 'treatment', whereas 
>> I do have this effect when I run independent models using lme4 for each 
>> of the 4 traits (and in that case the nesting syntax is more 
>> straightforward).
>>
>> ##############################
>>
>> I also tried this:
>>
>> prior<-list(R=list(V=diag(16),nu=16),G=list(G1=list(V=diag(8),nu=8))) #  
>> 4traits x 2 treatments x 4 replicates // 4 traits x 2 treatments
>>
>> model <- MCMCglmm(
>>
>> fixed = cbind(varA, varB, varC, varD) ~ -1 + treatment + cov1 + cov2,
>>
>> random = ~ us(trait:treatment):replicate,
>>
>> rcov = ~ us(trait:treatment:replicate):units,
>>
>> prior = prior,
>> family = c("gaussian", "gaussian", "gaussian", "gaussian"), nitt = 
>> 100000, burnin = 5000,
>> thin = 25, data = dataset)
>>
>> but I get the following error message:
>>
>> Error in `[.data.frame`(data, , components[[1]]) :
>> undefined columns selected
>>
>> ##############################
>>
>> So, my first question is: What would be the syntax for this kind of 
>> nested structure (replicate nested within treatment) in MCMCglmm?
>>
>> My second question is: What would be the syntax if there was another 
>> level of variance related to this nested structure? I'm thinking of the 
>> situation where there would be a sort of treatment nesting within each 
>> replicate (e.g. if each of the four pairs of replicates was kept in a 
>> different incubator, or if each pair had been done at a different moment 
>> in time, etc.). Would you just add a 'replicate' random effect in the 
>> model? How?
>>
>> Any help with this would be fantastic.
>>
>> Thanks a lot in advance.
>>
>> Al
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
> --Allan Debelle - PhDResearch Fellow at University of Sussex
> Phone: +44 (0)748 095 3874
> Web: sites.google.com/site/allandebelle
> LinkedIn: fr.linkedin.com/in/allandebelle
>
>

-- 
Paul V. Debes
DFG Research Fellow
University of Turku
Department of Biology
Finland

Email: paul.debes at utu.fi