Message-ID: <4254f651-7026-d1c0-6f8c-5b56768652e6@ed.ac.uk>
Date: 2021-04-20T07:39:02Z
From: Jarrod Hadfield
Subject: Bivariate MCMCglmm - single value threshold response variable
In-Reply-To: <AM6PR0302MB3335CF9C3845E846DB063C8BF7499@AM6PR0302MB3335.eurprd03.prod.outlook.com>
Hi,
It's hard to diagnose without the full code and summary of the output.
Could you also give some indication of the sample sizes and data structure?
Cheers,
Jarrod
On 19/04/2021 22:51, 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.
> 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):
> <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
> <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.
[[alternative HTML version deleted]]