1
$\begingroup$

The goal

I want to run dominance analysis on a mixed-effects beta model, to approximate the relative importance of a set of predictors (2 factors, 1 scaled continuous, 1 continuous with splines). The random part of the beta model would be random intercepts, with a nested structure: (1|random1/random2).

The mixed-effects beta would be fitted with glmmTMB::glmmTMB(), the metric of model quality would be marginal R2 (MuMIn:.r.squaredGLMM()[1]), and dominance analysis would be with the domir::domin() package.

The obstacle

I get errors when trying to run domir::domin() on the mixed-effects beta model, with two different errors messages, one without nested effects: Error: Invalid grouping factor specification; another when there is nesting: “Error: unparseable formula for grouping factor”. A model without random effects works fine. I get warnings from MuMIn about small mu and that the effects of zero-inflation and dispersion were ignored, but I don't think these are the real obstacle.

A reproducible example

Building on a subset of the real data throwing the errors, I try to provide a reproducible example. A question on how to run dominance analyses with a mixed-effects model has been asked before, but I haven’t found a solution there. I hope this catches the attention of Josep Luchman, author of the domir package, who I have seen to be very active in this forum :) Many thanks in advance, and thank you for the amazing package!

Cheers!

#### Load libraries

library(domir)
#> Warning: package 'domir' was built under R version 4.2.3
library(glmmTMB)
library(MuMIn)
#> Warning: package 'MuMIn' was built under R version 4.2.3
library(splines)


