Facoltà di Medicina e Chirurgia a.a. 2008-2009
Corso di Laurea Magistrale in Biotecnologie Mediche
Corso Integrato di Innovazione Biotecnologica
Analisi Multivariata dei Dati Sperimentali (Borelli)
L'analisi della varianza ("Anova")
scegliere un modello per un outcome continuo normale vs. covariate di tipo factor
|
un esempio di Crawley www <- "http://www.dmi.units.it/~borelli/s3v/08borelli_anova/growth.txt"
weights <- read.table( www , header = TRUE ) attach(weights) tratto dai Seminari per il Dipartimento di Scienze della Vita |
i dataset
library(HSAUR)
weightgain
source
type weightgain
1 Beef Low 90
2 Beef Low 76
...
19 Beef High 117
20 Beef High 111
21 Cereal Low 107
22 Cereal Low 95
...
38 Cereal High 77
39 Cereal High 86
40 Cereal High 92
1 Beef Low 90
2 Beef Low 76
...
19 Beef High 117
20 Beef High 111
21 Cereal Low 107
22 Cereal Low 95
...
38 Cereal High 77
39 Cereal High 86
40 Cereal High 92
foster
litgen
motgen weight
1 A A 61.5
2 A A 68.2
...
8 A B 60.2
9 A I 52.5
...
17 A J 39.6
18 B A 60.3
...
60 J J 42.0
61 J J 54.0
1 A A 61.5
2 A A 68.2
...
8 A B 60.2
9 A I 52.5
...
17 A J 39.6
18 B A 60.3
...
60 J J 42.0
61 J J 54.0
skulls
epoch mb bh bl nh
1 c4000BC 131 138 89 49
2 c4000BC 125 131 92 48
3 c4000BC 131 132 99 50
4 c4000BC 119 132 96 44
5 c4000BC 136 143 100 54
6 c4000BC 138 137 89 56
...
148 cAD150 140 135 103 48
149 cAD150 147 129 87 48
150 cAD150 136 133 97 51
1 c4000BC 131 138 89 49
2 c4000BC 125 131 92 48
3 c4000BC 131 132 99 50
4 c4000BC 119 132 96 44
5 c4000BC 136 143 100 54
6 c4000BC 138 137 89 56
...
148 cAD150 140 135 103 48
149 cAD150 147 129 87 48
150 cAD150 136 133 97 51
perché si dice "analisi della varianza"?
cosa si deve ricordare (le parole di Everitt)
- anova / manova; i disegni fattoriali
- bilanciato / non bilanciato
- il modello matematico della Anova
- la notazione di Wilkinson e Rogers
- le ipotesi del test F
- indipendendenza delle osservazioni
- time plot
- normalità dei dati
- qqplot / shapiro.test
- omoschedasticità
- bartlett.test
- il modello matematico della Manova
- attenzione alla Manova con >2 livelli
- Hotelling - Lawley
- Wilks
- Roy
- Pillai
esempio numero 1
data("weightgain", package = "HSAUR")
names(weightgain)
tapply(weightgain$weightgain, list(weightgain$source, weightgain$type), mean)
tapply(weightgain$weightgain, list(weightgain$source, weightgain$type), sd)
plot.design(weightgain)
verifichiamo le ipotesi dell'Anova
plot(1:length(weightgain$weightgain), weightgain$weightgain)
qqnorm(weightgain$weightgain)
shapiro.test(weightgain$weightgain)
bartlett.test(weightgain$weightgain ~ weightgain$source)
bartlett.test(weightgain$weightgain ~ weightgain$type)
vediamo la tavola Anova
modello <- aov (weightgain$weightgain ~ weightgain$source * weightgain$type)
summary(modello)
# cosa significa? ce lo dice Everitt
vediamo le stime su mu, gamma, beta e gammabeta
coef(modello)
coef(modello)[[1]]
coef(modello)[[1]] + coef(modello)[[2]]
coef(modello)[[1]] + coef(modello)[[3]]
coef(modello)[[1]] + coef(modello)[[2]] + coef(modello)[[3]] + coef(modello)[[4]]
tapply(weightgain$weightgain, list(weightgain$source, weightgain$type), mean)
capiamo cosa fa il termine di interazione
interaction.plot (weightgain$type, weightgain$source, weightgain$weightgain)
la diagnostica del modello
plot (modello)
esempio numero 2
data("foster", package = "HSAUR")
names(foster)
design is unbalanced = non orthogonality
uni <- rep(1,length(foster$weight))
contingenza <- xtabs(uni ~ litgen + motgen, data = foster)
contingenza
verifichiamo le ipotesi dell'Anova (esercizio)
vediamo le tavole Anova
summary(aov (weight ~ litgen * motgen, data = foster))
summary(aov (weight ~ motgen * litgen, data = foster))
procedura di multiple comparison
analisi <- aov (weight ~ litgen * motgen, data = foster)
onesto <- TukeyHSD (analisi, "motgen" )
onesto
plot(onesto)
la diagnostica del modello (esercizio)
esempio numero 3
data("skulls", package = "HSAUR")
names(skulls)
medie <- aggregate(skulls [ , c("mb", "bh", "bl", "nh")] , list (epoch = skulls$epoch) , mean)
medie
verifichiamo le ipotesi dell'Anova (esercizio)
vediamo le tavole Manova - analisi multivariata
analisi <- manova (cbind(mb, bh, bl, nh) ~ epoch, data = skulls)
summary(analisi)
summary(analisi, test = "Wilks")
summary(analisi, test = "Hotelling-Lawley")
summary(analisi, test = "Roy")
vediamo le tavole univariate, F test
summary.aov(manova (cbind(mb, bh, bl, nh) ~ epoch, data = skulls))
effettuiamo dei test multivariati accoppiati
summary(manova (cbind(mb, bh, bl, nh) ~ epoch, data = skulls), subset = epoch %in% c("c4000BC", "c3300BC"))
summary(manova (cbind(mb, bh, bl, nh) ~ epoch, data = skulls), subset = epoch %in% c("c4000BC", "c1850BC"))
summary(manova (cbind(mb, bh, bl, nh) ~ epoch, data = skulls), subset = epoch %in% c("c4000BC", "c200BC"))
summary(manova (cbind(mb, bh, bl, nh) ~ epoch, data = skulls), subset = epoch %in% c("c4000BC", "cAD150"))
- la raccomandazione di Stevens (2001): alfa = 0.15 / m
la diagnostica del modello (?)
attenzione a non cadere in
trappola 4
dati <- c(13, 22, 25, 17, 20, 18, 20, 14, 20, 18, 23, 15, 22, 18, 24, 23, 16, 19, 19, 23, 21, 18, 21, 19, 19, 18, 19, 14, 24, 13, 22, 23, 17, 19, 17, 16, 23, 18, 16, 18, 19, 18, 17, 19, 19, 18, 25, 19, 24, 23, 16, 19, 16, 16, 17, 23, 21, 15, 16, 21, 18, 17, 17, 18, 17, 17, 19, 17, 17, 18, 17, 18, 17, 16, 18, 21, 22, 20, 23, 22, 22, 22, 21, 21, 22, 23, 22, 22, 20, 22)
fattore <- c( rep(1, 30), rep(2,30), rep(3,30) )
# inganno: sembra che non ci sia differenza
summary.lm ( lm (dati ~ factor(fattore) ))
summary.aov ( lm(dati ~ factor(fattore) ))
modello <- aov (dati ~ factor(fattore))
TukeyHSD(modello)
plot(TukeyHSD(modello))
mean(dati[1:30])
mean(dati[31:60])
mean(dati[61:90])
# invece!!! la verità ..
par (mfrow = c(1,3))
hist(dati[1:30])
hist(dati[31:60])
hist(dati[61:90])
bartlett.test(dati ~ factor(fattore))
un esercizio finale
Ispirandoci a un precedente seminario, importate il dataset growth.txt di M. Crawley,
www <- "http://www.dmi.units.it/~borelli/s3v/08borelli_anova/growth.txt"
weights <- read.table( www , header = TRUE )
attach(weights)
e seguite le istruzioni .