Soluzione esame del 18 Novembre 2025

Statistics III - CdL SSE

Autore/Autrice

Affiliazione

Tommaso Rigon

Università degli Studi di Milano-Bicocca

Homepage

Parte I: analisi dei dati

library(faraway)
data(beetle)
head(beetle, 4) # Mostro le prime 4 osservazioni

  conc affected exposed
1 24.8       23      30
2 24.6       30      30
3 23.0       29      31
4 21.0       22      30

# Punto (a)
beetle$prop <- beetle$affected / beetle$exposed # Proporzioni
plot(beetle$conc, beetle$prop, pch = 16, xlab = "Concentrazione", ylab = "Proporzione di decessi")

Il grafico mostra una chiara dipendenza positiva tra la concentrazione di ossido di etilene e la proporzione di decessi. Nota: gli abbellimenti estetici, come ad esempio i nomi corretti negli assi o la creazione della nuova variabile prop, erano apprezzati ma non richiesti.

# Punto (b)
m1 <- glm(cbind(affected, exposed - affected) ~ conc, data = beetle, family = binomial(link = "logit"))

La scelta del GLM appropriato in questo caso era una regressione per dati binomiali per la variabile affected (S_i), in cui exposed (m_i) rappresenta il numero di tentativi e conc (x_i) l’unica variabile esplicativa, per i = 1,\dots,10. La variabile risposta Y_i era pertanto pari alla proporzione Y_i = S_i / m_i, la cui distribuzione è stata trattata a lezione.

Un modello di regressione lineare semplice, applicato direttamente alle proporzioni Y_i è un approccio gravemente errato e corrisponde al “linear probability model” menzionato in classe, che è molto problematico per tutte le ragioni discusse a lezione e nel materiale online.

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

(Intercept)        conc 
 -6.0104702   0.3412697 

Otteniamo quindi che l’equazione stimata è \hat{\mathbb{E}}(Y_i) = \hat{\pi}_i = \frac{e^{-6.010 + 0.341 x_i}}{1 + e^{-6.010 + 0.341 x_i}}, \qquad i=1,\dots,10.

# Punto (b.ii)
m0 <- glm(cbind(affected, exposed - affected) ~ 1, data = beetle, family = binomial(link = "logit")) # Modello nullo

# 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 conc. 
# 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 = 67.522 (Wald) e W = 100.3 / \phi = 100.3 (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 8.217, 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. Questo significa che la probabilità di decesso varia con la concentrazione di etilene (conc).

# Punto (c)
newdata <- data.frame(conc = seq(from = 10, to = 25, length = 200))
pred_m1 <- predict(m1, newdata = newdata, type = "response")
plot(beetle$conc, beetle$prop, pch = 16, xlab = "Concentrazione", ylab = "Proporzione di decessi", xlim = c(10, 25), ylim = c(0, 1))
lines(newdata$conc, pred_m1)

# Vi erano innumerevoli sintassi equivalenti per produrre lo stesso risultato. Ecco un'alternativa:
curve(plogis(coef(m1)[1] + coef(m1)[2] * x), from = 10, to = 25, ylab = "Estimated probability", ylim = c(0, 1))
points(beetle$conc, beetle$prop, pch = 16)

Le probabilità previste dal modello sono compatibili con l’andamento medio dei dati, anche se si osserva una variabilità eccessiva rispetto alle aspettative. Questo è un primo segnale, informale, della possibile presenza di sovradispersione.

# Punto (d)
cooks.distance(m1) # Valori influenti (leva + anomali)

         1          2          3          4          5          6          7 
1.50778167 0.43184011 0.22981702 0.01590161 0.30091250 0.87479492 0.03666182 
         8          9         10 
1.15503547 0.11313239 0.47712746 

La 1a osservazione e l’8a osservazione hanno una distanza di Cooks elevata, superiore a 1, mentre la sesta osservazione è borderline (ulteriori analisi mostrano che la 6a osservazione ha residuo molto alto ma valore leva relativamente basso). Valori alti della distanza di Cooks identificano osservazioni influenti, cioè al contempo punti leva e potenzialmente anomale.

