bibliografia
i dataset del corso
studentiannoscorso.txt
vaevt.txt
swelling.txt
tubi.txt
Facoltà di Medicina e Chirurgia a.a. 2009-2010
Corso di Laurea Magistrale in Biotecnologie Mediche
Corso Integrato di Innovazione Biotecnologica
Analisi Multivariata dei Dati Sperimentali
Ancova = Regressione + Anova
scegliere e diagnosticare la bontà di un modello per un outcome continuo vs. covariate continue e covariate factor, con possibile interazione
- perché nella Ancova è necessario saper scegliere il modello giusto
www <- "http://www.dmi.units.it/~borelli/medici/studentiannoscorso.txt"
dataset <- read.table( www , header = TRUE )
attach(dataset)
m1 <- lm (peso ~ statura * genere)
m2 <- lm (peso ~ statura + genere)
m3 <- lm (peso ~ statura : genere)
m4 <- lm (peso ~ genere)
m5 <- lm (peso ~ statura)
m6 <- lm (peso ~ 1)
- attenzione: qui (purtroppo) conta l'ordine della regressione
lm (peso ~ statura * genere)
lm (peso ~ genere * statura )
summary.aov(lm (peso ~ statura * genere))
summary.aov(lm (peso ~ genere * statura ))
- partiamo dal modello massimale per giungere al modello minimale adeguato
summary(m1)
- nota: cosa significano i coefficienti? Verifichiamolo con:
lm
(peso[genere == "f"] ~ statura[genere=="f"] )
lm (peso[genere == "m"] ~ statura[genere=="m"] )
lm (peso[genere == "m"] ~ statura[genere=="m"] )
m2 <- lm (peso ~ statura + genere)
summary(m2)
anova(m1, m2)
m3 <- lm (peso ~ statura : genere)
summary(m3)
anova(m3, m2)
m4 <- lm (peso ~ statura + genere + 0 )
summary(m4)
anova(m4, m2)
plot(m2)
- esercizio: provate a scoprire se anche l' anno , lo sport e il fumo sono predittori che interagiscono con la statura e/o il numero di scarpe.
EX Terzo "compito per casa" . La dottoressa Elisa D'Este ha studiato le curve di crescita delle cellule dell'osteosarcoma umano U2OS vs. un clone cellulare denominato H2 in relazione ad una possibile sovraespressione proteica influenzata da un certo ormone steroideo (il Ponasterone A). Quale è il modello minimale adeguato per descrivere il dataset? Interpretate l'output di R ed effettuate le diagnostica del modello.
|
Commento finale. Il "take home message" di queste tre lezioni sul modello lineare si può riassumere in questi passi, incidentalmente posti sia in ordine cronologico che di importanza:
- ipotizzare un modello
- semplificare il modello
- diagnosticare
il modello
importiamo il dataset di Elisa D'Este
www <- "http://www.dmi.units.it/~borelli/biotec/elisadeste.txt"
elisa <- read.table( www , header = TRUE )
attach(elisa)
facciamoci un'idea di come varia il numero di cellule al variare dei giorni, condizionandoci sia alla linea cellulare che alla somministrazione dell'ormone
coplot (log(numerocellule) ~ giorni | linea : ponasterone )
sembrerebbe che un modello lineare possa descrivere il fenomeno. Proviamo a mettere pendenze e quote indipendenti tra loro (modello massimale):
mo1 <- lm (log(numerocellule) ~ linea * ponasterone * giorni)
summary(mo1)
Le diverse pendenze sono significative, ma non le quote. Riduciamo i parametri del modello:
mo2 <- lm (log(numerocellule) ~ linea * ponasterone * giorni - 1 )
summary(mo2)
Il ruolo del ponasterone appare ancora discutibile, soprattutto per quanto riguarda l'osteosarcoma. Ce ne possiamo rendere conto anche graficamente:
coplot (log(numerocellule) ~ giorni | linea )
Creiamo dunque un nuovo fattore a tre livelli
which(linea == "u2os")
which(linea == "h2" & ponasterone =="no")
nuovofattore <- c( rep( "h2_ponSI", each = 7) , rep("h2_ponNO", each = 7) , rep("u2os", each = 14) )
nuovofattore <- as.factor(nuovofattore)
nuovofattore
Proviamo ora a raffinare il modello:
mo3 <- lm (log(numerocellule) ~ nuovofattore * giorni - 1 )
summary(mo3)
controlliamo di "non aver buttato via il bambino con l'acqua":
anova (mo2, mo3)
effettuiamo ora la diagnostica del modello
plot(mo3)
Notiamo che c'è eteroschedastiucità dei residui, e che i punti numeri 21 e 19 potrebbero essere potenzialmente influenti. Proviamo innanzitutto ad escluderli dall'analisi
elisa[19,]
elisa[21,]
mo3ridotto <- update( mo3, subset = ( (numerocellule != 57147863) & (numerocellule != 946237161) ) )
summary(mo3ridotto)
Le stime e le incertezze sui coefficienti sono cambiate di molto poco. Si notino i gradi di libertà del modello, per assicurarci di aver eliminato solo due righe del dataset.
Verifichiamo ora l'eteroschedasticità dei residui:
plot(mo3ridotto)
La situazione è migliorata notevolmente.