Regressione logistica in R: glm() con family = binomial

Quando l'esito è sì/no

La regressione lineare prevede un numero. Ma molte delle domande che vale la pena modellare sono binarie: il cliente abbandona, il paziente guarisce, l'email viene cliccata. Adattare una retta a un esito 0/1 si rompe subito: la retta prevede tranquillamente probabilità di −0.3 o 1.4, che non hanno senso.

La regressione logistica risolve il problema modellando la probabilità dell'esito attraverso la trasformazione in log-odds (logit): log(p / (1 − p)) = intercept + slope × x. La scala dei log-odds copre l'intera retta dei numeri, quindi un'equazione lineare ci sta in modo naturale, e tornando indietro ogni previsione viene compressa in (0, 1) lungo la familiare curva a S. Il prezzo del trucco: i coefficienti vivono sulla scala dei log-odds, e tutta l'arte di leggere una regressione logistica sta nel ritradurli in qualcosa che le persone capiscono.

Stimare il modello: glm() con family = binomial

glm() (modello lineare generalizzato) è il fratello maggiore di lm(); family = binomial seleziona la regressione logistica. In mtcars, am registra il tipo di cambio (1 = manuale, 0 = automatico): le auto che consumano poco tendono ad avere il cambio manuale?

fit <- glm(am ~ mpg, data = mtcars, family = binomial)
summary(fit)

Due cose prima di leggere l'output. Primo, family = binomial non è facoltativo: se lo ometti, glm() stima in silenzio i minimi quadrati ordinari. Secondo, l'esito deve essere binario: 0/1, logico o un factor a due livelli (R modella la probabilità del secondo livello).

Ora il summary, blocco per blocco:

• Coefficients: l'Estimate di mpg è circa 0.31, ed è una pendenza in log-odds: ogni mpg in più aggiunge 0.31 ai log-odds di avere il cambio manuale. Positivo significa "alza la probabilità", negativo significa "la abbassa"; oltre al segno, l'intuizione di nessuno funziona su questa scala, ed è per questo che esiste la prossima sezione.
• z value e Pr(>|z|): stessa logica dei t-test della regressione (Estimate ÷ Std. Error, poi un p-value per "potrebbe essere zero?"), solo con un'approssimazione normale, da cui z invece di t. Qui p ≈ 0.011: è improbabile che l'associazione tra mpg e tipo di cambio sia rumore.
• Null deviance e Residual deviance: la devianza è la misura di cattivo adattamento nel mondo dei glm (più piccola è meglio). La null deviance (43.2 su 31 df) è il modello con la sola intercetta; la residual deviance (29.7 su 30 df) è il tuo. Il calo di circa 13.6 per il costo di un predittore in df è l'analogo glm di "l'R quadro è salito".
• AIC: un punteggio per confrontare modelli che bilancia adattamento e complessità; vince il più basso. Da solo non significa nulla, è utile tra modelli candidati sugli stessi dati.

Dai log-odds agli odds ratio: exp(coef())

L'esponenziale sposta i coefficienti dalla scala additiva dei log-odds alla scala moltiplicativa degli odds:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)

print(exp(coef(fit)))
print(exp(confint(fit)))

exp(0.307) ≈ 1.36, ed ecco la frase onesta da memorizzare: "ogni mpg in più moltiplica gli odds di un cambio manuale per circa 1.36." Un odds ratio sopra 1 alza gli odds, sotto 1 li abbassa, esattamente 1 significa nessun effetto; ecco perché il verdetto dell'intervallo di confidenza per gli odds ratio è "l'intervallo esclude 1?" (non zero; lo zero era il confine sulla scala dei log-odds).

Attenzione al linguaggio: gli odds non sono probabilità. Odds = p / (1 − p), quindi una probabilità di 0.75 corrisponde a odds di 3. Moltiplicare gli odds per 1.36 non equivale a moltiplicare la probabilità per 1.36, e quando l'esito è frequente la differenza è grande. Un odds ratio di 2 per un esito raro si comporta come "circa il doppio del rischio"; per un esito che si verifica nel 50% dei casi, decisamente no. Non riportare mai un odds ratio con le parole del rischio relativo ("1.36 volte più probabile") a meno che l'esito non sia raro.

Probabilità previste: il tranello di type = "response"

