1
$\begingroup$

I am looking to fit a glmm for a project dealing with the presence and absence of a species of fish. I have been reading up on mixed modeling with hierarchical nesting and thought I had a decent grasp on what I was doing but it seems as though I was wrong. (I am not an expert at statistical by any means)

My DV is the presence/absence of the fish (0/1) and my predictor variables are dissolved oxygen (DO), percent cover (per_cover), percent pebbles between 6-11mm (per_pebble_611), and number of impoundments in the sub-watershed (Num_imp).

I have attempted to fit the model using several different versions of code but I have issues with singular fit, non-convergence, and questions on the random term notation. I have checked for issues with multicollinearity and no variables seems to be an issue.

The samples are stream sites nested within sub-watersheds (HUC12) nested within larger watersheds (HUC4).

The first formula I tried I included the random effect (1|HUC4_f/HUC12_f/SiteID). I have not been able to find any information about using three tiered nested random effects or how to notate them in the code but I attempted this anyway. My model is as follows:

glmm <- glmer(Pres ~ DO + per_cover + per_pebble_611 + Num_imp + (1|HUC4_f/HUC12_f/SiteID), family = "binomial")
Random effects:
 Groups                  Name        Variance  Std.Dev. 
 SiteID:(HUC12_f:HUC4_f) (Intercept) 1.018e-09 3.191e-05
 HUC12_f:HUC4_f          (Intercept) 2.203e+00 1.484e+00
 HUC4_f                  (Intercept) 1.030e-09 3.209e-05
Number of obs: 115, groups:  SiteID:(HUC12_f:HUC4_f), 115; HUC12_f:HUC4_f, 42; HUC4_f, 2

Fixed effects:
               Estimate Std. Error z value Pr(>|z|)   
(Intercept)    -7.58519    2.31424  -3.278  0.00105 **
DO              0.59430    0.23873   2.489  0.01280 * 
per_cover       2.94201    1.41749   2.076  0.03794 * 
per_pebble_611  0.08019    0.02937   2.731  0.00632 **
Num_imp        -0.01747    0.01143  -1.528  0.12646   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) DO     pr_cvr p__611
DO          -0.870                     
per_cover   -0.401  0.032              
pr_pbbl_611 -0.477  0.362  0.078       
Num_imp     -0.027 -0.174 -0.013 -0.231
convergence code: 0
boundary (singular) fit: see ?isSingular

This code results in singular fit of the model and I am also not sure if the random effect term specified is even a valid notation. In an attempt to simplify the model and use a notation similar to others I have found online I took out the SiteID term and used this code.

glmm <- glmer(Pres ~ DO + per_cover + per_pebble_611 + Num_imp + (1|HUC4_f/HUC12_f), family = "binomial")
Random effects:
 Groups         Name        Variance  Std.Dev.
 HUC12_f:HUC4_f (Intercept) 2.202e+00 1.483955
 HUC4_f         (Intercept) 1.666e-05 0.004081
Number of obs: 115, groups:  HUC12_f:HUC4_f, 42; HUC4_f, 2

Fixed effects:
               Estimate Std. Error z value Pr(>|z|)   
(Intercept)    -7.58463    2.31402  -3.278  0.00105 **
DO              0.59427    0.23871   2.489  0.01279 * 
per_cover       2.94168    1.41739   2.075  0.03795 * 
per_pebble_611  0.08019    0.02936   2.731  0.00632 **
Num_imp        -0.01747    0.01143  -1.528  0.12642   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) DO     pr_cvr p__611
DO          -0.870                     
per_cover   -0.401  0.032              
pr_pbbl_611 -0.477  0.362  0.078       
Num_imp     -0.027 -0.174 -0.013 -0.231
convergence code: 0
Model failed to converge with max|grad| = 0.0156959 (tol = 0.002, component 1)

However, when running this model I get errors saying the model failed to converge. I have seen other suggestions on similar posts that recomend setting the nAGQ argument equal to 0. When I do this the model does converge however I am not sure what this is actually doing to my model. So I then tried removing the largest level of nesting instead and used the code below.

