7
$\begingroup$

I study a colonially-nesting bird species. I am trying to perform an AICc evaluation of GLMMs for a nest site selection study. I collected data at nest sites and paired random sites. I want to evaluate habitat characteristics (fixed effect) while accounting for colony sites (random effect) to predict nest site (response variable). Nest site is binary (nest site / paired random location). All data are standardized so as to have a normal distribution ((value-mean)/sd).

I have 5 colonies that have between 4-9 nests in them, which comes out to 38 nests and 38 random pairs that I use in the analysis. When I run the GLMMs in r with lme4 (glmer function), I get a result saying that I have 0 variance for each of 60 models. This is problematic because this doesn't happen with I remove the random effect, but I want to keep it to account for spatial autocorrelation in colonies.

I have read a few articles on this site and others about how the low number of groups may be suppressing my variance, however, when I run the models slightly differently (incorrectly structured, but basically lumping all the paired sites together to form one "colony" with 38 observations and comparing it to the other small colonies of 4-9 nests), it DOES give me a variance. So I don't think the number of groups itself is to blame.. Maybe it's the combination between small samples in group and small number of groups? Still though.. Any help would be great. If I can't figure it out, I think I'd have to resort to a glm fixed effect model without colony as a random effect. It would be incorrect, but it's a start? What would you do if you couldn't get this to work. Thank you.

Something I'd like to stress: I am not a statistician or a programmer. I catch birds. I have taken stats classes and I have a working grasp of some things, but let's face it, sometimes it's fleeting. If answers could please not be too esoteric and focus more on the PRACTICAL, as in, "You should do this. You should do that", I would be much obliged. I'm trying to get this analysis done with ASAP. Thank you again.

Some links I have consulted:

Random effect equal to 0 in generalized linear mixed model

Why do I get zero variance of a random effect in my mixed model, despite some variation in the data?

https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#singular-models-random-effect-variances-estimated-as-zero-or-correlations-estimated-as---1

http://rpubs.com/bbolker/4187

A sample of my code below.

helpmeobiwan <-list(NestPlot = c(1, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0),NumDeadJun = c( 0.1409216, -0.1932639,-0.5274494,-0.5274494, 0.1409216, -0.5274494, -0.5274494 , 0.4751071, -0.5274494 , 2.1460347 ,-0.5274494, -0.1932639, 0.8092926, -0.5274494, -0.5274494 ,-0.5274494 ,-0.1932639, 0.1409216, -0.5274494, -0.5274494 ,-0.5274494, -0.5274494 ,-0.5274494,  0.1409216,-0.5274494, -0.5274494 ,-0.5274494,  0.1409216, -0.5274494,  0.1409216, -0.5274494, -0.5274494, -0.5274494, -0.1932639, -0.1932639, -0.5274494,  0.4751071 , 0.1409216 ,-0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.1932639, -0.5274494, -0.5274494 ,-0.5274494 ,-0.5274494,  0.1409216, -0.5274494, -0.5274494, -0.1932639, -0.5274494, -0.5274494, -0.5274494,  0.1409216, -0.5274494, -0.5274494  ,3.1485912 , 2.4802202,  1.4776637, -0.5274494 , 2.8144057, -0.5274494, -0.5274494,  1.1434781,  3.8169623,  3.8169623 ,-0.1932639, -0.5274494  ,1.4776637 , 1.8118492, -0.5274494),RandomPair = c(  "Madera2" ,  "Starfire1", "Madera2" ,  "Madera3" ,  "Starfire1" ,"Starfire1", "Starfire2", "Madera1" , "Madera3"   ,"Starfire2" ,"Starfire2", "Madera1",   "Madera2",  "Starfire1", "Starfire1" ,"Starfire1", "Madera1",   "Madera2" ,  "Starfire1", "Starfire1", "Starfire1", "Madera1" ,  "Starfire1", "Starfire1", "Madera1",   "Madera1" , "Starfire1", "Madera2" ,  "Madera1",   "Madera2" ,  "Madera1" ,  "Madera1"  , "Starfire1" ,"Starfire1", "Starfire1" ,"Starfire1" ,"Madera2"  , "Madera2",   "Starfire2" ,"Starfire2", "Starfire2" ,"Madera3" ,  "Madera3" ,  "Madera3" ,  "Madera3" ,  "Madera3" ,  "Starfire2", "Starfire2", "Starfire2", "Starfire2" ,"Starfire2", "Madera3",  "Madera3" ,  "Starfire2", "Madera3" ,  "Madera1"  , "Starfire2" ,"Starfire1", "Madera2" ,  "Madera3" ,  "Madera3"  , "Madera2"  , "Madera3"   ,"Starfire2", "Madera3",   "Starfire1", "Madera3"  , "Starfire2", "Starfire1", "Madera3",   "Starfire1", "Starfire2" ,"Madera1" ,  "Starfire2", "Starfire2", "Madera1"  ))




