# %%%%%%%%%%%%%%%%%%%%%%
# ANALISI DATI SUICIDIO
# %%%%%%%%%%%%%%%%%%%%%%
# Regressione logistica 
# VD = variabili suicidarie, PENSIERO SUICIDARIO-PIANIFICAZIONE SUICIDARIA-TENTATO SUICIDIO
# Predittori di interesse: dipendenza da sostanze per ciascuna sostanza
# 2008
# IRSEX      =   SEX (1-M; 2-F)
# MHSUITHK   =   PENSIERO SUICIDARIO
# MHSUIPLN   =   PIANIFICAZIONE SUICIDARIA
# MHSUITRY   =   Tentamen
# NDSSDNSP   =   Nicotine Dependence Syndrome Scale (NDSS)
# DEPNDALC   =   ALCOHOL
# DEPNDCOC   =   COCAINE
# DEPNDHER   =   HEROINE
# DEPNDANL   =   PAIN RELIEVERS
# DEPNDSED   =   SEDATIVES
# DEPNDSTM   =   STIMULANTS
# DEPNDHAL   =   HALLUCINOGENS
# DEPNDINH   =   INHALANTS
# DEPNDMRJ   =   MARIJUANA
# DEPNDTRN   =   TRANQUILIZERS

## Ricaviamo la directory locale della sessione corrente R
## e mettiamo qui il file RData:
getwd()

## Carichiamo il dataset
load("dataset08.RData")

# %%%%%%%%%%%%%%%%%%%%%%
# Vecchi Preamboli...
## Convertiamo IRSEX in fattore: 1-Maschio; 2-Femmina
dataset08_reduced$IRSEX<-factor(dataset08_reduced$IRSEX,labels=names(attributes(dataset08_reduced$IRSEX)$labels))
## ATTENZIONE: Convertiamo "MISSING = -9" in SPSS --> NA in R
dataset08_reduced$MHSUITHK<-ifelse(dataset08_reduced$MHSUITHK=="-9",NA,dataset08_reduced$MHSUITHK)
dataset08_reduced$MHSUIPLN<-ifelse(dataset08_reduced$MHSUIPLN=="-9",NA,dataset08_reduced$MHSUIPLN)
dataset08_reduced$MHSUITRY<-ifelse(dataset08_reduced$MHSUITRY=="-9",NA,dataset08_reduced$MHSUITRY)
## ATTENZIONE: Rimuoviamo i missing data nelle variabili di "SUICIDAL OUTCOME" [MHSUITHK, MHSUIPLN, MHSUITRY]
dataset08_reduced=na.omit(dataset08_reduced)
summary(dataset08_reduced)
# %%%%%%%%%%%%%%%%%%%%%%

## Regressione logistica: PENSIERO SUICIDARIO
Mod.PS_08<-glm(MHSUITHK ~ IRSEX + NDSSDNSP + DEPNDALC + DEPNDCOC + DEPNDHER + DEPNDANL + DEPNDSED + DEPNDSTM + 
                          DEPNDHAL + DEPNDINH + DEPNDMRJ + DEPNDTRN, data = dataset08_reduced, family=binomial)
summary(Mod.PS_08)
## Regressione logistica: PIANIFICAZIONE SUICIDARIA
Mod.PianS_08<-glm(MHSUIPLN ~ IRSEX + NDSSDNSP + DEPNDALC + DEPNDCOC + DEPNDHER + DEPNDANL + DEPNDSED + DEPNDSTM + 
                             DEPNDHAL + DEPNDINH + DEPNDMRJ + DEPNDTRN, data = dataset08_reduced, family=binomial)
summary(Mod.PianS_08)
## Regressione logistica: TENTATO SUICIDIO
Mod.Tent_08<-glm(MHSUITRY ~ IRSEX + NDSSDNSP + DEPNDALC + DEPNDCOC + DEPNDHER + DEPNDANL + DEPNDSED + DEPNDSTM + 
                            DEPNDHAL + DEPNDINH + DEPNDMRJ + DEPNDTRN, data = dataset08_reduced, family=binomial)
summary(Mod.Tent_08)

## probabilità previste dal modello
prob = predict(Mod.PS_08, dataset08_reduced,type="response")
prob[1:6]
length(prob)

## criterio decisionale per la probablità soglia "y=1"
k = 0.10085327

## assegnamo a ciascuna osservazione "y=1" se prob > k
Classificazione.binaria<-ifelse(prob>k,1,0)

## costruiamo la tabella stimolo-risposta
matching.st.risp <-rep(NA,length(prob))
for(i in 1:length(matching.st.risp)){
if(dataset08_reduced$MHSUITHK[i]==0 & Classificazione.binaria[i]==0) matching.st.risp[i]<-"CorrRej"
if(dataset08_reduced$MHSUITHK[i]==0 & Classificazione.binaria[i]==1) matching.st.risp[i]<-"F"
if(dataset08_reduced$MHSUITHK[i]==1 & Classificazione.binaria[i]==1) matching.st.risp[i]<-"H"
if(dataset08_reduced$MHSUITHK[i]==1 & Classificazione.binaria[i]==0) matching.st.risp[i]<-"Miss"
                                     }