##### This dataset is an "anonymised" subset of my data, which looks ok to reproduce the errors
df <- structure(list(response = c(0.21, 0.151, 0.256, 0.011, 0.26, 
                                  0.176, 0.087, 0.064, 0.123, 0.312, 0.148, 0.001, 0.099, 0.265, 
                                  0.222, 0.07, 0.171, 0.108, 0.115, 0.119, 0.2, 0.01, 0.092, 0.042, 
                                  0.067, 0.389, 0.035, 0.046, 0.041, 0.143, 0.042, 0.154, 0.351, 
                                  0.047, 0.068, 0.1, 0.062, 0.009, 0.224, 0.114, 0.063, 0.399, 
                                  0.026, 0.047, 0.07, 0.205, 0.31, 0.204, 0.092, 0.133, 0.209, 
                                  0.012, 0.022, 0.207, 0.233, 0.103, 0.072, 0.344, 0.063, 0.216, 
                                  0.302, 0.43, 0.046, 0.205, 0.082, 0.007, 0.027, 0.022, 0.036, 
                                  0.095, 0.192, 0.146, 0.184, 0.072, 0.122, 0.029, 0.043, 0.032, 
                                  0.312, 0.058), cont1 = c(121, 261, 140, 162, 47, 118, 47, 74, 
                                                           184, 162, 141, 184, 37, 184, 123, 146, 43, 220, 176, 107, 171, 
                                                           82, 171, 49, 162, 151, 40, 170, 111, 47, 61, 74, 177, 130, 147, 
                                                           125, 28, 35, 158, 51, 184, 88, 110, 184, 68, 78, 62, 46, 138, 
                                                           162, 116, 51, 127, 199, 51, 52, 114, 91, 184, 244, 243, 115, 
                                                           72, 117, 151, 39, 56, 30, 70, 41, 79, 171, 126, 114, 34, 184, 
                                                           62, 82, 60, 34), cont2 = c(0.52, 0.42, 0.301, 0.515, 0.271, 0.517, 
                                                                                      0.319, 0.375, 0.446, 0.406, 0.53, 0.581, 0.245, 0.313, 0.295, 
                                                                                      0.21, 0.35, 0.501, 0.434, 0.495, 0.312, 0.437, 0.473, 0.563, 
                                                                                      0.442, 0.268, 0.546, 0.543, 0.579, 0.482, 0.378, 0.521, 0.259, 
                                                                                      0.586, 0.555, 0.474, 0.404, 0.587, 0.392, 0.389, 0.423, 0.31, 
                                                                                      0.533, 0.475, 0.41, 0.329, 0.389, 0.45, 0.502, 0.41, 0.428, 0.563, 
                                                                                      0.56, 0.542, 0.498, 0.554, 0.558, 0.331, 0.423, 0.406, 0.319, 
                                                                                      0.257, 0.562, 0.453, 0.598, 0.436, 0.468, 0.49, 0.254, 0.493, 
                                                                                      0.495, 0.346, 0.364, 0.558, 0.29, 0.531, 0.382, 0.402, 0.389, 
                                                                                      0.42), factor1 = structure(c(2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 
                                                                                                                   2L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 
                                                                                                                   1L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 
                                                                                                                   2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 
                                                                                                                   1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
                                                                                                                   2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L), levels = c("1", "2"), class = "factor"), 
                     factor2 = structure(c(2L, 3L, 3L, 2L, 3L, 1L, 3L, 2L, 1L, 
                                           2L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 2L, 1L, 1L, 1L, 2L, 3L, 1L, 
                                           2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 
                                           3L, 1L, 1L, 1L, 2L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 
                                           1L, 2L, 2L, 2L, 1L, 2L, 3L, 3L, 3L, 1L, 1L, 3L, 2L, 3L, 2L, 
                                           1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 2L, 2L), levels = c("1", 
                                                                                                   "2", "3"), class = "factor"), rand1 = structure(c(33L, 19L, 
                                                                                                                                                     2L, 24L, 11L, 24L, 16L, 15L, 13L, 7L, 35L, 24L, 7L, 5L, 4L, 
                                                                                                                                                     31L, 29L, 7L, 7L, 8L, 18L, 12L, 13L, 8L, 27L, 8L, 1L, 24L, 
                                                                                                                                                     9L, 6L, 22L, 24L, 17L, 17L, 7L, 10L, 34L, 21L, 8L, 23L, 13L, 
                                                                                                                                                     8L, 9L, 30L, 25L, 29L, 13L, 8L, 26L, 24L, 24L, 35L, 7L, 24L, 
                                                                                                                                                     8L, 8L, 33L, 8L, 13L, 24L, 32L, 14L, 20L, 7L, 8L, 2L, 24L, 
                                                                                                                                                     29L, 24L, 24L, 8L, 8L, 24L, 33L, 8L, 28L, 3L, 8L, 8L, 24L
                                                                                                   ), levels = c("01", "02", "03", "04", "05", "06", "07", "08", 
                                                                                                                 "09", "10", "11", "12", "13", "14", "15", "16", "17", "18", 
                                                                                                                 "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", 
                                                                                                                 "29", "30", "31", "32", "33", "34", "35"), class = "factor"), 
                     rand2 = structure(c(7L, 38L, 28L, 59L, 51L, 61L, 10L, 53L, 
                                         52L, 62L, 37L, 58L, 42L, 12L, 47L, 41L, 35L, 24L, 15L, 9L, 
                                         55L, 8L, 1L, 21L, 6L, 57L, 20L, 3L, 5L, 18L, 19L, 50L, 48L, 
                                         39L, 15L, 43L, 23L, 11L, 29L, 31L, 52L, 60L, 5L, 34L, 45L, 
                                         35L, 25L, 40L, 56L, 16L, 59L, 37L, 15L, 59L, 21L, 54L, 7L, 
                                         40L, 52L, 33L, 46L, 4L, 26L, 15L, 2L, 30L, 44L, 14L, 3L, 
                                         61L, 60L, 27L, 36L, 7L, 32L, 13L, 49L, 17L, 60L, 22L), levels = c("01", 
                                                                                                           "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", 
                                                                                                           "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", 
                                                                                                           "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", 
                                                                                                           "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", 
                                                                                                           "42", "43", "44", "45", "46", "47", "48", "49", "50", "51", 
                                                                                                           "52", "53", "54", "55", "56", "57", "58", "59", "60", "61", 
                                                                                                           "62"), class = "factor")), class = "data.frame", row.names = c(NA, 
                                                                                                                                                                          -80L))
