bbolker / broom.mixed

tidy methods for mixed models in R
227 stars 23 forks source link

VarStruct information from `lme` tidy missing? #96

Open SchmidtPaul opened 4 years ago

SchmidtPaul commented 4 years ago

Hi there, great work so far! Maybe I am missing something here, but I get the impression that when fitting a model with heteroscedascity in the error term via nlme, the tidy(effects="ran_pars") does not extract all the information as it only provides the sigma estimate.

library(tidyverse)
library(nlme)

dat <- agridat::mcconway.turnip %>%
  mutate(densf = density %>% as.factor)

mod <- nlme::lme(fixed   = yield ~ gen * date * densf, 
                 random  = ~ 1|block,
                 weights = varIdent(form=~1|date),
                 data    = dat)

what tidy gives me

broom.mixed::tidy(mod, effects="ran_pars")
#> Registered S3 method overwritten by 'broom.mixed':
#>   method      from 
#>   tidy.gamlss broom
#> # A tibble: 2 x 4
#>   effect   group    term           estimate
#>   <chr>    <chr>    <chr>             <dbl>
#> 1 ran_pars block    sd_(Intercept)     1.26
#> 2 ran_pars Residual sd_Observation     2.08

what I want (ignoring the block variance for now)

mod$modelStruct$varStruct %>% 
  coef(unconstrained=FALSE, allCoef=TRUE) %>% 
  tibble::enframe(name="grp", value="varStruct") %>% 
  mutate(sigma         = mod$sigma) %>% 
  mutate(StandardError = sigma*varStruct) %>% 
  mutate(Variance      = StandardError^2)
#> # A tibble: 2 x 5
#>   grp       varStruct sigma StandardError Variance
#>   <chr>         <dbl> <dbl>         <dbl>    <dbl>
#> 1 21Aug1990      1     2.08          2.08     4.31
#> 2 28Aug1990      1.90  2.08          3.94    15.5

Created on 2020-12-08 by the reprex package (v0.3.0.9001) broom.mixed 0.2.6 nlme 3.1-149

Thanks in advance!

bbolker commented 4 years ago

I'll take a look when I get a chance (pull requests welcome!)

daranzolin commented 3 years ago

Ah, I was running around looking for these values too!

SchmidtPaul commented 3 years ago

I found out about this alternative solution since then:

# ...continuing reprex from above

# remotes::install_github('m-clark/mixedup')
mod %>% mixedup::extract_het_var()
#>   X21Aug1990 X28Aug1990
#> 1      4.306     15.505

Here's the permalink to the code.