Tab.st.risp<-table(matching.st.risp)
Tab.st.risp
H<-Tab.st.risp[3]/(Tab.st.risp[3]+ Tab.st.risp[4])
F<-Tab.st.risp[2]/(Tab.st.risp[2]+ Tab.st.risp[1])
dprime=qnorm(H)-qnorm(F);
names(dprime)="d'"
## visualizziamo i risultati
dprime;F;H


# %%%%%%%%%%%%%%%%%%%%%%%%%%%
# Curve ROC per modelli Logit
# con pacchetto pROC:
install.packages("pROC",dep=TRUE)
library(pROC)
# %%%%%%%%%%%%%%%%%%%%%%%%%%%

ROC_<-roc(dataset08_reduced$MHSUITHK,prob,smooth=FALSE,ci=TRUE)
## Il comando roc() genera diversi valori probablità soglia "k"
ROC_$thresholds
## dai quali determina poi analogamente a quanto visto sopra, coppie (F,H)
## denominate secondo la nomenclatura medica H = sensitivity, F = 1-specificity
## Hit (o "sensitivity"):
ROC_$sensitivities
## Falsi allarmi (o 1-specificity):
1-ROC_$specificities

## scegliamo ad esempio il criterio k=0.10085327
k=0.10085327
which(round(ROC_$thresholds,4)==round(k,4))
# [1] 19
## visualizziamo la coppia (F;H) relativa a questo criterio che sarà naturalmente uguale a quella calcolate "a mano"
ROC_$sensitivities[19]
1-ROC_$specificities[19]
## d-prime
qnorm(ROC_$sensitivities[19])-qnorm(1-ROC_$specificities[19])
## uguale al valore calcolato
dprime

## Scegliendo diversi valori per k, si generano diverse coppie (F,H)
## Il valore d' associato a queste coppie non è stabile, come ci si aspetterebbe dalla teoria SDT. 
## Le coppie (F,H) si dispongono su ROC empiriche solo approssimativamente curvilinee

## Scegliamo ad esempio il criterio 0.17938000
which(round(ROC_$thresholds,4)==round(0.17938000,4))
# [1] 62
ROC_$sensitivities[62]
1-ROC_$specificities[62]
## d-prime leggermente diverso...
qnorm(ROC_$sensitivities[62])-qnorm(1-ROC_$specificities[62])

## scegliamo ad esempio il criterio 0.38964615
which(round(ROC_$thresholds,4)==round(0.38964615,4))
# [1] 136
ROC_$sensitivities[136]
1-ROC_$specificities[136]
## d-prime ancora leggermente diverso...
qnorm(ROC_$sensitivities[136])-qnorm(1-ROC_$specificities[136])


## Tutti d' diversi!
plot(x=ROC_$thresholds,
     y=qnorm(ROC_$sensitivities)-qnorm(1-ROC_$specificities),
     ylim=c(0,4),xlab="soglia di probabilità k",ylab="d'",bty="n")
dprime<-qnorm(ROC_$sensitivities)-qnorm(1-ROC_$specificities)
summary(dprime) 
mean_dprime<-mean(dprime[!is.na(dprime) & !is.infinite(dprime)])
mean_dprime 
abline(h=mean_dprime,col="red",lwd=2)

## Ne consegue che il locus delle coppie (H,F) generate dal modello e basate su diversi valori k
## sarà solo approssimativamente una curva ROC, che chiameremo curva ROC empirica:
par(pty = "s") # area di disegno "quadrata"
plot(x=1-ROC_$specificities,y=ROC_$sensitivities,col="red",type="l",lwd=2)
abline(a=0,b=1,lty="dashed")

## sovrapponiamo la cruva ROC teorica per d' = 0.7619664
F=seq(0,1,by=0.01)
H=pnorm( 0.7619664 + qnorm(F) )
points(x=F,y=H,col="red",type="l",lty="dotted")

## sovrapponiamo la cruva ROC per un modello "ridotto", contenente solo Alcool
Mod.PS_08_RID<-glm(MHSUITHK ~ IRSEX + DEPNDALC, data = dataset08_reduced, family=binomial)
prob_RID = predict(Mod.PS_08_RID, dataset08_reduced,type="response")
ROC_RID<-roc(dataset08_reduced$MHSUITHK,prob_RID,smooth=FALSE,ci=TRUE)
points(x=1-ROC_RID$specificities,y=ROC_RID$sensitivities,col="orange",type="l",lwd=2)

## Test per confrontare i valori AUC nei modelli "nested", completo e ridotto:
roc.test(ROC_,ROC_RID,paired=TRUE)
roc.test(ROC_,ROC_RID,paired=TRUE,method="bootstrap")