[R-sig-ME] Bivariate MCMCglmm - single value threshold response variable

Jarrod Hadfield j@h@d||e|d @end|ng |rom ed@@c@uk
Tue Apr 20 09:39:02 CEST 2021


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 using ed.ac.uk>
> *Sent:* 08 April 2021 14:28
> *To:* Tara Cox [RPG] <bstlc using leeds.ac.uk>; 
> r-sig-mixed-models using r-project.org <r-sig-mixed-models using 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 using 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]]



More information about the R-sig-mixed-models mailing list