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

[Jose Valdebenito Chavez]

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