#Lineární model s metodou nejmenších čtverců
x<-c(0:20)
y<-0+2*x+rnorm(21,0,3)

par(mar=c(3,2,1,1))
plot(x,y,pch=19,col="blue",lwd=3)

model.lm1<-lm(y~x)
lines(x,model.lm1$coefficients[1]+model.lm1$coefficients[2]*x,
      col="lightblue",lwd=3)

# Otestování předpokladů Gauss-Markovovy věty
hist(model.lm1$residuals)
mean(model.lm1$residuals)
cor.test(model.lm1$residuals,y)

# Nyní vygenerujeme data, aby předpoklady určitě neplatily
x<-c(0:20)
y<-0+2*x+rlnorm(21,0,2)

plot(x,y,pch=19,col="blue")

model.lm2<-lm(y~x)
lines(x,model.lm2$coefficients[1]+model.lm2$coefficients[2]*x,
      col="navyblue",lwd=3)

hist(model.lm2$residuals)
mean(model.lm2$residuals)
cor.test(model.lm2$residuals,y)

# Zkusíme zlogaritmovat - obezlička, kterou jde občas obejít problém, ale dopouštíme se nepřesnosti
ylog<-log10(y)

plot(x,ylog,pch=19,col="blue")
model.log<-lm(ylog~x)
lines(x,model.log$coefficients[1]+model.log$coefficients[2]*x,
      col="cornflowerblue",lwd=3)
hist(model.log$residuals)
mean(model.log$residuals)
cor.test(model.log$residuals,y)


# lepší je použít nepramaterickou metodu Theil-Senův estimátor (je třeba do mblm() zadat argument repeated=FALSE)
install.packages("mblm")
library("mblm")

model.mblm<-mblm(y~x,repeated=FALSE)
lines(x,model.mblm$coefficients[1]+model.mblm$coefficients[2]*x,
      col="black",lwd=3)


# GLM - exponenciální quasipoissonovská regrese - linkovací funkce je log() (použijeme proto inverzní: exp())
# často využívaná pro proměnné ohraničené z jedné strany (koncentrace, hlučnost, počty jedinců nebo jednotek)
x<-c(0:20)
y<-exp(2+0.4*x+rnorm(21,0,1)) # intercept 2, sklon 0.4

plot(x,y,pch=19,col="blue")

modelglm1<-glm(y~x,family=poisson(link="log"))
lines(x,exp(modelglm1$coefficients[1]+modelglm1$coefficients[2]*x),
      col="red",lwd=3)


# GLM - logistická quasibinomiální regrese - linkovací funkce je logit() (použijeme proto inverzní: expit() neboli exp(x)/(1+exp(x)))
# často využívaná pro proměnné ohraničené z obou stran (podíly, procentuální hodnoty, pH, pravděpodobnost)
x<-c(0:20)
y<-1-1/(1+exp(-5+0.5*x+rnorm(21,0,1))) # intercept -5, sklon 0.5

plot(x,y,pch=19,col="blue")

modelglm2<-glm(y~x,family=quasibinomial(link="logit"))
lines(exp(modelglm2$coefficients[1]+modelglm2$coefficients[2]*x)/(1+exp(modelglm2$coefficients[1]+modelglm2$coefficients[2]*x)),
      col="forestgreen",lwd=3)