m1 <- glmer(NestPair ~ NumDeadJun+ (1|RandomPair), family=binomial, data=helpmeobiwan)
Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ Can you explain a bit more the connection between your nest sites and the paired random sites? Why are you considering paired random sites in the first place? How many paired random sites are you using per nest site and how do you pick them? Does any pair of random sites include a nest site? Any other details you can provide in an addendum to your original post would be helpful. $\endgroup$
    – Isabella Ghement
    Commented Jun 20, 2019 at 8:45
  • 1
    $\begingroup$ I'm comparing the habitat at nest sites to non-nest sites to see if these birds have a preference for certain habitat characteristics (canopy area, number of trees, etc). If nest sites differ from random, it suggests the birds have a preference. Pairing is a way to eliminate autocorrelation of grouped nest sites. So instead of comparing all nests together against total random points, you compare colony nests together against paired nearby sites and you see how different a nest site is from the places closeby. $\endgroup$
    – PIJA
    Commented Jun 20, 2019 at 14:05
  • $\begingroup$ However I didnt know to collect data in this way when I began, so I have autocorrelated nest sites spatially grouped into colonies, and completely randomized points chosen from across the habitat. I wanted to maintain my colony groups because those nests will be similar to each other, but all I have left are imperfect ways to compare the random sites. Options: 1) Continue as is. Keep colony groups with paired number of random points chosen from across the landscape 2) Lump all random points together as 1 group and compare to each colony group 3) Consider all points independent $\endgroup$
    – PIJA
    Commented Jun 20, 2019 at 14:08

2 Answers 2

Reset to default
3
$\begingroup$

I tried with the glmmADMB package, an alternative of lme4 for linear mixed modelling. You can install this package with this code:

install.packages("R2admb")
install.packages("glmmADMB", 
repos=c("http://glmmadmb.r-forge.r-project.org/repos",
        getOption("repos")),
type="source")

Then you go:

library(glmmADMB)
helpmeobiwan <-list(NestPlot = c(1, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0),NumDeadJun = c( 0.1409216, -0.1932639,-0.5274494,-0.5274494, 0.1409216, -0.5274494, -0.5274494 , 0.4751071, -0.5274494 , 2.1460347 ,-0.5274494, -0.1932639, 0.8092926, -0.5274494, -0.5274494 ,-0.5274494 ,-0.1932639, 0.1409216, -0.5274494, -0.5274494 ,-0.5274494, -0.5274494 ,-0.5274494,  0.1409216,-0.5274494, -0.5274494 ,-0.5274494,  0.1409216, -0.5274494,  0.1409216, -0.5274494, -0.5274494, -0.5274494, -0.1932639, -0.1932639, -0.5274494,  0.4751071 , 0.1409216 ,-0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.1932639, -0.5274494, -0.5274494 ,-0.5274494 ,-0.5274494,  0.1409216, -0.5274494, -0.5274494, -0.1932639, -0.5274494, -0.5274494, -0.5274494,  0.1409216, -0.5274494, -0.5274494  ,3.1485912 , 2.4802202,  1.4776637, -0.5274494 , 2.8144057, -0.5274494, -0.5274494,  1.1434781,  3.8169623,  3.8169623 ,-0.1932639, -0.5274494  ,1.4776637 , 1.8118492, -0.5274494),RandomPair = c(  "Madera2" ,  "Starfire1", "Madera2" ,  "Madera3" ,  "Starfire1" ,"Starfire1", "Starfire2", "Madera1" , "Madera3"   ,"Starfire2" ,"Starfire2", "Madera1",   "Madera2",  "Starfire1", "Starfire1" ,"Starfire1", "Madera1",   "Madera2" ,  "Starfire1", "Starfire1", "Starfire1", "Madera1" ,  "Starfire1", "Starfire1", "Madera1",   "Madera1" , "Starfire1", "Madera2" ,  "Madera1",   "Madera2" ,  "Madera1" ,  "Madera1"  , "Starfire1" ,"Starfire1", "Starfire1" ,"Starfire1" ,"Madera2"  , "Madera2",   "Starfire2" ,"Starfire2", "Starfire2" ,"Madera3" ,  "Madera3" ,  "Madera3" ,  "Madera3" ,  "Madera3" ,  "Starfire2", "Starfire2", "Starfire2", "Starfire2" ,"Starfire2", "Madera3",  "Madera3" ,  "Starfire2", "Madera3" ,  "Madera1"  , "Starfire2" ,"Starfire1", "Madera2" ,  "Madera3" ,  "Madera3"  , "Madera2"  , "Madera3"   ,"Starfire2", "Madera3",   "Starfire1", "Madera3"  , "Starfire2", "Starfire1", "Madera3",   "Starfire1", "Starfire2" ,"Madera1" ,  "Starfire2", "Starfire2", "Madera1"  ))
dontworryLeia <- helpmeobiwan
dontworryLeia$RandomPair <- as.factor(dontworryLeia$RandomPair)
attach(dontworryLeia)
mod <- glmmadmb(NestPlot ~ NumDeadJun + (1|RandomPair), family='binomial', data=dontworryLeia)
mod
summary(mod)
drop1(mod)

