ANOVA 2

Monica Marabelli

15 Gennaio 2016

ANOVA a una via

Nella scorsa esercitazione abbiamo visto che nell’Analisi dellaVarianza (ANOVA) si considerano le medie di una variabiledipendente (quantitativa) negli strati indotti dalle modalitá diuna seconda variabile (qualitativa), detta variabile esplicativa ofattore.

Si vuole verificare l’ipotesi nulla che tutte le medie siano ugualicontro l’ipotesi alternativa che almeno una coppia di mediepresenti una differenza statisticamente significativa.

Esercizio 1E’ stato condotto uno studio per valutare la durata della vita tra imembri di sovranitá, aristocrazia e nobiltá in Inghilterra.

I dati sono nel file age.xls

setwd("X:/")age <- read.table (file="age.csv",

header=T,sep=",",dec=".")

str(age)

’data.frame’: 700 obs. of 2 variables:$ class : Factor w/ 3 levels "aris","gent",..: 1 3 2 1 1 2 1 3 3 1 ...$ age_death: int 81 43 70 54 90 56 50 63 38 56 ...

Visualizziamo graficamente i dati

boxplot(age_death~class, data=age)

ANOVA

modello <- aov (lm(age_death~class, data=age))summary(modello)

Df Sum Sq Mean Sq F value Pr(>F)class 2 2452 1226.2 4.073 0.0174 *

Residuals 697 209849 301.1– – –Signif. codes: 0 `***´ 0.001 `**´ 0.01 `*´ 0.05 `.´ 0.1 ` ´ 1

Pvalue < 0.05: c’é almeno una coppia di medie la cui differenzaé statisticamente significativa (per α = 0.05).

Test post hoc

Per individuare quali sono le medie significativamente diversetra loro, applichiamo il test HSD (honest significant difference)di Tukey

TukeyHSD (modello)

Tukey multiple comparisons of means95% family-wise confidence level

Fit: aov(formula = lm(age_death ~ class, data = age))

$classdiff lwr upr p adj

gent-aris -2.376541 -5.944069 1.1909858 0.2617559sovr-aris -4.787597 -8.752292 -0.8229026 0.0130255sovr-gent -2.411056 -6.369596 1.5474843 0.3257161

Verifica degli assunti

Prima di procedere all’analisi ANOVA, occorre tuttavia verificarei seguenti assunti:

I normalitá dei residui entro gruppiI omoschedasticitá dei residui entro gruppi (uguali varianze)I indipendenza dei residui entro gruppi

Normalitá dei residui

shapiro.test (resid(aov(lm(age_death~class, data=age)))[age$class == "gent"])

Shapiro-Wilk normality testdata: resid(aov(lm(age_death ~class, data = age)))[age$class == "gent"]W = 0.961, p-value = 1.597e-06

shapiro.test (resid(aov(lm(age_death~class, data=age)))[age$class == "sovr"])

Shapiro-Wilk normality testdata: resid(aov(lm(age_death ~class, data = age)))[age$class == "sovr"]W = 0.9756, p-value = 0.003201

Normalitá dei residui

shapiro.test (resid(aov(lm(age_death~class, data=age)))[age$class == "aris"])

Shapiro-Wilk normality testdata: resid(aov(lm(age_death ~class, data = age)))[age$class == "aris"]W = 0.9778, p-value = 0.0004311

Pvalue < 0.05 per tutte e tre le classi. Rifiutiamo l’ipotesi nullache i residui siano distribuiti normalmente.

Omoschedasticitá dei residui

bartlett.test (age_death~class, data=age)

Bartlett test of homogeneity of variances

data: age_death by classBartlett’s K-squared = 3.3476, df = 2, p-value = 0.1875

Pvalue > 0.05: non rifiutiamo l’ipotesi nulla che le varianze sianoomogenee.

Indipendenza dei residui

Installiamo in R il pacchetto lmtestI PackagesI Install package(s)

Indipendenza dei residui

I Scegliere uno dei CRAN Mirrors italiani (Milano)I Cercare il pacchetto lmtest

Indipendenza dei residui

library(lmtest) # Rendiamo disponibili le funzioni del pacchetto

Applichiamo il test di Durbin–Watson

dwtest (aov(lm(age_death~class, data=age)))

Durbin-Watson test

data: aov(lm(age_death ~ class, data = age))DW = 2.0574, p-value = 0.7765alternative hypothesis: true autocorrelation is greater than 0

Pvalue > 0.05: non rifiuto l’ipotesi nulla di indipendenza deiresidui.

Analisi non parametrica

Siccome l’ipotesi di normalitá non é stata confermata,effettuiamo un’analisi non parametrica: applichiamo il test diKruskal-Wallis.

kruskal.test (age_death~class, data=age)

Kruskal-Wallis rank sum test

data: age_death by classKruskal-Wallis chi-squared = 8.3987, df = 2, p-value = 0.01501

Pvalue < 0.05. Rifiuto l’ipotesi nulla: l’appartenenza alle treclassi sociali ha un effetto significativamente diverso sulladurata della vita.

ANOVA a blocchi randomizzati

Nell’analisi a blocchi randomizzati, gli individui sono suddivisi ingruppi omogenei, detti blocchi; gli individui di ogni bloccosaranno poi attribuiti in modo casuale (random) ai trattamenti.

Esercizio 2

Alcuni maialini d’India sono sottoposti a 4 diete differenti. 5gruppi di 4 animali ciascuno sono tenuti in identiche condizioniambientali. Entro blocco le condizioni ambientali sonoidentiche, tra blocchi le condizioni ambientali possono risultarediverse. Gli animali di ciascun gruppo sono assegnati a casoad ognuna delle 4 diete.

Il nostro scopo é studiare l’effetto del tipo di dieta sul peso deimaialini d’India, cioé verificare se esiste una differenza diincremento di peso tra le diverse diete somministrate.

Leggiamo i dati

I dati sono contenuti nel file maialini.csv

setwd("Y:/STATISTICA")dati <- read.table (file="maialini.csv",

header=T,sep=",",dec=".")

head(dati)

Blocco Dieta Peso1 1 1 1.52 1 2 2.73 1 3 2.14 1 4 1.35 2 1 1.46 2 2 2.9

Trasformiamo i dati nel formato corretto

str(dati)

’data.frame’: 20 obs. of 3 variables:$ Blocco: int 1 1 1 1 2 2 2 2 3 3 ...$ Dieta : int 1 2 3 4 1 2 3 4 1 2 ...$ Peso : num 1.5 2.7 2.1 1.3 1.4 2.9 2.2 1 1.4 2.1 ...

Per poter effettuare l’analisi, le variabili trattamento e bloccodevono essere fattori.Ricordarsi di trasformare i dati prima di effettuare l’analisi!

dati$Blocco <- factor(dati$Blocco)dati$Dieta <- factor(dati$Dieta)

ANOVA a blocchi randomizzati

Effettuiamo l’analisi della varianza

modello2 <- aov(lm(Peso ~ Dieta + Blocco, data=dati))summary(modello2)

Df Sum Sq Mean Sq F value Pr(>F)Dieta 3 8.154 2.7178 41.866 1.24e-06 ***

Blocco 4 0.393 0.0982 1.513 0.26Residuals 12 0.779 0.0649

– – –Signif. codes: 0 `***´ 0.001 `**´ 0.01 `*´ 0.05 `.´ 0.1 ` ´ 1

