# Stima modelli di regressione logistica in R - Esempi

# 2.3 Neonati sotto peso
# Il file Neonati.txt contiene i dati su:
# indicatore nato sottopeso (>= 2500 gr = 0, < 2500 gr = 1)
# Madre fumatrice (0 = no, 1 = si')
# Peso della madre inizio gravidanza (in pounds, 1 pound circa 0.5 kg)

# par(mfrow = c(1,1))
neonati = read.table("neonati.txt", header = TRUE)
head(neonati)

# Dimensioni dati (n.ro casi e variabili)
dim(neonati)

# Obiettivo: sottopeso del neonato dipende dalle caratteristiche della madre ?

str(neonati)
table(neonati$Sottopeso)

# Analisi 'univariata'

table(neonati$Sottopeso,neonati$Madre_Fumatrice)
par(mfrow = c(1,2))
plot(neonati$Peso_della_Madre, neonati$Sottopeso)
boxplot(neonati$Peso_della_Madre~neonati$Sottopeso)

# Modello di regressione logistica per Yi (sottopeso) 
# e x1 = madre fumatrice (1 = Si'), x2= peso madre 
# pi= (P(Yi=1)), logit(pi)=log(pi/1-pi)=beta0+beta1xi1+beta2xi2
# Stima massima verosimiglianza del modello:

neonati.fit = glm(Sottopeso ~ ., data = neonati, family = binomial)
summary(neonati.fit)

anova(neonati.fit)

# Differenza devianza modello nullo e modello completo
LRT= 234.67-224.34
LRT

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

# Calcolo ddds ratios (OR=exp(betaxk)
OR=exp(neonati.fit$coefficients)
OR

# Confronto Y osservati e attesi secondo il modello stimato
table((predict(neonati.fit, type = "response") > 0.5), neonati$Sottopeso)

# Grafico con curve stimate (cambiando il simbolo grafico)
# Nb.:  valori asse y!
plot(neonati$Peso_della_Madre, neonati$Sottopeso, pch = "|", xlab = "Peso della madre", ylab = "Neonato sottopeso",
yaxt= "n")
axis(2, at = c(0, 1), las = 1)
curve(predict(neonati.fit, newdata = data.frame(Madre_Fumatrice = 0, Peso_della_Madre=x),
type = "response"), add = TRUE)
curve(predict(neonati.fit, newdata = data.frame(Madre_Fumatrice = 1, Peso_della_Madre=x),
type = "response"), add = TRUE, lty = 2)
legend(68, 1, lty = c(1, 2), legend = c("Fumatrice No", "Fumatrice Sì"), bty = "n")

# 2.4 Dati su immatricolati ai corsi della ex_Facoltà di Economia, a.a. 2009/2010 (situazione all'inizio II anno)
# Variabili rilevate
# corso:
# EC01/EC01E = Economia aziendale, EC11/EC11E = Economia, commercio internazionale e mercati finanziari
# EC21 = Statistica e informatica per l'azienda, la finanza e l'assicurazione
# sesso: M, F
# anno: anno nascita (ultime 2 cifre)
# provincia residenza
# scuola: 1-8 = magistrale/tecnico prof.le, 9-11 = liceo, 12 = titolo straniero
# voto: voto maturità 
# test ingresso (non obbligatorio e non selettivo): 0 = non superato, 1 = superato (6/16), 
# 2 = non sostenuto
# voto test 
# numEsamiTotali
# creditiTotali
# numEsami (no idoneità)
# crediti (no CFU idoneità)
# mediaEsami: voto medio 
# reiscrizioni: 1 = iscritto 2^ anno 2010/11, 0 = altrimenti
# matematica: 1 superato esame Matematica 1^ anno

s=read.table("studenti.csv", header = TRUE, sep = ";")
head(s)
str(s)

# Dati mancanti 
# Dati mancanti in R codificati con NA
# Se dati importati da file esterno, dati mancanti in genere altre codifiche (*, cella vuota, 9/99/999). 
# Se codici vari (plausibili), meglio sostuire tramite na.strings = codice per dato mancante dentro read.table, 
# oppure mediante <dataframe>[<dataframe> == 999] <- NA
# Per individuare dati mancanti
is.na(s)
is.na(s$corso)
summary(s)

# Alcune funzioni R prevedono argomenti specifici per escludere gli NA
# Esempio: media voto test
mean(s$votoTest)

# In alternativa
mean(s$votoTest, na.rm = TRUE)

# Tuttavia, non tutte le funzioni prevedono simili argomenti per cui è consigliabile 
# avere "pieno controllo" (identificazione e gestione) degli NA

