docs/tutorials/formula_interface/formula_interface.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Formula Interface Tutorial: Revisiting French Motor Third-Party Liability Claims\n",
    "\n",
    "\n",
    "**Intro**\n",
    "\n",
    "This tutorial showcases the formula interface of `glum`. It allows for the specification of the design matrix and the response variable using so-called [Wilkinson-formulas](https://www.jstor.org/stable/2346786) instead of constructing it by hand. This kind of model specification should be familiar to R users or those who have used the `statsmodels` or `linearmodels` Python packages before. This tutorial aims to introduce the basics of working with formulas to other users, as well as highlighting some important differences between `glum`s and other packages' formula implementations.\n",
    "\n",
    "For a more in-depth look at how formulas work, please take a look at the [documentation of `formulaic`](https://matthewwardrop.github.io/formulaic/), the package on which `glum`'s formula interface is based.\n",
    "\n",
    "\n",
    "**Background**\n",
    "\n",
    "This tutorial reimplements and extends the combined frequency-severity model from Chapter 4 of the [GLM tutorial](tutorials/glm_french_motor_tutorial/glm_french_motor.html). If you would like to know more about the setting, the data, or GLM modeling in general, please check that out first.\n",
    "\n",
    "**Sneak Peak**\n",
    "\n",
    "Formulas can provide a concise and convenient way to specify many of the usual pre-processing steps, such as converting to categorical types, creating interactions, applying transformations, or even spline interpolation. As an example, consider the following formula:\n",
    "\n",
    "```\n",
    "{ClaimAmountCut / Exposure} ~ C(DrivAge, missing_method='convert') * C(VehPower, missing_method=\"zero\") + bs(BonusMalus, 3)\n",
    "```\n",
    "\n",
    "Despite its brevity, it describes all of the following:\n",
    " - The outcome variable is the ratio of `ClaimAmountCut` and `Exposure`.\n",
    " - The predictors should include the interactions of the categorical variables `DrivAge` and `VehPower`, as well as those two variables themselves. (Even though they behave as such, neither the individual variables nor their interaction will be dummy-encoded by glum. For categoricals with many levels, this can lead to a substantial performance improvement over dummy encoding, especially for the interaction.)\n",
    " - If there are missing values in `DrivAge`, they should be treated as a separate category.\n",
    " - On the other hand, missing values in `VehPower` should be treated as all-zero indicators.\n",
    " - The predictors should also include a third degree B-spline interpolation of `BonusMalus`.\n",
    "\n",
    "The following chapters demonstrate each of these features in some detail, as well as some additional advantages of using the formula interface."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Contents\n",
    "* [1. Load and Prepare Datasets from Openml](#1.-Load-and-Prepare-Datasets-from-Openml)\n",
    "* [2. Reproducing the model from the GLM Tutorial](#2.-Reproducing-the-model-from-the-GLM-Tutorial)\n",
    "* [3. Categorical Variables](#3.-Categorical-Variables)\n",
    "* [4. Interactions and Structural Full-rankness](#4.-Interactions-and-Structural-Full-rankness)\n",
    "* [5. Fun with Functions](#5.-Fun-with-Functions)\n",
    "* [6. Miscellaneous Features](#6.-Miscellaneous-Features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pytest\n",
    "import scipy.optimize as optimize\n",
    "import scipy.stats\n",
    "from dask_ml.preprocessing import Categorizer\n",
    "from sklearn.metrics import mean_absolute_error\n",
    "from sklearn.model_selection import ShuffleSplit\n",
    "from glum import GeneralizedLinearRegressor\n",
    "from glum import TweedieDistribution\n",
    "\n",
    "from load_transform_formula import load_transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Load and Prepare Datasets from Openml<a class=\"anchor\"></a>\n",
    "[back to table of contents](#Table-of-Contents)\n",
    "\n",
    "First, we load in our [dataset from openML](\"https://www.openml.org/d/41214\") and apply several transformations. In the interest of simplicity, we do not include the data loading and preparation code in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ClaimNb</th>\n",
       "      <th>Exposure</th>\n",
       "      <th>Area</th>\n",
       "      <th>VehPower</th>\n",
       "      <th>VehAge</th>\n",
       "      <th>DrivAge</th>\n",
       "      <th>BonusMalus</th>\n",
       "      <th>VehBrand</th>\n",
       "      <th>VehGas</th>\n",
       "      <th>Density</th>\n",
       "      <th>Region</th>\n",
       "      <th>ClaimAmount</th>\n",
       "      <th>ClaimAmountCut</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IDpol</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>0.10000</td>\n",
       "      <td>D</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>50</td>\n",
       "      <td>B12</td>\n",
       "      <td>Regular</td>\n",
       "      <td>1217</td>\n",
       "      <td>R82</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>0.77000</td>\n",
       "      <td>D</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>50</td>\n",
       "      <td>B12</td>\n",
       "      <td>Regular</td>\n",
       "      <td>1217</td>\n",
       "      <td>R82</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0</td>\n",
       "      <td>0.75000</td>\n",
       "      <td>B</td>\n",
       "      <td>6</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>50</td>\n",
       "      <td>B12</td>\n",
       "      <td>Diesel</td>\n",
       "      <td>54</td>\n",
       "      <td>R22</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0</td>\n",
       "      <td>0.09000</td>\n",
       "      <td>B</td>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>50</td>\n",
       "      <td>B12</td>\n",
       "      <td>Diesel</td>\n",
       "      <td>76</td>\n",
       "      <td>R72</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>0</td>\n",
       "      <td>0.84000</td>\n",
       "      <td>B</td>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>50</td>\n",
       "      <td>B12</td>\n",
       "      <td>Diesel</td>\n",
       "      <td>76</td>\n",
       "      <td>R72</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6114326</th>\n",
       "      <td>0</td>\n",
       "      <td>0.00274</td>\n",
       "      <td>E</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>50</td>\n",
       "      <td>B12</td>\n",
       "      <td>Regular</td>\n",
       "      <td>3317</td>\n",
       "      <td>R93</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6114327</th>\n",
       "      <td>0</td>\n",
       "      <td>0.00274</td>\n",
       "      <td>E</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>95</td>\n",
       "      <td>B12</td>\n",
       "      <td>Regular</td>\n",
       "      <td>9850</td>\n",
       "      <td>R11</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6114328</th>\n",
       "      <td>0</td>\n",
       "      <td>0.00274</td>\n",
       "      <td>D</td>\n",
       "      <td>6</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>50</td>\n",
       "      <td>B12</td>\n",
       "      <td>Diesel</td>\n",
       "      <td>1323</td>\n",
       "      <td>R82</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6114329</th>\n",
       "      <td>0</td>\n",
       "      <td>0.00274</td>\n",
       "      <td>B</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>50</td>\n",
       "      <td>B12</td>\n",
       "      <td>Regular</td>\n",
       "      <td>95</td>\n",
       "      <td>R26</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6114330</th>\n",
       "      <td>0</td>\n",
       "      <td>0.00274</td>\n",
       "      <td>B</td>\n",
       "      <td>7</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>54</td>\n",
       "      <td>B12</td>\n",
       "      <td>Diesel</td>\n",
       "      <td>65</td>\n",
       "      <td>R72</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>678013 rows × 13 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         ClaimNb  Exposure Area  VehPower  VehAge  DrivAge  BonusMalus  \\\n",
       "IDpol                                                                    \n",
       "1              0   0.10000    D         5       0        5          50   \n",
       "3              0   0.77000    D         5       0        5          50   \n",
       "5              0   0.75000    B         6       1        5          50   \n",
       "10             0   0.09000    B         7       0        4          50   \n",
       "11             0   0.84000    B         7       0        4          50   \n",
       "...          ...       ...  ...       ...     ...      ...         ...   \n",
       "6114326        0   0.00274    E         4       0        5          50   \n",
       "6114327        0   0.00274    E         4       0        4          95   \n",
       "6114328        0   0.00274    D         6       1        4          50   \n",
       "6114329        0   0.00274    B         4       0        5          50   \n",
       "6114330        0   0.00274    B         7       1        2          54   \n",
       "\n",
       "        VehBrand   VehGas  Density Region  ClaimAmount  ClaimAmountCut  \n",
       "IDpol                                                                   \n",
       "1            B12  Regular     1217    R82          0.0             0.0  \n",
       "3            B12  Regular     1217    R82          0.0             0.0  \n",
       "5            B12   Diesel       54    R22          0.0             0.0  \n",
       "10           B12   Diesel       76    R72          0.0             0.0  \n",
       "11           B12   Diesel       76    R72          0.0             0.0  \n",
       "...          ...      ...      ...    ...          ...             ...  \n",
       "6114326      B12  Regular     3317    R93          0.0             0.0  \n",
       "6114327      B12  Regular     9850    R11          0.0             0.0  \n",
       "6114328      B12   Diesel     1323    R82          0.0             0.0  \n",
       "6114329      B12  Regular       95    R26          0.0             0.0  \n",
       "6114330      B12   Diesel       65    R72          0.0             0.0  \n",
       "\n",
       "[678013 rows x 13 columns]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = load_transform()\n",
    "with pd.option_context('display.max_rows', 10):\n",
    "    display(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Reproducing the Model From the GLM Turorial<a class=\"anchor\"></a>\n",
    "\n",
    "Now, let us start by fitting a very simple model. As usual, let's divide our samples into a training and a test set so that we get valid out-of-sample goodness-of-fit measures. Perhaps less usually, we do not create separate `y` and `X` data frames for our label and features – the formula will take care of that for us.\n",
    "\n",
    "We still have some preprocessing to do:\n",
    " - Many of the ordinal or nominal variables are encoded as integers, instead of as categoricals. We will need to convert these so that `glum` will know to estimate a separate coefficient for each of their levels.\n",
    " - The outcome variable is a transformation of other columns. We need to create it first.\n",
    "\n",
    "As we will see later on, these steps can be incorporated into the formula itself, but let's not overcomplicate things at first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "ss = ShuffleSplit(n_splits=1, test_size=0.1, random_state=42)\n",
    "train, test = next(ss.split(df))\n",
    "\n",
    "df = df.assign(PurePremium=lambda x: x[\"ClaimAmountCut\"] / x[\"Exposure\"])\n",
    "\n",
    "glm_categorizer = Categorizer(\n",
    "    columns=[\"VehBrand\", \"VehGas\", \"Region\", \"Area\", \"DrivAge\", \"VehAge\", \"VehPower\"]\n",
    ")\n",
    "df_train = glm_categorizer.fit_transform(df.iloc[train])\n",
    "df_test = glm_categorizer.transform(df.iloc[test])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example demonstrates the basic idea behind formulas: the outcome variable and the predictors are separated by a tilde (`~`), and different prefictors are separated by plus signs (`+`). Thus, formulas provide a concise way of specifying a model without the need to create dataframes by hand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>intercept</th>\n",
       "      <th>VehBrand[B1]</th>\n",
       "      <th>VehBrand[B10]</th>\n",
       "      <th>VehBrand[B11]</th>\n",
       "      <th>VehBrand[B12]</th>\n",
       "      <th>VehBrand[B13]</th>\n",
       "      <th>VehBrand[B14]</th>\n",
       "      <th>VehBrand[B2]</th>\n",
       "      <th>VehBrand[B3]</th>\n",
       "      <th>VehBrand[B4]</th>\n",
       "      <th>...</th>\n",
       "      <th>VehAge[1]</th>\n",
       "      <th>VehAge[2]</th>\n",
       "      <th>VehPower[4]</th>\n",
       "      <th>VehPower[5]</th>\n",
       "      <th>VehPower[6]</th>\n",
       "      <th>VehPower[7]</th>\n",
       "      <th>VehPower[8]</th>\n",
       "      <th>VehPower[9]</th>\n",
       "      <th>BonusMalus</th>\n",
       "      <th>Density</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>coefficient</th>\n",
       "      <td>2.88667</td>\n",
       "      <td>-0.064157</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.231868</td>\n",
       "      <td>-0.211061</td>\n",
       "      <td>0.054979</td>\n",
       "      <td>-0.270346</td>\n",
       "      <td>-0.071453</td>\n",
       "      <td>0.00291</td>\n",
       "      <td>0.059324</td>\n",
       "      <td>...</td>\n",
       "      <td>0.008117</td>\n",
       "      <td>-0.229906</td>\n",
       "      <td>-0.111796</td>\n",
       "      <td>-0.123388</td>\n",
       "      <td>0.060757</td>\n",
       "      <td>0.005179</td>\n",
       "      <td>-0.021832</td>\n",
       "      <td>0.208158</td>\n",
       "      <td>0.032508</td>\n",
       "      <td>0.000002</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 60 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "             intercept  VehBrand[B1]  VehBrand[B10]  VehBrand[B11]  \\\n",
       "coefficient    2.88667     -0.064157            0.0       0.231868   \n",
       "\n",
       "             VehBrand[B12]  VehBrand[B13]  VehBrand[B14]  VehBrand[B2]  \\\n",
       "coefficient      -0.211061       0.054979      -0.270346     -0.071453   \n",
       "\n",
       "             VehBrand[B3]  VehBrand[B4]  ...  VehAge[1]  VehAge[2]  \\\n",
       "coefficient       0.00291      0.059324  ...   0.008117  -0.229906   \n",
       "\n",
       "             VehPower[4]  VehPower[5]  VehPower[6]  VehPower[7]  VehPower[8]  \\\n",
       "coefficient    -0.111796    -0.123388     0.060757     0.005179    -0.021832   \n",
       "\n",
       "             VehPower[9]  BonusMalus   Density  \n",
       "coefficient     0.208158    0.032508  0.000002  \n",
       "\n",
       "[1 rows x 60 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula = \"PurePremium ~ VehBrand + VehGas + Region + Area + DrivAge + VehAge + VehPower + BonusMalus + Density\"\n",
    "\n",
    "TweedieDist = TweedieDistribution(1.5)\n",
    "t_glm1 = GeneralizedLinearRegressor(\n",
    "    family=TweedieDist,\n",
    "    alpha_search=True,\n",
    "    l1_ratio=1,\n",
    "    fit_intercept=True,\n",
    "    formula=formula,\n",
    ")\n",
    "t_glm1.fit(df_train, sample_weight=df[\"Exposure\"].values[train])\n",
    "\n",
    "pd.DataFrame(\n",
    "    {\"coefficient\": np.concatenate(([t_glm1.intercept_], t_glm1.coef_))},\n",
    "    index=[\"intercept\"] + t_glm1.feature_names_,\n",
    ").T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Categorical Variables<a class=\"anchor\"></a>\n",
    "\n",
    "`glum` also provides extensive support for categorical variables. The main function one needs to be aware of in the context of categoricals is simply called `C()`. A variable placed within it is always converted to a categorical, regardless of its type.\n",
    "\n",
    "A huge part of tabmat's/glum's performance advantage is that categoricals need not be one-hot encoded, but are treated as if they were. For this reason, we do not support using other coding schemes within the formula interface. If one needs to use other categorical encodings than one-hot, they can always do so manually (or even using `formulaic` directly) before the estimation.\n",
    "\n",
    "Let's try it out on our dataset!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ClaimNb             int64\n",
       "Exposure          float64\n",
       "Area               object\n",
       "VehPower            int64\n",
       "VehAge              int64\n",
       "DrivAge             int64\n",
       "BonusMalus          int64\n",
       "VehBrand           object\n",
       "VehGas             object\n",
       "Density             int64\n",
       "Region             object\n",
       "ClaimAmount       float64\n",
       "ClaimAmountCut    float64\n",
       "PurePremium       float64\n",
       "dtype: object"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_train_noncat = df.iloc[train]\n",
    "df_test_noncat = df.iloc[test]\n",
    "\n",
    "df_train_noncat.dtypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even though some of the variables are integers in this dataset, they are handled as categoricals thanks to the `C()` function. Strings, such as `VehBrand` or `VehGas` would have been handled as categorical by default anyway, but using the `C()` function never hurts: if applied to something that is already a caetgorical variable, it does not have any effect outside of the feature name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>intercept</th>\n",
       "      <th>C(VehBrand)[B1]</th>\n",
       "      <th>C(VehBrand)[B10]</th>\n",
       "      <th>C(VehBrand)[B11]</th>\n",
       "      <th>C(VehBrand)[B12]</th>\n",
       "      <th>C(VehBrand)[B13]</th>\n",
       "      <th>C(VehBrand)[B14]</th>\n",
       "      <th>C(VehBrand)[B2]</th>\n",
       "      <th>C(VehBrand)[B3]</th>\n",
       "      <th>C(VehBrand)[B4]</th>\n",
       "      <th>...</th>\n",
       "      <th>C(VehAge)[1]</th>\n",
       "      <th>C(VehAge)[2]</th>\n",
       "      <th>C(VehPower)[4]</th>\n",
       "      <th>C(VehPower)[5]</th>\n",
       "      <th>C(VehPower)[6]</th>\n",
       "      <th>C(VehPower)[7]</th>\n",
       "      <th>C(VehPower)[8]</th>\n",
       "      <th>C(VehPower)[9]</th>\n",
       "      <th>BonusMalus</th>\n",
       "      <th>Density</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>coefficient</th>\n",
       "      <td>2.88667</td>\n",
       "      <td>-0.064157</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.231868</td>\n",
       "      <td>-0.211061</td>\n",
       "      <td>0.054979</td>\n",
       "      <td>-0.270346</td>\n",
       "      <td>-0.071453</td>\n",
       "      <td>0.00291</td>\n",
       "      <td>0.059324</td>\n",
       "      <td>...</td>\n",
       "      <td>0.008117</td>\n",
       "      <td>-0.229906</td>\n",
       "      <td>-0.111796</td>\n",
       "      <td>-0.123388</td>\n",
       "      <td>0.060757</td>\n",
       "      <td>0.005179</td>\n",
       "      <td>-0.021832</td>\n",
       "      <td>0.208158</td>\n",
       "      <td>0.032508</td>\n",
       "      <td>0.000002</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 60 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "             intercept  C(VehBrand)[B1]  C(VehBrand)[B10]  C(VehBrand)[B11]  \\\n",
       "coefficient    2.88667        -0.064157               0.0          0.231868   \n",
       "\n",
       "             C(VehBrand)[B12]  C(VehBrand)[B13]  C(VehBrand)[B14]  \\\n",
       "coefficient         -0.211061          0.054979         -0.270346   \n",
       "\n",
       "             C(VehBrand)[B2]  C(VehBrand)[B3]  C(VehBrand)[B4]  ...  \\\n",
       "coefficient        -0.071453          0.00291         0.059324  ...   \n",
       "\n",
       "             C(VehAge)[1]  C(VehAge)[2]  C(VehPower)[4]  C(VehPower)[5]  \\\n",
       "coefficient      0.008117     -0.229906       -0.111796       -0.123388   \n",
       "\n",
       "             C(VehPower)[6]  C(VehPower)[7]  C(VehPower)[8]  C(VehPower)[9]  \\\n",
       "coefficient        0.060757        0.005179       -0.021832        0.208158   \n",
       "\n",
       "             BonusMalus   Density  \n",
       "coefficient    0.032508  0.000002  \n",
       "\n",
       "[1 rows x 60 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula_cat = (\n",
    "    \"PurePremium ~ C(VehBrand) + C(VehGas) + C(Region) + C(Area) \"\n",
    "    \"+ C(DrivAge) + C(VehAge) + C(VehPower) + BonusMalus + Density\"\n",
    ")\n",
    "\n",
    "t_glm3 = GeneralizedLinearRegressor(\n",
    "    family=TweedieDist,\n",
    "    alpha_search=True,\n",
    "    l1_ratio=1,\n",
    "    fit_intercept=True,\n",
    "    formula=formula_cat,\n",
    ")\n",
    "t_glm3.fit(df_train_noncat, sample_weight=df[\"Exposure\"].values[train])\n",
    "\n",
    "pd.DataFrame(\n",
    "    {\"coefficient\": np.concatenate(([t_glm3.intercept_], t_glm3.coef_))},\n",
    "    index=[\"intercept\"] + t_glm3.feature_names_,\n",
    ").T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, prediction works as expected with categorical variables. `glum` keeps track of the levels present in the training dataset, and makes sure that categorical variables in unseen datasets are also properly aligned, even if they have missing or unknown levels.<sup>3</sup> Therefore, one can simply use predict, and `glum` does The Right Thing™ by default.\n",
    "\n",
    "<sup>3</sup>: This is made possible due to `glum` saving a [`ModelSpec` object](https://matthewwardrop.github.io/formulaic/guides/model_specs/), which contains any information necessary for reapplying the transitions that were done during the formula materialization process. It is especially relevant in the case of [stateful transforms](https://matthewwardrop.github.io/formulaic/guides/transforms/), such as creating categorical variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([303.77443311, 548.47789523, 244.34438579, ..., 109.81572865,\n",
       "        67.98332028, 297.21717383])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_glm3.predict(df_test_noncat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Interactions and Structural Full-Rankness<a class=\"anchor\"></a>\n",
    "\n",
    "One of the biggest strengths of Wilkinson-formuals lie in their ability of concisely specifying interactions between terms. `glum` implements this as well, and in a very efficient way: the interactions of categorical features are encoded as a new categorical feature, making it possible to interact high-cardinality categoricals with each other. If this is not possible, because, for example, a categorical is interacted with a numeric variable, sparse representations are used when appropriate. In general, just as with `glum`'s categorical handling in general, you can be assured that `glum` you don't have to worry too much about the actual implementation, and can expect that `glum` will do the most efficient thing behind the scenes.\n",
    "\n",
    "Let's see how that looks like on the insurance example! Suppose that we expect `VehPower` to have a different effect depending on `DrivAge` (e.g. performance cars might not be great for new drivers, but may be less problematic for more experienced ones). We can include the interaction of these variables as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>intercept</th>\n",
       "      <th>C(VehBrand)[B1]</th>\n",
       "      <th>C(VehBrand)[B10]</th>\n",
       "      <th>C(VehBrand)[B11]</th>\n",
       "      <th>C(VehBrand)[B12]</th>\n",
       "      <th>C(VehBrand)[B13]</th>\n",
       "      <th>C(VehBrand)[B14]</th>\n",
       "      <th>C(VehBrand)[B2]</th>\n",
       "      <th>C(VehBrand)[B3]</th>\n",
       "      <th>C(VehBrand)[B4]</th>\n",
       "      <th>...</th>\n",
       "      <th>C(DrivAge)[4]:C(VehPower)[8]</th>\n",
       "      <th>C(DrivAge)[5]:C(VehPower)[8]</th>\n",
       "      <th>C(DrivAge)[6]:C(VehPower)[8]</th>\n",
       "      <th>C(DrivAge)[0]:C(VehPower)[9]</th>\n",
       "      <th>C(DrivAge)[1]:C(VehPower)[9]</th>\n",
       "      <th>C(DrivAge)[2]:C(VehPower)[9]</th>\n",
       "      <th>C(DrivAge)[3]:C(VehPower)[9]</th>\n",
       "      <th>C(DrivAge)[4]:C(VehPower)[9]</th>\n",
       "      <th>C(DrivAge)[5]:C(VehPower)[9]</th>\n",
       "      <th>C(DrivAge)[6]:C(VehPower)[9]</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>coefficient</th>\n",
       "      <td>2.88023</td>\n",
       "      <td>-0.069076</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.221037</td>\n",
       "      <td>-0.211854</td>\n",
       "      <td>0.052355</td>\n",
       "      <td>-0.272058</td>\n",
       "      <td>-0.074836</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.052523</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.147844</td>\n",
       "      <td>-0.03567</td>\n",
       "      <td>0.504407</td>\n",
       "      <td>0.682528</td>\n",
       "      <td>-0.106569</td>\n",
       "      <td>-0.308257</td>\n",
       "      <td>0.173206</td>\n",
       "      <td>0.010684</td>\n",
       "      <td>-0.220273</td>\n",
       "      <td>0.070334</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 102 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "             intercept  C(VehBrand)[B1]  C(VehBrand)[B10]  C(VehBrand)[B11]  \\\n",
       "coefficient    2.88023        -0.069076               0.0          0.221037   \n",
       "\n",
       "             C(VehBrand)[B12]  C(VehBrand)[B13]  C(VehBrand)[B14]  \\\n",
       "coefficient         -0.211854          0.052355         -0.272058   \n",
       "\n",
       "             C(VehBrand)[B2]  C(VehBrand)[B3]  C(VehBrand)[B4]  ...  \\\n",
       "coefficient        -0.074836              0.0         0.052523  ...   \n",
       "\n",
       "             C(DrivAge)[4]:C(VehPower)[8]  C(DrivAge)[5]:C(VehPower)[8]  \\\n",
       "coefficient                     -0.147844                      -0.03567   \n",
       "\n",
       "             C(DrivAge)[6]:C(VehPower)[8]  C(DrivAge)[0]:C(VehPower)[9]  \\\n",
       "coefficient                      0.504407                      0.682528   \n",
       "\n",
       "             C(DrivAge)[1]:C(VehPower)[9]  C(DrivAge)[2]:C(VehPower)[9]  \\\n",
       "coefficient                     -0.106569                     -0.308257   \n",
       "\n",
       "             C(DrivAge)[3]:C(VehPower)[9]  C(DrivAge)[4]:C(VehPower)[9]  \\\n",
       "coefficient                      0.173206                      0.010684   \n",
       "\n",
       "             C(DrivAge)[5]:C(VehPower)[9]  C(DrivAge)[6]:C(VehPower)[9]  \n",
       "coefficient                     -0.220273                      0.070334  \n",
       "\n",
       "[1 rows x 102 columns]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula_int = (\n",
    "    \"PurePremium ~ C(VehBrand) + C(VehGas) + C(Region) + C(Area)\"\n",
    "    \" + C(DrivAge) * C(VehPower) + C(VehAge) + BonusMalus + Density\"\n",
    ")\n",
    "\n",
    "t_glm4 = GeneralizedLinearRegressor(\n",
    "    family=TweedieDist,\n",
    "    alpha_search=True,\n",
    "    l1_ratio=1,\n",
    "    fit_intercept=True,\n",
    "    formula=formula_int,\n",
    ")\n",
    "t_glm4.fit(df_train, sample_weight=df[\"Exposure\"].values[train])\n",
    "\n",
    "pd.DataFrame(\n",
    "    {\"coefficient\": np.concatenate(([t_glm4.intercept_], t_glm4.coef_))},\n",
    "    index=[\"intercept\"] + t_glm4.feature_names_,\n",
    ").T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that, in addition to the interactions, the non-interacted variants of `DrivAge` and `VehPower` are also included in the model. This is a result of using the `*` operator to interact the variables. Using `:` instead would only include the interactions, and not the marginals. (In short, `a * b` is equivalent to `a + b + a:b`.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['C(VehPower)[4]',\n",
       " 'C(VehPower)[5]',\n",
       " 'C(VehPower)[6]',\n",
       " 'C(VehPower)[7]',\n",
       " 'C(VehPower)[8]',\n",
       " 'C(VehPower)[9]',\n",
       " 'C(DrivAge)[0]:C(VehPower)[4]',\n",
       " 'C(DrivAge)[1]:C(VehPower)[4]',\n",
       " 'C(DrivAge)[2]:C(VehPower)[4]',\n",
       " 'C(DrivAge)[3]:C(VehPower)[4]',\n",
       " 'C(DrivAge)[4]:C(VehPower)[4]',\n",
       " 'C(DrivAge)[5]:C(VehPower)[4]',\n",
       " 'C(DrivAge)[6]:C(VehPower)[4]',\n",
       " 'C(DrivAge)[0]:C(VehPower)[5]',\n",
       " 'C(DrivAge)[1]:C(VehPower)[5]',\n",
       " 'C(DrivAge)[2]:C(VehPower)[5]',\n",
       " 'C(DrivAge)[3]:C(VehPower)[5]',\n",
       " 'C(DrivAge)[4]:C(VehPower)[5]',\n",
       " 'C(DrivAge)[5]:C(VehPower)[5]',\n",
       " 'C(DrivAge)[6]:C(VehPower)[5]',\n",
       " 'C(DrivAge)[0]:C(VehPower)[6]',\n",
       " 'C(DrivAge)[1]:C(VehPower)[6]',\n",
       " 'C(DrivAge)[2]:C(VehPower)[6]',\n",
       " 'C(DrivAge)[3]:C(VehPower)[6]',\n",
       " 'C(DrivAge)[4]:C(VehPower)[6]',\n",
       " 'C(DrivAge)[5]:C(VehPower)[6]',\n",
       " 'C(DrivAge)[6]:C(VehPower)[6]',\n",
       " 'C(DrivAge)[0]:C(VehPower)[7]',\n",
       " 'C(DrivAge)[1]:C(VehPower)[7]',\n",
       " 'C(DrivAge)[2]:C(VehPower)[7]',\n",
       " 'C(DrivAge)[3]:C(VehPower)[7]',\n",
       " 'C(DrivAge)[4]:C(VehPower)[7]',\n",
       " 'C(DrivAge)[5]:C(VehPower)[7]',\n",
       " 'C(DrivAge)[6]:C(VehPower)[7]',\n",
       " 'C(DrivAge)[0]:C(VehPower)[8]',\n",
       " 'C(DrivAge)[1]:C(VehPower)[8]',\n",
       " 'C(DrivAge)[2]:C(VehPower)[8]',\n",
       " 'C(DrivAge)[3]:C(VehPower)[8]',\n",
       " 'C(DrivAge)[4]:C(VehPower)[8]',\n",
       " 'C(DrivAge)[5]:C(VehPower)[8]',\n",
       " 'C(DrivAge)[6]:C(VehPower)[8]',\n",
       " 'C(DrivAge)[0]:C(VehPower)[9]',\n",
       " 'C(DrivAge)[1]:C(VehPower)[9]',\n",
       " 'C(DrivAge)[2]:C(VehPower)[9]',\n",
       " 'C(DrivAge)[3]:C(VehPower)[9]',\n",
       " 'C(DrivAge)[4]:C(VehPower)[9]',\n",
       " 'C(DrivAge)[5]:C(VehPower)[9]',\n",
       " 'C(DrivAge)[6]:C(VehPower)[9]']"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[name for name in t_glm4.feature_names_ if \"VehPower\" in name]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The attentive reader might have also noticed that the first level of each categorical variable is omitted from the model. This is a manifestation of the more general concept of [ensuring structural full-rankedness](https://matthewwardrop.github.io/formulaic/guides/contrasts/#guaranteeing-structural-full-rankness)<sup>4</sup>. By default, `glum` and `formulaic` will try to make sure that one does not fall into the [Dummy Variable Trap](https://en.wikipedia.org/wiki/Dummy_variable_(statistics)). Moreover, it even does it in the case of (possibly multi-way) interactions involving categorical variables. It will always drop the necessary number of levels, and no more. If you want to opt out of this behavior (for example because you would like to penalize all levels equally), simply set the `drop_first` parameter during model initialization to `False`. If one only aims to include all levels of a certain variable, and not others, it is possible to do so by using the `spans_intercept` parameter (e.g. `C(VehPower, spans_intercept=False)` would include all levels of `VehPower` even if `drop_first` is set to `True`).\n",
    "\n",
    "<sup>4</sup>: Note, that it does not guarantee that the design matrix is actually full rank. For example, two identical numerical variables will still lead to a rank-deficient design matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Fun with Functions<a class=\"anchor\"></a>\n",
    "\n",
    "The previous example is only scratching the surface of what formulas are capable of. For example, they are capable of evaluating arbitrary Python expressions, which act as if they saw the columns of the input data frame as local variables (`pandas.Series`). The way to tell `glum` that a part of the formula should be evaluated as a Python expression before applying the formula grammar to it is to enclose it in curly braces. As an example, we can easily do the following within the formula itself:\n",
    "\n",
    " 1. Create the outcome variable on the fly instead of doing it beforehand.\n",
    " 2. Include the logarithm of a certain variable in the model.\n",
    " 3. Include a basis spline interpolation of a variable to capture non-linearities in its effect.\n",
    "\n",
    "1\\. works because because formulas can contain [Python operations](https://matthewwardrop.github.io/formulaic/guides/grammar/). 2. and 3. work because formulas are evaluated within a context that is aware of a number of [transforms](https://matthewwardrop.github.io/formulaic/guides/transforms/). To be precise, 2. is a regular transform and 3. is a stateful transform.\n",
    "\n",
    "Let's try it out!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>intercept</th>\n",
       "      <th>VehBrand[B1]</th>\n",
       "      <th>VehBrand[B10]</th>\n",
       "      <th>VehBrand[B11]</th>\n",
       "      <th>VehBrand[B12]</th>\n",
       "      <th>VehBrand[B13]</th>\n",
       "      <th>VehBrand[B14]</th>\n",
       "      <th>VehBrand[B2]</th>\n",
       "      <th>VehBrand[B3]</th>\n",
       "      <th>VehBrand[B4]</th>\n",
       "      <th>...</th>\n",
       "      <th>VehPower[4]</th>\n",
       "      <th>VehPower[5]</th>\n",
       "      <th>VehPower[6]</th>\n",
       "      <th>VehPower[7]</th>\n",
       "      <th>VehPower[8]</th>\n",
       "      <th>VehPower[9]</th>\n",
       "      <th>bs(BonusMalus, 3)[1]</th>\n",
       "      <th>bs(BonusMalus, 3)[2]</th>\n",
       "      <th>bs(BonusMalus, 3)[3]</th>\n",
       "      <th>np.log(Density)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>coefficient</th>\n",
       "      <td>3.808829</td>\n",
       "      <td>-0.060201</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.242194</td>\n",
       "      <td>-0.202517</td>\n",
       "      <td>0.063471</td>\n",
       "      <td>-0.345415</td>\n",
       "      <td>-0.072546</td>\n",
       "      <td>0.00777</td>\n",
       "      <td>0.079391</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.113038</td>\n",
       "      <td>-0.127255</td>\n",
       "      <td>0.060209</td>\n",
       "      <td>0.005577</td>\n",
       "      <td>-0.032114</td>\n",
       "      <td>0.207355</td>\n",
       "      <td>3.178178</td>\n",
       "      <td>0.361951</td>\n",
       "      <td>8.231846</td>\n",
       "      <td>0.121944</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 62 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "             intercept  VehBrand[B1]  VehBrand[B10]  VehBrand[B11]  \\\n",
       "coefficient   3.808829     -0.060201            0.0       0.242194   \n",
       "\n",
       "             VehBrand[B12]  VehBrand[B13]  VehBrand[B14]  VehBrand[B2]  \\\n",
       "coefficient      -0.202517       0.063471      -0.345415     -0.072546   \n",
       "\n",
       "             VehBrand[B3]  VehBrand[B4]  ...  VehPower[4]  VehPower[5]  \\\n",
       "coefficient       0.00777      0.079391  ...    -0.113038    -0.127255   \n",
       "\n",
       "             VehPower[6]  VehPower[7]  VehPower[8]  VehPower[9]  \\\n",
       "coefficient     0.060209     0.005577    -0.032114     0.207355   \n",
       "\n",
       "             bs(BonusMalus, 3)[1]  bs(BonusMalus, 3)[2]  bs(BonusMalus, 3)[3]  \\\n",
       "coefficient              3.178178              0.361951              8.231846   \n",
       "\n",
       "             np.log(Density)  \n",
       "coefficient         0.121944  \n",
       "\n",
       "[1 rows x 62 columns]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula_fun = (\n",
    "    \"{ClaimAmountCut / Exposure} ~ VehBrand + VehGas + Region + Area\"\n",
    "    \" + DrivAge + VehAge + VehPower + bs(BonusMalus, 3) + np.log(Density)\"\n",
    ")\n",
    "\n",
    "t_glm5 = GeneralizedLinearRegressor(\n",
    "    family=TweedieDist,\n",
    "    alpha_search=True,\n",
    "    l1_ratio=1,\n",
    "    fit_intercept=True,\n",
    "    formula=formula_fun,\n",
    ")\n",
    "t_glm5.fit(df_train, sample_weight=df[\"Exposure\"].values[train])\n",
    "\n",
    "pd.DataFrame(\n",
    "    {\"coefficient\": np.concatenate(([t_glm5.intercept_], t_glm5.coef_))},\n",
    "    index=[\"intercept\"] + t_glm5.feature_names_,\n",
    ").T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To allow for even more flexibility, you can add custom transformations that are defined in the context from which the call is made. E.g., we can define a transformation that takes the logarithm of ``VehAge + 1`` after casting it to numeric. To make the formula recognize this transform, you need to explicitly set ``context=0`` when calling the fit method (note that this differs from ``formulaic``'s default, which is already ``context=0``)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>intercept</th>\n",
       "      <th>_log_plus_one(VehAge)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>coefficient</th>\n",
       "      <td>5.046712</td>\n",
       "      <td>-0.151043</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             intercept  _log_plus_one(VehAge)\n",
       "coefficient   5.046712              -0.151043"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def _log_plus_one(x):\n",
    "    return np.log(pd.to_numeric(x) + 1)\n",
    "\n",
    "formula_custom_fun = (\n",
    "    \"{ClaimAmountCut / Exposure} ~ _log_plus_one(VehAge)\"\n",
    ")\n",
    "\n",
    "t_glm6 = GeneralizedLinearRegressor(\n",
    "    family=TweedieDist,\n",
    "    alpha_search=True,\n",
    "    l1_ratio=1,\n",
    "    fit_intercept=True,\n",
    "    formula=formula_custom_fun,\n",
    ")\n",
    "t_glm6.fit(df_train, sample_weight=df[\"Exposure\"].values[train], context=0)\n",
    "\n",
    "pd.DataFrame(\n",
    "    {\"coefficient\": np.concatenate(([t_glm6.intercept_], t_glm6.coef_))},\n",
    "    index=[\"intercept\"] + t_glm6.feature_names_,\n",
    ").T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Miscellaneous Features<a class=\"anchor\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variable Names\n",
    "\n",
    "`glum`'s formula interface provides a lot of control over how the resulting features are named. By default, it follows `formulaic`'s standards, but it can be customized by setting the `interaction_separator` and `categorical_format` paremeters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>intercept</th>\n",
       "      <th>DrivAge__0</th>\n",
       "      <th>DrivAge__1</th>\n",
       "      <th>DrivAge__2</th>\n",
       "      <th>DrivAge__3</th>\n",
       "      <th>DrivAge__4</th>\n",
       "      <th>DrivAge__5</th>\n",
       "      <th>DrivAge__6</th>\n",
       "      <th>VehPower__4</th>\n",
       "      <th>VehPower__5</th>\n",
       "      <th>...</th>\n",
       "      <th>DrivAge__4__x__VehPower__8</th>\n",
       "      <th>DrivAge__5__x__VehPower__8</th>\n",
       "      <th>DrivAge__6__x__VehPower__8</th>\n",
       "      <th>DrivAge__0__x__VehPower__9</th>\n",
       "      <th>DrivAge__1__x__VehPower__9</th>\n",
       "      <th>DrivAge__2__x__VehPower__9</th>\n",
       "      <th>DrivAge__3__x__VehPower__9</th>\n",
       "      <th>DrivAge__4__x__VehPower__9</th>\n",
       "      <th>DrivAge__5__x__VehPower__9</th>\n",
       "      <th>DrivAge__6__x__VehPower__9</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>coefficient</th>\n",
       "      <td>5.007277</td>\n",
       "      <td>1.497079</td>\n",
       "      <td>0.53565</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.152974</td>\n",
       "      <td>-0.210998</td>\n",
       "      <td>-0.205689</td>\n",
       "      <td>0.017896</td>\n",
       "      <td>-0.096153</td>\n",
       "      <td>-0.05484</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.143822</td>\n",
       "      <td>-0.002094</td>\n",
       "      <td>0.512258</td>\n",
       "      <td>0.730534</td>\n",
       "      <td>-0.280869</td>\n",
       "      <td>-0.367669</td>\n",
       "      <td>0.171063</td>\n",
       "      <td>0.022052</td>\n",
       "      <td>-0.270456</td>\n",
       "      <td>0.119634</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 56 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "             intercept  DrivAge__0  DrivAge__1  DrivAge__2  DrivAge__3  \\\n",
       "coefficient   5.007277    1.497079     0.53565         0.0   -0.152974   \n",
       "\n",
       "             DrivAge__4  DrivAge__5  DrivAge__6  VehPower__4  VehPower__5  \\\n",
       "coefficient   -0.210998   -0.205689    0.017896    -0.096153     -0.05484   \n",
       "\n",
       "             ...  DrivAge__4__x__VehPower__8  DrivAge__5__x__VehPower__8  \\\n",
       "coefficient  ...                   -0.143822                   -0.002094   \n",
       "\n",
       "             DrivAge__6__x__VehPower__8  DrivAge__0__x__VehPower__9  \\\n",
       "coefficient                    0.512258                    0.730534   \n",
       "\n",
       "             DrivAge__1__x__VehPower__9  DrivAge__2__x__VehPower__9  \\\n",
       "coefficient                   -0.280869                   -0.367669   \n",
       "\n",
       "             DrivAge__3__x__VehPower__9  DrivAge__4__x__VehPower__9  \\\n",
       "coefficient                    0.171063                    0.022052   \n",
       "\n",
       "             DrivAge__5__x__VehPower__9  DrivAge__6__x__VehPower__9  \n",
       "coefficient                   -0.270456                    0.119634  \n",
       "\n",
       "[1 rows x 56 columns]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula_name = \"PurePremium ~ DrivAge * VehPower\"\n",
    "\n",
    "t_glm7 = GeneralizedLinearRegressor(\n",
    "    family=TweedieDist,\n",
    "    alpha_search=True,\n",
    "    l1_ratio=1,\n",
    "    fit_intercept=True,\n",
    "    formula=formula_name,\n",
    "    interaction_separator=\"__x__\",\n",
    "    categorical_format=\"{name}__{category}\",\n",
    ")\n",
    "t_glm7.fit(df_train, sample_weight=df[\"Exposure\"].values[train])\n",
    "\n",
    "pd.DataFrame(\n",
    "    {\"coefficient\": np.concatenate(([t_glm7.intercept_], t_glm7.coef_))},\n",
    "    index=[\"intercept\"] + t_glm7.feature_names_,\n",
    ").T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Intercept Term\n",
    "\n",
    "Just like in the case of the non-formula interface, the presence of an intercept is determined by the `fit_intercept` argument. In case that the formula specifies a different behavior (e.g., adding `+0` or `-1` while `fit_intercept=True`), an error will be raised."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "formula_noint = \"PurePremium ~ DrivAge * VehPower - 1\"\n",
    "\n",
    "with pytest.raises(ValueError, match=\"The formula sets the intercept to False\"):\n",
    "    t_glm8 = GeneralizedLinearRegressor(\n",
    "        family=TweedieDist,\n",
    "        alpha_search=True,\n",
    "        l1_ratio=1,\n",
    "        fit_intercept=True,\n",
    "        formula=formula_noint,\n",
    "        interaction_separator=\"__x__\",\n",
    "        categorical_format=\"{name}__{category}\",\n",
    "    ).fit(df_train, sample_weight=df[\"Exposure\"].values[train])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### One-Sided Formulas\n",
    "\n",
    "Even when using formulas, the outcome variable can be specified as a vector, as in the interface without formulas. In that case the supplied formula should be one-sided (not contain a `~`), and only describe the right-hand side of the regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>intercept</th>\n",
       "      <th>DrivAge__0</th>\n",
       "      <th>DrivAge__1</th>\n",
       "      <th>DrivAge__2</th>\n",
       "      <th>DrivAge__3</th>\n",
       "      <th>DrivAge__4</th>\n",
       "      <th>DrivAge__5</th>\n",
       "      <th>DrivAge__6</th>\n",
       "      <th>VehPower__4</th>\n",
       "      <th>VehPower__5</th>\n",
       "      <th>...</th>\n",
       "      <th>DrivAge__4__x__VehPower__8</th>\n",
       "      <th>DrivAge__5__x__VehPower__8</th>\n",
       "      <th>DrivAge__6__x__VehPower__8</th>\n",
       "      <th>DrivAge__0__x__VehPower__9</th>\n",
       "      <th>DrivAge__1__x__VehPower__9</th>\n",
       "      <th>DrivAge__2__x__VehPower__9</th>\n",
       "      <th>DrivAge__3__x__VehPower__9</th>\n",
       "      <th>DrivAge__4__x__VehPower__9</th>\n",
       "      <th>DrivAge__5__x__VehPower__9</th>\n",
       "      <th>DrivAge__6__x__VehPower__9</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>coefficient</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.713298</td>\n",
       "      <td>0.783505</td>\n",
       "      <td>0.205914</td>\n",
       "      <td>0.016085</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000094</td>\n",
       "      <td>0.223685</td>\n",
       "      <td>4.66123</td>\n",
       "      <td>4.736272</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.144927</td>\n",
       "      <td>0.001657</td>\n",
       "      <td>0.515373</td>\n",
       "      <td>0.714834</td>\n",
       "      <td>-0.325666</td>\n",
       "      <td>-0.370935</td>\n",
       "      <td>0.20417</td>\n",
       "      <td>0.013222</td>\n",
       "      <td>-0.273913</td>\n",
       "      <td>0.115693</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 56 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "             intercept  DrivAge__0  DrivAge__1  DrivAge__2  DrivAge__3  \\\n",
       "coefficient        0.0    1.713298    0.783505    0.205914    0.016085   \n",
       "\n",
       "             DrivAge__4  DrivAge__5  DrivAge__6  VehPower__4  VehPower__5  \\\n",
       "coefficient         0.0    0.000094    0.223685      4.66123     4.736272   \n",
       "\n",
       "             ...  DrivAge__4__x__VehPower__8  DrivAge__5__x__VehPower__8  \\\n",
       "coefficient  ...                   -0.144927                    0.001657   \n",
       "\n",
       "             DrivAge__6__x__VehPower__8  DrivAge__0__x__VehPower__9  \\\n",
       "coefficient                    0.515373                    0.714834   \n",
       "\n",
       "             DrivAge__1__x__VehPower__9  DrivAge__2__x__VehPower__9  \\\n",
       "coefficient                   -0.325666                   -0.370935   \n",
       "\n",
       "             DrivAge__3__x__VehPower__9  DrivAge__4__x__VehPower__9  \\\n",
       "coefficient                     0.20417                    0.013222   \n",
       "\n",
       "             DrivAge__5__x__VehPower__9  DrivAge__6__x__VehPower__9  \n",
       "coefficient                   -0.273913                    0.115693  \n",
       "\n",
       "[1 rows x 56 columns]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula_onesie = \"DrivAge * VehPower\"\n",
    "\n",
    "t_glm8 = GeneralizedLinearRegressor(\n",
    "    family=TweedieDist,\n",
    "    alpha_search=True,\n",
    "    l1_ratio=1,\n",
    "    fit_intercept=False,\n",
    "    formula=formula_onesie,\n",
    "    interaction_separator=\"__x__\",\n",
    "    categorical_format=\"{name}__{category}\",\n",
    ")\n",
    "t_glm8.fit(\n",
    "    X=df_train, y=df_train[\"PurePremium\"], sample_weight=df[\"Exposure\"].values[train]\n",
    ")\n",
    "\n",
    "pd.DataFrame(\n",
    "    {\"coefficient\": np.concatenate(([t_glm8.intercept_], t_glm8.coef_))},\n",
    "    index=[\"intercept\"] + t_glm8.feature_names_,\n",
    ").T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Missing Values in Categorical Columns\n",
    "\n",
    "By default, `glum` raises a `ValueError` when it encounters a missing value in a categorical variable (`\"raise\"` option). However, there are two other options for handling these cases. They can also be treated as if they represented all-zeros indicators (`\"zero\"` option, which is also the way `pandas.get_dummies` works) or missing values can be treated as their own separate category (`\"convert\"` option).\n",
    "\n",
    "Similarly to the non-formula-based interface, `glum`'s behavior can be set globally using the `cat_missing_method` parameter during model initialization. However, formulas provide some additional flexibility: the `C` function has a `missing_method` parameter, with which users can select an option on a column-by-column basis. Here is an example of doing that (although our dataset does not have any missing values, so these options have no actual effect in this case):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>intercept</th>\n",
       "      <th>C(DrivAge, missing_method='zero')[0]</th>\n",
       "      <th>C(DrivAge, missing_method='zero')[1]</th>\n",
       "      <th>C(DrivAge, missing_method='zero')[2]</th>\n",
       "      <th>C(DrivAge, missing_method='zero')[3]</th>\n",
       "      <th>C(DrivAge, missing_method='zero')[4]</th>\n",
       "      <th>C(DrivAge, missing_method='zero')[5]</th>\n",
       "      <th>C(DrivAge, missing_method='zero')[6]</th>\n",
       "      <th>C(VehPower, missing_method='convert')[4]</th>\n",
       "      <th>C(VehPower, missing_method='convert')[5]</th>\n",
       "      <th>C(VehPower, missing_method='convert')[6]</th>\n",
       "      <th>C(VehPower, missing_method='convert')[7]</th>\n",
       "      <th>C(VehPower, missing_method='convert')[8]</th>\n",
       "      <th>C(VehPower, missing_method='convert')[9]</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>coefficient</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.786703</td>\n",
       "      <td>0.742765</td>\n",
       "      <td>0.239528</td>\n",
       "      <td>0.096531</td>\n",
       "      <td>0.071118</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.201078</td>\n",
       "      <td>4.637267</td>\n",
       "      <td>4.679391</td>\n",
       "      <td>4.863387</td>\n",
       "      <td>4.77263</td>\n",
       "      <td>4.749673</td>\n",
       "      <td>4.970188</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             intercept  C(DrivAge, missing_method='zero')[0]  \\\n",
       "coefficient        0.0                              1.786703   \n",
       "\n",
       "             C(DrivAge, missing_method='zero')[1]  \\\n",
       "coefficient                              0.742765   \n",
       "\n",
       "             C(DrivAge, missing_method='zero')[2]  \\\n",
       "coefficient                              0.239528   \n",
       "\n",
       "             C(DrivAge, missing_method='zero')[3]  \\\n",
       "coefficient                              0.096531   \n",
       "\n",
       "             C(DrivAge, missing_method='zero')[4]  \\\n",
       "coefficient                              0.071118   \n",
       "\n",
       "             C(DrivAge, missing_method='zero')[5]  \\\n",
       "coefficient                                   0.0   \n",
       "\n",
       "             C(DrivAge, missing_method='zero')[6]  \\\n",
       "coefficient                              0.201078   \n",
       "\n",
       "             C(VehPower, missing_method='convert')[4]  \\\n",
       "coefficient                                  4.637267   \n",
       "\n",
       "             C(VehPower, missing_method='convert')[5]  \\\n",
       "coefficient                                  4.679391   \n",
       "\n",
       "             C(VehPower, missing_method='convert')[6]  \\\n",
       "coefficient                                  4.863387   \n",
       "\n",
       "             C(VehPower, missing_method='convert')[7]  \\\n",
       "coefficient                                   4.77263   \n",
       "\n",
       "             C(VehPower, missing_method='convert')[8]  \\\n",
       "coefficient                                  4.749673   \n",
       "\n",
       "             C(VehPower, missing_method='convert')[9]  \n",
       "coefficient                                  4.970188  "
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula_missing = \"C(DrivAge, missing_method='zero') + C(VehPower, missing_method='convert')\"\n",
    "\n",
    "t_glm9 = GeneralizedLinearRegressor(\n",
    "    family=TweedieDist,\n",
    "    alpha_search=True,\n",
    "    l1_ratio=1,\n",
    "    fit_intercept=False,\n",
    "    formula=formula_missing,\n",
    "\n",
    ")\n",
    "t_glm9.fit(\n",
    "    X=df_train, y=df_train[\"PurePremium\"], sample_weight=df[\"Exposure\"].values[train]\n",
    ")\n",
    "\n",
    "pd.DataFrame(\n",
    "    {\"coefficient\": np.concatenate(([t_glm9.intercept_], t_glm9.coef_))},\n",
    "    index=[\"intercept\"] + t_glm9.feature_names_,\n",
    ").T"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "glum-dev",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.2"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}