#########################
# 5. I modelli lineari generalizzati
#########################




##  raul = read.csv(file.choose(), header = TRUE)

www = "http://www.biostatisticaumg.it/dataset/raul.csv"
raul = read.csv(www, header = TRUE)
attach(raul)
head(raul)



### 5.1    I dettagli da conoscere

### Vocabolario 5.1 — modello lineare generalizzato.

### 5.1.1    La funzione di collegamento:


### Esercizio 5.1 Disegnate una sigmoide.  
# Vedete che la ascissa y (e sottolineo ascissa, ossia,l’asse orizzontale) 
# è una variabile numerica continua, mentre  sigmoide  diventa un valore di probabilità, 
# sempre compreso tra 0 ed 1 (sull’asse verticale delle ordinate):

y = seq(-6, 6, 0.1)
p = exp(y) / (1 + exp(y))
plot(y, p)




## 5.1.3    Interpretare una regressione logistica
### A. il caso di un fattore


Ca125High = ( Ca125 > median(Ca125) )

table(Esito, Ca125High)


### odds ratio

(818 * 38) / (750 * 4)


### log - odds ratio:

log((818 * 38) / ( 750 * 4))

### sorpresa: nel summary ...


indicedirischio3 = Esito ~  Ca125High
logistico3 = glm(indicedirischio3, family = binomial)
summary(logistico3)



exp(2.338) # odds ratio

### altra sorpresa ...

p = 4 / (818 + 4)
p
exp(-5.3206) / (1 + exp(-5.3206))

log(p/(1-p))

exp(-5.3206) / (1 + exp(-5.3206))



### infine


-5.3206 + 2.3381  # [1] -2.9825
exp(-2.9825) / (1 + exp(-2.9825))
38/(750+38)





### B. il caso di una variabile numerica



indicedirischio2 = Esito ~  Ca125
modello = glm(indicedirischio2, family = binomial)
summary(modello)

## Per esempio,ipotizziamo Ca125 = 132... 


-4.132 + 0.014 * 132
 exp(-2.284) / (1+exp(-2.284))
 
 
 ??? siamo sicuri di essere in presenza di un modello predittivo ??? 
 

####  Esercizio 5.2 Commentate l’andamento di questa sigmoide


plot(Ca125, (as.numeric(Esito)-1))
xx = seq(0,300)
pl = -4.132 + 0.014 * xx
yy = exp(pl)/(1+exp(pl))
lines(xx,yy, col = "orange")


### Esercizio 5.3 caccia all'errore:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3594101/




###  5.1.4    Problemi con lo standard error


www = "http://www.biostatisticaumg.it/dataset/tossicologia.csv"
tossicologia = read.csv(www, header = TRUE)
attach(tossicologia)
table(causa)

alcolemia

prova = glm(alcolemia ~ genere + eta + causa + tossicologico,
            family = binomial)
            
summary(prova)  ## con "poca" info non si fa statistica



detach(tossicologia) ## importante detach per non fare caos



www = "http://www.biostatisticaumg.it/dataset/epilessie.csv"
epilessie = read.csv( www , header = TRUE )
attach(epilessie)
head(epilessie)


## collinearità:

table(eoneuro, causa)


prova = glm(prognosi ~ primoEEG + causa + eoneuro + svilpsicom+ etaallesordio, family = binomial)
summary(prova)

detach(epilessie) ## importante detach per non fare caos




###5.1.5    La sovradispersione


www = "http://www.biostatisticaumg.it/dataset/gdm.csv"
gdm = read.csv(www , header = TRUE)
attach(gdm)
head(gdm)




prova = glm( factor(GDM) ~ Ngrav  + BMIpreGrav , family = binomial)
summary(prova)

prova2 = glm( factor(GDM) ~ Ngrav  + BMIpreGrav , family = quasibinomial)
summary(prova2)


detach(gdm)




####################

5.2    La meta finale: valutare l’accuratezza del modello logistico

####################



www = "http://www.biostatisticaumg.it/dataset/roma.csv"
roma = read.csv(www, header = TRUE)
# roma = read.csv(file.choose(), header = TRUE)
attach(roma)
head(roma)
tail(roma)





moorerelation = Histology ~ Menopause * logHE4 + Menopause * logCA125
mooremodel = glm(moorerelation, family = binomial)
summary(mooremodel) ### 117.49




maximalrelation = Histology ~ logHE4 + logCA125 + logCA19.9 + logCEA + AgePatient + Menopause

maximalmodel = glm(maximalrelation, family = binomial)
step(maximalmodel)


attempt0 = Histology ~ Menopause + logHE4 + logCA125
attempt1 = Histology ~ Menopause * logHE4 + logCA125
attempt2 = Histology ~ logHE4 + Menopause * logCA125
attempt3 = Histology ~ Menopause + logHE4 * logCA125
moorerelation = Histology ~ Menopause * logHE4 + Menopause * logCA125

modelattempt0 = glm(attempt0, family = binomial)
modelattempt1 = glm(attempt1, family = binomial)
modelattempt2 = glm(attempt2, family = binomial)
modelattempt3 = glm(attempt3, family = binomial)
mooremodel = glm(moorerelation, family = binomial)


AIC(modelattempt0, modelattempt1, modelattempt2, modelattempt3, mooremodel) 



attempt4 = Histology ~ Menopause + logHE4 + logCA125 + I(logCA125^2)
attempt5 = Histology ~ Menopause + logHE4 + I(logHE4^2) + logCA125 

modelattempt4 = glm(attempt4, family = binomial)
modelattempt5 = glm(attempt5, family = binomial)

AIC(modelattempt0, modelattempt4, modelattempt5, mooremodel) 





modelattempt0
summary(modelattempt0)





## curva ROC

modelattempt0


#  (Intercept)  Menopausepost         logHE4       logCA125  
# -14.3770         0.9378         2.3382         0.6845  


###  Esercizio 5.3

(-1 + as.numeric(Menopause)) 



Menopause01 = -1 + as.numeric(Menopause)

iSN =  -14.38 + 0.94 * Menopause01 + 2.34 * logHE4 + 0.68 * logCA125

boxplot(iSN ~ Histology)




install.packages("pROC")
library(pROC)



iSNroc = roc(Histology ~ iSN, auc=TRUE)
iSNroc



plot(iSNroc)



iSNroc$response
iSNroc$predictor
iSNroc$sensitivities
iSNroc$specificities
iSNroc$thresholds


iSNroc$sensitivities^2 + iSNroc$specificities^2

plot(iSNroc$thresholds, iSNroc$sensitivities^2 + iSNroc$specificities^2, "l")


### esercizi


names(raul)
head(raul)

umgrisk = Ldhtre + 24 / Ldhuno

annalisa = roc(Esito ~ umgrisk, auc=TRUE)
annalisa
plot(annalisa)