Extract heterogeneous variances

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.

Examples

library(nlme)
#> 
#> 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)
#> # A tibble: 2 × 2
#>   group  variance
#>   <chr>     <dbl>
#> 1 Male       3.10
#> 2 Female     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: Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')

extract_het_var(model)
#> # A tibble: 2 × 3
#>   group     name      variance
#>   <chr>     <chr>        <dbl>
#> 1 Subject.1 SexMale       2.07
#> 2 Subject.1 SexFemale     2.96

# compare with
model <-
  glmmTMB(distance ~ age + Sex + (1 | Subject), dispformula = ~ Sex,
          data = Orthodont)

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 × 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)
}

Arguments

Value

Details

See also

Examples