# Provincia di residenza: veri dati mancanti ?

s$provincia=factor(s$provincia, levels=c(levels(s$provincia),  "NAP"))
s$provincia[is.na(s$provincia)] = "NAP"

# Esempio, selezionare solo casi validi su una variabile con is.na()
s1 <- s[!is.na(s$votoTest), ]
head(s1)

# o casi validi su più variabili
s2=s[!is.na(s$votoTest) & !is.na(s$scuola), ]

# o casi validi su tutto data frame (attenzione a riduzione numerosità)

s3=na.omit(s)
head(s3)

# Controlli sui dati (data preparation/data editing) 
# Filtri per selezionare sottoinsiemi di casi/variabili rispetto a specifiche caratteristiche
# Selezione anno nascita per i maschi
s[s$sesso == "M", "anno"]
s$anno[s$sesso == "M"]

# Selezione sesso per i nati <= '90 
s$sesso[s$anno <= 90]

# Selezione dei valori della variabile sesso per cui anno<=90 e provincia=TS:
s$sesso[s$anno <= 90 & s$provincia == "TS"]

# Selezione di sottoinsiemi di casi

sf=s[s$sesso== "F", ]
head(sf)

# Creazione nuova variabile da variabili nel data frame 
# Numero di idoneita' sostenute

idoneita=s$numEsamiTot-s$numEsami
idoneita

# Aggiungere nuova variabile al data frame
s=cbind(s, idoneita)
head(s)

# Valore predittivo del test ingresso sui risultati conseguiti al 1^ anno, controllanodo 
# per caratteristiche dello studente

# Obiettivo: alcuni risultati di carriera sono "prevedibili" a partire dall'esito del test?
# Risultati 1^ anno:
# 1. conseguimento 12 CFU (no idoneità) - soglia richiesta per superamento OFA per coloro che non hanno 
# sostenuto/superato il test
# 2. conseguimento di almeno 30 CFU (no idoneità) - proxy di apprendimento e capacità organizzativa 
# raggiunte dallo studente

# Modello logit per P(Y = 1[>= 12 CFU]) e P(Y = 1 [>= 30 CFU])
# logit(pi)=log(pi/1-pi)=sum(betakxik)

# xik: sesso, voto maturità, esito test, residenza (TS/fuoriTS), scuola (Liceo/NoLiceo)

# Creazione variabili dipendenti

crediti12=as.numeric(s$crediti>=12)
crediti12

crediti30=as.numeric(s$crediti>=30)
crediti30

# Controllo 
table(crediti12,s$crediti)
table(crediti30,s$crediti)
s=cbind(s, crediti12, crediti30)
head(s)

summary(s[,17:18])

# Codifica variabili esplicative

# Residenza TS/FuoriTS

residenza.ric=s$provincia

# Creazione dei livelli variabile residenza.ric
l <- levels(residenza.ric)
# Individuazione livelli non TS
ind <- c(which(l != "TS"))
# Attribuzione label "FuoriTS" ai livelli non TS 
levels(residenza.ric)[ind] <- "FuoriTS"
s <- cbind(s, residenza.ric)
head(s)

# Scuola Liceo (codici 9-11)/NoLiceo
Liceo=as.numeric(s$scuola>=9 & s$scuola <=11)
Liceo
s=cbind(s, Liceo)
head(s)

# Analisi esplorativa
# Y1(crediti12), Y2(crediti30), sesso, voto maturità, esito test, residenza (TS/fuoriTS), 
# scuola (Liceo/NoLiceo)

table(crediti12)
table(crediti30)
table(s$sesso)
voto_alto=as.numeric(s$voto>=80)
table(voto_alto)
table(s$test)
table(s$residenza.ric)

par(mfrow = c(1, 2))
plot(s$voto, s$crediti12, yaxt="n")
axis(2, at = c(0, 1), las = 1)
boxplot(s$voto~s$crediti12)
plot(s$voto, s$crediti30, yaxt="n")
axis(2, at = c(0, 1), las = 1)
boxplot(s$voto~s$crediti30)

# Analisi (Modelli) univariate
s.sesso = glm(s$crediti12 ~ s$sesso, family = binomial(link = logit))
summary(s.sesso)

s.voto = glm(crediti12 ~ s$voto, family = binomial(link = logit))
summary(s.voto)

s.test = glm(crediti12 ~ as.factor(s$test), family = binomial(link = logit))
summary(s.test)

s.residenza = glm(crediti12 ~ s$residenza.ric, family = binomial(link = logit))
summary(s.residenza)

