# Stima modelli di regressione logistica in R - Esempi

# 1. “Costruzione” modello
# (“migliore” modello nello specifico contesto applicativo)

# In generale inclusione di:

# variabili scientificamente rilevanti in relazione al contesto (problema/disciplina)
# con controllo di tutti i fattori di confondimento, evitando overfitting
# ma anche specificazioni parsimoniose (risultati più stabili e maggiore generalizzazione)
# buon adattamento ai dati

# 1. Rappresentazione grafica di logit(pi) = beta0 + beta1xi
# grafico 1: beta0 = 1, beta1 = .5

beta0=1
beta1=.5
beta0
beta1

# visualizzazione grafici accostati (mfrow = 'make frame by row')
par(mfrow=c(1,2))
# finestra del grafico suddivisa in 1 riga e 2 colonne

# grafico 1
curve(expr=exp(beta0+beta1*x)/(1+exp(beta0+beta1*x)), xlim=c(-15,15), ylab="p")

# grafico 2: beta0 = 1, beta1 = -.5
beta1= -.5
curve(expr=exp(beta0+beta1*x)/(1+exp(beta0+beta1*x)), xlim=c(-15,15), ylab="p")

# 0 < p < 1
# beta1 > 0: relazione positiva tra x e p (negativa se beta1 < 0)
# pendenza della curva dipende dai valori x (derivata di p rispetto a x = beta1p(1-p))
#  p >.5: immagine speculare rispetto a p < .5

par(mfrow=c(1,1))

# Cosa succede al crescere(decrescere) di beta1 ?

# 2. Alcuni esempi
# 2.1 Challenger 

# Il file challenger.txt contiene le registrazioni dei guasti occorsi alle guarnizioni 
# ad anelli in 23 lanci dello space shuttle e le temperature atmosferiche al momento 
# del lancio (space shuttle Challenger distrutto a 73 secondi dal lancio il 28 gennaio 1986). 
# La variabile failures vale 1 se si è verificato un guasto e vale 0 altrimenti. 
# La variabile temp è la temperatuta (gradi Fahrenheit) --50 °F = 10 °C, 80 °F = 27 °C--

chl=read.table("challenger.txt", TRUE)
head(chl)
str(chl)
table(chl$failures)
table(chl$temp)

# Obiettivo: la probabilità di un guasto dipende dalla temperatura?

plot(chl$temp, chl$failures)
# sono tutti ?

par(mfrow = c(1, 2))
plot(chl$temp, chl$failures)
plot(jitter(chl$temp), chl$failures)

boxplot(chl$temp~chl$failures)

# Modello di regressione logistica per Yi (failures) e xi(temp)
# pi= (P(Yi=1)), logit(pi)=log(pi/1-pi)=beta0+beta1xi
# Stima massima verosimiglianza del modello:

mod.chl = glm(formula = failures ~ temp, family = binomial(link = logit), data=chl)

# Argomenti (principali) di glm:

# 'formula' specifica il modello: a sx di ~ var. risposta, a dx ~ var. esplicative (puo' essere omessa)

# 'family =' tipo di modello che deve essere stimato in relazione a var. risposta 
# (Binomiale generalizzazione di Beornulli - distribuzione per Y), 'link' indica il 
# legame(collegamento) tra E(Y) e var. esplicative (in modello logistico il collegamento 
# avviene tramite trasformazione logit)

#'data =' nome del data frame 

# Contenuto dell'oggetto mod.chl
mod.chl

# Modello stimato: logit(p.hat)=15.0429-0.2322x
# Probabilità (stimata) di guasto decresce all'aumentare della temperatura
# Anche (molte) altre informazioni (componenti) contenute: 
names(mod.chl)

# Coefficienti beta0 e beta1
mod.chl$coefficients


# Valori Y
mod.chl$y
# Probabilità stimate 
mod.chl$fitted.values

#Valori Y osservati e probabilità 
confronto= cbind(chl$failures,mod.chl$fitted.values)
confronto


# Che dire dei coefficienti?
summary(mod.chl)

beta1 = mod.chl$coef[2]
beta1se=summary(mod.chl)$coeff[2, 2]

# Calcolo del test di Wald per beta1
Wald = beta1/beta1se
Wald

#Calcolo p-value = 2(1-PHI(|Wald|))
pvalue=2 * (1 - pnorm(abs(Wald)))
pvalue

# Determinare probabilità di guasto per una data temperatura.
# P0 = exp(beta0+beta1*x0)/(1+exp(beta0+beta1*x0))

# Calcolo predittore lineare beta0+beta1*x0 per x0=55

predlin0 = mod.chl$coef[1] + mod.chl$coef[2]*55

# Calcolo probabilità P0
P0 = exp(predlin0)/(1 + exp(predlin0))
P0

# Calcolo predittore lineare per tutti i valori x(=temp) osservati
predlin = mod.chl$coef[1] + mod.chl$coef[2]*chl$temp
predlin

# Calcolo probabilità P
P = exp(predlin)/(1 + exp(predlin))
P

