---
title: "custom transformers + <br/> patsy & statsmodels"
subtitle: "Lecture 12"
author: "Dr. Colin Rundel"
footer: "Sta 663 - Spring 2025"
format:
  revealjs:
    theme: slides.scss
    transition: fade
    slide-number: true
    self-contained: true
execute: 
  echo: true
---

```{python setup}
#| include: false
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

import sklearn
sklearn.set_config(display="text")

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn.model_selection import GridSearchCV, KFold, StratifiedKFold, train_test_split
from sklearn.metrics import classification_report

plt.rcParams['figure.dpi'] = 200

np.set_printoptions(
  edgeitems=30, linewidth=75,
  precision = 4, suppress=True
  #formatter=dict(float=lambda x: "%.5g" % x)
)

pd.set_option("display.width", 100)
pd.set_option("display.max_columns", 10)
pd.set_option("display.precision", 4)
```

```{r r_setup}
#| include: false
knitr::opts_chunk$set(
  fig.align="center",
  cache=FALSE
)

library(lme4)


local({
  hook_err_old <- knitr::knit_hooks$get("error")  # save the old hook
  knitr::knit_hooks$set(error = function(x, options) {
    # now do whatever you want to do with x, and pass
    # the new x to the old hook
    x = sub("## \n## Detailed traceback:\n.*$", "", x)
    x = sub("Error in py_call_impl\\(.*?\\)\\: ", "", x)
    hook_err_old(x, options)
  })
  
  hook_warn_old <- knitr::knit_hooks$get("warning")  # save the old hook
  knitr::knit_hooks$set(warning = function(x, options) {
    x = sub("<string>:1: ", "", x)
    hook_warn_old(x, options)
  })
})
```


# Custom sklearn transformers

## FunctionTransformer

The simplest way to create a new transformer is to use `FunctionTransformer()` from the preprocessing submodule which allows for converting a Python function into a transformer.

::: {.xsmall}
```{python}
from sklearn.preprocessing import FunctionTransformer
X = pd.DataFrame({"x1": range(1,6), "x2": range(5, 0, -1)})
```
:::

:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
log_transform = FunctionTransformer(np.log)
lt = log_transform.fit(X)
lt
lt.transform(X)
```
:::

::: {.column width='50%' .fragment}
```{python}
lt.get_params()
dir(lt)
```
:::
::::


## Input types

::: {.xsmall}
```{python, error=TRUE}
def interact(X, y = None):
  return np.c_[X, X[:,0] * X[:,1]]
X = pd.DataFrame({"x1": range(1,6), "x2": range(5, 0, -1)})
Z = np.array(X)
```
:::

. . .

:::: {.columns .xsmall}
::: {.column width='50%'}
```{python, error=TRUE}
FunctionTransformer(
  interact
).fit_transform(X)
```
:::

::: {.column width='50%'}
```{python, error=TRUE}
FunctionTransformer(
  interact
).fit_transform(Z)
```
:::
::::

. . .

:::: {.columns .xsmall}
::: {.column width='50%'}
```{python, error=TRUE}
FunctionTransformer(
  interact, validate=True
).fit_transform(X)
```
:::

::: {.column width='50%'}
```{python, error=TRUE}
FunctionTransformer(
  interact, validate=True
).fit_transform(Z)
```
:::
::::

::: {.aside}
The `validate` argument both checks that `X` is 2d as well as converts it to a np.array
:::


## Build your own transformer

For a more full featured transformer, it is possible to construct it as a class that inherits from `BaseEstimator` and `TransformerMixin` classes from the `base` submodule.

::: {.small}
```{python}
from sklearn.base import BaseEstimator, TransformerMixin

class scaler(BaseEstimator, TransformerMixin):
  def __init__(self, m = 1, b = 0):
    self.m = m
    self.b = b
  
  def fit(self, X, y=None):
    return self
  
  def transform(self, X, y=None):
    return X*self.m + self.b
