Grab other model attributes from `fit()`

[AmeliaMN]

< -- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->

Feature

I use the Stat2 textbook in my regression course, and one homework problem asks students to do a randomization test for a multiple linear regression model. But rather than collecting the coefficient estimates from each randomization sample, the question asks students to find some summary statistic about the model, e.g. $R^2$. In other words, instead of running

null_fits <- gss %>%
  specify(hours ~ age + college) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 100, type = "permute") %>%
  fit()

and getting a bunch of null slopes for age and college, I would like to get a bunch of null $R^2$ values. I know I can do this with purrr,

null_fits <- gss %>%
  specify(hours ~ age + college) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 100, type = "permute")  %>%
  split(.$replicate) %>%
  map(~lm(hours~age+college, data = .x)) %>%
  map(summary) %>%
  map_dbl("r.squared")

But I wish it was better-facilitated by infer itself.

[simonpcouch]

I appreciate you bringing this up, @AmeliaMN!

For now, fit.infer just applies tidy() to the output of each glm() fit. We could maybe open up an interface here for passing an arbitrary model summary function? For the output you mention, folks could pass glance, but I could also imagine folks passing augment, summary, or identity. With glance, that interface could feel something like:

library(tidymodels)

null_dist <- 
   gss %>%
   specify(hours ~ age + college) %>%
   hypothesize(null = "independence") %>%
   generate(reps = 2, type = "permute")

# continue to `tidy() %>% select()` by default
null_dist %>%
   fit()
#> # A tibble: 6 × 3
#> # Groups:   replicate [2]
#>   replicate term          estimate
#>       <int> <chr>            <dbl>
#> 1         1 intercept     41.7    
#> 2         1 age            0.00508
#> 3         1 collegedegree -1.43   
#> 4         2 intercept     40.4    
#> 5         2 age            0.0214 
#> 6         2 collegedegree  0.328

# optional `summary` (or otherwise named) argument
null_dist %>%
   fit(summary = glance)
#> # A tibble: 2 × 13
#> # Groups:   replicate [2]
#>   replicate r.squared adj.r.squ…¹ sigma stati…² p.value    df logLik   AIC   BIC
#>       <int>     <dbl>       <dbl> <dbl>   <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl>
#> 1         1   0.00121    -0.00281  14.8   0.302   0.740     2 -2057. 4121. 4138.
#> 2         2   0.00275    -0.00126  14.8   0.686   0.504     2 -2056. 4121. 4137.
#> # … with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>, and
#> #   abbreviated variable names ¹​adj.r.squared, ²​statistic

Actual code for this output:

null_dist %>% 
   nest() %>% 
   rowwise() %>% 
   mutate(mod = list(lm(hours ~ age + college, data = data)), 
          sum = list(glance(mod))) %>% 
   select(-data, -mod) %>% 
   unnest(cols = sum)

Any thoughts on this interface? This feels like a nice way of integrating broom and functions-as-arguments into teaching these pipelines, though I'm not sure if that's Astronaut Stuff Too Early.

One technical hiccup here is that glance() doesn't output $R^2$ for glm() output, though glm() is what's used under the hood in fit.infer() right now.

Also, looping in @mine-cetinkaya-rundel. :)