tidiers for gamlss models

[IndrajeetPatil]

There are a few issues with the tidier function for gamlss models:

1. Although confidence intervals are available for estimate, the tidier doesn't display them.
2. glance shows that the degress of freedom are 0, which doesn't match with the degrees of freedom found with summary function. (This might be confusion on my part about how df is computed here, so ignore this if I am wrong.)
3. Output is not a tibble, like the rest of the tidiers.

# needed libraries
library(gamlss)
#> Loading required package: splines
#> Loading required package: gamlss.data
#> Loading required package: gamlss.dist
#> Loading required package: MASS
#> Loading required package: nlme
#> Loading required package: parallel
#>  **********   GAMLSS Version 5.1-2  **********
#> For more on GAMLSS look at http://www.gamlss.org/
#> Type gamlssNews() to see new features/changes/bug fixes.
library(tidyverse)

# model
g <- gamlss(
  y ~ pb(x),
  sigma.fo = ~ pb(x),
  family = BCT,
  data = abdom,
  method = mixed(1, 20)
)
#> GAMLSS-RS iteration 1: Global Deviance = 4771.925 
#> GAMLSS-CG iteration 1: Global Deviance = 4771.013 
#> GAMLSS-CG iteration 2: Global Deviance = 4770.994 
#> GAMLSS-CG iteration 3: Global Deviance = 4770.994

# model summary
summary(g)
#> ******************************************************************
#> Family:  c("BCT", "Box-Cox t") 
#> 
#> Call:  gamlss(formula = y ~ pb(x), sigma.formula = ~pb(x),  
#>     family = BCT, data = abdom, method = mixed(1, 20)) 
#> 
#> Fitting method: mixed(1, 20) 
#> 
#> ------------------------------------------------------------------
#> Mu link function:  identity
#> Mu Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) -64.44299    1.33397  -48.31   <2e-16 ***
#> pb(x)        10.69464    0.05787  184.80   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------------------------------
#> Sigma link function:  log
#> Sigma Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) -2.65041    0.10859 -24.407  < 2e-16 ***
#> pb(x)       -0.01002    0.00380  -2.638  0.00855 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------------------------------
#> Nu link function:  identity 
#> Nu Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)
#> (Intercept)  -0.1072     0.6296   -0.17    0.865
#> 
#> ------------------------------------------------------------------
#> Tau link function:  log 
#> Tau Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   2.4948     0.4261   5.855 7.86e-09 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------------------------------
#> NOTE: Additive smoothing terms exist in the formulas: 
#>  i) Std. Error for smoothers are for the linear effect only. 
#> ii) Std. Error for the linear terms maybe are not accurate. 
#> ------------------------------------------------------------------
#> No. of observations in the fit:  610 
#> Degrees of Freedom for the fit:  11.7603
#>       Residual Deg. of Freedom:  598.2397 
#>                       at cycle:  3 
#>  
#> Global Deviance:     4770.994 
#>             AIC:     4794.515 
#>             SBC:     4846.419 
#> ******************************************************************

# getting confidence intervals
confint(g)
#> Warning in vcov.gamlss(object, robust = robust): Additive terms exists in the  mu formula. 
#>   Standard errors for the linear terms maybe are not appropriate
#> Warning in vcov.gamlss(object, robust = robust): Additive terms exists in the  sigma formula. 
#>   Standard errors for the linear terms maybe are not appropriate
#>                 2.5 %    97.5 %
#> (Intercept) -67.05752 -61.82847
#> pb(x)        10.58121  10.80806

# tidying with broom
broom::tidy(g, conf.int = TRUE)
#> # A tibble: 6 x 6
#>   parameter term        estimate std.error statistic   p.value
#>   <chr>     <chr>          <dbl>     <dbl>     <dbl>     <dbl>
#> 1 mu        (Intercept) -64.4      1.33      -48.5   1.89e-210
#> 2 mu        pb(x)        10.7      0.0578    185.    0.       
#> 3 sigma     (Intercept)  -2.65     0.108     -24.5   8.09e- 93
#> 4 sigma     pb(x)        -0.0100   0.00378    -2.65  8.29e-  3
#> 5 nu        (Intercept)  -0.107    0.557      -0.192 8.48e-  1
#> 6 tau       (Intercept)   2.49     0.301       8.28  7.77e- 16

# glance with broom
broom::glance(g)
#> # A tibble: 1 x 6
#>      df logLik   AIC   BIC deviance df.residual
#>   <int>  <dbl> <dbl> <dbl>    <dbl>       <dbl>
#> 1     0 -2385. 4795. 4846.    4771.        598.

^{Created on 2018-10-07 by the reprex package (v0.2.1)}

[IndrajeetPatil]

Hmm, this tidier has moved to broom.mixed, right?

Tagging @bbolker here. Not sure how I can transfer this issue to that repo.

[alexpghayes]

I would open a new issue in broom.mixed and link it to this one.

[IndrajeetPatil]

@simonpcouch Do you think the tidier for gamlss can be removed from broom for the next release?

Currently, it is both in broom and broom.mixed and gives the following warning if both are attached:

Registered S3 method overwritten by 'broom.mixed':
  method      from 
  tidy.gamlss broom

[simonpcouch]

Hi @IndrajeetPatil! Thanks for checking in on this. Will run reverse dependency checks with galmss tidiers removed and see how many new issues we run into, but I'll plan on removing them.🙂

[IndrajeetPatil]

@simonpcouch Have you give this further thought? Are they going to be retained in broom or will they be retired in favor of broom.mixed?

[simonpcouch]

tidy.gamlss is now soft-deprecated! Thanks for the nudge, @IndrajeetPatil.

[IndrajeetPatil]

Given that gamlss tidiers were now removed from broom.mixed as well. This change in broom will have to be reverted :(
bbolker/broom.mixed#64

[bbolker]

Sorry about this, I've reverted. gamlss will continue to live in broom.mixed

[github-actions]

This issue has been automatically locked. If you believe you have found a related problem, please file a new issue (with a reprex: https://reprex.tidyverse.org) and link to this issue.