################################################################################
# Exercise 2
# Additive Mixed Model 
################################################################################

library(mgcv)
library(lattice)
library(nlme)

###### (2a)
?Soybean
str(Soybean)
plot(Soybean)
plot(Soybean, outer= ~ Year*Variety)
#plot(Soybean, outer= ~ Variety*Year)

###### (2b)
Soybean$logweight <- log(Soybean$weight)
xyplot(logweight ~ Time | Variety, Soybean, groups = Plot, type = "l") #Spaghetti-Plot per Variety

# xyplot(logweight ~ Time | Variety + Year, Soybean, groups = Plot, type = "l") #Spaghetti-Plot per Variety and Year
# xyplot(logweight ~ Time | Plot, Soybean, groups = Variety, type="l") #Spaghetti-Plot per Plot


###### (2c)
# # LMM with random intercepts per Plot
# # the same model with two different specifications of the random effect
# l1 <- lme(logweight ~ Year + Variety +  Time,
#           random = ~1|Plot, data = Soybean, method = "REML")
# l1 <- lme(logweight ~ Year + Variety +  Time,
#           random = list(Plot = ~ 1), data = Soybean, method = "REML")


# modell with Year, Variety,  smooth function over time (the same form for all plots)
# and random intercepts for each Plot
# argument method: the smoothing parameter estimation method
# use gamm() instead of gam() for random intercept
g1 <- gamm(logweight ~ Year + Variety +  s(Time, bs = "ps", m = c(2,1), k = 20),
           random = list(Plot = ~ 1), data = Soybean, method = "REML") 
g1
plot(g1$gam)

summary(g1$lme)
summary(g1$gam)


######## ALTERNATIVE for the estimation of the random intercepts
g1gam <- gam(logweight ~ Year + Variety +  s(Time, bs = "ps", m = c(2,1), k = 20) + 
               s(Plot, bs = "re"), data = Soybean, method = "REML") 
summary(g1gam)
## get variance components of random efffects
gam.vcomp(g1gam)
## predict each component of the linear predictor seperately
# column s(Plot) contains predicitons of random effects
pred <- predict(g1gam, type="terms")
# pred[1:20,] ## look at the first 20 predictions


###### (2d)
# modell with Year, Variety, smooth function over time (the same form for all plots)
# and random intercepts and slopes for each Plot
g2 <- gamm(logweight ~ Year + Variety +  s(Time, bs = "ps", m = c(2,1), k = 20),
           random = list(Plot = ~ 1 + Time), data = Soybean, method = "REML")
g2
plot(g2$gam)


###### (2e)
# modell with Year, smooth function over time per Variety (the same form for all plots)
# and random intercepts for each Plot
g3 <- gamm(logweight ~ Year +  s(Time, by = Variety, bs = "ps", m = c(2,1), k = 20),
           random = list(Plot = ~ 1), data = Soybean, method = "REML")
g3
par(mfrow=c(1,2))
plot(g3$gam)