---
title: "Generalized Linear Regression"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Generalized Linear Regression}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

source("_threads.R")
library(dplyr)
library(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](#the-probit-link).

## How it works

```{r}
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](https://en.wikipedia.org/wiki/Sigmoid_function) equation is applied. This means that the target SQL database type will need to support the exponent function.

```{r}
library(tidypredict)
tidypredict_sql(model, dbplyr::simulate_mssql())
```

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

```{r}
df %>%
  tidypredict_to_column(model) %>%
  head(10)
```

## 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:

```{r}
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)
))
```

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:

```{r}
tidypredict_test(probit_model)
```

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.

```{r}
pm <- parse_model(model)
str(pm, 2)
```

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

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()`

```{r}
tidypredict_test(model)
```

## parsnip

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

```{r}
library(parsnip)

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

tidypredict_fit(parsnip_model)
```

## 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.

```{r, eval = requireNamespace("LiblineaR", quietly = TRUE)}
liblinear_model <- logistic_reg(penalty = 0.1) %>%
  set_engine("LiblineaR") %>%
  fit(factor(am) ~ mpg + cyl, data = mtcars)

tidypredict_fit(liblinear_model)
```