Extract heterogeneous variances for nlme, glmmTMB, and brms models.

extract_het_var(model, digits = 3, scale = "var", ...)

# S3 method for lme
extract_het_var(model, digits = 3, scale = "var", ...)

# S3 method for glmmTMB
extract_het_var(model, digits = 3, scale = "var", ...)

# S3 method for brmsfit
extract_het_var(
  model,
  digits = 3,
  scale = "var",
  ...,
  ci_level = 0.95,
  return_all = FALSE
)

Arguments

model

An appropriate mixed model.

digits

Rounding. Default is 3.

scale

Return result on original standard deviation scale ('sd') or as variance ('var'), the default.

...

Other arguments specific to the method. Unused at present.

ci_level

For brms objects, confidence level < 1, typically above 0.90. A value of 0 will not report it. Default is .95.

return_all

For brms class objects, return all fitted values (TRUE) or only the distinct ones. Default is FALSE.

Value

A vector of the estimates on the variance scale

Details

For nlme models with heterogeneous variance, i.e. that contain something like varIdent(form = ~1|Group), this returns a more presentable version the estimates. Only tested with the varIdent case.

For glmmTMB, this serves as a wrapper for extract_cor_structure for models with for the diag function as part of the formula. See that function for details. For distributional models where the dispersion is explicitly modeled separately via disp = ~ , use the component argument of the other functions in this package.

For brms distributional models with a sigma ~ . formula, this produces the (unique) fitted values for the dispersion part of the model. As this is often just a single grouping variable to allow variance to vary over the group levels, only the distinct fitted values, which would be one value per group, are returned. If all fitted values are desired, set return_all to TRUE.

This function has not been tested except in the more simple model settings. It's unclear how well it will work with other model complications added.

See also

Examples

#> #> Attaching package: ‘nlme’
#> The following object is masked from ‘package:lme4’: #> #> lmList
library(mixedup) model <- lme( distance ~ age + Sex, data = Orthodont, random = ~ 1|Subject, weights = varIdent(form = ~ 1 | Sex) ) summary(model)
#> Linear mixed-effects model fit by REML #> Data: Orthodont #> AIC BIC logLik #> 432.3567 448.2805 -210.1784 #> #> Random effects: #> Formula: ~1 | Subject #> (Intercept) Residual #> StdDev: 1.839463 1.761758 #> #> Variance function: #> Structure: Different standard deviations per stratum #> Formula: ~1 | Sex #> Parameter estimates: #> Male Female #> 1.0000000 0.4541324 #> Fixed effects: distance ~ age + Sex #> Value Std.Error DF t-value p-value #> (Intercept) 18.919992 0.7285568 80 25.969138 0.0000 #> age 0.549887 0.0473096 80 11.623168 0.0000 #> SexFemale -2.321023 0.7629703 25 -3.042088 0.0055 #> Correlation: #> (Intr) age #> age -0.714 #> SexFemale -0.468 0.000 #> #> Standardized Within-Group Residuals: #> Min Q1 Med Q3 Max #> -3.25494115 -0.53013324 -0.02554914 0.51114273 3.03915818 #> #> Number of Observations: 108 #> Number of Groups: 27
extract_het_var(model)
#> Male Female #> 1 3.104 0.64
library(glmmTMB) # does not get the same estimates as nlme, but would get similar if modeled # using dispersion approach. model <- glmmTMB(distance ~ age + Sex + (1 | Subject) + diag(Sex + 0 | Subject), data = Orthodont)
#> Warning: 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
#> Warning: 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
#> Warning: 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
#> Warning: Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
extract_het_var(model)
#> group SexMale SexFemale #> 1 Subject.1 2.072 2.964
# compare with model <- glmmTMB(distance ~ age + Sex + (1 | Subject), dispformula = ~ Sex, data = Orthodont)
#> Warning: 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
#> Warning: 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
#> Warning: 'giveCsparse' has been deprecated; setting 'repr = "T"' for you
extract_fixed_effects(model, component = 'disp', exponentiate = TRUE)
#> Warning: Intervals could not be computed for fixed effects. Returning CI based on SE. #> If Fixed Effect SE is NaN, check random effects for zero variance estimates.
#> # A tibble: 2 x 7 #> term value se z p_value lower_2.5 upper_97.5 #> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> #> 1 Intercept 3.09 0.639 5.45 0 2.06 4.63 #> 2 SexFemale 0.202 0.069 -4.68 0 0.103 0.395
if (FALSE) { library(brms) model <- brm(bf(distance ~ age + Sex + (1 | Subject), sigma ~ Sex), data = Orthodont) extract_het_var(model) }