# Funzione 'predict()'
# Calcolo probabilità P per tutti i valori di x 
predict(mod.chl, type = "response")

# Calcolo predittore lineare per tutti i valori x(=temp) osservati 
# (default su scale lineare)
predict(mod.chl)

# Calcolo probabilità P per dati valori di x
predict(mod.chl, newdata=data.frame(temp=c(55,80)), type='response')

# Y stimati 
ystima=predict(mod.chl, newdata=chl, type="response")>0.5
ystima
as.numeric(ystima)
table(chl$failures,ystima)

# Aggiungere al grafico probabilità stimate
plot(chl$temp, chl$failures)
curve(predict(mod.chl, newdata=data.frame(temp=x), type = "response"), add = TRUE)

# Scelta Y=1(0) 
nofail=1-chl$failures
nofail

mod.chl$coefficients
newmod.chl = glm(nofail ~ temp, family = binomial(link = logit), data=chl)
summary(newmod.chl)

# Commenti ?


# 2.2 Sopravvivenza pazienti in reparto di cura intensiva
# Il file ICUdata.txt contiene i dati su 200 pazienti ammessi in un reparto di cura 
# intensiva per i quali sono state osservate:

# 􏰆 stato: stato del paziente (0 = vivo, 1 = deceduto)
# 􏰆 eta: età in anni del paziente
# 􏰆 causa: causa dell'ammissione (1= programmata, 2 = emergenza)
#  coscienza al momento dell'ammissione (1= no coma, 2= coma)�

# Obiettivo: quali caratteristiche osservate influenzano la probabilità di non sopravvivenza 
# e previsione probabilità di sopravvivenza  

ICU=read.table("ICUdata.txt", TRUE)
head(ICU)
str(ICU)

# Alcune variabili esplicative sono categoriali
# Ricodifica come fattori di var esplicative cat. con nomi (etichette) alle modalità (livelli)

ICU$causa = factor(ICU$causa)
levels(ICU$causa) = c("P", "E")
ICU$coscienza = factor(ICU$coscienza)
levels(ICU$coscienza) = c("SI", "NO")
str(ICU)

table(ICU$causa,ICU$stato)
table(ICU$coscienza, ICU$stato)
plot(ICU$eta, ICU$stato)
boxplot(ICU$eta~ICU$stato)

# Modello di regressione logistica per Yi (non sopravvivenza) 
# e x1 = età, x2= causa ammissione (= 1 se non programmata),
# x3= coscienza (=1 non cosciente)
# pi= (P(Yi=1)), logit(pi)=log(pi/1-pi)=beta0+beta1xi1+beta2xi2+beta3x3
# Stima massima verosimiglianza del modello:

ICU.fit = glm(stato ~ ., family = binomial, data = ICU)
summary(ICU.fit)

# Probabilità di non sopravvivenza aumenta al crescere dell'età, se l'ammissione avviene in una situazione 
# di emergenza e se il paziente è in stato di incoscienza

# Test significatività modello complessivo (nullità singoli coefficienti non accettata)
# Confronto, mediante LRT, modello nullo (solo intercetta) e modello completo

# Stima modello nullo

ICU.fit0 = glm(stato ~ 1, family = binomial, data = ICU)
summary(ICU.fit0)

# Probabilità di non sopravvivenza

#Probabilità stimata di non sopravvivenza nel modello nullo
# è pari a:
p0 = exp(coef(ICU.fit0))/(1 + exp(coef(ICU.fit0)))
p0

# pari alla frequenza di decessi osservata (40/200)
predict(ICU.fit0, type = "response")

anova(ICU.fit)

# Differenza devianza modello nullo e modello completo
LRT= 200.16-145.31
LRT

# P-value distribuzione test
pval=1-pchisq(LRT, 3)
pval

# Possiamo aggiungere al grafico (rispetto all'età) le curve stimate in corrispondenza delle altre esplicative
plot(ICU$eta,ICU$stato, pch = "|", ylim = c(0, 1))

# Differenze rispetto al precedente ?

curve(predict(ICU.fit, newdata = data.frame(eta = x, causa = "P", coscienza = "SI"), 
type = "response"), add = TRUE, lty = 1)
curve(predict(ICU.fit, newdata = data.frame(eta = x, causa = "E", coscienza = "SI"),
type = "response"), add = TRUE, lty = 2)
curve(predict(ICU.fit, newdata = data.frame(eta = x, causa = "P", coscienza = "NO"),
type = "response"), add = TRUE, lty = 3)
curve(predict(ICU.fit, newdata = data.frame(eta = x, causa = "E", coscienza = "NO"),
type = "response"), add = TRUE, lty = 4)
legend(18, 1, lty = c(1, 2, 3, 4), legend = c("P-C", "E-C", "P-N", "E-N"), bty = "n")

# Confrontoi Y osservati e attesi secondo il modello stimato
ystima=predict(ICU.fit, type="response")>0.5
ystima
as.numeric(ystima)
table(ICU$stato,ystima)