Bivariate MCMCglmm - single value threshold response variable

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

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