Generalized Linear Regression • tidypredict

Highlights & Limitations

Defaults to 0-to-1 predictions for binomial family models. That is akin to running predict(model, type = "response")
Only treatment contrasts (contr.treatment) are supported.
offset is supported
Categorical variables are supported
In-line functions in the formulas are not supported:
- OK - wt ~ mpg + am
- OK - mutate(mtcars, newam = paste0(am)) and then wt ~ mpg + newam
- Not OK - wt ~ mpg + as.factor(am)
- Not OK - wt ~ mpg + as.character(am)
Interval functions are not supported: tidypredict_interval() & tidypredict_sql_interval()
The probit link is approximated rather than reproduced exactly. See The probit link.

How it works

library(tidypredict)
library(dplyr)

df <- mtcars %>%
  mutate(char_cyl = paste0("cyl", cyl)) %>%
  select(wt, char_cyl, am)

model <- glm(am ~ wt + char_cyl, data = df, family = "binomial")

It returns a SQL query that contains the coefficients (model) evaluated against the correct variable or categorical variable value. In most cases the resulting SQL is one short CASE WHEN statement per coefficient. It appends the offset field or value, if one is provided.

For binomial models, the sigmoid equation is applied. This means that the target SQL database type will need to support the exponent function.

library(tidypredict)
tidypredict_sql(model, dbplyr::simulate_mssql())
#> <SQL> 1.0 / (1.0 + EXP(-(((20.8527831345691 + ([wt] * -7.85934263583836)) + (CASE WHEN ([char_cyl] = 'cyl6') THEN 1.0 WHEN NOT ([char_cyl] = 'cyl6') THEN 0.0 END * 3.10462643177453)) + (CASE WHEN ([char_cyl] = 'cyl8') THEN 1.0 WHEN NOT ([char_cyl] = 'cyl8') THEN 0.0 END * 5.37942092366098))))

Alternatively, use tidypredict_to_column() if the results are to be used or previewed in dplyr.

df %>%
  tidypredict_to_column(model) %>%
  head(10)
#>                      wt char_cyl am        fit
#> Mazda RX4         2.620     cyl6  1 0.96662269
#> Mazda RX4 Wag     2.875     cyl6  1 0.79605201
#> Datsun 710        2.320     cyl4  1 0.93208127
#> Hornet 4 Drive    3.215     cyl6  0 0.21242376
#> Hornet Sportabout 3.440     cyl8  0 0.30918450
#> Valiant           3.460     cyl6  0 0.03783629
#> Duster 360        3.570     cyl8  0 0.13875740
#> Merc 240D         3.190     cyl4  0 0.01450687
#> Merc 230          3.150     cyl4  0 0.01975984
#> Merc 280          3.440     cyl6  0 0.04399324

The probit link

Every inverse link tidypredict writes is exact, with one exception: probit. The probit inverse link is the standard normal CDF, pnorm(), and no SQL backend has one, so it is written as the Bowling et al. logistic approximation instead:

$\frac{1}{1 + \exp(-0.07056\, x^3 - 1.5976\, x)}$

The same expression is used on both the R and the SQL paths, so a probit model is the one place where tidypredict_fit() does not reproduce predict() to floating-point precision. The approximation’s error is about 0.014% of the probability, which works out to roughly 1e-4:

probit_model <- glm(
  am ~ wt + mpg,
  data = mtcars,
  family = binomial(link = "probit")
)

max(abs(
  predict(probit_model, mtcars, type = "response") -
    rlang::eval_tidy(tidypredict_fit(probit_model), mtcars)
))
#> [1] 0.0001389629

That is four orders of magnitude larger than the disagreement any other link produces, and larger than tidypredict_test()’s default threshold, so a probit model will be reported as failing:

tidypredict_test(probit_model)
#> tidypredict test results
#> Difference threshold: 1e-12
#> 
#> Fitted records above the threshold: 31
#> 
#> Max difference: 0.00013896287832874

The difference is the approximation, not a defect in the parsed model. If the exact probabilities matter more than a portable formula, pass a threshold that reflects the approximation’s error, or use predict() directly.

Under the hood

The parser reads several parts of the glm object to tabulate all of the needed variables. One entry per coefficient is added to the final table. Other variables are added at the end. Some variables are not required for every parsed model. For example, offset is listed because it’s part of the formula (call) of the model, if there were no offset in a given model, that line would not exist.

pm <- parse_model(model)
str(pm, 2)
#> List of 2
#>  $ general:List of 8
#>   ..$ model   : chr "glm"
#>   ..$ version : num 2
#>   ..$ type    : chr "regression"
#>   ..$ residual: int 28
#>   ..$ sigma2  : num 1
#>   ..$ family  : chr "binomial"
#>   ..$ link    : chr "logit"
#>   ..$ is_glm  : num 1
#>  $ terms  :List of 4
#>   ..$ :List of 5
#>   ..$ :List of 5
#>   ..$ :List of 5
#>   ..$ :List of 5
#>  - attr(*, "class")= chr [1:3] "parsed_model" "pm_regression" "list"

The output from parse_model() is transformed into a dplyr, a.k.a. Tidy Eval, formula. All categorical variables are evaluated using if_else().

tidypredict_fit(model)
#> 1/(1 + exp(-(20.8527831345691 + (wt * -7.85934263583836) + (ifelse(char_cyl == 
#>     "cyl6", 1, 0) * 3.10462643177453) + (ifelse(char_cyl == "cyl8", 
#>     1, 0) * 5.37942092366098))))

From there, the Tidy Eval formula can be used anywhere it can be evaluated. tidypredict provides three paths:

Use directly inside dplyr, mutate(df, !! tidypredict_fit(model))
Use tidypredict_to_column(model) to add it to a piped command set
Use tidypredict_sql(model, con) to retrieve the SQL statement

Prediction intervals are not available for glm models, so tidypredict_interval() and tidypredict_sql_interval() have no glm counterpart.

How it performs

Testing the tidypredict results is easy. The tidypredict_test() function automatically uses the glm model object’s data frame to compare tidypredict_fit() to the results given by predict()

tidypredict_test(model)
#> tidypredict test results
#> Difference threshold: 1e-12
#> 
#>  All results are within the difference threshold

parsnip

tidypredict also supports glm() model objects fitted via the parsnip package, using linear_reg() with the "glm" engine.

library(parsnip)

parsnip_model <- linear_reg() %>%
  set_engine("glm") %>%
  fit(am ~ wt + cyl, data = mtcars)

tidypredict_fit(parsnip_model)
#> 1.5203311478662 + (wt * -0.3729886164835) + (cyl * 0.0138854914772273)

LiblineaR

Binary logistic regression models fitted with LiblineaR::LiblineaR() are also supported, including logistic_reg() models fitted via parsnip with the "LiblineaR" engine. As with glm() binomial models, predictions are on the 0-to-1 probability scale for the second factor level of the outcome.

liblinear_model <- logistic_reg(penalty = 0.1) %>%
  set_engine("LiblineaR") %>%
  fit(factor(am) ~ mpg + cyl, data = mtcars)

tidypredict_fit(liblinear_model)
#> 1/(1 + exp(-(-1.78560849993736 + (mpg * 0.166363458887312) + 
#>     (cyl * -0.324861381084944))))