First RandomPair was not considered as a factor, which explains this transformation. I guess you meant NestPlot and not NestPair in your m1 model. Anyway this should work!

Improve this answer
$\endgroup$
7
$\begingroup$

Function glmer() uses by default the Laplace approximation, which is not optimal for dichotomous data. A better alternative is the adaptive Gaussian quadrature. You can use this method by setting argument nAGQ of glmer() to a higher number (e.g., 11 or 15) or alternatively using the GLMMadaptive package. In your example, it gives:

library("GLMMadaptive")
helpmeobiwan <- list(NestPlot = c(1, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0),NumDeadJun = c( 0.1409216, -0.1932639,-0.5274494,-0.5274494, 0.1409216, -0.5274494, -0.5274494 , 0.4751071, -0.5274494 , 2.1460347 ,-0.5274494, -0.1932639, 0.8092926, -0.5274494, -0.5274494 ,-0.5274494 ,-0.1932639, 0.1409216, -0.5274494, -0.5274494 ,-0.5274494, -0.5274494 ,-0.5274494,  0.1409216,-0.5274494, -0.5274494 ,-0.5274494,  0.1409216, -0.5274494,  0.1409216, -0.5274494, -0.5274494, -0.5274494, -0.1932639, -0.1932639, -0.5274494,  0.4751071 , 0.1409216 ,-0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.5274494, -0.1932639, -0.5274494, -0.5274494 ,-0.5274494 ,-0.5274494,  0.1409216, -0.5274494, -0.5274494, -0.1932639, -0.5274494, -0.5274494, -0.5274494,  0.1409216, -0.5274494, -0.5274494  ,3.1485912 , 2.4802202,  1.4776637, -0.5274494 , 2.8144057, -0.5274494, -0.5274494,  1.1434781,  3.8169623,  3.8169623 ,-0.1932639, -0.5274494  ,1.4776637 , 1.8118492, -0.5274494),RandomPair = c(  "Madera2" ,  "Starfire1", "Madera2" ,  "Madera3" ,  "Starfire1" ,"Starfire1", "Starfire2", "Madera1" , "Madera3"   ,"Starfire2" ,"Starfire2", "Madera1",   "Madera2",  "Starfire1", "Starfire1" ,"Starfire1", "Madera1",   "Madera2" ,  "Starfire1", "Starfire1", "Starfire1", "Madera1" ,  "Starfire1", "Starfire1", "Madera1",   "Madera1" , "Starfire1", "Madera2" ,  "Madera1",   "Madera2" ,  "Madera1" ,  "Madera1"  , "Starfire1" ,"Starfire1", "Starfire1" ,"Starfire1" ,"Madera2"  , "Madera2",   "Starfire2" ,"Starfire2", "Starfire2" ,"Madera3" ,  "Madera3" ,  "Madera3" ,  "Madera3" ,  "Madera3" ,  "Starfire2", "Starfire2", "Starfire2", "Starfire2" ,"Starfire2", "Madera3",  "Madera3" ,  "Starfire2", "Madera3" ,  "Madera1"  , "Starfire2" ,"Starfire1", "Madera2" ,  "Madera3" ,  "Madera3"  , "Madera2"  , "Madera3"   ,"Starfire2", "Madera3",   "Starfire1", "Madera3"  , "Starfire2", "Starfire1", "Madera3",   "Starfire1", "Starfire2" ,"Madera1" ,  "Starfire2", "Starfire2", "Madera1"  ))
helpmeobiwan <- as.data.frame(helpmeobiwan)

fm <- mixed_model(NestPlot ~ NumDeadJun, random = ~ 1 | RandomPair, 
                  family = binomial(), data = helpmeobiwan)

summary(fm)
#> 
#> Call:
#> mixed_model(fixed = NestPlot ~ NumDeadJun, random = ~1 | RandomPair, 
#>     data = helpmeobiwan, family = binomial())
#> 
#> Data Descriptives:
#> Number of Observations: 76
#> Number of Groups: 5 
#> 
#> Model:
#>  family: binomial
#>  link: logit 
#> 
#> Fit statistics:
#>   log.Lik      AIC      BIC
#>  -46.2248 98.44959 97.27791
#> 
#> Random effects covariance matrix:
#>                StdDev
#> (Intercept) 0.0477673
#> 
#> Fixed effects:
#>             Estimate Std.Err z-value  p-value
#> (Intercept)  -0.1568  0.2829 -0.5544 0.579304
#> NumDeadJun   -1.2274  0.4917 -2.4961 0.012558
#> 
#> Integration:
#> method: adaptive Gauss-Hermite quadrature rule
#> quadrature points: 11
#> 
#> Optimization:
#> method: hybrid EM and quasi-Newton
#> converged: TRUE
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.