# Using mixed effects models or clogit on paired data

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.

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),
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?

• can you post your data? Mar 11 at 19:11
• @GeorgeSavva, I have a portion of the data available here (everything that is needed for the model is included): github.com/rlumkes/bedsite-data Mar 11 at 20:19
• This isn't a programming problem, the models are behaving as expected. Since your data is balanced (2 controls and 1 case per stratum) then the mixed model will be singular (every bedsiteid has the same probability being a case so there is no 'random effect of site'). Also I think the clogit doesn't have enough data to separate (eg) the Cropland from the Forest habitat. Try asking on Cross-Validated for analysis advice. Mar 11 at 20:47
• Thanks that's good to know about the models at least, and I'll post over there! Mar 11 at 20:55