Soluzione esame del 10 Febbraio 2026

Statistics III - CdL SSE

Autore/Autrice

Affiliazione

Tommaso Rigon

Università degli Studi di Milano-Bicocca

Homepage

Parte I: analisi dei dati

library(ISLR2)
data(Default)
head(Default, 4)

  default student   balance   income
1      No      No  729.5265 44361.63
2      No     Yes  817.1804 12106.13
3      No      No 1073.5492 31767.14
4      No      No  529.2506 35704.49

boxplot(Default$balance ~ Default$default)

boxplot(Default$balance ~ Default$student)

Sia default che student risultano essere “correlate” con la variabile balance. Sia gli studenti che i clienti insolventi hanno un maggior debito sulla carta di credito. In particolare, nel prevedere la variabile default tramite le restanti variabili, bisognerà tenere a mente della dipendenza esistente tra i predittori balance e student.

# Punto (b)
m1 <- glm(default ~ student, data = Default, family = "binomial")

# Punto (b.i)
coef(m1) # In alternativa, anche summary(m1) mostrava le stime dei coefficienti

(Intercept)  studentYes 
 -3.5041278   0.4048871

Otteniamo quindi che l’equazione stimata è \hat{\mathbb{E}}(Y_i) = \hat{\pi}_i =\frac{e^{\hat{\beta}_1 + \hat{\beta}_2 d_i}}{1 + e^{\hat{\beta}_1 + \hat{\beta}_2 d_i}} =\frac{e^{-3.504 + 0.405 d_i}}{1 + e^{-3.504 + 0.405 d_i}}, \qquad i=1,\dots,10000, dove d_i è una variabile dummy che vale 1 quando l’i-esimo cliente è uno studente e 0 altrimenti.

# Punto (b.ii)
m0 <- glm(default ~ 1, data = Default, family = "binomial")

# Statistiche test, p-value e gradi di libertà
lmtest::waldtest(m0, m1, test = "Chisq") # Test di Wald
anova(m0, m1, test = "LRT") # Test log-rapporto di verosimiglianza

# Il test di Wald, più semplicemente, corrisponde allo z-value relativo al coefficiente studentYes. 
# Si ricordi che z-value^2 = Chisq. 
summary(m1)

Il valore delle statistiche test per i test di Wald e log-rapporto di verosimiglianza sono rispettivamente pari a: W_e = 12.392 (Wald) e W = 11.967 / \phi = 11.967 (log-rapporto di verosimiglianza), poichè \phi = 1 nel modello binomiale. Nel caso del test di Wald, era ugualmente corretto riportare la statistica test z-value presente nel summary, pari a 3.52, il cui quadrato coincide con W_e a meno di piccole discrepanze numeriche.

Come riportato nelle tabelle, i gradi di liberà di entrambi i test è q = 1, perché è uno solo il parametro vincolato nel sistema d’ipotesi. Quest’ultimo è infatti pari a H_0: \beta_2 = 0 e H_1: \beta_2 \neq 0.

Considerati i gradi di libertà dei test W_e e W e la loro distribuzione asintotica sotto H_0, si ottiene in entrambi i casi un p-value quasi pari a 0. In altri termini, si rigetta l’ipotesi nulla che \beta_2 = 0.

# Punto (b.iii)
exp(coef(m1)[2]) # Rapporto tra quote (odds ratio)

studentYes 
  1.499133

Il rapporto tra le quote nel caso in cui il cliente sia studente e nel caso in cui non lo sia è pari a 1.499. Questo significa che, marginalmente rispetto alle altre covariate (ci torniamo), essere studente è associato a una probabilità di insolvenza maggiore. In termini di rapporto tra quote, la quota di solvenza per uno studente è circa una volta e mezzo la quota per un non-studente.

Un errore ricorrente (che commettono anche i giornalisti del Daily Mail)

Ho notato un errore ricorrente tra gli studenti relativo all’interpretazione degli odds ratio, in questo caso pari a 1.499. Questo indica una variazione percentuale del 49.9\% delle quote, e non delle probabilità. Come raccontato anche a Statistica I, si tratta di quantità ben distinte.

probs <- c(plogis(sum(coef(m1))), plogis(coef(m1)[1])) # Alternativamente: predict(m1, data.frame(student = c("Yes", "No")), type = "response")
probs

            (Intercept) 
 0.04313859  0.02919501

È quindi chiaro che la probabilità di insolvenza di uno studente è maggiore di quella di un cliente non studente. Il rapporto tra probabilità (rischio relativo) è quindi pari a probs[1] / probs[2] = 1.477601, ovvero una variazione percentuale delle probabilità del 47\%.

Se invece si considerano le quote:

odds <- probs / (1 - probs)
odds

            (Intercept) 
 0.04508342  0.03007299

Il rapporto tra quote (odds ratio) è quindi pari a odds[1] / odds[2] = 1.499133, cioè una variazione percentuale del 49.9\%, che è appunto pari al valore e^{\hat{\beta}_2} calcolato prima.