```
:::

##

::: {.xsmall}
```{python}
X = pd.DataFrame({
  "x1": range(1,6), 
  "x2": range(5, 0, -1)}
); X
```
:::

::: {.columns .xsmall}
::: {.column}
```{python}
double = scaler(2)
double.get_params()
double.fit_transform(X)
```
:::
::: {.column .fragment}
```{python}
double.set_params(b=-3).fit_transform(X)
```
:::
:::




## What else do we get?

```{python}
#| include: false
np.set_printoptions(
  linewidth=80
)
```

::: {.small}
```{python}
print(
  np.array(dir(double))
)
```
:::

# Demo - Interaction Transformer


## Useful methods

We employed a couple of special methods that are worth mentioning in a little more detail.

* `validate_data()` & `_check_feature_names()` are functions from `sklearn.base.validate` they are responsible for setting and checking the `n_features_in_` and the `feature_names_in_` attributes respectively.

* In general one or both is run during `fit()` with `reset=True` in which case the respective attribute will be set.

* Later, in `transform()` one or both will again be called with `reset=False` and the properties of `X` will be checked against the values in the attribute.

* These are worth using as they promote an interface consistent with sklearn and also provide convenient error checking with useful warning / error messages.



## `check_is_fitted()`

This is another useful helper function from `sklearn.utils` - it is fairly simplistic in that it checks for the existence of a specified attribute. If no attribute is given then it checks for any attributes ending in `_` that do not begin with `__`.

Again this is useful for providing a consistent interface and useful error / warning messages.

See also the other `check*()` functions in `sklearn.utils`.


## Other custom estimators

If you want to implement your own custom modeling function it is possible, there are different Mixin base classes in `sklearn.base` that provide the common core interface.

| Class                       | Description                               |
|-----------------------------|-------------------------------------------|
| `base.BiclusterMixin`       | Mixin class for all bicluster estimators  |
| `base.ClassifierMixin`      | Mixin class for all classifiers           |
| `base.ClusterMixin`         | Mixin class for all cluster estimators    |
| `base.DensityMixin`         | Mixin class for all density estimators    |
| `base.RegressorMixin`       | Mixin class for all regression estimators |
| `base.TransformerMixin`     | Mixin class for all transformers          |
| `base.OneToOneFeatureMixin` | Provides get_feature_names_out for simple transformers |


# patsy

## patsy

> `patsy` is a Python package for describing statistical models (especially linear models, or models that have a linear component) and building design matrices. It is closely inspired by and compatible with the formula mini-language used in R and S.
> 
>
>
> ...
>
>
>
> Patsy’s goal is to become the standard high-level interface to describing statistical models in Python, regardless of what particular model or library is being used underneath.


## Formulas

::: {.small}
```{python}
from patsy import ModelDesc
```

```{python}
ModelDesc.from_formula("y ~ a + a:b + np.log(x)")
ModelDesc.from_formula("y ~ a*b + np.log(x) - 1")
```
:::

::: {.aside}
This general syntax is known as [Wilkinson Notation](https://www.mathworks.com/help/stats/wilkinson-notation.html) and come from Wilkinson and Rogers. "Symbolic description of factorial models for analysis of variance." J. Royal Statistics Society 22, pp. 392–399, 1973.
:::


## Model matrix


::: {.xsmall}
```{python}
from patsy import demo_data, dmatrix, dmatrices
```
:::

:::: {.columns .xxsmall}
::: {.column width='35%'}
```{python}
data = demo_data(
  "y", "a", "b", "x1", "x2"
)
pd.DataFrame(data)
```
:::

::: {.column width='65%' .fragment}
```{python}
dmatrix("a + a:b + np.exp(x1)", data)
```
:::
::::

::: {.aside}
Note the `T.` in `a[T.a2]` is there to indicate treatment coding (i.e. typical dummy coding)
:::


## Model matrices

::: {.xsmall}
```{python}
y, X  = dmatrices("y ~ a + a:b + np.exp(x1)", data)
```
:::

:::: {.columns .xsmall}
::: {.column width='33%'}
```{python}
y
```
:::

::: {.column width='66%'}
```{python}
X
```
:::
::::


## as DataFrames

::: {.small}
```{python}
dmatrix("a + a:b + np.exp(x1)", data, return_type='dataframe')
```
:::


## Formula Syntax

<br/>

::: {.small}
| Code     | Description                                       | Example                |
|:--------:|:--------------------------------------------------|:-------------------------------|
| `+`      | unions terms on the left and right                | `a+a` ⇒ `a` |
| `-`      | removes terms on the right from terms on the left | `a+b-a` ⇒ `b` |
|`:`       | constructs interactions between each term on the<br/>left and right | `(a+b):c` ⇒  `a:c + b:c`|
| `*`      | short-hand for terms and their interactions       | `a*b` ⇒ `a + b + a:b` |
| `/`      | short-hand for left terms and their interactions with<br/>right terms | `a/b` ⇒ `a + a:b` |
| `I()`    | used for arithmetic calculations                  | `I(x1 + x2)` |
| `Q()`    | used to quote column names, e.g. columns with<br/>spaces or symbols | `Q('bad name!')` |
| `C()`    | used for categorical data coding                  |  `C(a, Treatment('a2'))` |
:::

## Examples

:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
dmatrix("x:y", demo_data("x","y","z"))
dmatrix("x*y", demo_data("x","y","z"))
```
:::

