insurance_gam.ebm

 98class InsuranceEBM:
 99    """
100    Explainable Boosting Machine configured for insurance pricing.
101
102    Wraps interpretML's ExplainableBoostingRegressor and adds:
103    - Exposure-aware fitting via log(exposure) init_score
104    - Polars DataFrame input support
105    - predict() returning expected values (frequency, severity, or pure premium)
106      rather than raw log scores
107    - Deviance scoring with the correct GLM family
108
109    Parameters
110    ----------
111    loss : str
112        Loss function. One of 'poisson' (default), 'tweedie', 'gamma', 'mse',
113        'mae', 'huber'. Poisson is standard for claim frequency; Tweedie or
114        Gamma for pure premium or severity.
115    variance_power : float
116        Tweedie variance power. Only used when loss='tweedie'. Set to 1.0 for
117        Poisson, 2.0 for Gamma, 1.5 for compound Poisson-Gamma (pure premium).
118    interactions : int or str
119        Number of pairwise interaction terms, or a multiplier string of the form
120        ``'Nx'`` (e.g. ``'3x'``, ``'1x'``) meaning ``N * n_features`` interaction
121        terms. ``'3x'`` is the default and usually works well. Use ``0`` or ``'0x'``
122        to disable interactions entirely (reduces memory and training time).
123        Passed through to ExplainableBoostingRegressor as ``interactions``.
124    monotone_constraints : dict, optional
125        Mapping of {feature_name: direction} where direction is +1 (increasing)
126        or -1 (decreasing). Applied at fit time via the interpretML API.
127    **ebm_kwargs
128        Any additional keyword arguments forwarded to ExplainableBoostingRegressor.
129        Common overrides: n_estimators, learning_rate, max_bins, min_samples_leaf,
130        random_state, n_jobs.
131
132    Notes
133    -----
134    interpretML uses a log link for Poisson, Tweedie, and Gamma families.
135    Exposure is incorporated as an offset: init_score = log(exposure). When no
136    exposure is provided, all exposures are treated as 1 (no offset needed).
137    """
138
139    def __init__(
140        self,
141        loss: str = "poisson",
142        variance_power: float = 1.5,
143        interactions: Union[int, str] = "3x",
144        monotone_constraints: Optional[dict] = None,
145        **ebm_kwargs,
146    ) -> None:
147        if loss not in _SUPPORTED_LOSSES:
148            raise ValueError(f"loss must be one of {_SUPPORTED_LOSSES}, got '{loss}'")
149        self.loss = loss
150        self.variance_power = variance_power
151        self.interactions = interactions
152        self.monotone_constraints = monotone_constraints or {}
153        self.ebm_kwargs = ebm_kwargs
154        self.ebm_: Optional[ExplainableBoostingRegressor] = None
155        self._feature_names: Optional[list] = None
156
157    def _build_ebm(self, feature_names: list) -> ExplainableBoostingRegressor:
158        """Instantiate the underlying EBM with appropriate settings."""
159        # Resolve interactions: 'Nx' means N * n_features (e.g. '3x', '1x', '0x')
160        if isinstance(self.interactions, str):
161            m = re.match(r'^(\d+)x$', self.interactions)
162            if m:
163                n_interactions = int(m.group(1)) * len(feature_names)
164            else:
165                raise ValueError(
166                    f"interactions string must be of the form 'Nx' (e.g. '3x', '1x'), "
167                    f"got {self.interactions!r}"
168                )
169        else:
170            n_interactions = int(self.interactions)
171
172        # Map loss string to interpretML objective.
173        # For Tweedie, the variance power must be embedded in the objective
174        # string as "tweedie_deviance:variance_power=<p>". Without it,
175        # interpretML falls back to a default power and the deviance is wrong.
176        if self.loss == "tweedie":
177            objective = f"tweedie_deviance:variance_power={self.variance_power}"
178        else:
179            objective_map = {
180                "poisson": "poisson_deviance",
181                "gamma": "gamma_deviance",
182                "mse": "rmse",
183                "mae": "mae",
184                "huber": "huber",
185            }
186            objective = objective_map[self.loss]
187
188        # tweedie_power kept for reference; power is now embedded in objective string
189        tweedie_power = self.variance_power if self.loss == "tweedie" else None
190
191        # Build monotone constraints list aligned to feature order
192        mc = []
193        for fn in feature_names:
194            mc.append(self.monotone_constraints.get(fn, 0))
195
196        kwargs = dict(self.ebm_kwargs)
197        kwargs.setdefault("random_state", 42)
198        kwargs.setdefault("n_jobs", -1)
199
200        ebm_args = {
201            "interactions": n_interactions,
202            "feature_names": feature_names,
203            "feature_types": None,  # let interpretML infer
204            "objective": objective,
205            "monotone_constraints": mc if any(v != 0 for v in mc) else None,
206        }
207        # tweedie_exp_target is not needed — power is embedded in the objective string above.
208        # No additional kwarg required for Tweedie.
209
210        ebm_args.update(kwargs)
211        return ExplainableBoostingRegressor(**ebm_args)
212
213    def fit(
214        self,
215        X: Union[pl.DataFrame, pd.DataFrame],
216        y: Union[np.ndarray, list, pl.Series, pd.Series],
217        exposure: Optional[Union[np.ndarray, list, pl.Series, pd.Series]] = None,
218        sample_weight: Optional[Union[np.ndarray, list, pl.Series, pd.Series]] = None,
219    ) -> "InsuranceEBM":
220        """
221        Fit the EBM to insurance training data.
222
223        Parameters
224        ----------
225        X : polars.DataFrame or pandas.DataFrame
226            Feature matrix. Column names are used as feature identifiers.
227        y : array-like
228            Target variable. For Poisson frequency models this is claim count.
229            For Tweedie pure premium models this is incurred loss amount (may
230            include zeros). For Gamma severity models this is average claim
231            cost on non-zero rows only.
232        exposure : array-like, optional
233            Exposure measure (e.g. earned car years). Incorporated as a log
234            offset: init_score = log(exposure). If omitted, all exposures are
235            assumed to be 1.
236        sample_weight : array-like, optional
237            Row-level observation weights. Applied in addition to exposure.
238            Useful for credibility weighting or re-balancing portfolio mix.
239
240        Returns
241        -------
242        self
243        """
244        X_pd = _to_pandas(X)
245        self._feature_names = list(X_pd.columns)
246
247        y_arr = _ensure_array(y)
248        exp_arr = _ensure_array(exposure)
249        sw_arr = _ensure_array(sample_weight)
250
251        self.ebm_ = self._build_ebm(self._feature_names)
252
253        fit_kwargs: dict = {}
254
255        if exp_arr is not None:
256            # Exposure enters as a log offset (init_score)
257            if np.any(exp_arr <= 0):
258                raise ValueError("All exposure values must be strictly positive.")
259            fit_kwargs["init_score"] = np.log(exp_arr)
260
261        if sw_arr is not None:
262            fit_kwargs["sample_weight"] = sw_arr
263
264        self.ebm_.fit(X_pd, y_arr, **fit_kwargs)
265        return self
266
267    def _check_fitted(self) -> None:
268        if self.ebm_ is None:
269            raise RuntimeError("Model has not been fitted yet. Call fit() first.")
270
271    def predict(
272        self,
273        X: Union[pl.DataFrame, pd.DataFrame],
274        exposure: Optional[Union[np.ndarray, list, pl.Series, pd.Series]] = None,
275    ) -> np.ndarray:
276        """
277        Predict expected values on the original (not log) scale.
278
279        For log-link families (Poisson, Tweedie, Gamma), the EBM produces additive
280        scores on the log scale. This method applies exp() and multiplies by
281        exposure to return predicted counts or amounts.
282
283        Parameters
284        ----------
285        X : polars.DataFrame or pandas.DataFrame
286            Feature matrix. Must have the same columns as the training data.
287        exposure : array-like, optional
288            Exposure values for prediction. If provided, predictions are scaled
289            by exposure (equivalent to adding log(exposure) to the log score).
290            For rate models, pass exposure=None and scale externally.
291
292        Returns
293        -------
294        numpy.ndarray
295            Predicted values on the response scale.
296        """
297        self._check_fitted()
298        X_pd = _to_pandas(X)
299        log_scores = self.ebm_.predict(X_pd)
300
301        _log_link_losses = {"poisson", "tweedie", "gamma"}
302        if self.loss in _log_link_losses:
303            predictions = np.exp(log_scores)
304            if exposure is not None:
305                exp_arr = _ensure_array(exposure)
306                if np.any(exp_arr <= 0):
307                    raise ValueError("All exposure values must be strictly positive.")
308                predictions = predictions * exp_arr
309        else:
310            predictions = log_scores
311
312        return predictions
313
314    def predict_log_score(
315        self,
316        X: Union[pl.DataFrame, pd.DataFrame],
317    ) -> np.ndarray:
318        """
319        Return raw additive scores on the log scale (before exp transformation).
320
321        Useful for combining predictions from separate frequency and severity
322        models on a common additive scale.
323
324        Parameters
325        ----------
326        X : polars.DataFrame or pandas.DataFrame
327            Feature matrix.
328
329        Returns
330        -------
331        numpy.ndarray
332            Log-scale scores.
333        """
334        self._check_fitted()
335        X_pd = _to_pandas(X)
336        return self.ebm_.predict(X_pd)
337
338    def score(
339        self,
340        X: Union[pl.DataFrame, pd.DataFrame],
341        y: Union[np.ndarray, list, pl.Series, pd.Series],
342        exposure: Optional[Union[np.ndarray, list, pl.Series, pd.Series]] = None,
343    ) -> float:
344        """
345        Compute mean deviance on test data (lower is better).
346
347        The deviance family matches the loss used at fit time. For Tweedie,
348        the same variance_power is used.
349
350        Parameters
351        ----------
352        X : polars.DataFrame or pandas.DataFrame
353            Feature matrix.
354        y : array-like
355            Observed values. For Poisson frequency models, this must be raw
356            claim counts (not rates). When exposure is provided, the model
357            computes expected counts (mu = exp(log_score) * exposure) and
358            compares to y; passing claim rates instead of counts will give
359            incorrect deviance values.
360        exposure : array-like, optional
361            Exposure values used to scale predictions.
362
363        Returns
364        -------
365        float
366            Mean deviance (negated so higher = better, consistent with sklearn).
367        """
368        self._check_fitted()
369        y_arr = _ensure_array(y)
370        exp_arr = _ensure_array(exposure)
371        y_pred = self.predict(X, exposure=exposure)
372        weights = exp_arr  # use exposure as weights in deviance calculation
373
374        if self.loss == "poisson":
375            dev = _deviance_poisson(y_arr, y_pred, weights)
376        elif self.loss == "gamma":
377            dev = _deviance_gamma(y_arr, y_pred, weights)
378        elif self.loss == "tweedie":
379            dev = _deviance_tweedie(y_arr, y_pred, self.variance_power, weights)
380        else:
381            # MSE for non-GLM losses
382            if weights is not None:
383                dev = float(np.average((y_arr - y_pred) ** 2, weights=weights))
384            else:
385                dev = float(np.mean((y_arr - y_pred) ** 2))
386
387        return -dev  # higher is better convention
388
389    @property
390    def feature_names(self) -> list:
391        """Feature names from the training data."""
392        self._check_fitted()
393        return self._feature_names
394
395    def __repr__(self) -> str:
396        fitted = self.ebm_ is not None
397        return (
398            f"InsuranceEBM(loss='{self.loss}', "
399            f"variance_power={self.variance_power}, "
400            f"interactions='{self.interactions}', "
401            f"fitted={fitted})"
402        )

 96class RelativitiesTable:
 97    """
 98    Extract and present insurance-standard relativity tables from a fitted InsuranceEBM.
 99
100    A relativity is defined as exp(score_bin - score_base) where the base bin is the
101    modal bin (the one with the highest training weight). A relativity of 1.0 means
102    the bin contributes the same as the base level; >1 means higher risk, <1 lower.
103
104    Parameters
105    ----------
106    model : InsuranceEBM
107        A fitted InsuranceEBM instance.
108
109    Examples
110    --------
111    >>> rt = RelativitiesTable(model)
112    >>> rt.table("driver_age")
113    shape: (10, 3)
114    ┌────────────────┬───────────┬────────────┐
115    │ bin_label      ┆ raw_score ┆ relativity │
116    ...
117    """
118
119    def __init__(self, model: InsuranceEBM) -> None:
120        model._check_fitted()
121        self.model = model
122        self._ebm = model.ebm_
123
124    def table(self, feature: str) -> pl.DataFrame:
125        """
126        Return relativity table for a single feature.
127
128        Parameters
129        ----------
130        feature : str
131            Name of the feature (must be a main effect, not an interaction).
132
133        Returns
134        -------
135        polars.DataFrame
136            Columns: bin_label (str), raw_score (float), relativity (float).
137            The base bin has relativity = 1.0.
138        """
139        names, scores = _get_ebm_shape(self._ebm, feature)
140        base_idx = _modal_bin_idx(self._ebm, feature)
141
142        # Clip base_idx to valid range (can happen if weights vector is shorter)
143        base_idx = min(base_idx, len(scores) - 1)
144        base_score = scores[base_idx]
145
146        relativities = np.exp(scores - base_score)
147
148        return pl.DataFrame(
149            {
150                "bin_label": names.tolist(),
151                "raw_score": scores.tolist(),
152                "relativity": relativities.tolist(),
153            }
154        )
155
156    def plot(
157        self,
158        feature: str,
159        kind: str = "bar",
160        ax=None,
161        title: Optional[str] = None,
162    ):
163        """
164        Bar or line chart of relativities for a feature.
165
166        Parameters
167        ----------
168        feature : str
169            Feature name.
170        kind : str
171            'bar' (default) or 'line'.
172        ax : matplotlib.axes.Axes, optional
173            Axes to draw on. Creates a new figure if None.
174        title : str, optional
175            Plot title. Defaults to the feature name.
176
177        Returns
178        -------
179        matplotlib.axes.Axes
180        """
181        import matplotlib.pyplot as plt
182
183        df = self.table(feature)
184        labels = df["bin_label"].to_list()
185        rels = df["relativity"].to_list()
186
187        if ax is None:
188            fig, ax = plt.subplots(figsize=(max(8, len(labels) * 0.6), 5))
189
190        if kind == "bar":
191            colours = ["#d62728" if r > 1.0 else "#2ca02c" for r in rels]
192            ax.bar(range(len(labels)), rels, color=colours, alpha=0.8, edgecolor="white")
193        elif kind == "line":
194            ax.plot(range(len(labels)), rels, marker="o", color="#1f77b4", linewidth=2)
195            ax.fill_between(range(len(labels)), 1.0, rels, alpha=0.15, color="#1f77b4")
196        else:
197            raise ValueError(f"kind must be 'bar' or 'line', got '{kind}'")
198
199        ax.axhline(1.0, color="black", linewidth=0.8, linestyle="--")
200        ax.set_xticks(range(len(labels)))
201        ax.set_xticklabels(labels, rotation=45, ha="right", fontsize=9)
202        ax.set_ylabel("Relativity")
203        ax.set_title(title or feature)
204        ax.grid(axis="y", alpha=0.3)
205
206        return ax
207
208    def export_excel(self, path: Union[str, Path]) -> None:
209        """
210        Export relativity tables for all main-effect features to an Excel workbook.
211
212        One sheet per feature. Interaction terms are skipped (they have 2-D score
213        arrays that don't translate cleanly to a flat relativity table).
214
215        Parameters
216        ----------
217        path : str or Path
218            Output path for the .xlsx file.
219
220        Notes
221        -----
222        Requires openpyxl. Install with: pip install insurance-gam[excel]
223        """
224        try:
225            import openpyxl  # noqa: F401
226        except ImportError:
227            raise ImportError(
228                "openpyxl is required for Excel export. "
229                "Install with: pip install insurance-gam[excel]"
230            )
231
232        path = Path(path)
233        feature_names = self.model.feature_names
234        frames: dict[str, pl.DataFrame] = {}
235
236        for feature in feature_names:
237            try:
238                df = self.table(feature)
239                frames[feature] = df
240            except ValueError:
241                # Skip interaction terms or features that can't be extracted
242                continue
243
244        if not frames:
245            raise RuntimeError("No main-effect features found to export.")
246
247        # Write via pandas ExcelWriter — polars converts to pandas per sheet
248        import pandas as pd
249        with pd.ExcelWriter(path, engine="openpyxl") as writer:
250            for sheet_name, df in frames.items():
251                # Excel sheet names are limited to 31 chars
252                safe_name = sheet_name[:31]
253                df.to_pandas().to_excel(writer, sheet_name=safe_name, index=False)
254
255    def summary(self) -> pl.DataFrame:
256        """
257        Summary table of all main-effect features.
258
259        Returns
260        -------
261        polars.DataFrame
262            Columns: feature (str), n_bins (int), min_relativity (float),
263            max_relativity (float), range (float).
264
265            range = max_relativity / min_relativity, a measure of how much
266            leverage the factor has over premium.
267        """
268        rows = []
269        for feature in self.model.feature_names:
270            try:
271                df = self.table(feature)
272                rels = df["relativity"].to_numpy()
273                rows.append(
274                    {
275                        "feature": feature,
276                        "n_bins": len(rels),
277                        "min_relativity": float(np.min(rels)),
278                        "max_relativity": float(np.max(rels)),
279                        "range": float(np.max(rels) / np.maximum(np.min(rels), 1e-10)),
280                    }
281                )
282            except ValueError:
283                continue
284
285        return pl.DataFrame(rows).sort("range", descending=True)

 96class MonotonicityEditor:
 97    """
 98    Post-fit monotonicity enforcement for EBM shape functions.
 99
100    Applies isotonic regression to bin scores for a named feature, adjusting
101    the stored term_scores_ in-place. The base model object (ebm_) is modified
102    directly, so predictions from InsuranceEBM.predict() will reflect the change.
103
104    Parameters
105    ----------
106    model : InsuranceEBM
107        A fitted InsuranceEBM instance. The underlying ebm_ will be modified in-place.
108
109    Notes
110    -----
111    The adjustment is not reversible once applied (the original scores are
112    overwritten). Keep a copy of the model if you need to compare before/after
113    predictions on a holdout set.
114    """
115
116    def __init__(self, model: InsuranceEBM) -> None:
117        model._check_fitted()
118        self.model = model
119        self._ebm = model.ebm_
120
121    def enforce(
122        self,
123        feature: str,
124        direction: Literal["increase", "decrease", "auto"] = "auto",
125    ) -> "MonotonicityEditor":
126        """
127        Enforce monotonicity on a feature's shape function via isotonic regression.
128
129        The bin scores (term_scores_) are updated in-place. After calling this,
130        model.predict() will reflect the monotone shape.
131
132        Parameters
133        ----------
134        feature : str
135            Feature name.
136        direction : str
137            'increase', 'decrease', or 'auto'. When 'auto', the dominant
138            direction is detected from the existing shape function.
139
140        Returns
141        -------
142        self
143            Returns self for method chaining.
144        """
145        term_idx, scores = _get_term_scores(self._ebm, feature)
146
147        if direction == "auto":
148            direction = _detect_direction(scores)
149
150        if direction not in ("increase", "decrease"):
151            raise ValueError(f"direction must be 'increase', 'decrease', or 'auto', got '{direction}'")
152
153        if _is_monotone(scores, direction):
154            # Already monotone — nothing to do
155            return self
156
157        # interpretML stores a special first bin for missing values (index 0)
158        # We preserve the missing-value score and only enforce on the main bins
159        if len(scores) > 1:
160            main_scores = scores[1:]
161            main_scores_mono = _isotonic_regression(main_scores, direction)
162            new_scores = np.concatenate([[scores[0]], main_scores_mono])
163        else:
164            new_scores = _isotonic_regression(scores, direction)
165
166        self._ebm.term_scores_[term_idx] = new_scores
167        return self
168
169    def check(self, feature: str, direction: Optional[str] = None) -> bool:
170        """
171        Check whether a feature's shape function is monotone.
172
173        Parameters
174        ----------
175        feature : str
176            Feature name.
177        direction : str, optional
178            'increase' or 'decrease'. If None, checks both directions (returns
179            True if monotone in either direction).
180
181        Returns
182        -------
183        bool
184            True if the shape function is monotone in the specified direction.
185        """
186        _, scores = _get_term_scores(self._ebm, feature)
187        main_scores = scores[1:] if len(scores) > 1 else scores
188
189        if direction is None:
190            return _is_monotone(main_scores, "increase") or _is_monotone(main_scores, "decrease")
191        return _is_monotone(main_scores, direction)
192
193    def plot_before_after(
194        self,
195        feature: str,
196        scores_before: np.ndarray,
197        direction: Optional[str] = None,
198        ax=None,
199    ):
200        """
201        Side-by-side comparison of shape function before and after monotonicity enforcement.
202
203        Because enforce() modifies scores in-place, you must capture the scores
204        before calling enforce(). Use MonotonicityEditor.get_scores() to obtain
205        a copy before enforcement.
206
207        Parameters
208        ----------
209        feature : str
210            Feature name.
211        scores_before : numpy.ndarray
212            Scores before enforcement (use get_scores() to capture these).
213        direction : str, optional
214            Direction label for the title.
215        ax : matplotlib.axes.Axes, optional
216            Two-element axes array. Creates a new figure if None.
217
218        Returns
219        -------
220        matplotlib.figure.Figure
221        """
222        import matplotlib.pyplot as plt
223
224        _, scores_after = _get_term_scores(self._ebm, feature)
225        main_before = scores_before[1:] if len(scores_before) > 1 else scores_before
226        main_after = scores_after[1:] if len(scores_after) > 1 else scores_after
227
228        n = max(len(main_before), len(main_after))
229        x = np.arange(n)
230
231        if ax is None:
232            fig, axes = plt.subplots(1, 2, figsize=(14, 5))
233        else:
234            axes = ax
235            fig = axes[0].get_figure()
236
237        for a, scores, title in zip(
238            axes,
239            [main_before, main_after],
240            ["Before enforcement", f"After enforcement ({direction or 'monotone'})"],
241        ):
242            a.plot(x[:len(scores)], scores, marker="o", color="#1f77b4", linewidth=2)
243            a.set_title(f"{feature} — {title}")
244            a.set_xlabel("Bin index")
245            a.set_ylabel("Score (log scale)")
246            a.grid(alpha=0.3)
247
248        fig.tight_layout()
249        return fig
250
251    def get_scores(self, feature: str) -> np.ndarray:
252        """
253        Return a copy of the current scores for a feature.
254
255        Use this to capture scores before calling enforce() if you want to
256        plot a before/after comparison.
257
258        Parameters
259        ----------
260        feature : str
261            Feature name.
262
263        Returns
264        -------
265        numpy.ndarray
266            Copy of term_scores_ for this feature.
267        """
268        _, scores = _get_term_scores(self._ebm, feature)
269        return scores.copy()