---
title: "patsy + <br/> statsmodels"
subtitle: "Lecture 18"
author: "Dr. Colin Rundel"
footer: "Sta 663 - Spring 2023"
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

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=80,
  precision = 5, suppress=True
  #formatter=dict(float=lambda x: "%.5g" % x)
)

pd.set_option("display.width", 150)
pd.set_option("display.max_columns", 10)
pd.set_option("display.precision", 6)
```

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

# 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

```{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")
```

## Model matrix


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

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

::: {.column width='50%' .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

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

:::: {.columns .small}
::: {.column width='50%'}
```{python}
y
```
:::

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


## as DataFrames

```{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 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 right terms | `a/b` ⇒ `a + a:b` |
| `I()`    | used for calculating arithmetic calculations      | `I(x1 + x2)` |
| `Q()`    | used to quote column names, e.g. columns with spaces or symbols | `Q('bad name!')` |
| `C()`    | used for categorical data coding                  |  `C(a, Treatment('a2'))` |
:::

## Examples

:::: {.columns .small}
::: {.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"))
dmatrix("x*(y+z)", demo_data("x","y","z"))
```
:::
::::

## Intercept Examples (-1)

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

```
:::

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

## Intercept Examples (0)

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

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


## Design Info

One of the keep 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).

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


## Stateful transforms

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

:::: {.columns .small}
::: {.column width='50%' .fragment}
```{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-lego PatsyTransformer

If you would like to use a Patsy formula in a scikitlearn pipeline, it is possible via the `PatsyTransformer` from the scikit-lego library (`sklego`).

::: {.small}
```{python}
from sklego.preprocessing import PatsyTransformer
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 .small}
::: {.column width='50%'}
```{python}
pt = PatsyTransformer("x*a + np.log(x)")
pt.fit_transform(X)
```
:::

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




## B-splines

Patsy also has support for B-splines and other related models,

:::: {.columns .small}
::: {.column width='50%'}
```{python}
d = pd.read_csv("data/d1.csv")
sns.relplot(x="x", y="y", data=d)
```
:::

::: {.column width='50%'}
```{python}
y, X = dmatrices("y ~ bs(x, df=6)", data=d)
X
```
:::
::::


## What is `bs(x)[i]`?

::: {.small}
```{python}
#| out-width: 75%
bs_df = ( 
  dmatrix("bs(x, df=6)", data=d, return_type="dataframe")
  .drop(["Intercept"], axis = 1)
  .assign(x = d["x"])
  .melt(id_vars="x")
)
sns.relplot(x="x", y="value", hue="variable", kind="line", data = bs_df, aspect=1.5)
```
:::


## Fitting a model

::: {.small}
```{python}
from sklearn.linear_model import LinearRegression
lm = LinearRegression(fit_intercept=False).fit(X,y)
lm.coef_
```
:::

. . .

::: {.small}
```{python}
#| out-width: 66%
plt.figure(layout="constrained")
sns.lineplot(x=d["x"], y=lm.predict(X).ravel(), color="r")
sns.scatterplot(x="x", y="y", data=d)
plt.show()
```
:::


## sklearn SplineTransformer

:::: {.columns .small}
::: {.column width='50%'}
```{python}
from sklearn.preprocessing import SplineTransformer

p = make_pipeline(
  SplineTransformer(
    n_knots=6, 
    degree=3, 
    include_bias=True
  ),
  LinearRegression(fit_intercept=False)
).fit(
  d[["x"]], d["y"]
)
```
:::

::: {.column width='50%'}
```{python}
plt.figure()
sns.lineplot(x=d["x"], y=p.predict(d[["x"]]).ravel(), color="k")
sns.scatterplot(x="x", y="y", data=d)
plt.show()
```
:::
::::


## Comparison

::: {.small}
```{python}
plt.figure()
sns.lineplot(x=d["x"], y=p.predict(d[["x"]]).ravel(), color="k", label = "sklearn")
sns.lineplot(x=d["x"], y=lm.predict(X).ravel(), color="r", label = "patsy")
sns.scatterplot(x="x", y="y", data=d)
plt.show()
```
:::



## Why different?

For patsy the number of splines is determined by `df` while for sklearn this is determined by `n_knots + degree - 1`.

::: {.small}
```{python}
#| out-width: 75%
p = p.set_params(splinetransformer__n_knots = 5).fit(d[["x"]], d["y"])

plt.figure(layout="constrained")
sns.lineplot(x=d["x"], y=p.predict(d[["x"]]).ravel(), color="k", label = "sklearn")
sns.lineplot(x=d["x"], y=lm.predict(X).ravel(), color="r", label = "patsy")
sns.scatterplot(x="x", y="y", data=d)
plt.show()
```
:::


##

but that is not the whole story, if we examine the bases we also see they differ slightly between implementations

::: {.small}
```{python}
#| out-width: 50%
bs_df = pd.DataFrame(
  SplineTransformer(n_knots=6, degree=3, include_bias=True).fit_transform(d[["x"]]),
  columns = ["bs["+ str(i) +"]" for i in range(8)]
).assign(
  x = d.x
).melt(
  id_vars = "x"
)
sns.relplot(x="x", y="value", hue="variable", kind="line", data = bs_df, aspect=1.5)
```
:::


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

<br/>

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

. . .

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

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

::: {.small}
```{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).


::: {.small}
```{python}
X_dc = pd.get_dummies(X)
model = sm.OLS(endog=y, exog=X_dc)
model
np.array(dir(model))
```
:::


## Fitting and summary

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


## Formula interface

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

::: {.small}
```{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 .small}
::: {.column width='50%'}
```{python}
res.params
res.bse
```
:::

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

## Diagnostic plots

:::: {.columns .small}
::: {.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

::: {.small}
```{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 .small}
::: {.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

> 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.
>
> `pop` - Population, either `Vic` (Victoria) or `other` (New South Wales or Queensland).


:::: {.columns .small}
::: {.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)

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

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

. . .

What is wrong now?

. . .

Behavior for dealing with missing data can be handled via `missing`, possible values are `"none"`, `"drop"`, and `"raise"`. 

```{python}
model = sm.GLM(y, X, family = sm.families.Binomial(), missing="drop")
```

## Fit and summary

::: {.small}
```{python}
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].

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

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


## Formula interface

::: {.small}
```{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.

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


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

##

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

## 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)
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)
g = g.map(sns.scatterplot, "Days", "Reaction")
g = g.map(sns.lineplot, "Days", "pred")
plt.show()
```

## Random intercept and slope model

::: {.small}
```{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

::: {.small}
```{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)
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.
:::

# Odds and ends

## t-test and z-test for equality of means

::: {.small}
```{python}
books = pd.read_csv("data/daag_books.csv")
cm = sm.stats.CompareMeans(
  sm.stats.DescrStatsW( books.weight[books.cover == "hb"] ),
  sm.stats.DescrStatsW( books.weight[books.cover == "pb"] )
)
```

```{python}
print(cm.summary())
print(cm.summary(use_t=False))
```
:::

##

::: {.small}
```{python}
print(cm.summary(usevar="unequal"))
```
:::


## Contigency tables

::: {.small}
Below are data from the GSS and a survery of Duke students in a intro stats class - the question asked about how concerned the respondent was about the effect of global warming on polar ice cap melt.
:::

:::: {.columns .small}
::: {.column width='50%'}
```{python}
gss = pd.DataFrame({"US": [454, 226], 
                    "Duke": [56,32]}, 
                   index=["A great deal", "Not a great deal"])
gss

```
:::

::: {.column width='50%'}
```{python}
tbl = sm.stats.Table2x2(gss.to_numpy())
print(tbl)
```
:::
::::

. . .

::: {.small}
```{python}
print(tbl.summary())
print(tbl.test_nominal_association()) # chi^2 test of independence
```
:::