::: {.column width='50%' .fragment}
```{python}
dmatrix("x/y", demo_data("x","y","z"))
```
:::
::::

##

::: {.xsmall}
```{python}
dmatrix("x*(y+z)", demo_data("x","y","z"))
```
:::

## Intercept Examples (-1)

:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
dmatrix("x", demo_data("x","y","z"))
dmatrix("x-1", demo_data("x","y","z"))

```
:::

::: {.column width='50%'}
```{python}
dmatrix("-1 + x", demo_data("x","y","z"))
```
:::
::::

## Intercept Examples (0)

:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
dmatrix("x+0", demo_data("x","y","z"))
dmatrix("x-0", demo_data("x","y","z"))
```
:::

::: {.column width='50%'}
```{python}
dmatrix("x - (-0)", demo_data("x","y","z"))
```
:::
::::


## Design Info

One of the key features of the design matrix object is that it retains all the necessary details (including stateful transforms) that are necessary to apply to new data inputs (e.g. for prediction).

::: {.xsmall}
```{python}
d = dmatrix("a + a:b + np.exp(x1)", data, return_type='dataframe')
d.design_info
```
:::


## Stateful transforms

::: {.xsmall}
```{python}
data = {"x1": np.random.normal(size=10)}
new_data = {"x1": np.random.normal(size=10)}
```
:::

:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
d = dmatrix("scale(x1)", data)
d
np.mean(d, axis=0)
```
:::

::: {.column width='50%'}
```{python}
pred = dmatrix(d.design_info, new_data)
pred
np.mean(pred, axis=0)
```
:::
::::


## scikit-learn + Patsy

The state of affairs here is a bit of a mess at the moment - previously the `sklego` package implemented a `PatsyTransformer` class that has since been deprecated in favor of the `FormulaicTransformer` which uses the `formulaic` package for formula handling.

. . .

::: {.center style="font-size: 200px"}
🤷
:::

## A PatsyTransformer

::: {.xsmall}
```{python}
from patsy import dmatrix, build_design_matrices
from sklearn.utils.validation import check_is_fitted
from sklearn.base import BaseEstimator, TransformerMixin