str(df)
#> 'data.frame':    80 obs. of  7 variables:
#>  $ response: num  0.21 0.151 0.256 0.011 0.26 0.176 0.087 0.064 0.123 0.312 ...
#>  $ cont1   : num  121 261 140 162 47 118 47 74 184 162 ...
#>  $ cont2   : num  0.52 0.42 0.301 0.515 0.271 0.517 0.319 0.375 0.446 0.406 ...
#>  $ factor1 : Factor w/ 2 levels "1","2": 2 2 2 1 2 1 2 2 2 1 ...
#>  $ factor2 : Factor w/ 3 levels "1","2","3": 2 3 3 2 3 1 3 2 1 2 ...
#>  $ rand1   : Factor w/ 35 levels "01","02","03",..: 33 19 2 24 11 24 16 15 13 7 ...
#>  $ rand2   : Factor w/ 62 levels "01","02","03",..: 7 38 28 59 51 61 10 53 52 62 ...

#### Test model-fitting functions and metric of model quality

mod <- glmmTMB::glmmTMB(response ~ scale(cont2) + bs(cont1, df = 3) + factor1 + factor2 + (1|rand1/rand2), dispformula = ~ 1, data = df, family = beta_family(link = "logit"))
summary(mod)
#>  Family: beta  ( logit )
#> Formula:          
#> response ~ scale(cont2) + bs(cont1, df = 3) + factor1 + factor2 +  
#>     (1 | rand1/rand2)
#> Data: df
#> 
#>      AIC      BIC   logLik deviance df.resid 
#>   -184.6   -158.4    103.3   -206.6       69 
#> 
#> Random effects:
#> 
#> Conditional model:
#>  Groups      Name        Variance Std.Dev.
#>  rand2:rand1 (Intercept) 0.13042  0.3611  
#>  rand1       (Intercept) 0.02887  0.1699  
#> Number of obs: 80, groups:  rand2:rand1, 62; rand1, 35
#> 
#> Dispersion parameter for beta family (): 22.1 
#> 
#> Conditional model:
#>                    Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)        -2.73257    0.36252  -7.538 4.78e-14 ***
#> scale(cont2)       -0.50095    0.08824  -5.677 1.37e-08 ***
#> bs(cont1, df = 3)1  2.02663    0.90326   2.244   0.0249 *  
#> bs(cont1, df = 3)2 -0.25883    0.67948  -0.381   0.7033    
#> bs(cont1, df = 3)3  1.53706    0.60982   2.521   0.0117 *  
#> factor12           -0.14068    0.22337  -0.630   0.5288    
#> factor22            0.06847    0.24910   0.275   0.7834    
#> factor23           -0.08611    0.33131  -0.260   0.7949    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
MuMIn::r.squaredGLMM(mod)
#> Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
#> Warning in r.squaredGLMM.glmmTMB(mod): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#>            R2m       R2c
#> [1,] 0.4398609 0.6686709


##### domir::domin

# Without any random effects
domir::domin(response ~ scale(cont2) + bs(cont1, df = 3) + factor1 + factor2,
             reg =  function(y)  glmmTMB::glmmTMB(formula = y, dispformula = ~ 1, data = df, family = beta_family(link = "logit")),
             fitstat = list(\(x) list(R2m = MuMIn::r.squaredGLMM(x)[[1]]), "R2m"))
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.2 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.2 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.2 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.2 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.2 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.2 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.2 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.2 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.2 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.2 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.2 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.2 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.2 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.2 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.2 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Overall Fit Statistic:      0.157678 
#> 
#> General Dominance Statistics:
#>                   General Dominance Standardized Ranks
#> scale(cont2)           0.1187231767  0.752947071     1
#> bs(cont1, df = 3)      0.0367193338  0.232875464     2
#> factor1                0.0014161674  0.008981390     3
#> factor2                0.0008193067  0.005196075     4
#> 
#> Conditional Dominance Statistics:
#>                         IVs: 1      IVs: 2       IVs: 3       IVs: 4
#> scale(cont2)      1.165115e-01 0.120016192 0.1196132751  0.118751758
#> bs(cont1, df = 3) 3.379900e-02 0.037772258 0.0378437056  0.037462371
#> factor1           6.118088e-06 0.002503043 0.0018023728  0.001353136
#> factor2           1.239085e-03 0.002567201 0.0006211075 -0.001150167
#> 
#> Complete Dominance Designations:
#>                           Dmnated?scale(cont2) Dmnated?bs(cont1, df = 3)
#> Dmnates?scale(cont2)                        NA                      TRUE
#> Dmnates?bs(cont1, df = 3)                FALSE                        NA
#> Dmnates?factor1                          FALSE                     FALSE
#> Dmnates?factor2                          FALSE                     FALSE
#>                           Dmnated?factor1 Dmnated?factor2
#> Dmnates?scale(cont2)                 TRUE            TRUE
#> Dmnates?bs(cont1, df = 3)            TRUE            TRUE
#> Dmnates?factor1                        NA              NA
#> Dmnates?factor2                        NA              NA

