# GLM
# overdispersion
# priklad vcely
rm(list=ls())
load("bees.RData")
head(bees)
attach(bees)

# vykresleni dat
x11(width=11,height=9)
plot(time,number, cex = 0.75, main="Bees activity")

# Poissonovska regrese
m1<-glm(number~ time + I(time^2) ,data=bees,family=poisson)
summary(m1)

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

# jaky bude stupen polynomu?
# mod.opt<-step(m1,scope = list(upper = ~ time+I(time^2)+I(time^3)+I(time^4)+I(time^5)+I(time^6), lower = ~ 1))
# summary(mod.opt)
# m1 <- mod.opt
# anova(mod.opt,m1,test="Chisq")

# vykresleni modelu
plot(time,number, cex = 0.75, main="Bees activity")
ab<-range(time)
xx<-seq(ab[1],ab[2],length.out=100)
# jeste vykreslime i interval spolehlivosti
predicted.log <- predict(m1,list(time=xx),type = "link", se.fit = T)
CI.L.log <- exp(predicted.log$fit - 1.96 * predicted.log$se.fit)
CI.H.log <- exp(predicted.log$fit + 1.96 * predicted.log$se.fit)
x<-c(xx,rev(xx))
y<-c(CI.L.log,rev(CI.H.log))
polygon(x,y,col="gray75",border="gray75")
yy<-predict(m1,list(time=xx),type="response")
lines(xx,yy,col="red")

m1a <- update(m1,family=quasipoisson)
# jeste vykreslime i interval spolehlivosti
predicted.log <- predict(m1a,list(time=xx),type = "link", se.fit = T)
CI.L.log <- exp(predicted.log$fit - 1.96 * predicted.log$se.fit)
CI.H.log <- exp(predicted.log$fit + 1.96 * predicted.log$se.fit)
x<-c(xx,rev(xx))
y<-c(CI.L.log,rev(CI.H.log))
polygon(x,y,col="gray75",border="gray75")
yy<-predict(m1a,list(time=xx),type="response")
lines(xx,yy,col="red")

detach(bees)

#-----------------------------------------------------------------------
# priklad melanoma
rm(list=ls())
load("melanoma.RData")
melanoma
attach(melanoma)

# Vykresleni dat - tzv. mosaic plot
plot(xtabs(Count ~ Type + Site),main="Mosaic plot")  # "xtabs" vytvori kontingencni tabulku

# da se testovat i klasickym Pearsonem
# chisq.test(xtabs(Count ~ Type + Site))

# model pro nezavisla data
m1 <- glm(Count ~ Type + Site, family=poisson, data=melanoma)
# model pro zavisla
m2 <- glm(Count ~ Type * Site, family=poisson, data=melanoma)

# stejne jako Pearson
#sum(residuals(m1, type = "pearson")^2)

# uzitecna tabulka pro male kontingencni tabulky
round(xtabs(residuals(m1) ~ Type + Site), 1)

anova(m1, m2, test = "Chisq")

# graficke porovnani
plot(xtabs(fitted(m1) ~ Type + Site, data=melanoma),main="Independent data")
plot(xtabs(fitted(m2) ~ Type + Site, data=melanoma),main="Full model")

detach(melanoma)

# ---------------------------------------------------------------------------