class PatsyTransformer(TransformerMixin, BaseEstimator):
    def __init__(self, formula):
        self.formula = formula

    def fit(self, X, y=None):
        m = dmatrix(self.formula, X)
        assert np.array(m).shape[0] == np.array(X).shape[0]
        self.design_info_ = m.design_info
        return self

    def transform(self, X):
        check_is_fitted(self, 'design_info_')
        return build_design_matrices([self.design_info_], X)[0]
```
:::


::: {.aside}
Based on [`koaning/scikit-lego`](https://github.com/koaning/scikit-lego/commit/85c5074d94b54b4948bb2a662ee80c87a8f80df0)
:::

##

::: {.xxsmall}
```{python}
df = pd.DataFrame({
  "y": [2, 2, 4, 4, 6], "x": [1, 2, 3, 4, 5],
  "a": ["yes", "yes", "no", "no", "yes"]
})
X, y = df[["x", "a"]], df[["y"]].values
```
:::

:::: {.columns .xxsmall}
::: {.column width='47%'}
```{python}
pt = PatsyTransformer("x*a + np.log(x)")
pt.fit_transform(X)
```
:::

::: {.column width='53%' .fragment}
```{python}
make_pipeline(
  PatsyTransformer("x*a + np.log(x)"),
  StandardScaler()
).fit_transform(X)
```
:::
::::



# statsmodels

## statsmodels

> statsmodels is a Python module that provides classes and functions for the estimation of many different statistical models, as well as for conducting statistical tests, and statistical data exploration. An extensive list of result statistics are available for each estimator. The results are tested against existing statistical packages to ensure that they are correct.

::: {.small}

```{python}
import statsmodels.api as sm
import statsmodels.formula.api as smf
import statsmodels.tsa.api as tsa
```

:::

. . .

`statsmodels` uses slightly different terminology for referring to `y` (dependent / response) and `x`  (independent / explanatory) variables. Specifically it uses `endog` and `exog` to refer to `y` and `x` variable(s) respectively.

This is particularly important when using the main API, less so when using the formula API.


## OpenIntro Loans data

::: {.small}
> This data set represents thousands of loans made through the Lending Club platform, which is a platform that allows individuals to lend to other individuals. Of course, not all loans are created equal. Someone who is a essentially a sure bet to pay back a loan will have an easier time getting a loan with a low interest rate than someone who appears to be riskier. And for people who are very risky? They may not even get a loan offer, or they may not have accepted the loan offer due to a high interest rate. It is important to keep that last part in mind, since this data set only represents loans actually made, i.e. do not mistake this data for loan applications!

For the full data dictionary see [here](https://www.openintro.org/data/index.php?data=loan50). We have removed some of the columns to make the data set more reasonably sized and also droped any rows with missing values.
:::

::: {.xsmall}
```{python}
loans = pd.read_csv("data/openintro_loans.csv")
loans
print(loans.columns)
```
:::


```{python echo=FALSE, out.width="66%"}
sns.pairplot(data = loans[["loan_amount","homeownership", "annual_income", "debt_to_income", "interest_rate", "public_record_bankrupt"]], hue="homeownership", corner=True)
```


## OLS

::: {.small}
```{python error=TRUE}
y = loans["loan_amount"]
X = loans[["homeownership", "annual_income", "debt_to_income", "interest_rate", "public_record_bankrupt"]]

model = sm.OLS(endog=y, exog=X)
```
:::

## What do you think the issue is here?

. . .

The error occurs because `X` contains mixed types - specifically we have categorical data columns which cannot be directly converted to a numeric dtype so we need to take care of the dummy coding for statsmodels (with this interface).

. . .

::: {.xsmall}
```{python}
#| error: true
X_dc = pd.get_dummies(X)
model = sm.OLS(endog=y, exog=X_dc)
```

```{python}
X_dc.dtypes
```
:::

. . .

::: {.xsmall}
```{python}
X_dc = pd.get_dummies(X, dtype='int')
model = sm.OLS(endog=y, exog=X_dc)
```
:::

##

::: {.xsmall}
```{python}
model
np.array(dir(model))
```
:::


## Fitting and summary

::: {.small}
```{python}
res = model.fit()
print(res.summary())
```
:::


## Formula interface

::: {.small}
Most of the modeling interfaces are also provided by `smf` (`statsmodels.formula.api`), in which case `patsy` is used to construct the model matrices.
:::

::: {.tiny}
```{python}
model = smf.ols(
  "loan_amount ~ homeownership + annual_income + debt_to_income + interest_rate + public_record_bankrupt",
  data = loans  
)
res = model.fit()
print(res.summary())
```
:::


## Result values and model parameters

:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
res.params
res.bse
```
:::