# With random effects, without nesting
domir::domin(response ~ scale(cont2) + bs(cont1, df = 3) + factor1 + factor2 + (1|rand1),
             reg =  function(y)  glmmTMB::glmmTMB(formula = y, dispformula = ~ 1, data = df, family = beta_family(link = "logit")),
             fitstat = list(\(x) list(R2m = MuMIn::r.squaredGLMM(x)[[1]]), "R2m"),
             consmodel = "(1|rand1)")
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in Ops.factor(1, rand1): '|' not meaningful for factors
#> Error: Invalid grouping factor specification, 1 | rand1

# With random effects, with nesting
domir::domin(response ~ scale(cont2) + bs(cont1, df = 3) + factor1 + factor2 + (1|rand1/rand2),
             reg =  function(y)  glmmTMB::glmmTMB(formula = y, dispformula = ~ 1, data = df, family = beta_family(link = "logit")),
             fitstat = list(\(x) list(R2m = MuMIn::r.squaredGLMM(x)[[1]]), "R2m"),
             consmodel = "(1|rand1/rand2)")
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
#> dispersion model are ignored
#> Warning: mu of 0.1 is too close to zero, estimate of random effect variances may
#>   be unreliable.
#> Error: unparseable formula for grouping factor

#### Session info
sessionInfo()
#> R version 4.2.0 (2022-04-22 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 22000)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=Spanish_Spain.utf8  LC_CTYPE=Spanish_Spain.utf8   
#> [3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C                  
#> [5] LC_TIME=Spanish_Spain.utf8    
#> 
#> attached base packages:
#> [1] splines   stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#> [1] MuMIn_1.47.5  glmmTMB_1.1.3 domir_1.0.1  
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.10           compiler_4.2.0        nloptr_2.0.1         
#>  [4] highr_0.9             TMB_1.8.1             tools_4.2.0          
#>  [7] boot_1.3-28           digest_0.6.29         lme4_1.1-29          
#> [10] evaluate_0.15         lifecycle_1.0.3       nlme_3.1-157         
#> [13] lattice_0.20-45       rlang_1.1.1           reprex_2.0.2         
#> [16] Matrix_1.4-1          cli_3.6.1             rstudioapi_0.13      
#> [19] yaml_2.3.5            mvtnorm_1.1-3         xfun_0.31            
#> [22] fastmap_1.1.0         coda_0.19-4           withr_2.5.0          
#> [25] stringr_1.5.0         knitr_1.39            fs_1.5.2             
#> [28] vctrs_0.6.3           stats4_4.2.0          grid_4.2.0           
#> [31] glue_1.6.2            survival_3.3-1        rmarkdown_2.14       
#> [34] multcomp_1.4-19       TH.data_1.1-1         minqa_1.2.4          
#> [37] magrittr_2.0.3        codetools_0.2-18      htmltools_0.5.2      
#> [40] emmeans_1.8.5-9000004 MASS_7.3-57           insight_0.19.1.3     
#> [43] xtable_1.8-4          numDeriv_2016.8-1.1   sandwich_3.0-1       
#> [46] stringi_1.7.6         estimability_1.4.1    zoo_1.8-10

Created on 2023-07-25 with reprex v2.0.2

EDIT (2 August 2023)

