[R-sig-ME] poor convergence due to problematic data? multinomial MCMCglmm model

[Walid Mawass]

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