First time poster but have been very grateful over the past couple months for this forum. First and foremost, I apologize in advance if I am not following the right procedures in asking a question.
I am searching for why my models come out singular.
I have 9 farms, and am treating farm as a random effect. I have 3 treatments in each farm with 2 replicates for each treatment (for most of the variables)
My response variable is count data on worms and due to overdispersion I am using a negative binomial distribution.
I have a number of explanatory variables that I reduced be testing for collinearity via spearman coefficients (>0.50) and making decisions on which to keep via biological importance, as well as a few removed that were clearly linearly dependent on one another. I also centered and scaled them.
From a full model of these variables, I built out approx. 120 models with glmer.nb() from lme4 and used model.sel() from the MuMIn package to determine the best model(s).
Unfortunately, many of my (better) models are returning the error "?isSingular" And in checking via
tt<-getME(mod1, "theta")
ll<-getME(mod1, "lower")
min(tt[ll==0])
or
theta <- getME(mod1,"theta")
diag.element <- getME(mod1,"lower")==0
any(theta[diag.element]<1e-5)
the models seem truly singular.
In reading through some stackexchange questions/answers and as well as some good websties (including the following https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#what-methods-are-available-to-fit-estimate-glmms) I've decided for my purposes, getting rid of the only random effect term does not make sense and I will not attempt to pursue a Bayesian approach.
But I would really like to understand why some of my models are singular and others are not.
What I've noted thus far:
Out of the 8 models with singularity problems, they all contain 4 of the same variables , one of these variables SNH_perc, has the same value throughout each farm due to it being measured at a kilometer radius from the center of the farm. I tried creating a little bit of noise in these values but this did not fix the singularity.
In using
VarCorr(mod1)
those with singularity problems have an extremely small value 1e-6 or so, while those without have about 0.1, indicating these models (often with more variables?) have less variation when more dimensions (or certain dimensions) are addedI got the eigenvalues of the Variance/Covariance Choleski Decomposition and they seem off but I am not sure how to interpret them either
an example for one model being
mod1 <- glmer.nb(abundance_SUM ~ SOM_perc + P2O5 + Cu + silt_percent + SNH_perc + (1|vineyard), data=firstrun.scale)
eigen() decomposition
$values
[1] 0.07609357, 0.07371038, 0.06832976, 0.06689690, 0.06391916, 0.06232956
$vectors
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] -0.2138548 0.2237181 0.64555958 0.77197345 1 0.922543154
[2,] -0.8997568 0.9746539 0.07141780 0.04855806 0 -0.007181145
[3,] -0.3617484 0.0000000 0.75959292 0.63379737 0 -0.380728292
[4,] 0.1176516 0.0000000 -0.01901948 0.00000000 0 -0.047982131
[5,] 0.0000000 0.0000000 0.02844624 0.00000000 0 0.034387977
[6,] 0.0000000 0.0000000 0.00000000 0.00000000 0 0.020584223
a sample of my dataset with the relevant explanatory variables (I'm not sure how much of it is useful to provide).
abundance_SUM farm treatment SOM_perc P2O5 Mg Cu silt_percent SNH_perc SNH_perc2
1 8 AT01 BG 7.598926 26.1835 13.40 1.66 70.65202 55.6 55.6
2 4 AT01 BG 7.964178 32.4680 14.00 1.68 71.18537 55.6 55.6
3 22 AT01 AC 7.848078 23.1355 13.40 1.66 70.92315 55.6 55.5
4 14 AT01 AC 8.192078 28.0170 15.95 1.51 71.15727 55.6 55.5
5 20 AT01 CC 8.058861 28.0720 13.65 1.46 67.84648 55.6 55.4
6 12 AT01 CC 8.239543 28.0400 13.25 2.06 70.80634 55.6 55.4
7 7 AT02 BG 6.622578 80.4045 13.35 7.04 57.07563 35.7 35.7
8 10 AT02 BG 6.962443 72.7000 14.30 4.89 57.49698 35.7 35.7
9 13 AT02 AC 7.263278 62.7370 13.95 6.71 61.04811 35.7 35.8
10 11 AT02 AC 6.953843 63.4000 12.15 7.68 61.04811 35.7 35.8
11 13 AT02 CC 7.258978 62.4965 14.40 54.15 65.65210 35.7 35.9
12 10 AT02 CC 7.074078 59.1835 12.75 52.07 65.65210 35.7 35.9
13 46 AT03 BG 4.194980 12.3990 20.05 10.71 46.27264 59.7 59.7
14 46 AT03 BG 4.704530 15.9010 20.50 10.11 46.27264 59.7 59.7
15 28 AT03 AC 5.749596 20.0230 17.70 8.04 47.82043 59.7 59.6
16 34 AT03 AC 5.874213 22.9345 18.70 10.09 47.82043 59.7 59.6
17 33 AT03 CC 5.169096 17.9850 21.50 9.59 46.08007 59.7 59.8
18 34 AT03 CC 5.267913 18.4005 23.25 10.52 46.08007 59.7 59.8
19 26 AT04 BG 5.865696 11.0740 10.15 15.10 42.56881 7.7 7.7
20 42 AT04 BG 5.753896 11.2875 10.15 12.89 42.56881 7.7 7.7
21 16 AT04 AC 3.861730 24.7900 13.45 29.98 34.74912 7.7 7.8
22 34 AT04 AC 3.268413 24.7310 11.65 30.02 34.74912 7.7 7.8
23 10 AT04 CC 4.936813 21.2350 10.70 21.65 34.88536 7.7 7.6
24 8 AT04 CC 4.807813 20.0090 11.20 21.28 34.88536 7.7 7.6
25 3 AT05 BG 4.395096 37.5385 13.70 78.91 59.70392 12.3 12.3
26 5 AT05 BG 4.863630 41.0965 13.75 75.80 59.70392 12.3 12.3
27 4 AT05 AC 5.839730 31.6870 14.70 88.75 62.55672 12.3 12.4
28 9 AT05 AC 5.237813 43.3220 12.60 80.68 62.55672 12.3 12.4
29 2 AT05 CC 5.710813 32.4700 14.95 106.10 58.84137 12.3 12.5
30 9 AT05 CC 4.885213 34.1055 12.00 99.67 58.84137 12.3 12.5
Thank you in advance for any hints in the right direction!
EDIT: Adding the summary output for both glm and glmm version of best model
Call:
MASS::glm.nb(formula = abundance_SUM ~ SOM_perc + P2O5 + Cu +
silt_percent + SNH_perc, data = firstrun.scale, init.theta = 6.692239628,
link = log)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.29715 -0.91392 -0.00525 0.60202 2.14098
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.75125 0.06393 43.038 < 2e-16 ***
SOM_perc -0.18855 0.07293 -2.585 0.009730 **
P2O5 -0.20610 0.06948 -2.966 0.003013 **
Cu -0.18915 0.08093 -2.337 0.019433 *
silt_percent -0.25421 0.07006 -3.628 0.000285 ***
SNH_perc 0.23310 0.07741 3.011 0.002603 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Negative Binomial(6.6922) family taken to be 1)
Null deviance: 108.025 on 53 degrees of freedom
Residual deviance: 55.154 on 48 degrees of freedom
AIC: 375.88
Number of Fisher Scoring iterations: 1
Theta: 6.69
Std. Err.: 1.85
2 x log-likelihood: -361.88
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(6.6922) ( log )
Formula: abundance_SUM ~ SOM_perc + P2O5 + Cu + silt_percent + SNH_perc + (1 | vineyard)
Data: firstrun.scale
AIC BIC logLik deviance df.resid
377.9 393.8 -180.9 361.9 46
Scaled residuals:
Min 1Q Median 3Q Max
-1.66497 -0.80737 -0.00517 0.65162 2.77370
Random effects:
Groups Name Variance Std.Dev.
vineyard (Intercept) 3.389e-11 5.821e-06
Number of obs: 54, groups: vineyard, 9
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.75126 0.06392 43.043 < 2e-16 ***
SOM_perc -0.18855 0.07374 -2.557 0.010564 *
P2O5 -0.20610 0.06699 -3.077 0.002093 **
Cu -0.18915 0.08376 -2.258 0.023925 *
silt_percent -0.25421 0.07260 -3.501 0.000463 ***
SNH_perc 0.23310 0.07526 3.097 0.001952 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) SOM_pr P2O5 Cu slt_pr
SOM_perc 0.030
P2O5 0.052 -0.006
Cu 0.068 -0.235 -0.332
silt_percnt 0.083 -0.174 0.272 0.022
SNH_perc -0.044 -0.406 -0.182 0.437 -0.084
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see ?isSingular
summary(mod1)
add the summary output from the modelglm.nb(abundance_SUM ~ SOM_perc + P2O5 + Cu + silt_percent + SNH_perc, data=firstrun.scale)
$\endgroup$