Many thanks to Joseph Luchman! Dominance analysis with a beta regression with random effects implemented with domir::domin() yields the same results as domir::domir(). See below the results using domir::domir() (reprex won't render it for some reason).

domir::domir(response ~ scale(cont2) + bs(cont1, df = 3) + factor1 + factor2,
         function(y)  {
           glmmTMB::glmmTMB(formula = update(y, . ~ . +(1|rand1/rand2)), dispformula = ~ 1, data = df, family = beta_family(link = "logit")) %>% MuMIn::r.squaredGLMM() %>% .[[1]]
           })

dominance analysis with domir::domir()

$\endgroup$

1 Answer 1

1
$\begingroup$

Riera,

A fairly easy fix. Take the (1|rand1/rand2) out of the formula and keep it only in consmodel argument like below:

> domir::domin(response ~ scale(cont2) + bs(cont1, df = 3) + factor1 + factor2,
+              reg =  function(y)  glmmTMB::glmmTMB(formula = y, dispformula = ~ 1, data = df, family = beta_family(link = "logit")),
+              fitstat = list(\(x) list(R2m = MuMIn::r.squaredGLMM(x)[[1]]), "R2m"),
+              consmodel = "(1|rand1/rand2)")
Overall Fit Statistic:      0.4398609 
Constant Model Fit Statistic:  0 

General Dominance Statistics:
                  General Dominance Standardized Ranks
scale(cont2)            0.321281139   0.73041529     1
bs(cont1, df = 3)       0.104204068   0.23690231     2
factor1                 0.005243451   0.01192070     4
factor2                 0.009132256   0.02076169     3

Conditional Dominance Statistics:
                        IVs: 1     IVs: 2      IVs: 3        IVs: 4
scale(cont2)      0.3433449821 0.33852835 0.315293963  0.2879572616
bs(cont1, df = 3) 0.1154408659 0.11804687 0.101928211  0.0814003273
factor1           0.0001384045 0.00932296 0.006256870  0.0052555700
factor2           0.0102260116 0.01748997 0.009185377 -0.0003723375

Complete Dominance Designations:
                          Dmnated?scale(cont2) Dmnated?bs(cont1, df = 3) Dmnated?factor1 Dmnated?factor2
Dmnates?scale(cont2)                        NA                      TRUE            TRUE            TRUE
Dmnates?bs(cont1, df = 3)                FALSE                        NA            TRUE            TRUE
Dmnates?factor1                          FALSE                     FALSE              NA              NA
Dmnates?factor2                          FALSE                     FALSE              NA              NA

There were 32 warnings (use warnings() to see them)

The domir::domin function was adding the (1|rand1/rand2) term to the formula which resulted in that term being included twice. It always adds the contents of consmodel to the formula argument when it estimates each sub-model and so, by including it in consmodel it will be included in each sub-model in the dominance analysis.

$\endgroup$
4
  • $\begingroup$ Thank you! This fixed the error. I want to run dominance analysis on another type of model (Dirichlet regression), but I've come across new issues. If you have time, could you like into it? Thank you very much! stats.stackexchange.com/questions/622422/… $\endgroup$
    – M. Riera
    Commented Jul 27, 2023 at 9:35
  • $\begingroup$ I would like to switch to domir(). How could I get a similar procedure to that performed by the "consmodel" argument from domin()? I've tried: .all = ~ (1|family/genus) and .adj = ~ (1|family/genus), but in both cases I get the error: "Error in eigen(h) : infinite or missing values in 'x'". Should I post this in a new question, over at Stack Overflow? $\endgroup$
    – M. Riera
    Commented Aug 1, 2023 at 8:53
  • 1
    $\begingroup$ .all and .adj parse using stats::terms() which will strip the parentheses from that term which causes it to fail as it's read by glmmTMB as 1|family/genus without parentheses. The fix is to add it using an update() to the formula in the function itself. If you can't get it to work recommend asking as a separate question and I can expand on this answer. $\endgroup$
    – jluchman
    Commented Aug 1, 2023 at 16:40
  • $\begingroup$ update() worked! I've updated my question to show the results, which are the same as those obtained with domin(). Thank you! $\endgroup$
    – M. Riera
    Commented Aug 2, 2023 at 7:11

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.