Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

VarStruct information from lme tidy missing? #96

Open
SchmidtPaul opened this issue May 15, 2020 · 3 comments
Open

VarStruct information from lme tidy missing? #96

SchmidtPaul opened this issue May 15, 2020 · 3 comments

Comments

@SchmidtPaul
Copy link

SchmidtPaul commented May 15, 2020

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
Copy link
Owner

bbolker commented May 17, 2020

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

@bbolker bbolker changed the title VarStruct information from lmer tidy missing? VarStruct information from lme tidy missing? Aug 28, 2020
@daranzolin
Copy link

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

@SchmidtPaul
Copy link
Author

SchmidtPaul commented Apr 6, 2021

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants