[R-sig-ME] poor convergence due to problematic data? multinomial MCMCglmm model
Walid Mawass
w@||dm@w@@@10 @end|ng |rom gm@||@com
Mon Jun 27 21:43:16 CEST 2022
One thing you can do is increase your burn-in time to at least 5*10^5.
Right now you have it at 2000 for a 5000000 iterations which is quite
short. And of course like you mentioned the issue with polygyny is the lack
of representation in the datatset
Walid Mawass
On Mon, Jun 27, 2022, 2:55 PM Jose Valdebenito Chavez via
R-sig-mixed-models <r-sig-mixed-models using r-project.org> wrote:
> Hi,
>
> I am having problems with MCMCglmm models trying to describe drug
> resistance in birds.
>
> Response variable: present/absence of drug resistance. That is
> cbind(all.present, all.absent) which will be e.g. Anas crecca = 4, 1.
> Explanatory variable: mating system, categorical with two levels
> (monogamy, polygyny)
> Number of species (and number of rows): 31
> Actual number of individuals birds examined: 711
>
>
> The Model:
>
> prior=list(R = list(V = diag(1), nu = 0.02),
> G = list(G1 = list(V = diag(1), nu = 2, alpha.mu = c(0),
> alpha.V = diag(1))))
>
> mp3 <- MCMCglmm(cbind(all.present, all.absent) ~ msys,
> random=~species,
> ginverse=list(species=inv.phylo$Ainv),
> prior=prior,
> family = "multinomial2",
> data = data,
> nitt=5002000,burnin=2000,thin=5000)
>
> The Outcome:
>
> Iterations = 2001:4997001
> Thinning interval = 5000
> Sample size = 1000
> DIC: 172.8442
>
> G-structure: ~species
> post.mean l-95% CI u-95% CI eff.samp
> species 2.065. 2.773e-06 8.973 1000
>
> R-structure: ~units
> post.mean. l-95% CI u-95% CI eff.samp
> units 4.189. 0.004191 12.83 1000
>
> Location effects: cbind(all.present, all.absent) ~ msys
> post.mean l-95% CI u-95% CI eff.samp pMCMC
> (Intercept) 4.166 2.174 6.947 939.42
> 0.002 **
> msyspolygyny 434.100 63.507 702.711 26.66 <0.001 ***
> ---
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>
> By the outcome and trace plots the model seems fine except for
> msyspolygyny, and I haven’t been able to correct the issue so far.
> Increasing nitt to 5002000 improved the autocorrelation a bit but still
> far from good.
>
> diag(autocorr(mp3$Sol)[2, , ])
> (Intercept) msyspolygyny
> 0.01511989 0.94800971
>
> Convergence for polygyny is, not surprisingly, bad.
>
> gelman.diag(mcmc.list(mp3$Sol, mp3b$Sol))
> Potential scale reduction factors:
>
> Point est. Upper C.I.
> (Intercept) 0.999 1.00
> msyspolygyny 1.357 2.11
>
> Multivariate psrf
> 1.23
>
> Then I checked and I realised that the problem is likely to originate from
> the data itself that has just 3 species with a polygynous mating system,
> out of 31.
>
> Is there a way to improve the autocorrelation/convergence? Run it for even
> longer? Maybe by specifying a prior able to handle it?
>
> Thanks,
> Jose
>
>
>
> -------------------------
>
> José O. Valdebenito
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list