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

[Jarrod Hadfield]

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