Large standard error of prediction from parsnip vs base R

user_9408126 0
Viewed 125

It seems like predict is producing a standard error that is too large. I get 0.820 with a parsnip model but 0.194 with a base R model. 0.194 for a standard error seems more reasonable since about 2*0.195 above and below my prediction are the ends of the confidence interval. What is my problem/misunderstanding?

library(parsnip)
library(dplyr)

# example data
mod_dat <- mtcars %>%
  as_tibble() %>%
  mutate(cyl_8 = as.numeric(cyl == 8)) %>%
  select(mpg, cyl_8)

parsnip_mod <- logistic_reg() %>%
  set_engine("glm") %>%
  fit(as.factor(cyl_8) ~ mpg, data = mod_dat)

base_mod <- glm(as.factor(cyl_8) ~ mpg, data = mod_dat, family = "binomial")

parsnip_pred <- tibble(mpg = 18) %>%
  bind_cols(predict(parsnip_mod, new_data = ., type = 'prob'),
            predict(parsnip_mod, new_data = ., type = 'conf_int', std_error = T)) %>%
  select(!ends_with("_0"))

base_pred <- predict(base_mod, tibble(mpg = 18), se.fit = T, type = "response") %>%
  unlist()

# these give the same prediction but different SE
parsnip_pred
#> # A tibble: 1 x 5
#>     mpg .pred_1 .pred_lower_1 .pred_upper_1 .std_error
#>   <dbl>   <dbl>         <dbl>         <dbl>      <dbl>
#> 1    18   0.614         0.230         0.895      0.820
base_pred
#>          fit.1       se.fit.1 residual.scale 
#>      0.6140551      0.1942435      1.0000000

Created on 2020-06-04 by the reprex package (v0.3.0)

--EDIT--

As @thelatemail and @Limey said, using type="link" for the base model will give the standard error on the logit scale (0.820). However, I want the standard error on the probability scale. Is there an option in the parsnip documentation that I'm missing? I would like to use parsnip.

tidymodels predict confidence-interval r
2 Answers

@thelatemail is correct. From the online doc for predict.glm:

type
the type of prediction required. The default is on the scale of the linear predictors; the alternative "response" is on the scale of the response variable. Thus for a default binomial model the default predictions are of log-odds (probabilities on logit scale) and type = "response" gives the predicted probabilities.

The default is to report using the logit scale,, 'response' requests results on the raw probability scale. It's not obvious from the parsnip::predict documentation that I found how that chooses the scale on which to return its results, but it's clear it's using the raw probability scale.

So both methods are returning correct answers, they're just using different scales.

I don't want to steal an accepted solution from @thelatemail, so invite them to post a similar answer to this.

user_134348710

As @thelatemail said, you can get the standard error on the probability scale with parsnip using the arguments: type="raw", opts=list(se.fit=TRUE, type="response"). But at that point, you might as well use a base model since the output is exactly the same. However, this is still useful if you are already using a parsnip model and you want the standard error output of a base model.

library(parsnip)
library(dplyr)

mod_dat <- mtcars %>%
  as_tibble() %>%
  mutate(cyl_8 = as.numeric(cyl == 8)) %>%
  select(mpg, cyl_8)

parsnip_mod <- logistic_reg() %>%
  set_engine("glm") %>%
  fit(as.factor(cyl_8) ~ mpg, data = mod_dat)

base_mod <- glm(as.factor(cyl_8) ~ mpg, data = mod_dat, family = "binomial")

predict(parsnip_mod, tibble(mpg = 18), type="raw",
        opts=list(se.fit=TRUE, type="response")) %>% 
  as_tibble()
#> # A tibble: 1 x 3
#>     fit se.fit residual.scale
#>   <dbl>  <dbl>          <dbl>
#> 1 0.614  0.194              1

predict.glm(base_mod, tibble(mpg = 18), se.fit = T, type="response") %>% 
  as_tibble()
#> # A tibble: 1 x 3
#>     fit se.fit residual.scale
#>   <dbl>  <dbl>          <dbl>
#> 1 0.614  0.194              1

Created on 2020-06-11 by the reprex package (v0.3.0)

user_94081260
Related