s.Liceo = glm(crediti12 ~ s$Liceo, family = binomial(link = logit))
summary(s.Liceo)

# Modello multivariato (tutte var esplicative)

s.multi = glm(crediti12 ~ s$sesso+s$voto+as.factor(s$test)+s$residenza.ric+s$Liceo, 
              family = binomial(link = logit))
summary(s.multi)

# Modello multivariato senza var sesso
s.multi1 = glm(crediti12 ~ s$voto+as.factor(s$test)+s$residenza.ric+s$Liceo, 
              family = binomial(link = logit))
summary(s.multi1)

# Modello multivariato senza test 
s.multi2 = glm(crediti12 ~ s$sesso+s$voto+s$residenza.ric+s$Liceo, 
               family = binomial(link = logit))
summary(s.multi2)

# Modello multivariato senza voto
s.multi3 = glm(crediti12 ~ s$sesso+as.factor(s$test)+s$residenza.ric+s$Liceo, 
               family = binomial(link = logit))
summary(s.multi3)
# test: 0 = non superato, 1 = superato (6/16), 2 = non sostenuto

# Modello multivariato con voto test 
s.multi4 = glm(crediti12 ~ s$sesso+s$voto+s$votoTest+s$residenza.ric+s$Liceo, 
               family = binomial(link = logit))
summary(s.multi4)

# Ricodifica test

test.ric=s$test>=2
test.ric
table(test.ric)

# Modello multivariato con test ricodificato  
s.multi5 = glm(crediti12 ~ s$sesso+s$voto+s$residenza.ric+s$Liceo+test.ric, 
               family = binomial(link = logit))
summary(s.multi5)

s$Liceo=factor(s$Liceo)
levels(s$Liceo) = c("No Liceo", "Liceo")
str(s$Liceo)

test.ric=factor(test.ric)
levels(test.ric) = c("Test sostenuto", "Test non sostenuto")
str(test.ric)

# Modello multivariato senza sesso e test ricodificato
s.multi6 = glm(crediti12 ~ voto+residenza.ric+Liceo+test.ric, data=s,
               family = binomial(link = logit))
summary(s.multi6)

# Grafico probabilità stimata Modello specificazione multi6

par(mfrow=c(1,1))
plot(s$voto, s$crediti12, pch = "|", xlab = "Voto maturità", ylab = "Prob. 12 CFU",yaxt = "n")
curve(predict(s.multi6, newdata = data.frame(voto = x, residenza.ric = "TS", test.ric = "Test sostenuto", Liceo = "Liceo"),
              type = "response"), add = TRUE, col = 2, lwd = 2, lty = 1)
curve(predict(s.multi6, newdata = data.frame(voto = x, residenza.ric = "TS", test.ric = "Test non sostenuto", Liceo = "Liceo"),
              type = "response"), add = TRUE, lty = 2)
curve(predict(s.multi6, newdata = data.frame(voto = x, residenza.ric = "TS", test.ric = "Test sostenuto", Liceo = "No Liceo"),
              type = "response"), add = TRUE, lty = 3)
curve(predict(s.multi6, newdata = data.frame(voto = x, residenza.ric = "TS", test.ric = "Test non sostenuto", Liceo = "No Liceo"),
              type = "response"), add = TRUE, lty = 4)
curve(predict(s.multi6, newdata = data.frame(voto = x, residenza.ric = "FuoriTS", test.ric = "Test sostenuto", Liceo = "Liceo"),
              type = "response"), add = TRUE, lty = 5)
curve(predict(s.multi6, newdata = data.frame(voto = x, residenza.ric = "FuoriTS", test.ric = "Test non sostenuto", Liceo = "Liceo"),
              type = "response"), add = TRUE, lty = 6)
curve(predict(s.multi6, newdata = data.frame(voto = x, residenza.ric = "FuoriTS", test.ric = "Test sostenuto", Liceo = "No Liceo"),
              type = "response"), add = TRUE, lty = 7)
curve(predict(s.multi6, newdata = data.frame(voto = x, residenza.ric = "FuoriTS", test.ric = "Test non sostenuto", Liceo = "No Liceo"),
              type = "response"), add = TRUE, lty = 8)
legend(55, 1, lty = c(1, 2, 3, 4, 5, 6, 7, 8), legend = c("TS-Test-L", "TS-NoTest-L", "TS-Test-NoL", "TS-NoTest-NoL", 
       "FuoriTS-Test-L", "FuoriTS-NoTest-L", "FuoriTS-Test-NoL", "FuoriTS-NoTest-NoL"), bty = "n")