persistant autocorrelation in binomial MCMCglmm
Hi Jarrod, I have an update on the model performance with a probit link - the autocorrelation is much better behaved with the "threshold" family. All ACF values are <0.1 for iteration and thinning #s that were resulting in autocorrelation using the logit link. Looking at the latent variables, though, a lot of the distributions included values below -7 (lowest was -10). All of the means where within the -7 to 7 range though, because only a few estimates per observation tended to reach very low negative values. Are *any* estimates outside the range considered problematic? I noticed on the forum a post where someone else had this issue ( https://stat.ethz.ch/pipermail/r-sig-mixed-models/2012q3/019067.html) and I also tried the chi square prior you suggested for that problem (V=1, nu=1000, alpha.mu=0, alpha.V=1) but the result was the same. In terms of the data and system, I would expect an extremely low, near zero probability of interaction for some of these dyads because they are not using similar areas and so are not physically able to interact. Is this signal perhaps too strong? If my goal is to weed out these improbable interactions, though, will the model not serve this purpose? Many thanks, Christina Christina M. Aiello Biologist- U.S. Geological Survey Las Vegas Field Station 160 N. Stephanie St. Henderson, NV 89074 (702) 481-3957 caiello at usgs.gov
On Thu, May 4, 2017 at 10:17 AM, Aiello, Christina <caiello at usgs.gov> wrote:
Hi Jarrod, Appreciate the quick response and thoughts 1) I thought I had checked the absolute value of the latent variables, but now that I look again, I must not have examined both ends of the distribution. There are 167 observations whose latent variable distributions dip below -20 (minimum was -32). In each of these cases, the left tail of the distribution includes very few estimates at such low values. Do you think this has to do with the rarity of a response of 1 in the dataset? Or might this be indicative of another problem? 2) I'll give the probit link a try today and see how the results compare 3) I can definitely let the chains run longer, and continue increasing thinning? I was hesitant to keep upping the values because I haven't seen many published analyses using iterations and intervals beyond what I've been trying. I was worried having to run the chain so long might be indicative of other underlying problems that I wasn't considering. Many thanks! Christina Christina M. Aiello Biologist- U.S. Geological Survey Las Vegas Field Station 160 N. Stephanie St. Henderson, NV 89074 (702) 481-3957 caiello at usgs.gov On Wed, May 3, 2017 at 9:31 PM, Jarrod Hadfield <j.hadfield at ed.ac.uk> wrote:
Hi Christina, 1/ The model syntax looks fine, it is just that MCMCglmm is not very efficient for this type of problem. You say there are no numerical issues because the latent variable is under 20. However, are they commonly under -20? 2/ Centering and scaling the covariates will not effect the mixing because MCMCglmm is block updating all location effects. Moving from a logistic model (family="categorical") to a probit model (family="threshold") will probably improve mixing, and the inferences will be pretty much the same. 3/ You could also up the number of iterations - but perhaps this takes too much time? Link and Eaton's recommendation not to thin is fine if you are not worried about filling your hard drive. You reduce the Monte Carlo error by saving additional correlated samples, but if the total number of samples you can store is limited you are better storing uncorrelated samples (obtained by thinning) because this reduces the Monte Carlo error more. Cheers, Jarrod On 04/05/2017 01:24, Aiello, Christina wrote:
Dear list, I'm very new to MCMCglmm but have done my best to read-up on Jarrod Hadfield's package documents, tutorials and various examples posted online and discussed on this forum. I'm having trouble fitting what I thought was a fairly simple binomial mixed effects model using MCMCglmm. I'll start by describing the data, then the model, then my problem and questions: My dataset is comprised of unique dyads - pairs of animals located at one of four sites (C1, C2, R1, R2). The response variable, 'contact' indicates that the dyad did (1) or did not (0) interact over the course of the study. The unique id of the members of the dyad are 'tort1' and 'tort2'. Because individuals appear in multiple dyads, I've included a random effect for tortID using the multiple membership function available in the package to account for the non-independence of observations and the fact that some individuals may have a tendency to contact more than others. For fixed effects, in this simplified model I only have one categorical variable, 'site' (which I would have entered as a random effect but I only have 4 levels) and one continuous variable, 'overlap' - which is an estimate of space-use similarity for each dyad. I centered and scaled this variable by the non-zero mean value and standard deviation (though I've also tried the model without centering). This may be relevant to my problem: 'overlap's distribution is highly skewed and mostly zero values - similarly, the response variable 'contact' is rare and characterized by mostly zeros. table(datafi$contact, datafi$site)
C1 C2 R1 R2 0 241 229 176 181 1 35 24 14 9 The model: pr<-list( R= list(V=1, n=0, fix=1), G= list(G1=list(V=1, n=0.002)) ) m1 <- MCMCglmm( fixed = contact ~ (1 + site + overlap ) , random = ~mm(tort1 +tort2), data = datafi, family = "categorical", verbose = FALSE, pr=TRUE, pl=TRUE, prior=pr, nitt=4100000 , thin=2000 , burnin= 100000 )
summary(m1)
Iterations = 100001:4098001
Thinning interval = 2000
Sample size = 2000
DIC: 207.4525
G-structure: ~mm(tort1 + tort2)
post.mean l-95% CI u-95% CI eff.samp
tort1+tort2 2.128 0.0002693 5.488 414.6
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 1 1 1 0
Location effects: contact ~ (1 + site + overlap)
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) -2.2102 -3.6055 -0.7437 1505.6 0.004 **
siteC2 -0.4143 -2.7572 1.4982 1808.4 0.708
siteR1 -1.2543 -4.0794 0.8424 1069.0 0.268
siteR2 -1.4753 -3.9300 0.9782 1348.9 0.205
overlap 3.9025 2.8273 5.1260 488.5 <5e-04 ***
As far as I can tell, the chains themselves look good and if I run
multiple
chains and run the Gelman-Rubin diagnostic, the PSRF values are all 1 or
1.01. The parameter estimates are consistent and make sense. The problem
lies in the autocorrelation - large amounts in the 'overlap' variable and
many of the random intercepts. Here's a sample of the autocorr results:
sort(diag(autocorr(m1$Sol)[2,,]))
##these are the worst offenders
(Intercept) tort1.4534 tort1.3719 tort1.33 tort1.3620
tort1.30 tort1.3045
0.0926484964 0.0938622549 0.1009204749 0.1049459123 0.1065261665
0.1090179501 0.1237642453
siteR2 tort1.3150 tort1.5579 siteR1 tort1.5473
tort1.2051 tort1.3092
0.1339370132 0.1359132027 0.1383816060 0.1506535457 0.1639062068
0.1682852625 0.1683907054
tort1.5044 tort1.804 tort1.5141 tort1.5103 tort1.4148
tort1.4678 tort1.4428
0.1752670493 0.1767909176 0.1865412328 0.1919929722 0.2257633018
0.2318115800 0.2521806794
tort1.3633 tort1.3335 tort1.5101 tort1.3043 tort1.26
tort1.2014 tort1.6
0.2577034325 0.2593673083 0.2602145001 0.2717718040 0.3487288823
0.3748047689 0.4478979043
overlap
0.5400325556
autocorr.diag(m1$VCV)
tort1+tort2 units Lag 0 1.00000000 NaN Lag 2000 0.58292962 NaN Lag 10000 0.12771910 NaN Lag 20000 0.05262786 NaN Lag 1e+05 0.01757316 NaN I've attempted to fit the model with both uninformative (shown above) and parameter expanded priors ( pr2<-list( R= list(V=1, n=0, fix=1), G=list(G1=list(V=1, nu=1, alpha.mu =0, alpha.V=1000)) )), with parameter expanded priors performing slightly worse. I've attempted incrementally larger iteration, thinning, and burn in values, increasing the thinning to as high as 2000 with a large burn-in (100000) in hopes of improving convergence and reducing autocorrelation. I've tried slice sampling and saw little improvement. Nothing I tried while retaining this model structure improved the acfs. I checked the latent variable estimates and all were under 20, with mean of -5. The only way I was able to reduce the autocorrelation was to fit a model without the random effect, which isn't ideal as I'm ignoring repeated measures of individuals among dyads. I've read on this forum that random effects in binomial models are notoriously hard to estimate with this package and I've also read that one should not just increase thinning to deal with the problem (MEE 2012 Link & Eaton <http://onlinelibrary.wiley.com/store/10.1111/j.2041-210X.20 11.00131.x/asset/j.2041-210X.2011.00131.x.pdf?v=1&t=j29nl332 &s=e0f97f28309122f2bbfa66bccb0cd445696e2f15>). Interestingly, I have count responses association with all interacting dyads and I can fit zero truncated models to those responses just fine with the same fixed and random effects. My questions are then: 1) Do you think there is something inherently wrong with the data or just problems with the mixing algorithms in the context of this data? 2) Are there any other changes to the MCMCglmm specification I might try to improve mixing? Or any problems with my current specification? 3) Any suggestions on where to go from here? I would greatly appreciate any insights and happy to provide further info as needed! Christina [[alternative HTML version deleted]]
_______________________________________________ R-sig-mixed-models at 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.