---
title: "A credit scoring case study"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{A credit scoring case study}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
has_rf <- requireNamespace("randomForest", quietly = TRUE)
has_igraph <- has_rf && requireNamespace("igraph", quietly = TRUE)
knitr::opts_chunk$set(eval = has_rf, purl = has_rf)
```

Credit scoring is the running example across `Proximum`, `rankimp` and
`e2tree`, for three reasons. The regulatory context makes explanation a
requirement rather than a nicety; the data are large enough that the $O(n^2)$
cost of a proximity matrix bites; and the predictors are correlated enough that
different importance measures genuinely disagree.

The portfolio here is synthetic. Public credit panels that can be
redistributed have had the protected attributes stripped out of them, and the
third question below cannot be asked without one. `?loans` records what was
built into the data and `inst/data-raw/loans.R` is the specification as code.
Read the numbers as a demonstration of what the diagnostics say, not as a fact
about consumer lending.

```{r setup}
library(Proximum)
library(randomForest)

str(loans, give.attr = FALSE)
```

Two columns are not predictors. `vintage` is the origination year, held back
for the second question, and `applicant_group` is a protected attribute that
the model is not allowed to see and that the third question is about.

```{r model}
predictors <- setdiff(names(loans), c("applicant_group", "vintage"))

set.seed(1)
rows <- sample(nrow(loans), 500)
book <- loans[rows, ]

rf <- randomForest(default ~ ., data = book[, predictors], ntree = 500)
px <- as_proximity(rf, newdata = book[, predictors])
summary(px)
```

## 1. Are there segments, and are they risk segments?

A lender's segments are drawn by hand: thin file, high utilisation, short
tenure. The forest draws its own, and they are in the proximity rather than in
any coefficient. `autoplot(type = "network")` keeps the pairs the forest puts
above a threshold and reads the communities of the graph that remain.

```{r network, eval = has_igraph, fig.width = 6.5, fig.height = 5}
set.seed(11)
picture <- autoplot(px, type = "network", threshold = 0.2)
picture
```

The communities are not labelled with anything the forest was told. Whether
they are risk segments is a question to be asked of them afterwards:

```{r communities, eval = has_igraph}
community <- picture$layers[[2]]$data$community
segments <- aggregate(
  cbind(default == "yes", score, utilisation) ~ community,
  data.frame(community, book), mean
)
names(segments) <- c("community", "default_rate", "score", "utilisation")
segments$n <- as.vector(table(community))
segments[order(-segments$default_rate), ]

mean(book$default == "yes")
```

Read the table rather than this paragraph for the figures, since a community
label is a clustering of a particular fit and the exact counts move with the
seed. What does not move is the shape of it. One community of about sixty
borrowers defaults at close to nine in ten, against a book rate of one in six.
It is neither the low-score community nor the high-utilisation community: its
mean score sits near 614 and its mean utilisation near 0.84, while the two
communities that have one of those without the other default at well under a
tenth. The forest has found an interaction, and the proximity shows it as a
block, because a block of observations the model keeps together is what an
interaction looks like once you stop reading coefficients.

Some communities hold a single borrower. Those are not segments, they are what
a threshold of 0.2 leaves isolated, and reporting them as segments would be
reading the threshold rather than the data.

## 2. Does the structure survive a shift in the window?

The portfolio spans four origination years and the relationship moves across
them: by construction, utilisation carries three times the weight in 2022 that
it carries in 2019. A lender who refits annually wants to know whether the
model has reorganised its view of the book or merely re-estimated it.

The question cannot be answered by comparing two proximity matrices on
different borrowers. It is answered by applying both forests to the *same*
borrowers and comparing what each makes of them:

```{r windows}
early <- loans$vintage %in% c("2019", "2020")
late  <- loans$vintage %in% c("2021", "2022")

set.seed(3)
held <- loans[late, ][sample(sum(late), 400), ]

forest_on <- function(subset, seed) {
  set.seed(seed)
  randomForest(default ~ ., data = loans[subset, predictors], ntree = 500)
}
represent <- function(fit) as_proximity(fit, newdata = held[, predictors])
```

A difference is only a difference against the noise floor, and the noise floor
here is how much two refits of the same window disagree:

```{r floor}
replicates <- lapply(1:4, function(i) represent(forest_on(early, 100 + i)))
stability(replicates)
```

Four refits of the early window agree with each other at a Mantel correlation
of 0.99, and the six pairwise comparisons span 0.990 to 0.991. That is the
floor. Now the two windows:

```{r shift}
mantel_test(represent(forest_on(early, 1)),
            represent(forest_on(late, 2)), n_perm = 999)
```

0.86, against a floor of 0.99. The two forests are far from independent, which
is the expected part: they are fitted to the same kind of borrower and agree on
the broad ordering. But the gap is many times the spread between refits, so the
answer to the question is that the forest does reorganise, and the monitoring
that would catch it is a Mantel correlation against last year's model rather
than a comparison of accuracies, which would have shown nothing.

## 3. What is the protected attribute worth to the model that never saw it?

`applicant_group` was not in the formula. It affects the default probability
nowhere in the generating process. It does shift the income and the score,
which are in the formula, so the question is whether the representation the
forest builds separates the groups anyway, and by how much.

[`permanova()`] partitions the variation in the proximity across the terms of
a design, so the outcome and the attribute can be entered together and the
attribute read after the outcome is accounted for:

```{r fairness}
permanova(px, ~ default + applicant_group, data = book, n_perm = 999)
```

The outcome accounts for 3.1 per cent of the variation in how the forest
arranges the borrowers. The attribute the model never received accounts for a
further 0.4 per cent, and that share is reliably not zero: across 999
permutations of the borrowers, at most one reproduced it.

Read the two numbers together and neither alone. 0.4 per cent is small, and a
diagnostic that called it alarming would be worthless the first time it met a
real portfolio. What it is not is zero, and it is measured on a model that was
never given the attribute, on data where the attribute causes nothing. This is
the leakage that comes free with correlated predictors, and its size is the
thing worth tracking between refits.

## What this does not show

The proximity is a description of the model, not of the borrowers. Everything
above is a statement about how this forest arranges this book, and a different
forest on the same book will arrange it differently, which is what the first
part of question 2 measures. None of it is a statement about who defaults, and
none of it is a fairness certificate: a `permanova()` that reports 0.4 per cent
has measured one thing about one representation, and the question of whether a
lending decision is fair is not settled by a percentage.