Using mixed effects models or clogit on 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.

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?