# GLM
# priklad Potemnik skladistni
rm(list=ls())
load("beetle.RData")
beetle
attach(beetle)
proportions <- killed/population

# Logisticka regrese
# Vykresleni dat
x11(width=12,height=9)
plot(dose,proportions,pch=20,xlab="CS2",ylab="umrtnost")

# definice GLM modelu
glmmod <- glm(formula = cbind(killed,population - killed) ~ dose, family = binomial(logit),data = beetle)
glmmod <- glm(formula = proportions ~ dose, family = binomial(logit), weights=population,data = beetle) # alternativa

summary(glmmod)

# vykresleni regresni krivky
t <- seq(from = min(dose), to = max(dose), length.out=100)
lines(t,predict.glm(glmmod,data.frame(dose=t),type="response"),col="red")

# overeni normality residui
plot(glmmod,which=2)
shapiro.test(residuals(glmmod))

# Probitovy model
# Vykresleni dat
plot(dose,proportions,pch=20,xlab="CS2",ylab="umrtnost")

# definice GLM modelu
glmmod <- glm(formula = cbind(killed, population - killed) ~ dose, family = binomial(probit),data = beetle)
summary(glmmod)

# vykresleni regresni krivky
t <- seq(from=min(dose),to=max(dose),length.out=100)
lines(t,predict.glm(glmmod,data.frame(dose=t),type="response"),col="red")

# overeni normality residui
plot(glmmod,which=2)
shapiro.test(residuals(glmmod))

# komplementarni log-log
# Vykresleni dat
plot(dose,proportions,pch=20,xlab="CS2",ylab="umrtnost")

# definice GLM modelu
glmmod <- glm(formula = cbind(killed, population - killed) ~ dose, family = binomial(link="cloglog"),data = beetle)
summary(glmmod)

# vykresleni regresni krivky
t <- seq(from=min(dose),to=max(dose),length.out=100)
lines(t,predict.glm(glmmod,data.frame(dose=t),type="response"),col="red")

# overeni normality residui
plot(glmmod,which=2)
shapiro.test(residuals(glmmod))
detach(beetle)

#--------------------------------------------------------------------------------
# priklad AIDS ve VB
rm(list=ls())
load("aids.RData")
head(aids)
attach(aids)

# Poissonovska regrese

time <- 1:36
# vykresleni dat
plot(time,number,type="p",pch=20,ylab="number of AIDS", xlab="time",main="Poisson regression")

# definice modelu
ind <- which(number!=0)
num <- number[ind]
tim <- time[ind]
glm1 <- glm(num ~ 0 + tim, family=poisson("identity"))
#glm1 <- lm(number ~ time)
summary(glm1)

glm2 <- glm(number ~ time, family=poisson("log"))
summary(glm2)
glm3 <- glm(number ~ time, family=poisson("sqrt"))
summary(glm3)

# vykresleni regresnich krivek
#points(time,fitted(glm1),col=2)
xx<-seq(0,37,length.out=100)
yA.logit<-predict(glm1,list(tim=xx),type="response")
lines(xx,yA.logit,col=2,lwd=2)
yB.logit<-predict(glm2,list(time=xx),type="response")
lines(xx,yB.logit,col=3,lwd=2)
yC.logit<-predict(glm3,list(time=xx),type="response")
lines(xx,yC.logit,col=4,lwd=2)
legend("topleft",col=c(2,3,4),lty=c(1,1,1),lwd=c(2,2,2),legend=c("linear","log-linear","square-root"))

# overeni normality residui
plot(glm1,which=2)
shapiro.test(residuals(glm1))
plot(glm2,which=2)
shapiro.test(residuals(glm2))
plot(glm3,which=2)
shapiro.test(residuals(glm3))

#  x<-residuals(glm3)
#   library(nortest)
#   shapiro.test(x)
#   ad.test(x)
#   cvm.test(x)
#   pearson.test(x)
#   sf.test(x)
#   lillie.test(x)

detach(aids)
# ---------------------------------------------------------------------------