---
title: "NaiveBayes models"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{NaiveBayes models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
if (
  requireNamespace("klaR", quietly = TRUE) &&
    requireNamespace("naivebayes", quietly = TRUE)
) {
  library(tidypredict)
  library(dplyr)
  eval_code <- TRUE
} else {
  eval_code <- FALSE
}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = eval_code
)
```

| Function                                                      |Works|
|---------------------------------------------------------------|-----|
|`tidypredict_fit()`, `tidypredict_sql()`, `parse_model()`      |  ✔  |
|`tidypredict_to_column()`                                      |  ✗  |
|`tidypredict_test()`                                           |  ✗  |
|`tidypredict_interval()`, `tidypredict_sql_interval()`         |  ✗  |
|`parsnip`                                                      |  ✔  |

`klaR::NaiveBayes()` and `naivebayes::naive_bayes()` fit naive Bayes
classifiers. Predicting with such a model
multiplies the prior probability of each class by one conditional density per
predictor, and then normalizes those products into posterior probabilities.
Working on the log scale turns the products into sums, which makes the posterior
probabilities the softmax of the summed log densities.

Numeric predictors contribute a normal log density, and categorical predictors
contribute a `case_when()` lookup of the conditional probability of the observed
level. `naivebayes::naive_bayes()` additionally supports Poisson densities for
integer predictors when fit with `usepoisson = TRUE`. Kernel density estimates
cannot be expressed this way, so models fit with `usekernel = TRUE` are not
supported.

Because these models predict one probability per outcome class,
`tidypredict_fit()` returns a *named list* of expressions, one for each class,
rather than a single expression. Since the output is a list,
`tidypredict_to_column()` and `tidypredict_test()` are not supported.

## `tidypredict_` functions

```{r}
model <- klaR::NaiveBayes(Species ~ ., data = iris)
```

- Create the R formulas, one per class
    ```{r}
fit <- tidypredict_fit(model)
names(fit)
fit[["setosa"]]
    ```

- Add the predictions to the original table
    ```{r}
library(dplyr)

iris %>%
  mutate(!!!tidypredict_fit(model)) %>%
  glimpse()
    ```

- Confirm that the results match the model's `predict()` results
    ```{r}
probs <- sapply(fit, \(f) rlang::eval_tidy(f, iris))
all.equal(unname(probs), unname(predict(model, iris)$posterior))
    ```

### Extreme values of a numeric predictor

There is one case where `tidypredict_fit()` does not reproduce `predict()`, and
it is deliberately left in place.

Both packages evaluate each normal density and then replace any density that
came out as exactly zero with their `threshold` argument, `0.001` by default. A
normal density only reaches zero by underflowing the smallest number a double
can hold, which takes a value roughly 38 standard deviations from that class's
mean. `tidypredict` works on the log scale throughout, where such a value is an
ordinary negative number and nothing underflows, so the substitution never
happens.

The substitution is drastic when it fires. Consider an iris with a
`Sepal.Length` of 20, well beyond the largest in the data:

```{r}
outlier <- iris[1, ]
outlier$Sepal.Length <- 20

predict(model, outlier)$posterior
```

`setosa` has the narrowest spread of the three classes, so its density is the
first to underflow. Replacing it with `0.001` makes it enormous next to the
other two, which are around `1e-162` and `1e-97` but were computed honestly, so
the class that fits the value worst is the one that wins. This happens silently
unless the value is extreme enough for *every* class to underflow, which is the
only case either package warns about. `tidypredict` instead continues on the
log scale and reports the ordering the fitted model implies:

```{r}
sapply(tidypredict_fit(model), \(f) rlang::eval_tidy(f, outlier))
```

Reproducing the substitution is not possible on the log scale in any case.
Testing `exp(log_density) == 0` instead of `density == 0` picks out a different
set of values: over 200,000 draws from the band where the underflow happens the
two tests disagree on about 1,000 of them, because the two routes round
differently in the denormal range. A SQL backend, working in its own floating
point, would disagree again in a third way.

## parsnip

`parsnip` fitted models are also supported by `tidypredict`:
```{r}
library(parsnip)
library(discrim)

p_model <- naive_Bayes() %>%
  set_engine("klaR", usekernel = FALSE) %>%
  fit(Species ~ ., data = iris)
```

```{r}
tidypredict_fit(p_model)[["virginica"]]
```

## `naivebayes::naive_bayes()`

The `naivebayes` package is supported in the same way:
```{r}
nb_model <- naivebayes::naive_bayes(Species ~ ., data = iris)

nb_fit <- tidypredict_fit(nb_model)
nb_fit[["setosa"]]
```

```{r}
nb_probs <- sapply(nb_fit, \(f) rlang::eval_tidy(f, iris))
all.equal(
  unname(nb_probs),
  unname(predict(nb_model, iris[names(nb_model$tables)], type = "prob"))
)
```

The `"naivebayes"` parsnip engine works too, as long as `usekernel = FALSE`:
```{r}
nb_p_model <- naive_Bayes() %>%
  set_engine("naivebayes", usekernel = FALSE) %>%
  fit(Species ~ ., data = iris)

tidypredict_fit(nb_p_model)[["virginica"]]
```

## Parse model spec

Here is an example of the model spec:
```{r}
pm <- parse_model(model)
str(pm, 2)
```