Forklarbar maskinlæring

En modell for matlevering predikerer 41 minutter for én bestilling. Koordinatoren spør hvorfor akkurat denne bestillingen tar så lang tid. Utvikleren spør hvilke features modellen lener seg mest på generelt. En revisor spør om noen bydeler systematisk får lengre estimater. Det er tre ulike spørsmål. Fortolkbarhet, hvor godt et menneske kan forstå sammenhengen mellom input og modellens oppførsel, må derfor alltid knyttes til et konkret spørsmål og en målgruppe.

En lokal forklaring gjelder én prediksjon eller et lite område rundt den: «denne bestillingen fikk 9 minutter ekstra på grunn av avstanden». En global forklaring oppsummerer oppførselen over mange observasjoner: «avstand er den viktigste featuren, og lange avstander gir lengre tid». For en lineær modell gir de to samme bilde av hvordan en feature virker, fordi helningen er lik overalt. For en ikke-lineær modell med vekselvirkninger kan samme feature trekke opp for én observasjon og ned for en annen.

Det får konsekvenser når lokale forklaringer slås sammen til et globalt bilde. Snittet med fortegn kan bli nær null selv om featuren er avgjørende overalt, fordi positive og negative bidrag kansellerer hverandre. Til viktighet bruker vi derfor snittet av absoluttverdiene, og til retning ser vi på hele fordelingen. En forsiktig global analyse oppsummerer lokale forklaringer over mange representative observasjoner i stedet for å la én typisk observasjon stå for alle.

Modellen under er $f=2x_1{x}_2$, der $x_1$ bare virker sammen med $x_2$.

import numpy as np

rng = np.random.default_rng(16)
X = rng.uniform(-1, 1, size=(500, 2))   # bakgrunnsdata med to features
# 1 = ren vekselvirkning 2·x1·x2, 0 = additiv modell 2·x1 + x2
samspill = 1.0  # @param 0..1 0.1

def f(X):
    return samspill * 2 * X[:, 0] * X[:, 1] + (1 - samspill) * (2 * X[:, 0] + X[:, 1])

def effekt_x1(X, h=1e-4):
    Xh = X.copy()
    Xh[:, 0] += h
    return (f(Xh) - f(X)) / h           # lokal helning for x1 i hvert punkt

A = np.array([[0.5, 0.8]])
B = np.array([[0.5, -0.8]])
print(f"lokal effekt av x1 i A: {effekt_x1(A)[0]:+.2f}")
print(f"lokal effekt av x1 i B: {effekt_x1(B)[0]:+.2f}")
alle = effekt_x1(X)
print(f"globalt snitt med fortegn: {alle.mean():+.2f}")
print(f"globalt snitt av |effekt|: {np.abs(alle).mean():.2f}")

En krum overflate ser flat ut hvis du zoomer langt nok inn, slik jorda virker flat når du står på den. LIME (Local Interpretable Model-agnostic Explanations) bygger på den idéen: rundt ett punkt kan selv en komplisert modell tilnærmes med en rett linje. Metoden trenger bare å kunne spørre modellen om prediksjoner, så den virker på alle modelltyper.

Oppskriften har fem trinn. Velg observasjonen $x$ som skal forklares. Lag mange perturbasjoner $z$ (perturbasjon), kunstige varianter av $x$ der featureverdiene er endret litt. Spør svart-boksen $f$ om prediksjonen for hver av dem. Vekt hver nabo etter nærhet til $x$ med ƒLIME-nærhetsvekt, slik at nære naboer teller mest. Tilpass til slutt en enkel, tolkbar modell $g$, for eksempel lineær regresjon med få features, med vektet minste kvadrat. Hele målet er ƒLIME-mål: $g$ skal ligne $f$ lokalt, og en straff $\Omega (g)$ holder forklaringen enkel.

Den enkle modellen er en surrogatmodell: den trenes på svart-boksens output, ikke på den ekte fasiten. Hvor godt den etterligner svart-boksen i det vektede nabolaget, kalles fidelity. Høy fidelity betyr bare at forklaringen beskriver modellen godt. Den sier ingenting om at modellen er riktig, rettferdig eller kausal.

LIME kan være ustabil. Kernelbredden ${\sigma }^2$ bestemmer hvor stort nabolaget er, og et for stort nabolag blander flere krumme mønstre. Et annet seed gir andre naboer, og naive perturbasjoner kan lage umulige kombinasjoner, som en 15-åring med 30 års arbeidserfaring. Kjør derfor metoden flere ganger og rapporter fidelity sammen med koeffisientene.

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(4172)
x0 = np.array([0.4, -0.7])
naboer = x0 + rng.normal(0, 0.55, size=(300, 2))   # perturbasjoner rundt x0
bredde = 0.45  # @param 0.05..5 0.05

