library(faraway)
library(lme4)

#ANOVA vs. ML estimator

#kontrasty pro neusporadane a usporadane faktory
op <- options(contrasts=c("contr.sum", "contr.poly"))
#zmeni a ulozi predchozi hodnoty
#contr.sum znamena, ze soucet parametru je 0, srov. s contr.treatment


#data - svetlost papiru v zavislosti na operatorovi
data(pulp)

#odhad ANOVA
lmod <- aov(bright ~ operator, pulp)
summary(lmod)
coef(lmod)
options(op)

#odhad rozptylu nahodneho efektu
(0.447-0.106)/5


#odhad metodou maximalni verohodnosti
#pevny intercept, nahodny intercept podle operatora
#odhad REML
mmod <- lmer(bright ~ 1+(1|operator), pulp)
summary(mmod)
# v tomto konkretnim pripade obdobne vysledky



#BLOKY JAKO NAHODNE EFEKTY

data(penicillin)
summary(penicillin)

op <- options(contrasts=c("contr.sum", "contr.poly"))

lmod <- aov(yield ~ blend + treat, penicillin)
summary(lmod)
coef(lmod)

mmod <- lmer(yield ~ treat + (1|blend), penicillin)
summary(mmod)
options(op)

#anova test - pevny efekt (ML ratio)
mod1 <- lmer(yield ~ treat + (1|blend), penicillin, REML=F)
mod2 <- lmer(yield ~ 1 + (1|blend), penicillin, REML=F)
anova(mod1,mod2)

## predikce

#random effects
ranef(mmod)$blend
#fixed effects
(cc <- model.tables(lmod))

#shrinkage
cc[[1]]$blend/ranef(mmod)$blend

#best linear unbiased predictor
fixef(mmod)[[1]]+ranef(mmod)$blend

#rezidua - srovnavani s BLUP
qqnorm(resid(mmod),main="")
plot(fitted(mmod),resid(mmod),xlab="Fitted",ylab="Residuals")
abline(0,0)


#VNORENE EFEKTY

data(eggs)
summary(eggs)
cmod <- lmer(Fat ~ 1 + (1|Lab) + (1|Lab:Technician) + (1|Lab:Technician:Sample), data=eggs)
summary(cmod)

cmodr <- lmer(Fat ~ 1 + (1|Lab) + (1|Lab:Technician), data=eggs)
anova(cmod,cmodr)
#model se da zjednodusit


#############################
# LONGITUDINALNI MODEL

#PANEL STUDY OF INCOME DYNAMICS
data(psid)
head(psid)

library(lattice)
xyplot(income ~ year | person, psid, type="l", subset=(person < 21),strip=FALSE)
xyplot(log(income+100) ~ year | sex, psid, type="l")

#obecne
lmod <- lm(log(income) ~ I(year-78), subset=(person==1), psid)
coef(lmod)

#pro jednotlive osoby
slopes <- numeric(85);intercepts <- numeric(85)
for(i in 1:85){
  lmod <- lm(log(income) ~ I(year-78), subset=(person==i), psid)
  intercepts[i] <- coef(lmod)[1]
  slopes[i] <- coef(lmod)[2]
}

plot(intercepts,slopes,xlab="Intercept",ylab="Slope")
psex <- psid$sex[match(1:85,psid$person)]
boxplot(split(slopes,psex))
t.test(slopes[psex=="M"],slopes[psex=="F"])
t.test(intercepts[psex=="M"],intercepts[psex=="F"])
#zeny vice rostou, muzi maji vyssi plat (v roce 78)

#########

psid$cyear <- psid$year-78

mmod <- lmer(log(income) ~ cyear*sex +age+educ+(cyear|person),psid)
summary(mmod)

# pro interpretaci pouzijeme normalni aproximaci t-statistiky


### nebo takto

library(nlme)

fit <- lme( log(income) ~ cyear*sex +age+educ,
	random = ~ cyear | person,
	data = psid)

summary(fit)