sda models

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

sda::sda() fits shrinkage discriminant analysis models, including the diagonal variant (diagonal = TRUE). Predicting with such a model is a softmax over one linear predictor per outcome class, so 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.

Only the features that survive shrinkage appear in the fitted model, and the generated expressions reference just those.

`tidypredict_` functions

model <- sda::sda(as.matrix(iris[1:4]), iris$Species, verbose = FALSE)

Create the R formulas, one per class

fit <- tidypredict_fit(model)
names(fit)
#> [1] "setosa"     "versicolor" "virginica"
fit[["setosa"]]
#> exp(-13.1507094677661 + (Sepal.Length * 5.73266109594551) + (Sepal.Width * 
#>     11.3604293693889) + (Petal.Length * -16.8084846238171) + 
#>     (Petal.Width * -16.6240117338795))/(exp(-13.1507094677661 + 
#>     (Sepal.Length * 5.73266109594551) + (Sepal.Width * 11.3604293693889) + 
#>     (Petal.Length * -16.8084846238171) + (Petal.Width * -16.6240117338795)) + 
#>     exp(-2.17403913501322 + (Sepal.Length * -1.38327577793292) + 
#>         (Sepal.Width * -4.20634050168124) + (Petal.Length * 4.57106120700443) + 
#>         (Petal.Width * 2.49705634843865)) + exp(-32.2372116751408 + 
#>     (Sepal.Length * -4.3493853180126) + (Sepal.Width * -7.15408886770766) + 
#>     (Petal.Length * 12.2374234168126) + (Petal.Width * 14.1269553854408)))

Add the predictions to the original table

library(dplyr)

iris %>%
  mutate(!!!tidypredict_fit(model)) %>%
  glimpse()
#> Rows: 150
#> Columns: 8
#> $ Sepal.Length <dbl> 5.1, 4.9, 4.7, 4.6, 5.0, 5.4, 4.6, 5.0, 4.4, 4.9…
#> $ Sepal.Width  <dbl> 3.5, 3.0, 3.2, 3.1, 3.6, 3.9, 3.4, 3.4, 2.9, 3.1…
#> $ Petal.Length <dbl> 1.4, 1.4, 1.3, 1.5, 1.4, 1.7, 1.4, 1.5, 1.4, 1.5…
#> $ Petal.Width  <dbl> 0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.3, 0.2, 0.2, 0.1…
#> $ Species      <fct> setosa, setosa, setosa, setosa, setosa, setosa, …
#> $ setosa       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
#> $ versicolor   <dbl> 1.009018e-20, 1.005223e-16, 2.186404e-18, 1.5199…
#> $ virginica    <dbl> 3.704434e-39, 2.916177e-34, 2.957542e-36, 1.7209…

Confirm that the results match the model’s predict() results

probs <- sapply(fit, \(f) rlang::eval_tidy(f, iris))
posterior <- sda::predict.sda(model, as.matrix(iris[1:4]), verbose = FALSE)
all.equal(unname(probs), unname(posterior$posterior))
#> [1] "Mean relative difference: 4.212717e-08"

sda() rounds its posterior probabilities with zapsmall(), so expect agreement to about seven decimal places rather than exactly.

parsnip

parsnip fitted models are also supported by tidypredict:

library(parsnip)
library(discrim)

p_model <- discrim_linear() %>%
  set_engine("sda") %>%
  fit(Species ~ ., data = iris)

tidypredict_fit(p_model)[["virginica"]]
#> exp(-32.2372116751408 + (Sepal.Length * -4.3493853180126) + (Sepal.Width * 
#>     -7.15408886770766) + (Petal.Length * 12.2374234168126) + 
#>     (Petal.Width * 14.1269553854408))/(exp(-13.1507094677661 + 
#>     (Sepal.Length * 5.73266109594551) + (Sepal.Width * 11.3604293693889) + 
#>     (Petal.Length * -16.8084846238171) + (Petal.Width * -16.6240117338795)) + 
#>     exp(-2.17403913501322 + (Sepal.Length * -1.38327577793292) + 
#>         (Sepal.Width * -4.20634050168124) + (Petal.Length * 4.57106120700443) + 
#>         (Petal.Width * 2.49705634843865)) + exp(-32.2372116751408 + 
#>     (Sepal.Length * -4.3493853180126) + (Sepal.Width * -7.15408886770766) + 
#>     (Petal.Length * 12.2374234168126) + (Petal.Width * 14.1269553854408)))

sda() is fit from a numeric matrix, so a model fit directly can only refer to the matrix columns it was given. Categorical predictors therefore work through the parsnip interface, which keeps the formula around and lets the dummy columns be written in terms of the original factors:

cars <- transform(mtcars, cyl = factor(cyl), gear = factor(gear))

c_model <- discrim_linear() %>%
  set_engine("sda") %>%
  fit(cyl ~ mpg + gear + disp, data = cars)

tidypredict_fit(c_model)[["8"]]
#> exp(-9.69338951858338 + (mpg * -0.300898442298615) + (ifelse(gear == 
#>     "4", 1, 0) * -0.288744001063274) + (ifelse(gear == "5", 1, 
#>     0) * 0.346792280462239) + (disp * 0.0491641125277206))/(exp(-5.90478040863942 + 
#>     (mpg * 0.450539182424721) + (ifelse(gear == "4", 1, 0) * 
#>     0.120334811858762) + (ifelse(gear == "5", 1, 0) * -0.326397688436066) + 
#>     (disp * -0.0361970033802603)) + exp(4.8234106925129 + (mpg * 
#>     -0.14546309773062) + (ifelse(gear == "4", 1, 0) * 0.18956123261121) + 
#>     (ifelse(gear == "5", 1, 0) * -0.0358060646840222) + (disp * 
#>     -0.0157055813201664)) + exp(-9.69338951858338 + (mpg * -0.300898442298615) + 
#>     (ifelse(gear == "4", 1, 0) * -0.288744001063274) + (ifelse(gear == 
#>     "5", 1, 0) * 0.346792280462239) + (disp * 0.0491641125277206)))

Parse model spec

Here is an example of the model spec:

pm <- parse_model(model)
str(pm, 2)
#> List of 3
#>  $ general    :List of 4
#>   ..$ model  : chr "sda"
#>   ..$ version: num 2
#>   ..$ type   : chr "multiclass_regression"
#>   ..$ family : chr "multinomial"
#>  $ classes    : chr [1:3] "setosa" "versicolor" "virginica"
#>  $ class_terms:List of 3
#>   ..$ :List of 5
#>   ..$ :List of 5
#>   ..$ :List of 5
#>  - attr(*, "class")= chr [1:3] "parsed_model" "pm_multiclass_regression" "list"

tidypredict_ functions

parsnip

Parse model spec

`tidypredict_` functions