Soluzione esame del 27 Gennaio 2026

Statistics III - CdL SSE

Autore/Autrice

Affiliazione

Tommaso Rigon

Università degli Studi di Milano-Bicocca

Homepage

Parte I: analisi dei dati

library(faraway)
data(anaesthetic)
head(anaesthetic, 4)

  breath tgrp
1      3    A
2      6    B
3      3    C
4      4    D

# Punto (a)
anaesthetic$breath[anaesthetic$breath == 0] <- 0.25

Questo passaggio era necessario per lo svolgimento delle domande successive, ma non conferiva di per sè alcun punteggio.

# Punto (b)
boxplot(anaesthetic$breath ~ anaesthetic$tgrp) # Boxplot (non mostrato nel seguito)

tapply(anaesthetic$breath, anaesthetic$tgrp, mean) # Media di ciascun gruppo

     A      B      C      D 
5.4000 3.2500 3.0625 4.3625 

tapply(anaesthetic$breath, anaesthetic$tgrp, sd) # Deviazione standard di ciascun gruppo

       A        B        C        D 
5.030434 2.953588 3.107583 3.064520 

Sulla base dei boxplot e delle statistiche descrittive calcolate, le differenze di efficacia (breath) tra i diversi metodi di anestesia non appaiono particolarmente marcate. Per confermare questa valutazione preliminare, è tuttavia necessario ricorrere a un test statistico formale.

# Punto (c)
m_lin <- lm(breath ~ tgrp, data = anaesthetic)

Si tratta di un modello di tipo ANOVA, un concetto che è stato trattato anche nel corso Analisi Statistica Multivariata. La risposta alla gran parte delle domande del punto (c) si poteva ottenere leggendo con attenzione i risultati del comando summary.

summary(m_lin)

Anzitutto, desumiamo che l’equazione che esprime la risposta media in funzione della variabile tgrp è la seguente \hat{\mathbb{E}}(Y_i) = \hat{\beta}_1 + \hat{\beta}_2 d_{i,B} + \hat{\beta}_3 d_{i, C} + \hat{\beta}_4 d_{i, D} = 5.4 - 2.15 \, d_{i,B} -2.34 \, d_{i, C} -1.04\, d_{i, D},\qquad i=1,\dots,80. dove d_{i, B}, d_{i, C} e d_{i, D} sono delle variabili dummy che identificano l’appartenenza dell’i-esima unità al gruppo B, C, D della variabile tgrp; il primo gruppo (A) è preso come categoria di riferimento. Di conseguenza, il valore \hat{\beta}_1 rappresenta la media del gruppo A, mentre i rimanenti coefficienti rappresentano la differenza tra la media del gruppo considerato e la media del gruppo A.

Di conseguenza, le previsioni per ciascuno dei quattro gruppi coincidono con le rispettive medie di gruppo. Non era necessario eseguire ulteriori comandi, poiché tali medie erano già state calcolate. In effetti, il comando predict produce proprio i valori delle medie.

predict(m_lin, newdata = data.frame(tgrp = c("A", "B", "C", "D")))

     1      2      3      4 
5.4000 3.2500 3.0625 4.3625 

Veniva infine richiesto di confrontare il modello stimato con il modello nullo, ossia di verificare l’ipotesi
H_0: \beta_2 = \beta_3 = \beta_4 = 0, contro l’alternativa H_1, secondo cui almeno uno dei coefficienti è diverso da zero. In questo contesto, tale sistema di ipotesi equivale a verificare se le medie dei gruppi siano uguali tra loro oppure no. Tradotto nel linguaggio del problema, si tratta di accertare, mediante un test statistico formale, se in media esistano differenze di efficacia tra le quattro differenti anestesie. Anche in questo caso, la risposta si otteneva tramite il summary, da cui si deduce in particolare che:

F-statistic: 1.774 on 3 and 76 DF,  p-value: 0.1592

In altri termini, la statistica test F è pari a 1.774 avente gradi di libertà pari a q = 3 (numero di parametri vincolati) e n - p = 76. Il p-value è pari a 0.1592, pertanto non c’è sufficiente evidenza contro l’ipotesi nulla. In altre parole, concludiamo che in media i trattamente anestetici sono indistinguibili tra loro e che eventuali differenze empiriche sono da attribuire all’effetto del caso. In alternativa, si potevano usare i seguenti comandi:

