m-clark / mixedup

An R package for extracting results from mixed models that are easy to use and viable for presentation.
http://m-clark.github.io/mixedup
MIT License
67 stars 2 forks source link

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

Closed SchmidtPaul closed 2 years ago

SchmidtPaul commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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

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