def svart_boks(X):
    score = 1.4 * X[:, 0] - 0.9 * X[:, 1] + 0.8 * X[:, 0] * X[:, 1]
    return 1 / (1 + np.exp(-score))

y = svart_boks(naboer)
vekt = np.exp(-np.sum((naboer - x0) ** 2, axis=1) / bredde)   # kernel: nære naboer teller mest
A = np.column_stack([np.ones(len(naboer)), naboer])
koef = np.linalg.lstsq(A * np.sqrt(vekt[:, None]), y * np.sqrt(vekt), rcond=None)[0]
rest = np.sum(vekt * (y - A @ koef) ** 2) / np.sum(vekt * (y - np.average(y, weights=vekt)) ** 2)
print(f"helninger {np.round(koef[1:], 3)}, lokal fidelity (vektet R²) {1 - rest:.3f}")

t = np.linspace(-1.2, 2.0, 100)                     # snitt langs feature 1 gjennom x0
snitt = np.column_stack([t, np.full_like(t, x0[1])])
plt.figure()
plt.plot(t, svart_boks(snitt), label="svart boks")
plt.plot(t, np.column_stack([np.ones_like(t), snitt]) @ koef, label="lokal surrogat")
plt.xlabel("feature 1 (feature 2 holdt fast i x0)")
plt.ylabel("modellens output")
plt.legend()
plt.tight_layout()
plt.show()

Shapley-verdier kommer fra samarbeidsspill. Spillere kan danne koalisjoner $S$, og $v(S)$ er verdien koalisjonen oppnår. Spørsmålet er hvordan den samlede gevinsten $v(N)-v(\varnothing )$ skal fordeles rettferdig. Det vanskelige er vekselvirkninger: en spiller kan bidra lite alene og mye sammen med en bestemt annen spiller. Hvor mye spilleren bidro, avhenger da av hvem som allerede var med.

Svaret er å ta snittet over alle rekkefølger spillerne kan komme inn i. For hver rekkefølge måler vi spillerens marginalbidrag: hvor mye verdien øker når spilleren slutter seg til dem som allerede er der. En shapley-verdi er gjennomsnittet over alle $n!$ rekkefølger, ƒshapley-verdi via rekkefølger. Den samme summen kan skrives over delmengder, ƒshapley-verdi, der vekten $\frac{|S|!(n-|S|-1)!}{n!}$ er andelen rekkefølger som har akkurat $S$ foran spilleren. For $n=3$ blir vektene $\frac{1}{3}$, $\frac{1}{6}$ og $\frac{1}{3}$ for $|S|=0$, 1 og 2.

Fordelingen er rettferdig i en presis forstand. Effektivitet, ƒeffektivitet: summen av Shapley-verdiene er nøyaktig $v(N)-v(\varnothing )$, så ingenting forsvinner og ingenting legges til. Symmetri: to spillere som bidrar likt i alle koalisjoner, får samme verdi. Nullspiller: en spiller som aldri endrer verdien, får 0. Effektiviteten er også den beste kontrollen når du regner for hånd.

I maskinlæring er spillerne features, og $v(S)$ er modellens output når bare featurene i $S$ har sine faktiske verdier, mens resten er erstattet av en referanse.

Figuren under er et lite spill for en modell som predikerer leveringstid: $v(S)$ er hvor mange minutter over referansen modellen gir når bare featurene i $S$ er med.

Koden regner de samme verdiene med rå kraft over alle rekkefølger.

from itertools import permutations
from math import factorial

spillere = ["Avstand", "Regn", "Ledige bud"]
v_avstand_regn = 16  # @param 0..30
# v(S): ekstra minutter når bare featurene i S har sin faktiske verdi
v = {(): 0, ("Avstand",): 10, ("Regn",): 4, ("Ledige bud",): -2,
     ("Avstand", "Regn"): v_avstand_regn, ("Avstand", "Ledige bud"): 6,
     ("Regn", "Ledige bud"): 2, ("Avstand", "Regn", "Ledige bud"): 12}

def verdi(S):
    return v[tuple(p for p in spillere if p in S)]

phi = {p: 0.0 for p in spillere}
for rekkefolge in permutations(spillere):   # alle 3! = 6 rekkefølger
    S = set()
    for p in rekkefolge:
        phi[p] += verdi(S | {p}) - verdi(S)  # marginalbidraget til p
        S.add(p)

