`extract_vc()` with `weights = varIdent()` in `lme()`

[SchmidtPaul]

Hi, thanks for my probably favourite non-CRAN-package.

I guess you may not be planning to make this possible, since extract_vc() "has functionality for simpler models", but I noticed that when fitting a heterogeneous error variance (i.e. diagonal variance structure on the R-side of a model) with lme(), I will not receive separate estimates for the error variance when using extract_vc(). Moreover, it seems like I get only one (the first?) of the multiple error variances.

Below is a reprex where I tried to extract the relevant info myself.

library(agridat)
library(glmmTMB)
library(mixedup)
library(nlme)
library(tidyverse)


# data --------------------------------------------------------------------
dat <- agridat::mcconway.turnip %>%
  mutate(unit = 1:n()) %>%
  mutate_at(vars(density, unit), as.factor)

# mod:lme -----------------------------------------------------------------
diag_lme <- lme(
  yield ~
    gen * date * density,
  random  = ~ 1 | block,
  weights = varIdent(form =  ~ 1 | date),
  data    = dat
)

# extract VC --------------------------------------------------------------
extract_vc(diag_lme)
#>          group    effect variance    sd sd_2.5 sd_97.5 var_prop
#> block    block Intercept    1.597 1.264  0.455   3.511     0.27
#> 1     Residual              4.306 2.075  1.548   2.782     0.73

# it should do something like this:
lme_Gside <- extract_vc(diag_lme) %>%
  filter(group != "Residual") %>% 
  as_tibble()

lme_Rside <- diag_lme$modelStruct$varStruct %>%
    coef(unconstrained = FALSE, allCoef = TRUE) %>%
    enframe(name = "group", value = "varStruct") %>%
    mutate(sigma         = diag_lme$sigma) %>%
    mutate(StandardError = sigma * varStruct) %>%
    mutate(variance      = StandardError ^ 2) %>%
    mutate(effect = "Residual") %>%
    select(effect, group, variance)

bind_rows(lme_Gside, lme_Rside)
#> # A tibble: 3 x 7
#>   group     effect    variance    sd sd_2.5 sd_97.5 var_prop
#>   <chr>     <chr>        <dbl> <dbl>  <dbl>   <dbl>    <dbl>
#> 1 block     Intercept     1.60  1.26  0.455    3.51     0.27
#> 2 21Aug1990 Residual      4.31 NA    NA       NA       NA   
#> 3 28Aug1990 Residual     15.5  NA    NA       NA       NA

^{Created on 2022-01-09 by the reprex package (v2.0.1)}

[SchmidtPaul]

My bad, I just realized I forgot about extract_het_var(). Nevertheless, my follow-up question is then: you are not planning on combining it with extract_vc() to obtain the table more or less like so?

lme_Gside <- extract_vc(diag_lme) %>%
  filter(group != "Residual") %>% 
  as_tibble()

lme_Rside <- extract_het_var(diag_lme) %>%
  pivot_longer(cols = everything(),
               names_to = "effect",
               values_to = "variance") %>% 
  mutate(group = "Residual")

bind_rows(lme_Gside, lme_Rside)
#> # A tibble: 3 x 7
#>   group    effect     variance    sd sd_2.5 sd_97.5 var_prop
#>   <chr>    <chr>         <dbl> <dbl>  <dbl>   <dbl>    <dbl>
#> 1 block    Intercept      1.60  1.26  0.455    3.51     0.27
#> 2 Residual X21Aug1990     4.31 NA    NA       NA       NA   
#> 3 Residual X28Aug1990    15.5  NA    NA       NA       NA

[m-clark]

Not at this time, but I am going to do some minor updates to this package here and there (hopefully soon), and will leave this open to consider it further. Seems feasible for those simple settings for which it's intended.

[m-clark]

@SchmidtPaul Much belatedly looking into this further, it looks like we have the following:

• nlme: this could maybe be done as you describe.
• glmmTMB: output already does this for the diag structure when extract_vc is called due to the way tmb works, though the residual variance isn't added to the group variances as is done in extract_het_var.
• brms: I think these aren't easily combinable, since one can model sigma without any random effects, leading to an error if extract_vc is called on such an object but which works fine with extract_het_var. Also, brms allows for a full modeling of sigma, which basically results in extracting fitted values for the sigma model.

Since they all potentially produce notably different output even for the same type of model, for now I've decided to add an include_het_var argument to extract_vc, which will return a list of extract_vc output and extract_het_var output if TRUE, with a default of FALSE. In addition, I thought it better to return a more 'tidy' result, so this should make it easy enough to just bind_rows on the nlme result as below.

> extract_vc(lme_het_var, scale = 'var', include_het_var = TRUE) %>% bind_rows()
# A tibble: 4 × 7
  group    effect    variance    sd sd_2.5 sd_97.5 var_prop
  <chr>    <chr>        <dbl> <dbl>  <dbl>   <dbl>    <dbl>
1 Subject  Intercept     3.38  1.84   1.35    2.51    0.522
2 Residual NA            3.10  1.76   1.44    2.16    0.478
3 Male     NA            3.10 NA     NA      NA      NA    
4 Female   NA            0.64 NA     NA      NA      NA