::: {.column width='50%'}
```{python}
res.rsquared
res.aic
res.bic
res.predict()
```
:::
::::

## Diagnostic plots

:::: {.columns .xsmall}
::: {.column width='50%'}
*QQ Plot*
```{python}
plt.figure()
sm.graphics.qqplot(res.resid, line="s")
plt.show()
```
:::

::: {.column width='50%'}
*Leverage plot*
```{python}
plt.figure()
sm.graphics.plot_leverage_resid2(res)
plt.show()
```
:::
::::


## Alternative model

::: {.tiny}
```{python}
res = smf.ols(
  "np.sqrt(loan_amount) ~ homeownership + annual_income + debt_to_income + interest_rate + public_record_bankrupt",
  data = loans  
).fit()
print(res.summary())
```
:::

##

:::: {.columns .xsmall}
::: {.column width='50%'}
```{python}
#| out-width: 80%
plt.figure()
sm.graphics.qqplot(res.resid, line="s")
plt.show()
```
:::

::: {.column width='50%'}
```{python}
#| out-width: 80%
plt.figure()
sm.graphics.plot_leverage_resid2(res)
plt.show()
```
:::
::::


## Bushtail Possums

::: {.small}
> Data representing possums in Australia and New Guinea. This is a copy of the data set by the same name in the DAAG package, however, the data set included here includes fewer variables.
:::

:::: {.columns .xsmall}
::: {.column width='66%'}

```{python}
possum = pd.read_csv("data/possum.csv")
possum
```
:::

::: {.column width='33%'}
```{r echo=FALSE, out.width="75%"}
knitr::include_graphics("imgs/possum.jpg")
```
:::
::::


## Logistic regression models (GLM)

::: {.small}
```{python error=TRUE}
y = pd.get_dummies( possum["pop"], drop_first = True, dtype="int")
X = pd.get_dummies( possum.drop(["site","pop"], axis=1), dtype="int")

model = sm.GLM(y, X, family = sm.families.Binomial())
```
:::

. . .

::: {.center}
What went wrong this time?
:::


## Missing values

Missing values can be handled via `missing` argument, possible values are `"none"`, `"drop"`, and `"raise"`. 

::: {.xsmall}
```{python}
#| eval: false
model = sm.GLM(y, X, family = sm.families.Binomial(), missing="drop")
res = model.fit()
print(res.summary())
```
:::


## Success vs failure

Note `endog` can be 1d or 2d for binomial models - in the case of the latter each row is interpreted as [success, failure].

::: {.tiny}
```{python error=TRUE}
y = pd.get_dummies( possum["pop"], dtype="int")
X = pd.get_dummies( possum.drop(["site","pop"], axis=1), dtype="int")

res = sm.GLM(y, X, family = sm.families.Binomial(), missing="drop").fit()
print(res.summary())
```
:::


## Formula interface

::: {.tiny}
```{python}
res = smf.glm(
  "pop ~ sex + age + head_l + skull_w + total_l + tail_l-1",
  data = possum, 
  family = sm.families.Binomial(), 
  missing="drop"
).fit()
print(res.summary())
```
:::


## sleepstudy data