n_fak = factorial(len(spillere))
for p in spillere:
    print(f"{p:11s} φ = {phi[p] / n_fak:6.2f}")
print(f"sum φ = {sum(phi.values()) / n_fak:.2f}, v(alle) - v(tom) = {verdi(spillere) - verdi([])}")

SHAP bruker Shapley-verdier til å forklare én prediksjon $f(x)$. Spillerne er featurene, og $v(S)$ er modellens forventede output når featurene i $S$ har verdiene fra $x$ og resten hentes fra et sett med bakgrunnsdata. Resultatet er en additiv forklaring, ƒadditiv forklaringsmodell:

Startpunktet ${\varphi }_0$ er baseline (baseline (SHAP)), typisk modellens gjennomsnittlige prediksjon over bakgrunnsdataene. Hver ${\varphi }_j$ er en featureattribusjon: hvor mye feature $j$ flytter denne prediksjonen bort fra baseline. Positiv verdi trekker opp, negativ trekker ned. Fortegnet er relativt til baseline. Det betyr ikke «bra» eller «dårlig», og det betyr ikke at featureverdien er lav. For klassifikatorer forklarer SHAP ofte log-odds og ikke sannsynlighet, så oppgi alltid hvilken output som forklares.

For en lineær modell med uavhengige features blir SHAP-verdiene eksakt ƒSHAP for lineær modell. En feature som har gjennomsnittsverdien, får bidrag 0, uansett hvor stor vekten er. Bidraget måler avvik fra referansen, ikke vekten alene.

Attribusjonene er heller ikke modellens interne vekter. I et nevralt nettverk går hver feature gjennom mange ikke-lineære veier, og ${\varphi }_j$ oppsummerer nettoeffekten for akkurat denne observasjonen.

Baseline er et valg. Bytter du bakgrunnsdata, endres ${\varphi }_0$ og dermed bidragene, men summen treffer fortsatt prediksjonen.

import itertools
import math
import numpy as np

x = np.array([1.2, -0.6, 0.8])
b0 = 0.0  # @param 0..1 0.1
b2 = 0.0  # @param 0..1 0.1
baseline = np.array([b0, 0.0, b2])            # bakgrunnsverdier for x0, x1, x2

def modell(v):
    return 0.25 + 0.7 * v[0] - 0.4 * v[1] + 0.3 * v[0] * v[2]

d = len(x)
phi = np.zeros(d)
for j in range(d):
    andre = [k for k in range(d) if k != j]
    for r in range(d):
        for S in itertools.combinations(andre, r):
            v = baseline.copy()
            v[list(S)] = x[list(S)]           # featurene i S får sin faktiske verdi
            uten = modell(v)
            v[j] = x[j]                       # ... og så legger vi til j
            vekt = math.factorial(r) * math.factorial(d - r - 1) / math.factorial(d)
            phi[j] += vekt * (modell(v) - uten)

print(f"baseline-output φ0 = {modell(baseline):.3f}")
print(f"bidrag φ = {np.round(phi, 3)}")
print(f"φ0 + sum φ = {modell(baseline) + phi.sum():.3f}, modell(x) = {modell(x):.3f}")

Ett SHAP-plott per prediksjon blir uoversiktlig når du har tusenvis av observasjoner. Et SHAP-sammendragsplott, ofte kalt beeswarm-plott, viser alle på én gang. Hver rad er en feature, og hver prikk er én observasjon. Tre regler er alt du trenger:

• •Rekkefølge: featurene er sortert etter gjennomsnittlig $|\varphi |$, med den viktigste øverst.
• •x-aksen: SHAP-verdien, altså hvor mye featuren flyttet den prediksjonen. Til høyre for null betyr opp, til venstre betyr ned. Aksen viser ikke featureverdien.
• •Farge: featureverdien, fra blå (lav) til rød (høy). Ligger de røde prikkene til høyre, trekker høye verdier opp. Ligger de til venstre, trekker høye verdier ned.

Formen sier også noe. En bred rad betyr at featuren har stor effekt for noen observasjoner. En lang hale på én side betyr sjeldne, men store effekter. Blandede farger på samme side av null betyr at sammenhengen ikke er monoton, eller at featuren virker sammen med andre.

Plottet over kommer fra en syntetisk modell for leveringstid. Avstand står øverst, med røde prikker til høyre: lange avstander gir mer tid. Hos Ledige bud ligger de røde prikkene litt til venstre og de blå langt til høyre: mange ledige bud trekker litt ned, få bud trekker kraftig opp. Høy verdi trekker altså ned. Regn har mest blå prikker nær null og noen få røde langt til høyre: regn betyr sjelden noe, men mye når det først regner. Klokkeslett har både røde og blå prikker til venstre og noen midtverdier til høyre, så effekten er ikke monoton. Ordrestørrelse ligger tett rundt null og er minst viktig.