H0 che le 4 diete abbiano lo stesso effetto: rifiutata.H0 che l’effetto sia uguale tra i 5 blocchi: non rifiutata.

ANOVA a due vie con replica

Si parla di ANOVA a due vie quando viene studiato l’effetto didue variabili esplicative e si valuta anche l’interazione tra diesse.

Esercizio 3

Vogliamo studiare l’effetto sulla capacitá vitale di: tipo di lavoro,etá e loro interazione.

Leggiamo i dati contenuti nel file capacity.csv

setwd("Y:/STATISTICA")hum_cap <- read.table (file="capacity.csv",

header=T,sep=",",dec=".")

Esercizio 3

str(hum_cap)

’data.frame’: 60 obs. of 3 variables:$ lavoro : Factor w/ 4 levels "A","B","C","D": 1 1 1 1 1 1 1 1 1 ...$ cl_eta : int 1 1 1 1 1 2 2 2 2 2 ...$ capacita: num 4.31 4.89 4.05 4.44 4.59 4.13 4.61 3.91 4.52 ...

Trasformiamo in fattore la variabile classi d’etá

hum_cap$cl_eta <- factor(hum_cap$cl_eta)

ANOVA a due vie con replica

modello3 <- aov(lm(capacita ~ lavoro*cl_eta, data=hum_cap))summary(modello3)

Df Sum Sq Mean Sq F value Pr(>F)lavoro 3 19.779 6.593 31.48 2.13e-11 ***cl_eta 2 12.309 6.154 29.38 4.65e-09 ***

lavoro:cl_eta 6 8.949 1.491 7.12 1.83e-05 ***Residuals 48 10.054 0.209

– – –Signif. codes: 0 `***´ 0.001 `**´ 0.01 `*´ 0.05 `.´ 0.1 ` ´ 1

I risultati indicano che sia il tipo di lavoro sia l’etá sial’interazione dei due fattori hanno un effetto significativo sullacapacitá vitale.