m_lin0 <- lm(breath ~ 1, data = anaesthetic)
anova(m_lin0, m_lin)

Analysis of Variance Table

Model 1: breath ~ 1
Model 2: breath ~ tgrp
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     79 1079.1                           
2     76 1008.5  3    70.628 1.7742 0.1592

che producono gli stessi identici risultati, con qualche dettaglio aggiuntivo.

# Punto (d)
par(mfrow = c(2, 2))
plot(m_lin, which = 1:4)

I grafici diagnostici non destano particolari preoccupazioni. Si osservano lievi differenze di variabilità tra i gruppi, ma, considerato il numero limitato di osservazioni, questa moderata eteroschedasticità è verosimilmente attribuibile all’effetto del caso.

Il grafico diagnostico Q–Q plot è quello che presenta le maggiori criticità: si evidenzia infatti uno scostamento piuttosto marcato dalla distribuzione normale. Ciò non risulta particolarmente sorprendente, innanzitutto perché la variabile risposta assume esclusivamente valori positivi e, di conseguenza, la distribuzione dei residui non può essere di tipo gaussiano. Tuttavia, come discusso a lezione, deviazioni dalla normalità non sempre hanno un impatto rilevante sulle conclusioni inferenziali; questo è proprio il caso in esame (si vedano i punti successivi).

In generale, le previsioni di un modello lineare possono assumere valori negativi, circostanza problematica quando la variabile risposta è limitata a valori positivi. Si osservi tuttavia che, in questo esempio specifico, le previsioni coincidono con le medie aritmetiche dei singoli gruppi; di conseguenza, esse non possono assumere valori negativi, come conseguenza della proprietà di internalità della media aritmetica.

Come già osservato, il modello lineare m_lin risulta complessivamente soddisfacente, pur essendo formalmente mal specificato. Un’alternativa naturale, che tiene conto del fatto che la variabile risposta assume valori positivi, è rappresentata da un modello Gamma. Come si vedrà, tale scelta risulta appropriata ma non modifica le conclusioni inferenziali raggiunte.

# Punto (e)
m_gamma <- glm(breath ~ tgrp, data = anaesthetic, family = "Gamma")

La risposta alla gran parte delle domande del punto (e) si poteva desumere dai risultati del comando summary. In alternativa, si potevano calcolare i coefficienti

coef(m_gamma)

(Intercept)       tgrpB       tgrpC       tgrpD 
 0.18518519  0.12250712  0.14134543  0.04404118 

da cui si ottiene che l’equazione che esprima la risposta media in funzione della variabile tgrp è \hat{\mathbb{E}}(Y_i) = \frac{1}{\hat{\beta}_1 + \hat{\beta}_2 d_{i,B} + \hat{\beta}_3 d_{i, C} + \hat{\beta}_4 d_{i, D}} = \frac{1}{0.185 + 0.123 \, d_{i,B} + 0.141 \, d_{i, C} + 0.044\, d_{i, D}}, \qquad i=1,\dots,80.

L’interpretazione dei coefficienti stimati in un modello gamma con legame canonico ha causato qualche difficoltà. La media della variabile aleatoria Y_A, cioè una generica variabile risposta appartenente al gruppo A, è pari a \hat{\mathbb{E}}(Y_A) = \frac{1}{\hat{\beta}_1} = \frac{1}{0.185} = 5.40. In altre parole, \hat{\beta}_1 rappresenta il reciproco della media del gruppo A. Inoltre, prendendo come esempio il gruppo B, la media di Y_B è pari a \hat{\mathbb{E}}(Y_B) = \frac{1}{\hat{\beta}_1 + \hat{\beta}_2} = \frac{1}{0.185 + 0.123} = 3.25. Pertanto, i valori \hat{\beta_2}, \hat{\beta}_3 e \hat{\beta}_4 corrispondono alla differenza tra il reciproco delle media del gruppo di riferimento ed il reciproco della media del gruppo A. Per esempio, si ha che \hat{\beta}_2 = \frac{1}{\mathbb{E}(Y_B)} - \frac{1}{\mathbb{E}(Y_A)}. Di consequenza, valori positivi di \hat{\beta}_2 implicano che \hat{\mathbb{E}}(Y_A) > \hat{\mathbb{E}}(Y_B); la disuguaglianza è invertita in caso di valori negativi di \hat{\beta}_2.

