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

[Tara Cox [RPG]]

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