mcmcglmm Priors: Auto correlation and extreme post mean values

Thanks for sharing the relevant data and model details and output. The
first and biggest concern I have is the use of Site variable for a Random
Intercept. There are only 175 observations, but you have 86 different
sites. On an average, without seeing the actual data, that means you have
~2 observations per site. There is not enough data to use Site for a Random
intercept as it exists.

You should use site only after you can reduce 86 sites into some
meaningful groupings based on criteria you use in Ecology (reduce to 3 or 4
groups). Additionally, the remaining factors create 2 * 5 * 4 = 40 cells
and you have only 175 observations ~ 4.4 observations per cell - this
number is also too low for you to use all of these factors in your
analysis.

All of the problems you are facing mostly likely are arising from the
wrong application of the techniques to your data. I would take the
following steps:
1. Run a bivariate means comparison using Anova after you declare your
Landuse, Human_presence etc as factors using ( aov( trappings ~ landuse)
2. In doing so you are going to change your assumptions of trappings from
count to continuous, but I think for exploration, it is okay.
3. Make sure you examine the means of the trappings for each level and
combine levels within a factor, when the sub-groups are not different from
each other.
4. This will reduce the number of groups you are working with since your
records are limited.

Only after undertaking the steps above should you consider, running any
models in lme4 or nlme. Ideally, if you could reduce those factors that
have 4 and 5 groups to having 2 each (if the differences are not
significant - see # 3) then you could do meaningful analyses. Similarly,
with Site, the 86 levels should be reduced to either 3 or 4 for you to use
them even for a descriptive analysis. Without adding in more data or
reducing your number of levels, you should not undertake any further
modeling.

Sree

On Fri, May 15, 2020 at 11:33 AM Alexander Botha <alexbotha555 at gmail.com>
wrote:

Hi Sree.
I have run multiple models with lmer and glmer using a combination of
predictor variables as either random or fixed effects. I was worried about
the collinearity in fixed effects (singular fit) when using lme4.
I hope these answer your questions. Please see below summaries of some of
the models I ran.
I have also attached the structure of my data in case it might help.
Thank you for your help.
Please let me know if you require anything else or if you have any
suggestions for data exploration.

*Data structure*
data.frame': 175 obs. of  8 variables:
 $ Site          : Factor w/ 86 levels "AA1","carcass",..: 66 66 66 68 68
68 21 21 21 32 ...
 $ Landuse       : Factor w/ 2 levels "farmland","reserve": 1 1 1 1 1 1 1
1 1 1 ...
 $ MP            : Factor w/ 4 levels "FM","FQ","NM",..: 3 1 4 3 2 4 3 2
1 3 ...
 $ Human_presence: Factor w/ 5 levels "0","1","2","3",..: 4 4 4 4 4 4 5 5
5 5 ...
 $ Jackal_trapped: int  1 0 0 0 1 0 0 0 0 0 ...

*model1<-lmer(Jackal_trapped~ Human_presence + (1|Site), data=data)*
REML criterion at convergence: 77.7

Scaled residuals:
    Min      1Q  Median      3Q     Max
-1.3261 -0.7568  0.0000  0.0000  2.6909

Random effects:
 Groups   Name        Variance Std.Dev.
 Site     (Intercept) 0.00000  0.00
 Residual             0.08413  0.29
Number of obs: 175, groups:  Site, 86

Fixed effects:
                Estimate Std. Error t value
(Intercept)      0.21951    0.04530   4.846
Human_presence1  0.16510    0.09232   1.788
Human_presence2 -0.21951    0.10230  -2.146
Human_presence3  0.08049    0.07911   1.017
Human_presence4 -0.21951    0.05456  -4.024

Correlation of Fixed Effects:
            (Intr) Hmn_p1 Hmn_p2 Hmn_p3
Humn_prsnc1 -0.491
Humn_prsnc2 -0.443  0.217
Humn_prsnc3 -0.573  0.281  0.254
Humn_prsnc4 -0.830  0.407  0.368  0.475
convergence code: 0

*model1.1<-lmer(Jackal_trapped~ Human_presence + MP + Landuse + (1|Site),
data=data)*
REML criterion at convergence: 88.5

Scaled residuals:
     Min       1Q   Median       3Q      Max
-1.36759 -0.73677 -0.00991  0.03404  2.68699

Random effects:
 Groups   Name        Variance Std.Dev.
 Site     (Intercept) 0.00000  0.0000
 Residual             0.08531  0.2921
Number of obs: 175, groups:  Site, 86

Fixed effects:
                 Estimate Std. Error t value
(Intercept)      0.344378   0.181977   1.892
Human_presence1  0.161919   0.094159   1.720
Human_presence2 -0.228924   0.103472  -2.212
Human_presence3  0.014851   0.154207   0.096
Human_presence4 -0.341485   0.179179  -1.906
MPFQ            -0.012836   0.060249  -0.213
MPNM             0.009494   0.060567   0.157
MPTQ            -0.122941   0.114023  -1.078
Landusereserve  -0.116350   0.169306  -0.687

Correlation of Fixed Effects:
            (Intr) Hmn_p1 Hmn_p2 Hmn_p3 Hmn_p4 MPFQ   MPNM   MPTQ
Humn_prsnc1 -0.159
Humn_prsnc2 -0.144  0.221
Humn_prsnc3 -0.868  0.159  0.138
Humn_prsnc4 -0.963  0.131  0.134  0.868
MPFQ        -0.248  0.158  0.052  0.094  0.052
MPNM        -0.178  0.082  0.047  0.003  0.002  0.594
MPTQ        -0.327  0.040  0.090  0.063  0.260  0.293  0.322
Landusersrv -0.929  0.004  0.017  0.838  0.947  0.021 -0.080  0.207
convergence code: 0
boundary (singular) fit: see ?isSingular

*model2<-lmer(Jackal_trapped~ Human_presence + Landuse + (1|Site),
data=data)*
Linear mixed model fit by REML ['lmerMod']
Formula: Jackal_trapped ~ Human_presence + Landuse + (1 | Site)
   Data: data

REML criterion at convergence: 79.4

Scaled residuals:
    Min      1Q  Median      3Q     Max
-1.3227 -0.7549  0.0000  0.0000  2.6842

Random effects:
 Groups   Name        Variance Std.Dev.
 Site     (Intercept) 0.00000  0.0000
 Residual             0.08455  0.2908
Number of obs: 175, groups:  Site, 86

Fixed effects:
                Estimate Std. Error t value
(Intercept)      0.28201    0.16877   1.671
Human_presence1  0.16510    0.09255   1.784
Human_presence2 -0.21951    0.10255  -2.140
Human_presence3  0.03049    0.15231   0.200
Human_presence4 -0.28201    0.17150  -1.644
Landusereserve  -0.06250    0.16255  -0.385

Correlation of Fixed Effects:
            (Intr) Hmn_p1 Hmn_p2 Hmn_p3 Hmn_p4
Humn_prsnc1 -0.132
Humn_prsnc2 -0.119  0.217
Humn_prsnc3 -0.902  0.146  0.132
Humn_prsnc4 -0.984  0.130  0.117  0.888
Landusersrv -0.963  0.000  0.000  0.854  0.948
convergence code: 0

*model3<-lmer(Jackal_trapped~ Human_presence + MP + (1|Landuse) (1|Site),
data=data)*
Formula: Jackal_trapped ~ Human_presence + MP + (1 | Landuse) + (1 | Site)
   Data: data

REML criterion at convergence: 87.2

Scaled residuals:
     Min       1Q   Median       3Q      Max
-1.35994 -0.74172 -0.01178  0.02922  2.68746

Random effects:
 Groups   Name        Variance Std.Dev.
 Site     (Intercept) 0.00000  0.0000
 Landuse  (Intercept) 0.00000  0.0000
 Residual             0.08504  0.2916
Number of obs: 175, groups:  Site, 86; Landuse, 2

Fixed effects:
                 Estimate Std. Error t value
(Intercept)      0.228251   0.067428   3.385
Human_presence1  0.162180   0.094010   1.725
Human_presence2 -0.227740   0.103294  -2.205
Human_presence3  0.103705   0.083907   1.236
Human_presence4 -0.224817   0.057214  -3.929
MPFQ            -0.011955   0.060140  -0.199
MPNM             0.006149   0.060276   0.102
MPTQ            -0.106692   0.111368  -0.958

Correlation of Fixed Effects:
            (Intr) Hmn_p1 Hmn_p2 Hmn_p3 Hmn_p4 MPFQ   MPNM
Humn_prsnc1 -0.418
Humn_prsnc2 -0.347  0.221
Humn_prsnc3 -0.443  0.285  0.227
Humn_prsnc4 -0.705  0.398  0.371  0.422
MPFQ        -0.614  0.158  0.052  0.139  0.100
MPNM        -0.682  0.083  0.049  0.130  0.245  0.598
MPTQ        -0.370  0.040  0.088 -0.208  0.202  0.295  0.347
convergence code: 0
boundary (singular) fit: see ?isSingular
*model4<-glmer(Jackal_trapped~Human_presence + MP + Landuse +(1|Site),
data = data, family = binomial)*

Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: Jackal_trapped ~ Human_presence + MP + Landuse + (1 | Site)
   Data: data

     AIC      BIC   logLik deviance df.resid
   104.1    135.7    -42.0     84.1      165

Scaled residuals:
    Min      1Q  Median      3Q     Max
-0.8317 -0.4895  0.0000  0.0000  2.0429

Random effects:
 Groups Name        Variance Std.Dev.
 Site   (Intercept) 0        0
Number of obs: 175, groups:  Site, 86

Fixed effects:
                  Estimate Std. Error z value Pr(>|z|)
(Intercept)     -5.455e-01  1.511e+00  -0.361    0.718
Human_presence1  7.428e-01  7.188e-01   1.033    0.301
Human_presence2 -3.302e+02  2.122e+07   0.000    1.000
Human_presence3  3.143e-02  1.238e+00   0.025    0.980
Human_presence4 -4.416e+01  7.035e+06   0.000    1.000
MPFQ            -2.413e-01  8.942e-01  -0.270    0.787
MPNM             7.619e-02  7.622e-01   0.100    0.920
MPTQ            -7.063e-01  1.039e+00  -0.680    0.497
Landusereserve  -6.420e-01  1.331e+00  -0.482    0.630

Correlation of Fixed Effects:
            (Intr) Hmn_p1 Hmn_p2 Hmn_p3 Hmn_p4 MPFQ   MPNM   MPTQ
Humn_prsnc1 -0.251
Humn_prsnc2  0.000  0.000
Humn_prsnc3 -0.861  0.210  0.000
Humn_prsnc4  0.000  0.000  0.000  0.000
MPFQ            -0.413  0.296  0.000  0.166  0.000
MPNM            -0.317  0.192  0.000  0.044  0.000  0.646
MPTQ             -0.383  0.116  0.000  0.083  0.000  0.430  0.459
Landusersrv -0.860  0.007  0.000  0.838  0.000  0.035 -0.115  0.158
convergence code: 0

*model5<-glmer(Jackal_trapped~Human_presence + MP + Landuse +(1|Site),
data = data, family = poisson)*
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: Jackal_trapped ~ Human_presence + MP + Landuse + (1 | Site)
   Data: data

     AIC      BIC   logLik deviance df.resid
   110.7    142.3    -45.3     90.7      165

Scaled residuals:
    Min      1Q  Median      3Q     Max
-0.6418 -0.4413  0.0000  0.0000  1.8248

Random effects:
 Groups Name        Variance Std.Dev.
 Site   (Intercept) 0        0
Number of obs: 175, groups:  Site, 86

Fixed effects:
                  Estimate Std. Error z value Pr(>|z|)
(Intercept)     -1.008e+00  1.278e+00  -0.789    0.430
Human_presence1  5.166e-01  5.851e-01   0.883    0.377
Human_presence2 -3.688e+01  2.122e+07   0.000    1.000
Human_presence3  3.101e-02  1.072e+00   0.029    0.977
Human_presence4 -3.153e+01  1.255e+06   0.000    1.000
MPFQ            -1.762e-01  7.501e-01  -0.235    0.814
MPNM             5.648e-02  6.197e-01   0.091    0.927
MPTQ            -5.145e-01  8.896e-01  -0.578    0.563
Landusereserve  -4.524e-01  1.134e+00  -0.399    0.690

Correlation of Fixed Effects:
            (Intr) Hmn_p1 Hmn_p2 Hmn_p3 Hmn_p4 MPFQ   MPNM   MPTQ
Humn_prsnc1 -0.267
Humn_prsnc2  0.000  0.000
Humn_prsnc3 -0.875  0.226  0.000
Humn_prsnc4  0.000  0.000  0.000  0.000
MPFQ        -0.392  0.302  0.000  0.160  0.000
MPNM        -0.299  0.185  0.000  0.039  0.000  0.618
MPTQ        -0.345  0.119  0.000  0.076  0.000  0.395  0.431
Landusersrv -0.865  0.010  0.000  0.848  0.000  0.033 -0.116  0.139
convergence code: 0

With kind regards,

Alexander Edward Botha
alexbotha555 at gmail.com 082 414 9030
PhD candidate
Mammal ecology


On Fri, May 15, 2020 at 4:25 PM sree datta <sreedta8 at gmail.com> wrote:

Hi Alexander,

Prior to running the model with "MCMCglmm", have you attempted to run a
log-linear model / a mixed-model with "nlme" or "lme4" / a decision tree? I
ask these questions to better understand what the simpler univariate and
multivariate approaches would reveal in terms of association between your
dependent variable and your predictor variables.

These are the exploratory models I would run to understand what would
represent reasonable priors to use for Bayesian based models. If you could
share more on some of your explorations with the data, that would be
helpful.

Sree

On Fri, May 15, 2020 at 9:29 AM Alexander Botha <alexbotha555 at gmail.com>
wrote:

Good day List,
My name is Alex, I am currently using the package mcmcglmm to determine
the
impact of lunar cycles, human presence and land use type (agricultural
vs
protected) on the trapping success of meso-predators for my PhD. I am
new
to MCMCglmm and I was wondering if you could assist with my problem.

Structure of data: I am testing if and how the change in lunar cycle
(factor) and human presence ( factor) impact trapping success. Trapping
success is count data but because we only trapped 1 individual at most,
the
data can also be considered binary. Human presence (HP) is split into 5
categories.

Model: I am running Human presence, Moon phase and Land use type as
fixed
effects and the site name as a random effect.

Problem: I have quantified human presence in various locations, with
areas
exposed to intense human pressures, there is a complete absence of any
successful trappings (HP2 and HP4 have no trapping success resulting in
only zeros). Post mcmcglmm, these parameters display auto correlation in
their graphs (if i set them as fixed or random effects) as well as when
using the geweke and gelman tests, but almost perfect mixing for the
others, especially when I run it with 5-20 million iterations. I also
see
poorly mixed VCV graphs if the degree of belief is less than 20. I have
used a variety of priors, see below, with no success. I have also
increased
the amount of itterations to 20 million, decreased the thinning factor,
increased the burnin and used poisson, zapoisson, zipoisson, ordinal and
threshold distributions, with threshold and ordinal having the best DIC
values. Together with this, I am also seeing extreme post mean values
for
the models that are displaying the lowest DIC values, which is not
something that I see in any of the literature (see below).

I have read Jarrod Hadfields awesome course and tutorial notes, as well
as
other literature in my field, online tutorials, information documents
and
the vignettes and help functions in R as well as the correspondence on
this
email list between Jarrod Hadfield and other users but I cannot seem to
figure out why the above is happening. My data set is quite small (176
entries in total, with about 40 entries per HP category) and according
to
what I have read, priors can have a large impact on small to moderately
sized data sets.

Questions:
1. Does the problem lie with my priors? And if so, do you have any tips
on
how to solve it?
2. Should I be using ordinal or threshold family? I have read that the
mixing of zi and za poisson models can be poor, and from my tries with
them, this seems to be the case.
3. Are either the extreme post mean values or the auto correlation in
certain parameters avoidable in this case?

*Examples of priors*
prior1<-list(R=list(V=diag(1)*1e-8,
nu=0.2),G=list(G1=list(V=diag(1)*1e-6,
nu=0.2))
prior2<-list(R=list(V=diag(1), nu=0.2),G=list(G1=list(V=diag(1), nu=2)))
prior3<-list(R=list(V=diag(2), nu=0, fix=2), G=list(G1=G1))
G1=list(V=diag(2), nu=2, alpha.mu=c(0,0), alpha.V=diag(2)*1000)
prior4:
A note: I have these with degree belief ranging from 0.0002 to 20, with
the
best mixing (apart from the auto correlation in certain parameters)
with nu
set at 20 for the R structure. If i set the G structures degree of
belief
to less than 20, it mixes the random effect extremely poorly. I have
also
used ranging values from 1e-1 to 1e-20 and fixed the variances at 1 and
2
using various different priors.
*Examples of model scripts:*
OM1.1<-MCMCglmm(Jackal_trapped ~ Human_presence +  MP + Landuse, random
=~Site, data = data, prior = prior1,
family = "ordinal", nitt = 5000000, thin = 500, burnin = 500000,
verbose =
FALSE, pl = TRUE, DIC=TRUE)
model1.1<-MCMCglmm(Jackal_trapped~trait,random=~us(trait):Human_presence
+

us(trait):Landuse,family="zapoisson",rcov=~idh(trait):units,data=data,prior=prior,nitt=500000,thin=500,burnin=200000,verbose=FALSE)

model1.2<-MCMCglmm(Jackal_trapped~trait,random=~idh(trait):Human_presence
+

idh(trait):Landuse,family="zapoisson",rcov=~idh(trait):units,data=data,prior=prior,nitt=500000,thin=500,burnin=200000,verbose=FALSE)

*Example of large post mean summary*
 family: I have used zapoisson, ordinal and zipoisson distributions
prior1<-list(R=list(V=diag(1)*1e-4,
nu=0.2),G=list(G1=list(V=diag(1)*1e-2,
nu=2)))
OM1.1<-MCMCglmm(Jackal_trapped ~ Human_presence +  MP + Landuse, random
=~Site, data = data, prior = prior1, family = "family", nitt = 5000000,
thin = 500, burnin = 500000, verbose = FALSE, pl = TRUE, DIC=TRUE)
 Iterations = 500001:4999501
 Thinning interval  = 500
 Sample size  = 9000

 DIC: 0.003093406

 G-structure:  ~Site

     post.mean l-95% CI u-95% CI eff.samp
Site    0.1943 0.001154   0.1913     9000

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units 1.755e+10 1.72e+09 3.97e+10     3632

 Location effects: Jackal_trapped ~ Human_presence + MP + Landuse
                              post.mean    l-95% CI  u-95% CI eff.samp
 pMCMC
(Intercept)                  -107789  -229797     9915     4666
0.0498 *
Human_presence1     49459   -12241      117809     6957         0.1002
Human_presence2    -49280  -153682      46279     8478          0.3184
Human_presence3     18642   -66603      104880     8529          0.6684
Human_presence4   -214029  -340008     -89259     6098          <1e-04
***
MPFQ                         -44627  -119345      24385       7328
0.1856
MPNM                          12362   -52404      72882       9000
0.6691
MPTQ                           -43907  -122103    22055     8357
 0.1942
Landusereserve            -51105  -138603    35033     9000
 0.2302

*Example of model outputs with auto correlation*
family: I have used zapoisson, ordinal and zipoisson distributions
prior2<-list(R=list(V=diag(1)*1e-6,
nu=20),G=list(G1=list(V=diag(1)*1e-6,
nu=20)))
OM2.1<-MCMCglmm(Jackal_trapped ~ Human_presence +  MP + Landuse, random
=~Site, data = data, prior = prior2,
family = "family", nitt = 5000000, thin = 500, burnin = 500000, verbose
=
FALSE, pl = TRUE, DIC=TRUE)
 Iterations = 500001:4999501
 Thinning interval  = 500
 Sample size  = 9000

 DIC: 182.6921

 G-structure:  ~Site
     post.mean  l-95% CI  u-95% CI eff.samp
Site 1.114e-06 5.171e-07 1.911e-06     8515

 R-structure:  ~units

      post.mean  l-95% CI u-95% CI eff.samp
units 1.106e-06 5.129e-07 1.87e-06     9000

 Location effects: Jackal_trapped ~ Human_presence + MP + Landuse

                             post.mean l-95% CI u-95% CI eff.samp
 pMCMC

(Intercept)               -0.70349 -1.20259 -0.28723    2.027     <
1e-04
***
Human_presence1   0.53752  0.41033  0.67240    6.030    < 1e-04 ***
Human_presence2  -0.44909 -0.65858 -0.03042    2.600    0.00289 **
Human_presence3  -0.11236 -0.48808  0.30079    1.683     0.53556
Human_presence4  -1.38806 -1.74144 -0.93892    2.829    < 1e-04 ***
MPFQ                        -0.38918 -0.57513 -0.23967    2.300    <
1e-04
***
MPNM                         0.14273 -0.05171  0.28830    1.625
0.14356

MPTQ                         -0.25338 -0.52899 -0.03162    1.477
0.01222
*
Landusereserve          -0.64469 -1.06522 -0.18880    2.239    < 1e-04
***

*Example: Using my variables as random effects*
*family: I have used zapoisson, ordinal and zipoisson distributions*
model1.1<-MCMCglmm(Jackal_trapped~trait,random=~us(trait):Human_presence
+

us(trait):Landuse,family="family",rcov=~idh(trait):units,data=data,prior=prior,nitt=500000,thin=500,burnin=200000,verbose=FALSE)
Iterations = 200001:499501
Thinning interval  = 500
Sample size  = 600

 DIC: 175.5933

 G-structure:  ~us(trait):Human_presence


        post.mean   l-95% CI  u-95% CI eff.samp
traitJackal_trapped:traitJackal_trapped.Human_presence 1.179e-06
5.584e-07
2.051e-06    472.7
traitza_Jackal_trapped:traitJackal_trapped.Human_presence    -1.243e-08
-5.964e-07 5.494e-07    600.0
traitJackal_trapped:traitza_Jackal_trapped.Human_presence    -1.243e-08
-5.964e-07 5.494e-07    600.0
traitza_Jackal_trapped:traitza_Jackal_trapped.Human_presence  1.172e-06
 5.252e-07 1.883e-06    507.3

 ~us(trait):Landuse

      post.mean   l-95% CI  u-95% CI eff.samp
traitJackal_trapped:traitJackal_trapped.Landuse       1.143e-06
4.876e-07
2.029e-06      600
traitza_Jackal_trapped:traitJackal_trapped.Landuse    7.237e-09
-5.758e-07
4.980e-07      600
traitJackal_trapped:traitza_Jackal_trapped.Landuse    7.237e-09
-5.758e-07
4.980e-07      600
traitza_Jackal_trapped:traitza_Jackal_trapped.Landuse 1.157e-06
4.896e-07
2.022e-06      600

 R-structure:  ~idh(trait):units

                             post.mean  l-95% CI  u-95% CI eff.samp
traitJackal_trapped.units    9.855e-07 4.689e-07 1.685e-06    704.6
traitza_Jackal_trapped.units 1.017e-06 4.490e-07 1.625e-06    600.0

 Location effects: Jackal_trapped ~ trait

                       post.mean l-95% CI u-95% CI eff.samp  pMCMC
(Intercept)               0.7789   0.7730   0.7875    8.419 <0.002 **
traitza_Jackal_trapped   -2.8911  -2.9052  -2.8741    6.982 <0.002 **

I hope I have shared enough info regarding my problem and that it makes
sense. I hope my post meets the requirements of the list and that it
does
not seem like I am making my problem somebody elses, I am honestly just
lost at the moment.
I welcome any constructive criticism and any other help you can provide.

Thank you for your help.
I look forward to your responses.

With kind regards,

Alexander Edward Botha
alexbotha555 at gmail.com 082 414 9030
PhD candidate
Mammal ecology

        [[alternative HTML version deleted]]