# GLM
# priklad motak zere tetreva
rm(list=ls())
load("motak.RData")

# budeme modelovat pomoci Gamma distribuce, ocistime od nulovych hodnot
cisty <- motak[motak$pomer.zkonzumovanych!=0,]
attach(cisty)

#Vykresleni dat
x11(16,9)
plot(tetrev.pocty,pomer.zkonzumovanych,pch=20,xlab="pocty tetreva",ylab="pomer zkonzumovanych tetrevu",main="Konzumace tetreva motakem")

# definice GLM modelu
glmmod <- glm(pomer.zkonzumovanych ~ I(tetrev.pocty^(-3)),data=cisty,family=Gamma(link="inverse"))
summary(glmmod)

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

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

# priklad toxic
rm(list=ls())
load("toxic.RData")
# posouzeni homogenity rozptylu
attach(toxic)
plot(METHOD,VOL,xlab="metoda",ylab="objem")  # graficky
library(car)
leveneTest(VOL~METHOD)        # levenuv test
bartlett.test(VOL~METHOD) # bartlettuv test
# testy nezamitaji homogenitu rozptylu => muzeme pouzit klasickou ANOVA

summary(mod.lm <- lm(VOL ~ METHOD, data = toxic))
anova(mod.lm)  # metoda ma vliv na objem (dal by se pouzit i t-test - stejne vysledky)

# hledani optimalniho LM modelu
lmmod.min <- lm(VOL~1,data=toxic)
lmmod.opt<-step(lmmod.min,scope = list(upper = ~(TEMP+CAT)*METHOD, lower = ~1))
summary(lmmod.opt)
# podle Step-wise procedury je nejvhodnejsi model VOL ~ METHOD + TEMP

mod1 <- lm(VOL~METHOD,data=toxic)
anova(lmmod.opt,mod1)
# podle anovy staci VOL ~ METHOD

# overeni normality residui
plot(lmmod.opt,which=2)
shapiro.test(resid(lmmod.opt))

plot(mod1,which=2)
shapiro.test(resid(mod1))

# Trocha te grafiky
matice <- data.frame(vol = VOL,meth = as.vector(METHOD, mode="numeric")-1,temp = TEMP)
attach(matice)
library(rgl) 
open3d()
scatter3d(vol ~ meth + temp, data = matice)
# model bez interakci

model1 <- lm(vol ~ (meth+temp)^2,data=matice)
b<-model1$coefficients
x <- c(0,1)
nn <- 10
y <- seq(from = min(temp), to=max(temp),length.out = nn)
X <- matrix(x,2,nn)
Y <- matrix(y,2,nn, byrow = T)
z <- b[1] + b[2]*X + b[3]*Y + b[4]*X*Y

library(rgl)
open3d()
plot3d(vol ~ meth+temp, type = "p",xlab = "metoda", ylab = "teplota", zlab = "objem")
persp3d(x,y,z, col = "blue", add = T, alpha = 0.5)
# model s interakcemi
detach(matice)

# hledani optimalniho GLM modelu
glmmod.min <- glm(VOL~1,data=toxic,family=gaussian) # stejne jako pro LM
glmmod.opt<-step(glmmod.min,scope = list(upper = ~METHOD*(TEMP+CAT), lower = ~1))
summary(glmmod.opt)
# overeni normality residui
plot(glmmod.opt,which=2)
shapiro.test(residuals(glmmod.opt))

# pro Gamma rozdeleni
glmmod.min <- glm(VOL~1,data=toxic,family=Gamma)
glmmod.opt<-step(glmmod.min,scope = list(upper = ~METHOD*(TEMP+CAT), lower = ~1))
summary(glmmod.opt)
# overeni normality residui
plot(glmmod.opt,which=2)
shapiro.test(residuals(glmmod.opt))

detach(toxic)