1 Analisi della varianza a due vie (modello lineare e GLM)
inserire formule
trattare un caso non 2x2 ma bilanciato
eventuale (futura) discussione sui casi sbilanciati.
1.1 Inizializzazione:
vengono caricati i packages necessari per realizzare questo documento
Data Frame Summary
X
Dimensions: 93 x 3Duplicates: 49
| No | Variable | Stats / Values | Freqs (% of Valid) | Graph | ||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | sfebbr [integer] | Mean (sd) : 5.5 (5.1) min < med < max: 0 < 4 < 30 IQR (CV) : 6 (0.9) | 20 distinct values |  |
||||||||
| 2 | trattamento [factor] | 1. B 2. C |
|
 |
||||||||
| 3 | livelloRischio [factor] | 1. alto 2. basso |
|
 |
Generated by summarytools 0.9.6 (R version 4.0.2)
2020-09-16
1.2 Il problema generale
Il data.frame trial1.new contiene parte dei risultati di un clinical trial.
In questa lezione non entrerò nel dettaglio di come sono stati raccolti i dati e come fosse stato pianificato l’esperimento, ma descrivo alcune caratteristiche dei dati che serviranno a delineare il problema statistico.
Le 93 unità sono i pazienti sottoposti ad una terapia antifungina. Sono suddivisi in pazienti ad alto e basso rischio (variabile livelloRischio).
Per ognuno dei livelli di tale variabile è stato assegnato a caso il trattamento sperimentale (variabile trattamento), per cui la suddivisione delle osservazioni nelle quattro caselle date dalle combinazioni dei due fattori non è esattamente bilanciata.
La variabile sfebbr rappresenta il numero di giorni trascorsi per lo sfebbramento dopo la somministrazione del trattamento
Si vuole studiare la dipendenza della variabile sfebbr (quantitativa) dalle variabili trattamento e livelloRischio (variabili qualitative con due livelli ciascuna).
Ovviamente impiegheremo la metodologia generale dei modelli lineari, ma in effetti lo scopo dell’analisi è vedere quale dei due trattamenti trattamento sia più efficace ai fini dello sfebbramento, ossia se il valore atteso della variabile sfebbramento differisce significativamente nei gruppi corrsiondenti ai due trattamenti.
L’analisi potrebbe essere complicata dal fatto che comunque si sa che il valore atteso di sfebbr ragionevolmente è influenzato dalla variabile livelloRischio, che potrebbe confondere il risultato.
Prima però effettuiamo qualche analisi preliminare
1.3 Analisi preliminare dei dati
Le osservazioni so raggruppate in quattro gruppi, poichè le due variabili di classificazione hanno due modalità sciascuna.
alto basso
B 29 15
C 32 17Intanto calcoliamo medie e dispersioni della variabile d’interesse livelloRischio in funzione delle altre due. Otterremmo una tavola di medie ed una di deviazione standard
livelloRischio
trattamento alto basso
B 6.275862 2.733333
C 6.718750 4.352941 livelloRischio
trattamento alto basso
B 5.344534 2.120198
C 6.118372 2.9142151.4 Omoschedasticità
Vediamo adesso se,almeno dal punto vista empirico, vi è una relazione fra dispersione e medie. E’ un tipo di rappresentazione sempre utile da fare tutte le volte che i dati sono suddivise in gruppi.
dal grafico appare evidente che c’è una relazione fra dispersione e medie.[^nota1] [^nota1]: In effetti in questo modulo del corso non affronto il problema di un vero e proprio test per saggiare l’ipotesi di uguaglianza fra le varianze all’interno dei gruppi (il cosiddetto test di Bartlett)
Che vi siano comunque forti dubbi sulla ragionevolezza dell’assunzione di uguaglianza delle varianze appare chiaro applicando un semplice modello additivo per la variabile sfebbr, e rappresentando graficamente la relazione fra i residui empirici e i valori calcolati.
L’aumento della dispersione appare chiaro dal tipico aspetto ad imbuto…
1.5 Una trasformata della variabile.
Sebbene non sia l’apptroccio metodologicamente più corretto, per ora affrontiamo in prima approssimazione il problema trasformando la variabile di risposta con \(log(y+\frac{1}{2})\). In questo modo eviteremo logaritmi di zeri e consideriamo anche il fatto che i gironi di sfebbramento in realtà sono stati rlievati troncando il numero di giorni trascorsi fra la somministrazione del trattamento e il primo momento in cui viene rilevato losfebbramento, per cui un valore di zero giorni significa in realtà che è passato meno di un giorno, per cui un valore riportato in sfebbr di \(y\) giorni in realtà corrisponde a una numero di giorni maggiore o uguale a zero e inferiore a uno, per cui il valore centrale della classe è \(y+\frac{1}{2}\), per cui applichiamo a questa la trasformata logaritmica.
- Riaffronteremo lo stesso problema con una metodologia più adeguata più avanti, quando introdurremo i cosiddetti modelli lineari generalizzati GLM *
Vediamo che intanto con questo cambio di variabile otteniamo uno smorzamento della relazione fra deviazione standard e medie
livelloRischio
trattamento alto basso
B 1.651687 0.9327282
C 1.683377 1.4010235 livelloRischio
trattamento alto basso
B 0.7262758 0.7959457
C 0.7777326 0.6287361
Gli ultimi due grafici differiscono solo per il cambio di scala in ordinata: il primo grafico potrebbe far credere ancora ad una relazione fra dispersioni e medie, ma in realtà questo è solo dovuto alla scala verticale (che comunque per default R calcola in modo da coprire quasi esattamente il range della variabile)
1.6 Analisi sulla variabile trasformata
Consideriamo due modelli, uno additivo
e l’altro moltiplicativo, ossia con interazioni
Il test per il confronto fra i due modelli (test sull’interazione), mediante anova(m1,m2), sostanzialmente saggia l’ipotesi che le interazioni siano nulle (la interazione, perchè in questo caso abbiamo un solo grado di libertà per le interazioni dato che si tratta di una tavola \(2 \times 2\))
Analysis of Variance Table
Model 1: xl ~ trattamento + livelloRischio
Model 2: xl ~ trattamento * livelloRischio
Res.Df RSS Df Sum of Sq F Pr(>F)
1 90 49.711
2 89 48.715 1 0.99687 1.8212 0.1806IL test risulta non significativo, per cui possiamo utilizzare il modello additivo m1
Call:
lm(formula = xl ~ trattamento + livelloRischio, data = X)
Residuals:
Min 1Q Median 3Q Max
-1.77809 -0.35043 -0.06888 0.49656 1.66300
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.5730 0.1249 12.589 < 2e-16 ***
trattamentoC 0.1818 0.1544 1.178 0.24206
livelloRischiobasso -0.4880 0.1622 -3.008 0.00341 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7432 on 90 degrees of freedom
Multiple R-squared: 0.1035, Adjusted R-squared: 0.08359
F-statistic: 5.196 on 2 and 90 DF, p-value: 0.007321.7 Perchè il modello con le interazioni è detto moltiplicativo?
codifica 0/1 delle due variabili esplicative l’interazione è un effetto moltiplicativo e coincide con l’effetto x12=x1*x2
xl x1 x2
[1,] 0.9162907 0 0
[2,] 1.8718022 1 0
[3,] 2.0149030 0 1
[4,] 1.7047481 1 0
[5,] 1.7047481 0 0
[6,] 2.5257286 1 1 xl x1 x2 x12
[1,] 0.9162907 0 0 0
[2,] 1.8718022 1 0 0
[3,] 2.0149030 0 1 0
[4,] 1.7047481 1 0 0
[5,] 1.7047481 0 0 0
[6,] 2.5257286 1 1 1
Call:
lm(formula = xl ~ x1 + x2 + x12)
Residuals:
Min 1Q Median 3Q Max
-1.62588 -0.48473 -0.01644 0.56791 1.73435
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.9327 0.1910 4.883 4.57e-06 ***
x1 0.4683 0.2621 1.787 0.07737 .
x2 0.7190 0.2353 3.056 0.00296 **
x12 -0.4366 0.3235 -1.350 0.18059
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7398 on 89 degrees of freedom
Multiple R-squared: 0.1215, Adjusted R-squared: 0.09188
F-statistic: 4.103 on 3 and 89 DF, p-value: 0.008921
Call:
lm(formula = xl ~ x1 + x2)
Residuals:
Min 1Q Median 3Q Max
-1.77809 -0.35043 -0.06888 0.49656 1.66300
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.0849 0.1549 7.005 4.33e-10 ***
x1 0.1818 0.1544 1.178 0.24206
x2 0.4880 0.1622 3.008 0.00341 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7432 on 90 degrees of freedom
Multiple R-squared: 0.1035, Adjusted R-squared: 0.08359
F-statistic: 5.196 on 2 and 90 DF, p-value: 0.00732
Call:
lm(formula = xl ~ trattamento + livelloRischio, data = X)
Residuals:
Min 1Q Median 3Q Max
-1.77809 -0.35043 -0.06888 0.49656 1.66300
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.5730 0.1249 12.589 < 2e-16 ***
trattamentoC 0.1818 0.1544 1.178 0.24206
livelloRischiobasso -0.4880 0.1622 -3.008 0.00341 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7432 on 90 degrees of freedom
Multiple R-squared: 0.1035, Adjusted R-squared: 0.08359
F-statistic: 5.196 on 2 and 90 DF, p-value: 0.00732
Call:
lm(formula = xl ~ x1 * x2)
Residuals:
Min 1Q Median 3Q Max
-1.62588 -0.48473 -0.01644 0.56791 1.73435
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.9327 0.1910 4.883 4.57e-06 ***
x1 0.4683 0.2621 1.787 0.07737 .
x2 0.7190 0.2353 3.056 0.00296 **
x1:x2 -0.4366 0.3235 -1.350 0.18059
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7398 on 89 degrees of freedom
Multiple R-squared: 0.1215, Adjusted R-squared: 0.09188
F-statistic: 4.103 on 3 and 89 DF, p-value: 0.0089211.8 Un’approccio più moderno e generalizzabile:i GLM
In effetti possiamo pensare che i tempi di sfebbramento seguano una distribuzione asimmetrica, approssimabile da una distribuzione gamma, con un valore atteso che potrebbe essere messo in relazione con un predittore lineare nelle stesse variabili che abbiamo considerato finora, ossia trattamento e livelloRischio.
Adattiamo un modello additivo glm1
Call:
glm(formula = y ~ trattamento + livelloRischio, family = Gamma(link = "log"),
data = X)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.5029 -0.5464 -0.3305 0.2335 1.7923
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.8538 0.1276 14.533 <2e-16 ***
trattamentoC 0.1801 0.1576 1.143 0.2560
livelloRischiobasso -0.5539 0.1656 -3.345 0.0012 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 0.5756454)
Null deviance: 54.283 on 92 degrees of freedom
Residual deviance: 47.855 on 90 degrees of freedom
AIC: 497.19
Number of Fisher Scoring iterations: 51.9 Confronto fra un GLM additivo e uno moltiplicativo
Adesso, adattiamo un modello moltiplicativo glm2,seguendo la stessa impostazione per il modello lineare impiegato per la variabile trasformata.
Call:
glm(formula = y ~ trattamento * livelloRischio, family = Gamma(link = "log"),
data = X)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.4292 -0.5974 -0.2467 0.3072 1.8889
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.91337 0.14123 13.548 < 2e-16 ***
trattamentoC 0.06332 0.19499 0.325 0.74616
livelloRischiobasso -0.73985 0.24188 -3.059 0.00294 **
trattamentoC:livelloRischiobasso 0.34276 0.33258 1.031 0.30552
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 0.5784149)
Null deviance: 54.283 on 92 degrees of freedom
Residual deviance: 47.244 on 89 degrees of freedom
AIC: 497.9
Number of Fisher Scoring iterations: 5Analysis of Deviance Table
Model 1: y ~ trattamento + livelloRischio
Model 2: y ~ trattamento * livelloRischio
Resid. Df Resid. Dev Df Deviance
1 90 47.855
2 89 47.244 1 0.61192L’interpretazione dei risultati è la stessa in fondo di quella che faremmo col modello lineare: possiamo utilizzare un modello addiditivo (glm1)
Call:
glm(formula = y ~ trattamento + livelloRischio, family = Gamma(link = "log"),
data = X)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.5029 -0.5464 -0.3305 0.2335 1.7923
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.8538 0.1276 14.533 <2e-16 ***
trattamentoC 0.1801 0.1576 1.143 0.2560
livelloRischiobasso -0.5539 0.1656 -3.345 0.0012 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 0.5756454)
Null deviance: 54.283 on 92 degrees of freedom
Residual deviance: 47.855 on 90 degrees of freedom
AIC: 497.19
Number of Fisher Scoring iterations: 5Da quest’ultimo summary ritroviamo ancora la stessa conclusione, ossia che il diverso trattamento non influenza il valore atteso dei tempi di sfebbramento. Mentre esiste un effetto significativo per il fattore livelloRischio