Koden gjør den samme lesingen automatisk for en liten SHAP-matrise: rangering etter snittet av $|\varphi |$ og retning ut fra fargen.

import numpy as np

navn = ["Avstand", "Ledige bud", "Regn", "Ordrestørrelse"]
# SHAP-verdier i minutter for seks leveranser (én kolonne per feature)
phi = np.array([[9.0, -4.0, 0.5, 1.0], [3.0, 6.0, 4.0, -0.5], [-5.0, -3.5, 0.0, 0.5],
                [7.5, 5.0, -1.0, 0.0], [-6.0, -2.0, 3.5, -1.0], [1.0, 4.5, 0.0, 0.5]])
# Normaliserte featureverdier, 0 = lav og 1 = høy (fargen i plottet)
verdi = np.array([[0.9, 0.7, 0.3, 0.8], [0.6, 0.1, 0.9, 0.2], [0.1, 0.8, 0.2, 0.6],
                  [0.8, 0.2, 0.1, 0.4], [0.0, 0.6, 0.8, 0.1], [0.5, 0.0, 0.2, 0.7]])
snu_regn = 0  # @param 0..1
if snu_regn:
    phi[:, 2] *= -1                     # bytt fortegn på hele Regn-kolonnen

viktighet = np.abs(phi).mean(axis=0)
for j in np.argsort(-viktighet):        # viktigst øverst, som i plottet
    r = np.corrcoef(verdi[:, j], phi[:, j])[0, 1]
    retning = "høy verdi trekker opp" if r > 0 else "høy verdi trekker ned"
    print(f"{navn[j]:15s} snitt |φ| = {viktighet[j]:4.2f}   {retning} (r = {r:+.2f})")

SHAP og LIME forklarer modellen, ikke verden. Et bidrag på +3 minutter fra «ordre siste time» betyr at modellens prediksjon steg fordi featuren hadde denne verdien. Det betyr ikke at leveransen ville gått 3 minutter fortere om vi kunne endret verdien. Modellen har lært samvariasjoner fra data, og en feature kan samvariere med utfallet fordi begge påvirkes av en tredje ting, eller fordi den er en proxy for noe annet. Postnummer kan for eksempel bære informasjon om inntekt eller boligtype. En presis forklaring av en slik modell viser at den bruker proxyen, men gjør ikke bruken riktig eller rettferdig.

Korrelerte features gir et ekstra problem. Når to features bærer nesten samme informasjon, kan æren fordeles på mange måter som gir like gode prediksjoner. Én modell kan bruke bare den første, en annen kan dele likt, og forklaringene blir helt forskjellige selv om prediksjonene er nesten like. Metoden må også bestemme hvordan en feature fjernes. Erstattes den uavhengig av de andre, kan det lage kombinasjoner som aldri finnes i virkeligheten. Betingede varianter respekterer sammenhengen bedre, men svarer på et litt annet spørsmål.

Hva gjør man da? Oppgi baseline, outputskala og bakgrunnsdata. Sjekk stabilitet, kombiner med domenekunnskap, og husk at årsaksspørsmål krever eksperimenter eller en kausal modell. En kontrafaktisk forklaring, den minste endringen som ville snudd beslutningen, svarer på enda et annet spørsmål.

Her er to lineære modeller som gir nesten like prediksjoner fordi x2 er en nesten-kopi av x1.

import numpy as np

rng = np.random.default_rng(3)
n = 200
stoy = 0.05  # @param 0..1 0.05
x1 = rng.normal(0, 1, n)            # ordre siste time (standardisert)
x2 = x1 + rng.normal(0, stoy, n)    # aktive brukere i appen: nesten en kopi
X = np.column_stack([x1, x2])
snitt = X.mean(axis=0)
x = X[0]                            # leveransen vi forklarer

for navn, w in [("A", np.array([2.0, 0.0])), ("B", np.array([1.0, 1.0]))]:
    phi = w * (x - snitt)           # lineær SHAP: w_j (x_j - snitt_j)
    print(f"modell {navn}: f(x) = {x @ w:+.2f}   φ = {np.round(phi, 2)}")

avvik = np.abs(X @ np.array([2.0, 0.0]) - X @ np.array([1.0, 1.0]))
print(f"største forskjell i prediksjon over datasettet: {avvik.max():.2f}")