# Punto (e.iii)
predict(m_gamma, newdata = data.frame(tgrp = c("A", "B", "C", "D")), type = "response")

     1      2      3      4 
5.4000 3.2500 3.0625 4.3625 

Le previsioni del modello gamma sono pari, ancora una volta, alle medie campionarie di ciascun gruppo. Non è un caso, come visto negli esercizi teorici. Di conseguenza, le previsioni di questo modello sono identiche a quelle del modello lineare.

# Punto (e.iv)
m_gamma0 <- glm(breath ~ 1, data = anaesthetic, family = "Gamma")
anova(m_gamma0, m_gamma, test = "LRT")

Analysis of Deviance Table

Model 1: breath ~ 1
Model 2: breath ~ tgrp
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1        79     63.146                     
2        76     58.884  3   4.2621   0.1511

Il sistema d’ipotesi del test log-rapporto di verosimiglianza è identico a quello del modello lineare, ossia di verificare l’ipotesi nulla
H_0: \beta_2 = \beta_3 = \beta_4 = 0, contro l’alternativa H_1, secondo cui almeno uno dei coefficienti è diverso da zero. Tale sistema di ipotesi equivale a verificare se le medie dei gruppi siano uguali tra loro oppure no.

Il default del comando anova non è il test LRT (log-rapporto di verosimiglianza), bensì il test F, che non abbiamo mai visto a lezione. I risultati sono corretti da un punto di vista statistico ma non sono quanto richiesto.

Il p-value associato al test W con q = 3 gradi di libertà (numero di vincoli, Df) è pari a 0.1511, che è molto simile a quello ottenuto nel caso del modello m_lin. Le conclusioni qualitative sono quindi le medesime, cioè che non v’è una differenza statisticamente significativa tra le medie dei gruppi.

Il valore del test W non è riportato tra le informazioni fornite dal comando anova. Il valore Deviance rappresenta la differenza tra le devianze dei due modelli, che tuttavia va riscalato per \hat{\phi} quando questo è diverso da 1. Di conseguenza:

deviance_diff <- deviance(m_gamma0) - deviance(m_gamma)
phi_hat <- summary(m_gamma)$dispersion
W <- deviance_diff / phi_hat
round(c(deviance_diff, phi_hat, W), 4)

[1] 4.2621 0.8042 5.2997

Pertanto si ottiene che la statistica test log-rapporto di verosimiglianza è W = 4.2621 / 0.8042 = 5.2997.

# Punto (f)
par(mfrow = c(2, 2))
plot(m_gamma, which = 1:4)

I grafici diagnostici non presentano alcun tipo di problematicità e sono migliorati rispetto a quelli del modello lineare.

# Punto (g)
fit_A <- predict(m_gamma, newdata = data.frame(tgrp = "A"), type = "response", se = TRUE)
z_test <- (fit_A$fit - 10) / fit_A$se.fit # statistica test
p_value <- 2 * (1 - pnorm(abs(z_test))) # p-value
round(c(z_test, p_value), 4)

      1       1 
-4.2481  0.0000 

Una statistica test basata sulla deviazione standard è la seguente z_\text{test} = \frac{\hat{\mu}_A - 10}{\hat{se}(\hat{\mu}_A)} = - 4.2481. Si può assumere, come per gli intervalli di confidenza, che la distribuzione di tale statistica test sia normale standard. Pertanto, il p-value associato è praticamente pari a 0. Di conseguenza, si rifiuta l’ipotesi nulla H_0: \mu_A = 10 in favore dell’ipotesti alternativa H_1 : \mu_A \neq 10.

# Punto (h)

# NON BISOGNAVA FARE NULLA

No, non ha senso verificare se nel modello m_gamma è presente sovradispersione. Il modello Gamma già prevede un parametro di dispersione, in questo caso pari a \hat{\phi} = 0.8042.