La domanda chiedeva se essere studente “aumenti” il rischio di insolvenza e questa era la parte più difficile e delicata, perché ci si chiede se esista o meno una relazione di causa-effetto. Fino a qui abbiamo registrato una significativa associazione tra rischio di insolvenza ed essere studenti. Come noto da Statistica I, associazione non implica causalità e in questo caso c’erano buone ragioni per ritenere che l’associazione fosse dovuta ad una terza variabile (balance), come vedremo nei punti successivi.

# Punto (c)
m2 <- glm(default ~ student + income + balance, data = Default, family = "binomial")

# Punto (c.i)
round(coef(m2), 3)

(Intercept)  studentYes      income     balance 
    -10.869      -0.647       0.000       0.006

Otteniamo quindi che l’equazione stimata è \hat{\mathbb{E}}(Y_i) = \hat{\pi}_i =\frac{e^{-10.869 -0.647 d_i + 0.000 x_{i1} + 0.006 x_{i2}}}{1 + e^{-10.869 -0.647 d_i + 0.000 x_{i1} + 0.006 x_{i2}}}, \qquad i=1,\dots,10000, dove x_{i1} e x_{i2} rappresentano i valori delle variabili income e balance, rispettivamente. L’esponenziale di ciascuno di questi coefficienti è interpretabile in termini di rapporti tra quote, tuttavia si rende necessaria una discussione dettagliata.

round(exp(coef(m2)[-1]), 3)

studentYes     income    balance 
     0.524      1.000      1.006

L’interpretazione dell’intercetta \hat{\beta}_1 = -10.869 non è di particolare interesse, perchè è associata la probabilità di insolvenza per un non-studente che non abbia nè reddito nè debiti. Quest’ultima risulta essere sostanzialmente zero (plogis(-10.869) è praticamente 0).
Il valore e^{\hat{\beta}_3} = 1, relativo alla variabile income, è interpretabile come rapporto tra quote ed è indistinguibile dal valore 1, come oltretutto confermato dal un test statistico formale fornito nel summary. In altri termini, il reddito non è associato a maggiore o minore rischio di insolvenza. Come vedremo in un punto successivo, questa variabile poteva essere omessa.
Il valore e^{\hat{\beta}_2} = 0.524, relativo alla variabile student, è interpretabile come rapporto tra quote ed è minore di 1 (nel modello m1 era invece maggiore di 1). Si noti inoltre che e^{\hat{\beta}_2} è significativamente minore di 1, come si desume dal p-value associato. Questo significa che, a parità di debito e reddito (balance e income), essere studente è associato a una probabilità di insolvenza minore rispetto ad un non studente. Questa conclusione è solo apparentemente in contraddizione con quanto illustrato dal modello m1 ed è un esempio del cosiddetto paradosso di Simpson, nel quale l’aggiunta di una covariata modifica la direzione o l’entità di un’associazione marginale. In altri termini, non è tanto la condizione di studente in quanto tale a comportare un rischio di insolvenza, quanto piuttosto il debito accumulato sulla carta di credito, che però a sua volta è, in media, maggiore tra gli studenti.
Il valore e^{\hat{\beta}_4} = 1.006, relativo alla variabile balance, è interpretabile come rapporto tra quote ed è maggiore di 1. In altre parole, a un aumento unitario (1\$) della variabile balance la quota di insolvenza è 1.006 volte la quota originaria. Considerata la scala della variabile di riferimento, questo è un aumento piuttosto rilevante. Per esempio, un aumento del debito medio di 500\$ comporta un rapporto tra quote pari a e^{500 \hat{\beta}_4} = 17.606, cioè una quota di insolvenza di circa 18 volte maggiore di quella originaria.

# Punto (d)
anova(m1, m2)

Analysis of Deviance Table

Model 1: default ~ student
Model 2: default ~ student + income + balance
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1      9998     2908.7                          
2      9996     1571.5  2   1337.1 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Poichè si tratta di modelli annidati è possibile condurre un test statistico formale, il quale rigetta l’ipotesi nulla. Criteri di informazione quali AIC e BIC erano forse strumenti meno ortodossi ma comunque conducevano alla medesima conclusione.

# Punto (e)
summary(m2)

Come evidente dal summary, la variabile income è rimovibile dal modello senza perdita informativa. Il p-value risultava infatti pari a 0.711.

Poichè il testo non lo specificava, per gli intervalli di confidenza c’erano due possibili soluzioni, entrambe corrette:

# Punto (f)
newdata <- data.frame(student = "Yes", income = mean(Default$income), balance = 100)

fit <- predict(m2, newdata = newdata, type = "response", se.fit = TRUE) 
fit$fit # Previsione

           1 
1.958948e-05

# Intervallo di confidenza alla Wald
fit$fit + c(-1, 1) * qnorm(0.975) * fit$se.fit

[1] 3.679883e-06 3.549907e-05

# Intervallo di confidenza "Wald trasformato"
fit_lin <- predict(m2, newdata = newdata, se.fit = TRUE) 
plogis(fit_lin$fit + c(-1, 1) * qnorm(0.975) * fit_lin$se.fit)

[1] 8.695777e-06 4.412974e-05

# Punto (g)

# NON BISOGNAVA FARE NULLA

No, non ha senso verificare se nel modello m2 è presente sovradispersione, perchè la variabile risposta è binaria e non può essere raggruppata in nessuna maniera a causa della presenza delle variabili continue balance e income.