Bivariate MCMCglmm - single value threshold response variable

Mon, Apr 19, 2021 2:51 PM

Hi Jarrod,

Thanks for getting back to me!

I've been attempting to run my models with your advised priors and it has improved convergence. However, for some of my models, the trace and density of plots of variableY (i.e. survival) still show poor mixing (see attached). The heidel.diag for at.level(variable, "X").ID:at.level(variable, "Y").ID also fails. I have tried playing around with the number of iterations, thin and burnin, but haven't had any luck. Is it possible to amend the prior to assist with convergence?

Further, as my hypotheses state that phenotypic trait X affects survival, I'd like to look at correlation between the two. I had considered using the below to extract these values, but am not sure what the best method would be. Would you be able to give any advice?

corr.calc <- model2$VCV[,"traitY:traitX.ID"]/(sqrt(model2$VCV[,"traitY:traitY.ID"])*sqrt(model2$VCV[,"traitX:traitX.ID"]))

Best wishes,
Tara

From: Jarrod Hadfield <j.hadfield at ed.ac.uk>
Sent: 08 April 2021 14:28
To: Tara Cox [RPG] <bstlc at leeds.ac.uk>; r-sig-mixed-models at r-project.org <r-sig-mixed-models at r-project.org>
Subject: Re: [R-sig-ME] Bivariate MCMCglmm - single value threshold response variable

Hi Tara,

Your second model is correct. You can ignore the warning (not an error)
that some fixed effects have been removed; often when you set up
equations involving at.level, base R's model.matrix will not
automatically delete extra parameters and so MCMCglmm checks up on this.

The model is not mixing because you are trying to estimate the residual
variance for the binary trait and this is not identifiable in the
likelihood. You should fix it at one which results in a probit model. To
do this have fix=2 when specifying the prior for R1.

Also, make sure you are using the latest version of MCMCglmm as there
was a bug with covu models thatmay be triggered under certain
circumstances and this has now been fixed.

Cheers,

Jarrod


On 08/04/2021 12:02, Tara Cox [RPG] wrote:
> This email was sent to you by someone outside the University.
> You should only click on links or attachments if you are certain that the email is genuine and the content is safe.
>
> Dear list,
>
> I am stuck attempting to run a bivariate MCMCglmm that uses response variables phenotypic trait 'x' and survival (yes/no) (hereafter referred to as 'x' and 'y', modelled Poisson and threshold, respectively). I have multiple repeats per individual for x, but a single value for y. From my research, there are two options for running this analysis:
>
> 1) As trait y possesses a single value per individual, there is no within-individual variance and so I fix the variance to 0.0001 in my prior, as per below:
>
> prior.chi = list(R = list(V = diag(c(1, 0.0001), 2, 2), nu = 2, fix=2),
>                           G = list(G1 = list(V = diag(2), nu = 1000, alpha.mu = c(0,0), alpha.V = diag(c(1, 1)),   #chi squared as this random effect (AnimalID) is specified to both response variables
>                                        G2 = list(V = diag(1), nu = 1,  alpha.mu = 0, alpha.V = diag(1000,1))))             #parameter expanded as this random effect (ObserverID) is specified only to response variable x, which is Poisson
>
> model1 <- MCMCglmm(cbind(x, y) ~ trait-1 +
>                            trait:Sex +
>                            at.level(trait,1):Age +
>                            at.level(trait,1):Age2 +
>                            at.level(trait,2):PopulationDensity +
>                            at.level(trait,2):LocalPopulationDensity +
>                            at.level(trait,2):FoodAbundance
>                          random =~ us(trait):AnimalID + idh(at.level(trait,1)):ObserverID,
>                          rcov =~ idh(trait):units,
>                          family = c("poisson","threshold"),
>                          prior = prior.chi,
>                          nitt=2250000,
>                          burnin=150000,
>                          thin=525,
>                          verbose = FALSE,
>                          pr=FALSE,
>                          data = data)
>
> 2) Use the stacked data/'covu' approach, where I get rid of the ID term for the threshold variable, and allow a covariance between the threshold residual and the Poisson ID term (as per supplementary material in Thomson et al. 2017: https://doi.org/10.1111%2Fevo.13169):
>
> prior.covu <- list(G = list(G1 = list(V = diag(1), nu = 1)),                                   # rand effect for ObserverID (fitted for x)
>                                R = list(R1 = list(V = diag(2), nu = 0.002, covu = TRUE),     # 2-way var-cov matrix of AnimalID for x, residual for y
>                                             R2 = list(V = diag(1), nu = 0.002)))                            # residual for x
>
> model2 <- MCMCglmm(x.y.data.stack ~ variable - 1 +
>                          trait:Sex +
>                            at.level(variable, "x"):Age
>                            at.level(variable, "x"):Age2 +
>                            at.level(variable, "y"):PopulationDensity +
>                            at.level(variable, "y"):LocalPopulationDensity +
>                            at.level(variable, "y"):FoodAbundance,
>                          random = ~ us(at.level(variable,"x")):ObserverID + us(at.level(variable, "x")):AnimalID,
>                          rcov = ~us(at.level(variable, "y")):AnimalID + idh(at.level(variable, "x")):units,
>                          family = NULL,            #specified already in the data
>                          prior = prior.covu,
>                          nitt=2250000,
>                          burnin=150000,
>                          thin=525,
>                          verbose = FALSE,
>                          pr=FALSE,
>                          data = data)
>
> My problem is that I cannot get my models for either method to run successfully.
>
> The convergence model 1 is poor, with trace+density plots showing poor mixing. I know that the method is viable, as I have successfully run multiple models without convergence issues using a full parameter expanded prior alongside Poisson and gaussian response variables. Therefore, I must be specifying my chi squared prior incorrectly. I've tried to read further into how to amend the prior, but I've reached my limit of understanding and keep getting stuck.
>
> When running model 2, I receive an error message stating some fixed effects are not estimable and have been removed, and that I should use an informative prior. I've tried specifying a parameter expanded prior by modifying the G structure to include 'alpha.mu = 0' and 'alpha.V = diag(1000,1)', but this doesn't help. When running the model with a single fixed effect (e.g. Age), I do not receive any error messages, but the model shows poor mixing.
>
> If anyone could provide any insight into which method might be better suited to my analysis, or how to improve the priors for either model, it would be much appreciated!
>
> Thanks a lot.
>
> Best wishes,
> Tara
>
>
>
>
>          [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336. Is e buidheann carthannais a th? ann an Oilthigh Dh?n ?ideann, cl?raichte an Alba, ?ireamh cl?raidh SC005336.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: Trace_density_variableSURVIVAL.png
Type: image/png
Size: 160002 bytes
Desc: Trace_density_variableSURVIVAL.png
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20210419/9349dee3/attachment-0001.png>

Bivariate MCMCglmm - single value threshold response variable

Thread (5 messages)