library(lme4)
library(ggplot2) 
library(dplyr)
#------------------------------------------------------------------------


#Example 1: Can time since calving predict protein content in cows, while accounting for individual cows?

#Get data
data("Milk")

#Test Assumptions

#Independence
#N/A

#Normality
hist(Milk$protein, col="cornsilk")

#Linearity & Homoscedasticity
plot(protein~Time, data=Milk) 
plot(lmer(protein~Time + (1|Cow), data=Milk))

#Model building
lm1_reg <- lm(protein~Time, data=Milk)
lm1_reg
summary(lm1_reg)

lm1_mix1 <- lmer(protein~Time + (1|Cow), REML=FALSE, data=Milk)
lm1_mix1
summary(lm1_mix1)

lm1_mix2 <- lmer(protein~Time + (Time|Cow), REML=FALSE, data=Milk)
lm1_mix2
summary(lm1_mix2)

#Model Comparison
AIC(lm1_reg)
AIC(lm1_mix1)
AIC(lm1_mix2) #lowest AIC

#Plots
Milk$pred2 <-predict(lm1_mix2, newdata=Milk)

ggplot(data=Milk, aes(x=Time, y=protein, col=Cow)) +
  geom_point(aes(color=Cow))+
  geom_smooth(method='lm', se=F, aes(colour=Cow), size=0.5)+
  geom_smooth(method="lm", se=T, aes(group=1), col="red", size=2)+
  xlab("Time") + ylab("protein")+
  theme(legend.position = "none")

ggplot(data=Milk, aes(x=Time, y=pred2, group=Cow)) +
  geom_point(aes(color=Cow))+
  geom_smooth(method='lm', se=F, aes(colour=Cow))+
  theme(legend.position = "None")+
  geom_smooth(method='lm', se=T, aes(group=1), col="red", size=2)+
  xlab("Time") + ylab("Predicted protein")+
  theme(legend.position = "None")
#------------------------------------------------------------------------


#Example 2: Can verbal IQ scores predict general IQ scores, while controlling for schools?

#Get data
data("bdf")

#subset for only schools with 24 students
bdf1 <-data.frame(IQ.verb=bdf$IQ.verb, IQ.perf=bdf$IQ.perf, schoolNR=factor(bdf$schoolNR))
bdf2 <- subset(bdf1, schoolNR %in% c(86,57,249,24,222,218,204,18,164,142))
head(bdf2)

#Test Assumptions

#Independence
#N/A

#Normality
hist(bdf2$IQ.perf, col="green")

#Linearity & Homoscedasticity
plot(IQ.perf~IQ.verb, data=bdf2)
plot(lmer(IQ.perf~IQ.verb +(1|schoolNR), data=bdf2))

#Model building
lm2_reg <- lm(IQ.perf~IQ.verb, data=bdf2)
lm2_reg
summary(lm2_reg)

lm2_mix1 <- lmer(IQ.perf~IQ.verb + (1|schoolNR), REML=FALSE, data=bdf2)
lm2_mix1
summary(lm2_mix1)

lm2_mix2 <- lmer(IQ.perf~IQ.verb + (IQ.verb|schoolNR), REML=FALSE, data=bdf2)
lm2_mix2
summary(lm2_mix2)

#Model Comparison
AIC(lm2_reg)  #lowest AIC, no but not much difference
AIC(lm2_mix1)
AIC(lm2_mix2)

#Plots
ggplot(data=bdf2, aes(x=IQ.verb, y=IQ.perf)) +
  geom_point(aes(color=schoolNR))+
  geom_smooth(method="lm", se=T, col="red", size=2)


bdf2$pred1 <-predict(lm2_mix1, newdata=bdf2)
bdf2$pred2 <-predict(lm2_mix2, newdata=bdf2)

ggplot(data=bdf2, aes(x=IQ.verb, y=pred1, col=schoolNR)) +
  geom_point(aes(color=schoolNR))+
  geom_smooth(method="lm", se=F, aes(colour=schoolNR))

ggplot(data=bdf2, aes(x=IQ.verb, y=pred2, col=schoolNR)) +
  geom_point(aes(color=schoolNR))+
  geom_smooth(method="lm", se=F, aes(colour=schoolNR))
#------------------------------------------------------------------------


#Example 3: Can time predict X-ray pixel intensities, while controlling for individual dogs?

#Get data
data("Pixel")

#Test Assumptions

#Independence
#N/A

#Normality
hist(Pixel$pixel, col="red")

#Linearity & Homoscedasticity
plot(pixel~day, data=Pixel) #may not be, especially at beginning
plot(lmer(pixel~day +(1|Dog), data=Pixel)) #looks like violates (certainly time frames)

#Model building
lm3_reg <- lm(pixel~day, data=Pixel)
lm3_reg
summary(lm3_reg)

lm3_mix1 <- lmer(pixel~day +(1|Dog), REML=FALSE, data=Pixel)
lm3_mix1
summary(lm3_mix1)

lm3_mix2 <- lmer(pixel~day +(day|Dog), REML=FALSE, data=Pixel)
lm3_mix2
summary(lm3_mix2)

#Model Comparison
AIC(lm3_reg)  
AIC(lm3_mix1) #lowest AIC, not significant
AIC(lm3_mix2)

#Plots
Pixel$pred <-predict(lm3_mix1, newdata=Pixel)

ggplot(data=Pixel, aes(x=day, y=pixel, col=Dog)) +
  geom_point(aes(color=Dog))+
  geom_smooth(method='lm', se=F, aes(colour=Dog), size=0.5)+
  geom_smooth(method="lm", se=T, aes(group=1), col="red", size=2)+
  xlab("Day") + ylab("Pixel intensity")

ggplot(data=Pixel, aes(x=day, y=pred, col=Dog)) +
  geom_point(aes(color=Dog))+
  geom_smooth(method='lm', se=F, aes(colour=Dog), size=0.5)+
  geom_smooth(method="lm", se=T, aes(group=1), col="red", size=2)+
  xlab("Day") + ylab("Pixel intensity")