Il bug di regressione logistica più comune in circolazione:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)
new_cars <- data.frame(mpg = c(15, 20, 25, 30))

print(predict(fit, newdata = new_cars))                      # log-odds!
print(predict(fit, newdata = new_cars, type = "response"))   # probabilities

La prima chiamata restituisce il default type = "link": previsioni sulla scala dei log-odds, valori negativi compresi. La seconda restituisce le probabilità vere e proprie. Se le tue "probabilità" escono negative o sopra 1, il motivo è questo. Esegui il blocco: un'auto da 15 mpg non ha praticamente nessuna probabilità di avere il cambio manuale, un'auto da 30 mpg molto probabilmente ce l'ha, e la curva a S si piega nel mezzo.

Classificazione: soglia e tabella di confusione

Le probabilità diventano classi previste scegliendo una soglia, con 0.5 come scelta di default, e la pagella onesta è una tabella di previsto contro reale:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)

probs <- predict(fit, type = "response")   # fitted probabilities, one per car
predicted <- probs > 0.5

print(table(actual = mtcars$am, predicted = predicted))

Le celle sulla diagonale sono le previsioni corrette; le due celle fuori diagonale sono i due errori diversi (prevedere manuale per un'automatica, e il contrario). L'accuratezza complessiva da sola può lusingare molto un modello: se il 95% dei clienti non abbandona, "prevedi che nessuno abbandoni" ottiene il 95% senza individuare nessuno di quelli che abbandonano, quindi guarda sempre entrambi i tipi di errore. E 0.5 è una convenzione, non una legge: quando i due errori hanno costi diversi, sposta la soglia di conseguenza.

Un'avvertenza onesta: questa tabella valuta il modello sugli stessi dati su cui è stato stimato, e questo lo lusinga. Una valutazione vera tiene da parte dati che il modello non ha mai visto.

Più predittori

Esattamente come con lm(): aggiungi termini con +, e ogni interpretazione acquista la precisazione "a parità degli altri":

fit2 <- glm(am ~ mpg + qsec, data = mtcars, family = binomial)
summary(fit2)
print(exp(coef(fit2)))

Ogni coefficiente elevato a esponente è ora il moltiplicatore degli odds per un aumento di un'unità di quel predittore tra auto uguali negli altri predittori. Il meccanismo si estende, ma si estendono anche le avvertenze della regressione lineare: predittori correlati si rimescolano a vicenda i coefficienti.

Avvertenze

• Separazione completa. Se un predittore divide perfettamente l'esito (ogni auto sopra un certo mpg è manuale, ogni auto sotto è automatica), il coefficiente di massima verosimiglianza tende all'infinito. R avvisa con glm.fit: fitted probabilities numerically 0 or 1 occurred e riporta coefficienti enormi con errori standard assurdi. Non usare quei numeri; semplifica il modello, raccogli più dati o usa un metodo penalizzato (i pacchetti brglm2 o logistf).
• Abbastanza eventi. Il vincolo che conta è il numero di casi dell'esito più raro, non il totale delle righe. Una vecchia regola pratica chiede nell'ordine di 10-15 eventi per predittore; gli esempi con 32 auto servono a insegnare il meccanismo, non sono un modello per campioni da pubblicazione.
• Gli odds ratio non sono rischi relativi quando l'esito è frequente: lo abbiamo visto sopra, lo ripetiamo perché i revisori se ne accorgeranno anche se tu non te ne accorgi.

Cosa ti porti a casa

• Esito binario → glm(y ~ x, data = df, family = binomial); non dimenticare mai la family.
• I coefficienti grezzi sono log-odds; exp(coef(fit)) dà gli odds ratio, e il valore nullo per i loro intervalli è 1.
• La frase tipo: "ogni aumento di un'unità di x moltiplica gli odds dell'esito per exp(b)."
• predict(..., type = "response") per le probabilità: il default restituisce log-odds, la confusione numero uno.
• Classifica con una soglia e giudica con una tabella di confusione; l'accuratezza da sola può mentire.
• Fai attenzione agli avvisi di separazione, conta i tuoi eventi e non travestire gli odds ratio da rischi relativi.

Prossimo passo: il meccanismo dietro ogni intervallo visto finora, cioè gli intervalli di confidenza con t.test(), confint() e prop.test().