bibliografia 





i dataset del corso


studentiannoscorso.txt

vaevt.txt

swelling.txt

tubi.txt





::biodata:: ::didattica:: ::in regalo:: ::utilità::
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



La regressione logistica: introduzione ai modelli lineari generalizzati
 scegliere e diagnosticare la bontà di un modello per un fattore (es. outcome  0/1)  vs. covariate continue e/o factor


 # Vogliamo scoprire se le variazioni a 48 ore del VaeVt e della PaFiO2 sono predittrici dell'outcome.

#  rm (list = ls() )
www <- "http://www.dmi.units.it/~borelli/biotec/vaevt.txt"
vaevt <- read.table( www , header = TRUE )
attach(vaevt)

deltavaevt <- quarantotto - basale
deltapafi <- pafi48 - pafi0

vaevt[1:10,]

# facciamo un po' di analisi esplorativa: sarà più "importante" deltapafi o deltavaevt? Vediamolo con i conditional density plot

layout( matrix(1:2, ncol = 2))
cdplot( factor(outcome) ~ deltavaevt)
cdplot( factor(outcome) ~ deltapafi)

# ancora un po' di analisi esplorativa: la probabilità di avere outcome infausto?

plot( deltapafi ~ deltavaevt, col = "white" )
points( deltapafi[outcome =="1"] ~ deltavaevt[outcome=="1"], col = "darkgreen", pch = 22 )
points( deltapafi[outcome =="0"] ~ deltavaevt[outcome=="0"], col = "red", pch = 19 )


# introduciamo un modello lineare generalizzato


modello1 <- glm( outcome ~ deltavaevt, family = binomial() )
summary(modello1)


# l'intercetta non è significativa, leviamola:

modello2 <- glm( outcome ~ deltavaevt - 1 , family = binomial() )
summary(modello2)


# vediamo se anche deltapafi è un valido predittore:

modello3 <- glm( outcome ~ deltavaevt + deltapafi - 1 , family = binomial() )
summary(modello3)


# a quanto pare, non serve aggiungere deltapafi nel modello. Lo controlliamo formalmente:


anova( modello2, modello3, test = "Chisq" )


# cosa significano i coefficienti di modello2? Vediamolo graficamente
 
modello2

plot( deltavaevt, outcome)

b <- modello2$coef[[1]]


probabilita <-  exp ( 0 + b * deltavaevt) /  ( 1 + exp( 0 + b * deltavaevt) )
probabilita

plot( deltavaevt, outcome )
points( deltavaevt, probabilita, col = "red")

min(deltavaevt)
max(deltavaevt)

x <- -285:405/100
lineacurva <-  exp ( 0 + b * x ) /  ( 1 + exp( 0 + b * x ) )
plot( deltavaevt , outcome)
points( deltavaevt , probabilita, col = "red")
points( x, lineacurva, col = "red", pch =".")
lines ( c( min(deltavaevt) , max(deltavaevt) ) , c(0.5, 0.5) , col = "blue" )

# un aumento di 1 unità di  deltavaevt  aumenta il log-odds di modello2$coef[[1]] . Se invece della stima puntuale desideriamo una stima intervallare?
 
confint(modello2, parm = "deltavaevt")

# passiamo dai logit alle probabilità

exp(coef(modello2)["deltavaevt"])
exp(confint(modello2, parm = "deltavaevt"))

# l'intervallo di fiducia al 95% è molto ampio, e noninformativo.
# Si noti però che il valore 1 non appartiene all'intervallo (descrizione equivalente di *)

# Visualizziamo con un bubble plot il modello2 raffrontandolo con il plot esplorativo di poco fa
 
prob2 <- predict (modello2, type = "response")
prob2

plot( deltapafi ~ deltavaevt, xlim = c(-0.3, 0.5), ylim = c(-200,300), col = "white" )
points( deltapafi[outcome =="1"] ~ deltavaevt[outcome=="1"], col = "darkgreen", pch = 22 )
points( deltapafi[outcome =="0"] ~ deltavaevt[outcome=="0"], col = "red", pch = 19 )
symbols( deltavaevt , deltapafi , circles = prob2, add = TRUE)

# si osservi infine il fenomeno della sovradispersione: la devianza residua è superiore ai gradi di libertà residui. Possiamo compensare questo fatto modificando il parametro di dispersione (cioè, la "forma" della variabile aleatoria che modella la componente casuale Y del glm.

 modello2bis <- glm( outcome ~ deltavaevt - 1 , family = quasibinomial() )
summary(modello2bis)

# attenzione. Avremmo sbagliato di grosso a tentare un modello lineare ..

modellosbagliato <- lm( outcome ~ deltavaevt)
summary(modellosbagliato)

plot(modellosbagliato)





Quarto "compito per casa" da inviarmi per posta.
Indicate nell'oggetto: Cognome  4  (es. Borelli 4)

(insensato scientificamente, ma didattico): provate a scoprire se il body mass index e il genere sono predittori del fumo nel dataset studentiannoscorso.