2
$\begingroup$

I'm new to the world of R and statistical modeling and struggling to find an appropriate way to visualize the results of a generalized linear mixed model with a negative binomial distribution (glmer.nb from the lme4 package). I'm using a dataset that looks something like this:

 $ Year             : Factor w/ 2 levels "2017","2018"
 $ Reampled         : Factor w/ 2 levels "No","Yes"
 $ Habitat          : Factor w/ 4 levels "MCF","MM","MMF","REGEN"
 $ Site             : Factor w/ 63 levels "MCF_001","MCF_002",...
 $ Disturbed        : Factor w/ 2 levels "Disturbed","Mature"
 $ Species          : Factor w/ 3 levels "MYLU","MYSE", "NoID"
 $ Count            : int  1 3 0 1 0 0 6 38 3 43 ...
 $ Bat.Survey.Nights: int  4 4 4 5 5 5 6 6 6 4 ...
 $ Avg.Snags        : num  -0.855 -0.855 -0.855 1.846 1.846 ...
 $ Avg.Understory   : num  -0.00715 -0.00715 -0.00715 -0.94871 -0.94871 ...
 $ Avg.Midstory     : num  -0.352 -0.352 -0.352 0.256 0.256 ...
 $ Avg.Canopy       : num  -1.061 -1.061 -1.061 0.695 0.695 ...
 $ Avg.Canopy.Cover : num  -0.831 -0.831 -0.831 0.506 0.506 ...
 $ Perc.Dec.Dom     : num  -0.493 -0.493 -0.493 -1.095 -1.095 ...

My model is based off of bat counts per site (with an offset to account for number of survey nights). Here, I am comparing vegetation between different bat species:

>nglmer.veg <- glmer.nb(Count ~ Avg.Snags + Avg.Understory*Species + 
                                  Avg.Midstory*Species + Avg.Canopy.Cover*Species + Perc.Dec.Dom +
                                  offset(log(Bat.Survey.Nights)) + (1|Site),
                                data = insect.data)
>summary(glmer.veg)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(1.0954)  ( log )
Formula: Count ~ Avg.Snags + Avg.Understory * Species + Avg.Midstory *  
    Species + Avg.Canopy.Cover * Species + Perc.Dec.Dom + offset(log(Bat.Survey.Nights)) +      (1 | Site)
   Data: insect.data
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))

     AIC      BIC   logLik deviance df.resid 
  1134.0   1187.4   -551.0   1102.0      191 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.0407 -0.6660 -0.2981  0.4499  3.9381 

Random effects:
 Groups Name        Variance Std.Dev.
 Site   (Intercept) 0.6728   0.8203  
Number of obs: 207, groups:  Site, 36

Fixed effects:
                             Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   -0.6367     0.2069  -3.077 0.002090 ** 
Avg.Snags                      0.7330     0.2026   3.617 0.000298 ***
Avg.Understory                 1.0821     0.2209   4.899 9.61e-07 ***
SpeciesMYSE                   -0.3875     0.2236  -1.733 0.083012 .  
SpeciesNoID                    1.4968     0.1998   7.490 6.90e-14 ***
Avg.Midstory                   0.5031     0.2037   2.470 0.013522 *  
Avg.Canopy.Cover              -0.3813     0.2338  -1.631 0.102914    
Perc.Dec.Dom                   0.9980     0.2163   4.614 3.96e-06 ***
Avg.Understory:SpeciesMYSE    -0.6613     0.2040  -3.241 0.001190 ** 
Avg.Understory:SpeciesNoID    -0.5353     0.1992  -2.687 0.007205 ** 
SpeciesMYSE:Avg.Midstory      -1.3015     0.2883  -4.514 6.36e-06 ***
SpeciesNoID:Avg.Midstory      -0.4005     0.1864  -2.149 0.031660 *  
SpeciesMYSE:Avg.Canopy.Cover   0.5597     0.2464   2.272 0.023114 *  
SpeciesNoID:Avg.Canopy.Cover   0.5246     0.2196   2.389 0.016916 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

What is the best way to visualize these results?

$\endgroup$
  • $\begingroup$ What's your goal for the visualization? Is this for yourself for checking assumptions / fit, or for presentation to others (eg, put in manuscript)? $\endgroup$ – gung Aug 5 at 18:06
  • $\begingroup$ @gung - I've checked residuals with DHARMa. Now i'm just trying to visualize results for a manuscript $\endgroup$ – dwash7 Aug 5 at 22:36
3
$\begingroup$

One good way to visualize the results of mixed models is via effect plots. These display the form and magnitude of the association between the expected outcome and some of the predictions in your model while keeping the rest of them at fixed values.

In R these are provided via, e.g., the effects package. You could fit the negative binomial mixed model with the adaptive Gaussian quadrature, which in general is considered to be better than the Laplace approximation using the GLMMadaptive package that I’ve written. Examples of effects plots with this package can be found here and here.

$\endgroup$
  • 2
    $\begingroup$ +1 for a great suggestion! However, you should probably mention somewhere in the answer that you are the author of the package you recommend. $\endgroup$ – Frans Rodenburg Aug 6 at 8:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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