heart=read.csv("heart.csv", header=TRUE)
heart
### ck= razina enzima "creatinine kinase"
### ha= broj pacijenata koji su u buducnosti dozivjeli 
##      srcani udar
### ok= broj onih koji nisu

### cilj je ispitati utjecaj ck na vjerojatnost pojave
##  srcanog udara

### koristimo GLM uz binomnu familiju razdioba
##  uz kanonsku funkciju veze
### odziv je proporcija ljudi sa srcanim udarom
##  za danu razinu ck (vidi predavanja), tj.
##  Y_i = B(m_i, p_i)/m_i uz 
##  p_i= exp(eta_i)/(1+exp(eta_i)) za eta_i= b_0 + b_1*ck_i
### parametar disperzije je phi_i=1/w_i za tezine w_i= m_i 
##  (dakle phi=1) -- vidi predavanja

### prvo crtamo proporcije u ovisnosti o ck
p <- heart$ha/(heart$ha+heart$ok) ## p_i je dakle y_i
plot(heart$ck,p,xlab="Creatinine kinase level",
     ylab="Proportion Heart Attack")



### priprema podataka za modeliranje
m=heart$ha + heart$ok ### ovo su m_i-evi
heart=cbind(heart, m=m, prop=p )
heart

### prilagodba modela koristeci funkciju glm
?glm

mod.0 <- glm(prop ~ ck, weights=m, 
             family=binomial(link=logit),
             data=heart)
mod.0

## Prvo provjeramo ima li naznaka da su narusene 
## pretpostavke naseg GLM-a
### To se radi analizom reziduala ("diagnostic plots") na
## slican nacin kao u normalnom lin. modelu
### umjesto klasicnih (standardiziranih) reziduala koriste
##  se tzv. reziduali devijance koji bi se, ako vrijede 
## pretpostavke GLM-a, trebali ponasati kao realizacije
## n N(0,phi) slucajnih varijabli

par(mfrow=c(2,2))
plot(mod.0)
### predicted values su zapravo linearni prediktori
### ipak, plot funkcija daje reziduale kao da se radi o 
##  normalnom linearnom modelu 

### napravit cemo gornji lijevi diagnosticki graf sami
par(mfrow=c(1,1))
### type="link" daje procijenjene linearne prediktore
##  hat{eta}_i, a residuals funkcija po defaultu
##  reziduale devijance
lin_predikt=predict(mod.0, type="link")
rez_dev=residuals(mod.0)
plot(lin_predikt, rez_dev,
     ylab="reziduali devijance", 
     xlab="linearni prediktor (eta_i)",
     pch=20)
abline(h=0, lty=2)
### kako bi provjerili ima li naznaka nekog trenda razlicitog
## od 0, tipicno radimo nekakvu neparametarsku regresiju
library(splines)
spl.0=smooth.spline(lin_predikt, rez_dev)
grid=seq(-2,12, by=0.1)
points(grid, predict(spl.0, x=grid)$y, type="l")
### trend implicira da imamo problem s pretpostavkom
## nezavisnosti odziva, a to je tipicno rezultat izostavljanja
## vazne kovarijate i/ili veza izmedu nekih kovarijata i odziva
## je npr. kvadratna ili kubicna umjesto linearne

### uocimo da graf procijenjenih vjerojatnosti u modelu
## u odnosu na ck, zajedno sa stvarnim proporcijama
## ne ukazuje na neke probleme u modelu

plot(heart$ck, heart$prop, xlab="Creatinine kinase level",
     ylab="Proportion Heart Attack")
grid2=seq(20,500, by=1)
lines(grid2, predict(mod.0, newdata=data.frame(ck=grid2), 
                     type="response"))


## --> diagnostika je bitna jer ako su pretpostavke GLM-a
## narusene, nema smisla raditi npr, testove usporedbe modela
## buduci da su one bazirane na tim pretpostavkama



### probajmo
mod.1 <- glm(prop ~ poly(ck, degree = 3), weights=m, 
             family=binomial(link=logit),
             data=heart)

lin_predikt=predict(mod.1, type="link")
rez_dev=residuals(mod.1)
plot(lin_predikt, rez_dev,
     ylab="reziduali devijance", 
     xlab="linearni prediktor (eta_i)",
     pch=20)
abline(h=0, lty=2)
spl.1=smooth.spline(lin_predikt, rez_dev)
grid=seq(-4,9, by=0.1)
points(grid, predict(spl.1, x=grid)$y, type="l")

### nema naznaka trenda

### fit se cini bolji nego za mod.0:
plot(heart$ck, heart$prop, xlab="Creatinine kinase level",
     ylab="Proportion Heart Attack")
lines(grid2, predict(mod.0, newdata=data.frame(ck=grid2), 
                     type="response"), col="red")
lines(grid2, predict(mod.1, newdata=data.frame(ck=grid2), 
                     type="response"))

### mozemo testirati je li dovoljan mod.0 u odnosu na mod.1
## jer su to ugnjezdeni modeli

### buduci da je phi=1, devijanca je 
##  2*(logvjer(satur. model) - logvjer(prilagodeni model))
### manja devijanca znaci bolju prilagodbu podacima,
### devijanca kod GLM igra istu ulogu kao RSS kod normalnog
##  linearnog modela
### razlika devijanci (tj. smanjenje devijance) 
##  kod usporedbe dva ugnjezdena modela 
##  u ovom slucaju je
##  2*(logvjer(veci model) - logvjer(manji model))
### velike razlike upucuju na to da je potreban veci model

anova(mod.0, mod.1, test="Chisq")
### Residual deviance = devijanca prilagodjenog modela,
## u ovom slucaju to su Model 1 i Model 2
### dakle p-vrijednost naseg testa je 
1-pchisq(32.676, df=2) ### df=2 jer smo dodali dvije kovarijate
# [1] 8.025956e-08
### dakle, ta p-vrijednost je dana u zadnjem stupcu

### mala p-vrijednost sugerira da odbacimo H_0 da je 
## mod.0 dovoljan u korist mod.1

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