glmm <- glmer(Pres ~ DO + per_cover + per_pebble_611 + Num_imp + (1|HUC12_f/SiteID), family = "binomial")
Random effects:
 Groups         Name        Variance  Std.Dev.
 SiteID:HUC12_f (Intercept) 1.535e-06 0.001239
 HUC12_f        (Intercept) 2.203e+00 1.484267
Number of obs: 115, groups:  SiteID:HUC12_f, 115; HUC12_f, 42

Fixed effects:
               Estimate Std. Error z value Pr(>|z|)   
(Intercept)    -7.58510    2.31425  -3.278  0.00105 **
DO              0.59429    0.23873   2.489  0.01280 * 
per_cover       2.94201    1.41750   2.075  0.03794 * 
per_pebble_611  0.08019    0.02937   2.731  0.00632 **
Num_imp        -0.01747    0.01143  -1.528  0.12646   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) DO     pr_cvr p__611
DO          -0.870                     
per_cover   -0.401  0.032              
pr_pbbl_611 -0.477  0.362  0.078       
Num_imp     -0.027 -0.174 -0.013 -0.231

This resulted in no errors with singular fit or non-convergence but I am unsure what is the correct formula/notation to run the analysis I am interested in interpreting.

Any insight would be greatly appreciated. Also I tried to include as much information as possible without making this post too lengthy but please let me know if any additional information would make giving advice easier/possible and I'll be happy to add it in.

$\endgroup$

1 Answer 1

2
$\begingroup$

Garrett, can you please clarify how many values of your response variable Pres you have available for each site identified by a SiteID? Do you have (i) multiple values of Pres or (ii) a single value of Pres for that site?

The summary for your most complex model includes this line:

Number of obs: 115, groups:  SiteID:(HUC12_f:HUC4_f), 115; HUC12_f:HUC4_f, 42; HUC4_f, 2

This seems to imply that you only have one value of Pres available per Site?

Because this is a mixed effects model with nested random grouping factors, you need to have multiple values of Pres for at least some of your sites. If this is not the case (i.e., all of your sites have exactly one value of Pres recorded), then you cannot treat site as a random grouping factor in your mixed effects model. You can however still include HUC4_f and HUC12_f as random grouoing factors in your model:

glmm <- glmer(Pres ~ DO + per_cover + per_pebble_611 + Num_imp +      
                     (1|HUC4_f/HUC12_f), 
                      family = "binomial")

A ruld of thumb states that a random grouping factor should have at least 5 levels/categories for it to be included in the model.

$\endgroup$
5
  • 1
    $\begingroup$ Each site has a single value for presence, either present (1) or absent (0). So then I would need to use the grouping factor (1|HUC4_f/HUC12_f)? $\endgroup$ Feb 20, 2021 at 3:37
  • $\begingroup$ That’s what I thought! You are on the right track with your comment - can you try fitting that model and posting the output in your original question? $\endgroup$ Feb 20, 2021 at 3:40
  • 1
    $\begingroup$ I actually included that in the original post as the second model and output listed. However, when I run that model I run into issues with the model failing to converge. But now that I am reading over your answer again, I think I would need to eliminate the HUC4 level in the grouping term as well because there are only 2 categories within that group. So that being said I guess I would just use (1|HUC12_f) as my random effect? $\endgroup$ Feb 20, 2021 at 3:51
  • $\begingroup$ Given that you only have 2 categories for HUC4, it makes sense that you would not include HUC4 as a random grouping factor in your model. Just as you suggested, include HUC12_f as the only random grouping factor. Then your data hierarchy will have two levels. The top level will consist of the various categories of HUC12_f. The bottom level will consist of the response observations - nested within each category of HUC12_f. You can include HUC4 as a predictor in your model if need be. $\endgroup$ Feb 20, 2021 at 4:11
  • $\begingroup$ With two levels for your data hierarchy, predictors can refer to either level of the hierarchy. For example, the HUC4 predictor (if you decide to include it) will refer to the top level of the hierarchy. Any predictor measured for each site will refer to the bottom level of the hierarchy. $\endgroup$ Feb 20, 2021 at 4:14

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.