---
title: "Robustly important variables in credit scoring"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Robustly important variables in credit scoring}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
ready <- all(vapply(
  c("randomForest", "ranger", "rsample"),
  requireNamespace, logical(1), quietly = TRUE
))
knitr::opts_chunk$set(eval = ready)
```

A lender must be able to say which characteristics drive a decision, and must be
able to defend the claim. "Our random forest's permutation importance ranked
income first" does not survive the obvious follow-up: *would a different method,
or a different seed, or a different fold, have said the same?*

This vignette answers that question end to end, using judges drawn along all
four axes at once, which are method, seed, fold and model family. It ends up
answering it wrongly, and being unable to tell that it has, which is the point.

## The portfolio

`applications` ships with the package. Eight hundred loan applications, seven
predictors, and a default indicator that comes out `"yes"` 23.0 per cent of the
time. It is generated rather than collected, so the ordering a ranking ought to
recover is known rather than argued over, and it is stored on the data frame.

```{r data}
library(rankimp)

truth <- sort(attr(applications, "effects"), decreasing = TRUE)
truth
```

Two features of it matter here, and both were built in on purpose. `income` and
`bureau_score` come from one latent creditworthiness and correlate at 0.84, so
they stand in for each other and the credit for the signal has to be divided
somehow. `prior_arrears` is the largest effect in the data and takes only six
distinct values, which is the case mean decrease in impurity handles badly. See
`?applications`.

## The panel

Two model families, two methods, three folds, two seeds: 24 judges.

```{r panel}
library(randomForest)
library(ranger)

set.seed(1)
rf_fit <- randomForest(default ~ ., data = applications, ntree = 300)
set.seed(1)
rgr_fit <- ranger(default ~ ., data = applications, num.trees = 300,
                  importance = "impurity", probability = TRUE)

J <- importance_judges(
  fit_list  = list(rf = rf_fit, ranger = rgr_fit),
  methods   = c("permutation", "mdi"),
  data      = applications,
  target    = "default",
  resamples = rsample::vfold_cv(applications, v = 3),
  seeds     = 1:2
)
J
```

Any one of those 24 rows, reported alone, would look like a result.

## What the panel agrees on

```{r consensus}
cr <- consensus_rank(J)
cr
```

Compare that against `truth` above. The consensus does **not** recover the
ordering. It promotes `debt_ratio` over `prior_arrears`, which is the largest
effect in the data, and it swaps `income` ahead of `bureau_score`, which is the
correlated pair being divided the wrong way round. Two errors, both in the top
four, and `tau_x` of 0.72 says the judges were a long way from unanimous while
saying nothing about which way they were wrong.

It gets the bottom of the table right. `employment_yrs` is fifth, where it
belongs, and the two predictors that enter the outcome nowhere come last.

## What survives resampling

```{r confsets}
cb <- rank_confsets(cr, n_boot = 500)
cb
```

Read the intervals against the truth rather than against the consensus.

`employment_yrs` has the interval `[5, 5]`, a point, and it is right.
`credit_lines` and `age` share `[6, 7]`, which is the procedure declining to
order two variables that have no order, and that is right too. So far the
uncertainty is doing what it was built to do.

The top four are the problem. `debt_ratio` gets `[1, 2]` and `prior_arrears`
gets `[1, 4]`, so resampling the judges is *more* confident about the variable
that is second in truth than about the one that is first.

```{r select}
rank_select(cb, threshold = 2)
```

That is the failure worth the whole vignette. Asked for the variables it can
place in the top two, the procedure returns one, confidently, and it is not the
largest effect in the data. Nothing in the output says so.

```{r select-more}
rank_select(cb, threshold = 4)
```

Widen the question and the answer becomes defensible: these four characteristics
drive default, ahead of the other three. That statement is true. Any statement
that orders the four is not, and the interval for `prior_arrears` is the only
thing hinting that the ordering is unsafe.

```{r plot, fig.width = 6, fig.height = 4}
autoplot(cb)
```

## Where the disagreement lives

A panel of 24 that disagrees is worth taking apart before it is averaged.

```{r clusters}
het <- judge_clusters(J)
het
```

Two groups, at `p = 0.005` against a single population. The useful question is
which axis they fall along, and the panel records the provenance of every judge:

```{r composition}
provenance <- attr(J, "provenance")
table(method = provenance$method, cluster = het$cluster)
table(model = provenance$model, cluster = het$cluster)
```

The seam is the **method**, and it is not approximate: all twelve MDI judges in
one group, all twelve permutation judges in the other, while randomForest and
ranger split six and six across both. Swapping the engine changes nothing.
Swapping how importance is defined changes the answer.

```{r centres}
het$centres
```

There is the mechanism. Permutation puts `prior_arrears` **first**, which is
where the truth puts it. Impurity puts it **fourth**. `prior_arrears` is a count
with six distinct values and mean decrease in impurity is biased against
predictors with few places to split, so the panel has rediscovered a known
property of the estimator without being told to look for it.

Averaging the two groups is what produced the consensus above. Twelve judges
were right about `prior_arrears` and twelve were wrong, and the median ranking
of the 24 landed between them, on second.

## What a consensus cannot do for you

Every judge here is a tree ensemble. Averaging over methods, seeds, folds and
two forest implementations measures how much the answer depends on those
choices, and nothing else. A bias that all 24 judges share, because they are all
built from trees, passes through the consensus untouched and comes out looking
like agreement.

That is not an abstract worry in this example, it is what happened. The
cardinality bias was visible only because the panel contained a method that does
not share it, and it was still strong enough to move the largest effect in the
data off the top of the consensus. Half the panel was enough to see the bias and
not enough to correct it.

Read the intervals as what they are. They describe how much the answer moves
when the judges are resampled, so they are narrow when the judges agree, and
judges agree when they share an assumption as readily as when they are right.
Widen the panel along the axis you are worried about, or the narrow confidence
set you get back is measuring the wrong thing.