Best analysis for site selection with paired data

Question

I'm running into a problem after trying what I thought would be a simple analysis. I have 47 sites where I measured a variety of habitat characteristics (canopy cover, habitat type, percent of bare ground, elevation, etc.). The habitat type consists of 4 categorical variables, canopy cover is continuous, and then there are the percentages. For each site, I also measured the same characteristics at two random sites, 50 m away. My goal is to see if/how these characteristics are informing the site selection of the original site. Because I'm looking at fine-scale selection, I want the two randoms to be paired to the site. A portion of my data is available at: github.com/rlumkes/bedsite-data

I originally tried a mixed effects model with SiteID as the random effect, but received a singular fit warning. Type refers to site (1) or random (0).

bedsites.random <- glmer(Type ~ Habitat + Canopy_Cover + 
                         X100cm_Cover + (1|BedsiteID), 
                         family = binomial(link = "logit"), 
                         data = bedsites)  

    boundary (singular) fit: see help('isSingular')

I surmised this was from only have one observation for each site, so I tried clogit for case control studies in R, only to end up with this warning and huge beta estimates:

bed.mod <- clogit(Type ~ Habitat + Canopy_Cover + X100cm_Cover + 
                  strata(BedsiteID), data = bedsite)

Warning message:
In coxexact.fit(X, Y, istrat, offset, init, control, weights = weights,  :
  Loglik converged before variable  1,2,4 ; beta may be infinite.

THEN someone told me to try a negative binomial, which resulted in this:

    summary(m1 <- glm.nb(Type ~ Habitat + Canopy_Cover, 
                         data = bedsite))

Call:
glm.nb(formula = Type ~ Habitat + Canopy_Cover, data = bedsite, 
    init.theta = 9353.492376, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.1609  -0.7308  -0.7308   0.5796   1.0996  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -1.343735   0.408254  -3.291 0.000997 ***
HabitatCRP        0.069801   0.578916   0.121 0.904031    
HabitatForest     0.079300   0.879014   0.090 0.928117    
HabitatGrassland  0.023234   0.464190   0.050 0.960081    
HabitatShrubs     0.702339   0.520194   1.350 0.176969    
Canopy_Cover      0.011405   0.006117   1.865 0.062243 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(9353.492) family taken to be 1)

    Null deviance: 103.266  on 140  degrees of freedom
Residual deviance:  95.877  on 135  degrees of freedom
AIC: 203.88

Number of Fisher Scoring iterations: 1


              Theta:  9353 
          Std. Err.:  115219 
Warning while fitting theta: iteration limit reached 

 2 x log-likelihood:  -189.882 
Warning messages:
1: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace >  :
  iteration limit reached
2: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace >  :
  iteration limit reached

I also tried running these models with only one set of randoms in case the two randoms per one observation was throwing it off, but I got the same warnings. I'm completely at a loss for how to analyze this. Any ideas?

Answer 1

With your data, the isSingular warning simply means that the data don't provide enough information to distinguish the variance among random intercepts from a value of 0. For evaluating whether there is some systematic difference between the chosen sites (Type=1) and the nearby "random" sites (Type=0), related to site characteristics, that probably doesn't matter a lot. You would get the same point estimates for coefficients whether or not you include the random intercepts.

A potentially bigger problem is in your modeling of the continuous predictors. For example, your model assumes a linear association between the log-odds of being Type=1 and the percentage of Canopy_Cover. Almost all of your Canopy_Cover values are close to 0, however, with less than 5% of observations at over 50% cover. Does the implicit simple linear association make sense in that context? Similarly, the X100cm_Cover values include a relatively small number of very large values.

A flexible fit of the continuous predictors would be highly preferable. The problem you face is that you only have 47 cases in the minority class, so you are already in risk of overfitting with your model that estimates 6 coefficients beyond the intercept. To avoid overfitting in a binary regression, you typically need about 15 members of the minority class per estimated coefficient.

Two alternative solutions come to mind. One is simply to evaluate the within-site paired differences between the Type=1 and the Type=0 locations for each characteristic, probably with a non-parametric approach. Another is to use a highly flexible model, with a slow learning rate to minimize overfitting, to identify the site characteristics that are most strongly associated with differences between Type=1 and Type=0. That could be, for example, a generalized additive model or a boosted tree.

Answer 2

Okay I was able to figure it out with help from a colleague. The package glmmTMB is a conditional model that was made for a dataset similar to this and makes sure to pair the sites while still looking at selection. Thank you for the feedback everyone!