Alcuni studenti hanno fatto uso di strumenti grafici diagnostici, come ad esempio plot(m1, which = 4), oppure si sono affidati a librerie esterne, come il comando influencePlot della libreria car. Si tratta di soluzioni corrette, per le quali non ho assegnato penalizzazioni, tuttavia nel caso del comando InfluencePlot vengono utilizzati i residui studentizzati, la cui definizione non è mai stata discussa a lezione nel caso dei GLM e non è nemmeno presente nei libri di testo. Ciò che la funzione rstudent sta calcolando nel caso dei GLM è la seguente “media pesata” tra residui di Pearson e residui di devianza standardizzati: \tilde{r}_{i,\text{stud}} = \text{sign}(y_i - \hat{\mu}_i) \sqrt{(1 - h_{i, W})\tilde{r}_{i, D}^2 + h_{i, W} \tilde{r}_{i, P}^2}, \qquad i=1,\dots,n. Pur fornendo una risposta ragionevole, fare uso di quantità e comandi R di cui non si conoscono i dettagli non è una prassi raccomandabile, soprattutto considerato che nulla di tutto ciò era richiesto dalla prova.

# Punto (e)
exp(coef(m1)[2])

    conc 
1.406733 

Il valore \exp(\hat{\beta}_2) = 1.4067 si può interpretare come rapporto tra quote (odds ratio) relativo all’aumento unitario della concentrazione di etilene. Tale rapporto è maggiore di 1 e abbastanza elevato, ad indicare un marcato aumento della probabilità di decesso al crescere della variabile conc.

# Punto (f)
new15 <- data.frame(conc = 15)
predict(m1, newdata = new15, type = "response") # Probabilità prevista: 0.2908159 

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

# Intervallo di confidenza alla Wald
fit15 <- predict(m1, newdata = new15, type = "response", se.fit = TRUE) 
fit15$fit + c(-1, 1) * qnorm(0.975) * fit15$se.fit

[1] 0.2086652 0.3729667

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

[1] 0.2158957 0.3791630

# Punto (g)
X2 <- sum(residuals(m1, type = "pearson")^2)
X2 # Statistica X2 di Pearson
1 - pchisq(X2, m1$df.residual) # p-value test bontà d'adattamento.
X2 / m1$df.residual # Stima del parametro di dispersione

È presente sovradispersione. Il valore di X^2 = 35.456 è elevato e al di fuori dei valori tipici, come confermato dal p-value associato, praticamente pari a 0. Si ottiene inoltre che \hat{\phi} = X^2 / (n - p) = 4.432.

mquasi <- glm(cbind(affected, exposed - affected) ~ conc, data = beetle, family = quasibinomial(link = "logit"))
summary(mquasi)

Lo stimatore di quasi-verosimiglianza produce le medesime stime per \beta ottenute in precedenza, ma corregge gli errori standard. Fortunatamente, questo non altera le conclusioni inferenziali a cui eravamo giunti allo step (b). Tuttavia, si noti che le osservazioni che in precedenza erano state registrate come “anomale” ed “influenti”, rientrano nella norma una volta che si tiene conto della sovradispersione:

cooks.distance(mquasi)

          1           2           3           4           5           6 
0.340202840 0.097436676 0.051853929 0.003587901 0.067895298 0.197381176 
          7           8           9          10 
0.008272056 0.260612233 0.025526217 0.107654922 

Non si trattava quindi di osservazioni influenti o anomale, bensì probabilmente di un modello non correttamente specificato.

# Punto (h)

# Intervallo di confidenza alla Wald
fit15_quasi <- predict(mquasi, newdata = new15, type = "response", se.fit = TRUE) 
fit15_quasi$fit + c(-1, 1) * qnorm(0.975) * fit15_quasi$se.fit

[1] 0.1178694 0.4637624

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

[1] 0.1505894 0.4867872

Anche in questo caso era possibile considerare uno dei due intervalli di confidenza alternativi (entrambi corretti). Confrontati con quelli ottenuti in precedenza, si nota che la loro ampiezza è aumentata come conseguenza della sovradispersione.