::: {.small}
> These data are from the study described in Belenky et al. (2003), for the most sleep-deprived group (3 hours time-in-bed) and for the first 10 days of the study, up to the recovery period. The original study analyzed speed (1/(reaction time)) and treated day as a categorical rather than a continuous predictor.

::: {.small}
```{python}
sleep = pd.read_csv("data/sleepstudy.csv")
sleep
```
:::
:::


::: {.aside}
These data come from the `sleepstudy` dataset in the `lme4` R package
:::

##

::: {.xsmall}
```{python}
g = sns.relplot(x="Days", y="Reaction", col="Subject", col_wrap=6, data=sleep, height=2)
```
:::

## Random intercept model

::: {.small}
```{python}
me_rand_int = smf.mixedlm(
  "Reaction ~ Days", data=sleep, groups=sleep["Subject"], 
  subset=sleep.Days >= 2
)
res_rand_int = me_rand_int.fit(method=["lbfgs"])
print(res_rand_int.summary())
```
:::


## lme4 version

::: {.small}
```{r}
summary(
  lmer(Reaction ~ Days + (1|Subject), data=sleepstudy)
)
```
:::



## Predictions

```{python echo=FALSE, out.width="90%"}
pred = sleep.assign(pred = res_rand_int.predict())
g = sns.FacetGrid(pred, col="Subject", col_wrap=6, height=2)
g = g.map(sns.scatterplot, "Days", "Reaction")
g = g.map(sns.lineplot, "Days", "pred")
plt.show()
```


::: {.aside}
The prediction is only taking into account the fixed effects here, not the group random effects.
:::


## Recovering random effects for prediction

::: {.small}
```{python}
# Multiply each RE by the random effects design matrix for each group
rex = [ 
  np.dot(
    me_rand_int.exog_re_li[j], 
    res_rand_int.random_effects[k]
  ) 
  for (j, k) in enumerate(me_rand_int.group_labels)
]
rex[0]

# Add the fixed and random terms to get the overall prediction
y_hat = res_rand_int.predict() + np.concatenate(rex)
```
:::

::: {.aside}
Based on code provide on [Stack Overflow](https://stats.stackexchange.com/questions/467543/including-random-effects-in-prediction-with-linear-mixed-model).
:::

##

```{python}
#| echo: false
pred = sleep.assign(pred = y_hat)
g = sns.FacetGrid(pred, col="Subject", col_wrap=6, height=2)
g = g.map(sns.scatterplot, "Days", "Reaction")
g = g.map(sns.lineplot, "Days", "pred")
plt.show()
```

## Random intercept and slope model

::: {.xsmall}
```{python}
me_rand_sl= smf.mixedlm(
  "Reaction ~ Days", data=sleep, groups=sleep["Subject"], 
  subset=sleep.Days >= 2,
  re_formula="~Days" 
)
res_rand_sl = me_rand_sl.fit(method=["lbfgs"])
print(res_rand_sl.summary())
```
:::

## lme4 version

::: {.xxsmall}
```{r}
summary(
  lmer(Reaction ~ Days + (Days|Subject), data=sleepstudy)
)
```
:::


## Prediction

```{python echo=FALSE, out.width="90%"}
# Dictionary of random effects estimates
re = res_rand_sl.random_effects

# Multiply each RE by the random effects design matrix for each group
rex = [me_rand_sl.exog_re_li[j] @ re[k] for (j, k) in enumerate(me_rand_sl.group_labels)]

# Add the fixed and random terms to get the overall prediction
rex = np.concatenate(rex)
y_hat = res_rand_sl.predict() + rex

pred = sleep.assign(pred = y_hat)
g = sns.FacetGrid(pred, col="Subject", col_wrap=6, height=2)
g = g.map(sns.scatterplot, "Days", "Reaction")
g = g.map(sns.lineplot, "Days", "pred")
plt.show()
```

::: {.aside}
We are using the same approach described previously to obtain the RE estimates and use them in the predictions.
:::


