[R-sig-ME] Bivariate MCMCglmm - single value threshold response variable
Tara Cox [RPG]
b@t|c @end|ng |rom |eed@@@c@uk
Wed Apr 28 00:28:20 CEST 2021
Hi Jarrod,
Below is a table detailing the number of individuals in my dataset and how many assays exists per individual for trait X. For the model I'd mentioned failing in my previous email, trait Y is decision to disperse Y/N (rather than survival):
Table showing the number of individuals that possess 1, 2, 3, 4, 5 or 6 assays for trait X.
Assay number
Total individuals
Total assays
1
2
3
4
5
6
Philopatric
15
2
1
0
0
0
18
22
Disperse
154
85
21
8
4
3
275
457
ALL
169
87
22
8
4
3
293
479
I assume the poor convergence is due to the small sample sizes for philopatric birds, but I was hoping there may be something that can be done. The model summary:
## Iterations = 60001:2199501
## Thinning interval = 500
## Sample size = 4280
##
## DIC: 2962.032
##
## G-structure: ~us(at.level(variable, "X")):ObserverID
##
## post.mean l-95% CI
## at.level(variable, "X"):at.level(variable, "X").ObserverID 0.2269 0.06213
## u-95% CI eff.samp
## at.level(variable, "X"):at.level(variable, "X").ObserverID 0.4546 4280
##
## ~us(at.level(variable, "X")):BirdID
##
## G-R structure below
##
## R-structure: ~us(at.level(variable, "DISP")):BirdID
##
## post.mean
## at.level(variable, "X").BirdID:at.level(variable, "X").BirdID 0.2458
## at.level(variable, "DISP").BirdID:at.level(variable, "X").BirdID -0.0168
## at.level(variable, "X").BirdID:at.level(variable, "DISP").BirdID -0.0168
## at.level(variable, "DISP").BirdID:at.level(variable, "DISP").BirdID 1.0000
## l-95% CI
## at.level(variable, "X").BirdID:at.level(variable, "X").BirdID 0.0005161
## at.level(variable, "DISP").BirdID:at.level(variable, "X").BirdID -0.3445199
## at.level(variable, "X").BirdID:at.level(variable, "DISP").BirdID -0.3445199
## at.level(variable, "DISP").BirdID:at.level(variable, "DISP").BirdID 1.0000000
## u-95% CI
## at.level(variable, "X").BirdID:at.level(variable, "X").BirdID 0.5498
## at.level(variable, "DISP").BirdID:at.level(variable, "X").BirdID 0.3423
## at.level(variable, "X").BirdID:at.level(variable, "DISP").BirdID 0.3423
## at.level(variable, "DISP").BirdID:at.level(variable, "DISP").BirdID 1.0000
## eff.samp
## at.level(variable, "X").BirdID:at.level(variable, "X").BirdID 1069
## at.level(variable, "DISP").BirdID:at.level(variable, "X").BirdID 1338
## at.level(variable, "X").BirdID:at.level(variable, "DISP").BirdID 1338
## at.level(variable, "DISP").BirdID:at.level(variable, "DISP").BirdID 0
##
## ~idh(at.level(variable, "X")):units
##
## post.mean l-95% CI u-95% CI eff.samp
## at.level(variable, "X").units 1.164 0.8185 1.505 1244
##
## Location effects: X.DISP.stack ~ variable - 1 + variable:Sex + at.level(variable, "X"):Age + at.level(variable, "X"):Age2 + at.level(variable, "X"):TentColour + at.level(variable, "X"):BranchOrientation + at.level(variable, "X"):BranchHeight + at.level(variable, "X"):TentPoles + at.level(variable, "X"):Assay.number + at.level(variable, "DISP"):PopulationDensity + at.level(variable, "DISP"):LocalPopDensity + at.level(variable, "DISP"):FoodAbundance
##
## post.mean l-95% CI u-95% CI
## variableX 1.551577 1.156333 1.920005
## variableDISP 2.012833 1.523589 2.491882
## variableX:Sex1 0.246286 0.002812 0.499340
## variableDISP:Sex1 0.438481 -0.138548 0.964070
## at.level(variable, "X"):Age 1.520044 0.795296 2.274839
## at.level(variable, "X"):Age2 -0.940711 -1.640312 -0.212713
## at.level(variable, "X"):TentColourG -0.522777 -0.880206 -0.168102
## at.level(variable, "X"):BranchOrientationP 0.266240 -0.045538 0.568197
## at.level(variable, "X"):BranchHeight1 -0.172651 -1.449845 1.128834
## at.level(variable, "X"):TentPoles1 -0.734384 -2.094550 0.526492
## at.level(variable, "X"):Assay.number 0.419480 0.297508 0.561766
## at.level(variable, "DISP"):PopulationDensity 0.508142 0.195178 0.799767
## at.level(variable, "DISP"):LocalPopDensity -0.350103 -0.582710 -0.092236
## at.level(variable, "DISP"):FoodAbundance 1.224312 0.412853 1.949664
## eff.samp pMCMC
## variableX 4280 < 2e-04 ***
## variableDISP 3946 < 2e-04 ***
## variableX:Sex1 4059 0.052336 .
## variableDISP:Sex1 4280 0.112150
## at.level(variable, "X"):Age 3982 < 2e-04 ***
## at.level(variable, "X"):Age2 4039 0.013084 *
## at.level(variable, "X"):TentColourG 4280 0.004206 **
## at.level(variable, "X"):BranchOrientationP 4280 0.096729 .
## at.level(variable, "X"):BranchHeight1 4051 0.772430
## at.level(variable, "X"):TentPoles1 3517 0.273832
## at.level(variable, "X"):Assay.number 4280 < 2e-04 ***
## at.level(variable, "DISP"):PopulationDensity 4070 < 2e-04 ***
## at.level(variable, "DISP"):LocalPopDensity 4280 0.005140 **
## at.level(variable, "DISP"):FoodAbundance 3493 0.000467 ***
Best wishes,
Tara
Tara Cox
Pronouns: she, her, hers
PhD researcher
Dugdale group, School of Biology
Faculty of Biological Sciences
University of Leeds
________________________________
From: Jarrod Hadfield <j.hadfield using ed.ac.uk>
Sent: 20 April 2021 08:39
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,
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><mailto:j.hadfield using ed.ac.uk>
Sent: 08 April 2021 14:28
To: Tara Cox [RPG] <bstlc using leeds.ac.uk><mailto:bstlc using leeds.ac.uk>; r-sig-mixed-models using r-project.org<mailto:r-sig-mixed-models using r-project.org> <r-sig-mixed-models using r-project.org><mailto: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):
>
> 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<mailto: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.
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list