# ANCOVA_fitness.R
# ANCOVA for T. dubius fitness

# Load necessary libraries
library(car)
library(ggplot2)
library(multcomp)
library(emmeans)

# Read in data set
fitdata <- read.table(header=T,colClasses=c(rep("numeric",2),"factor"),text="
eggs length treat
290  5.7   DietIG
 99  5.2   DietIG
340  5.5   DietIG
271  4.8   DietIG
200  5.2   DietIG
405  5.2   DietIG
178  5.1   DietIG
 48  5.0   DietIG
146  4.8   DietIG
184  4.9   DietIG
 66  4.6   DietCPW
 93  5.0   DietCPW
  9  5.4   DietCPW
404  5.4   DietCPW
244  5.1   DietCPW
195  4.9   DietCPW
343  5.2   DietCPW
516  5.0   DietCPW
215  4.6   DietCPW
412  5.6   DietCPW
167  4.8   DietCPW
316  5.2   DietCPW
334  5.3   DietCPW
 62  4.7   WildCPW
290  5.0   WildCPW
488  5.8   WildCPW
336  5.2   WildCPW
337  5.8   WildCPW
230  5.0   WildCPW
381  5.3   WildCPW
192  5.5   WildCPW
186  5.3   WildCPW
467  5.2   WildCPW
 59  4.9   WildCPW
323  5.7   WildCPW
291  4.9   WildCPW
164  5.3   WildCPW
142  5.3   WildCPW
269  5.0   WildCPW
329  5.4   WildCPW
235  5.0   WildCPW
")

# Choose y and x variables
fitdata <- transform(fitdata,y=eggs)
fitdata <- transform(fitdata,x=length)

# Print data set
print(fitdata)

# Graphics using ggplot2
print(ggplot(fitdata,aes(x,y,group=treat,color=treat))+
geom_point()+
stat_smooth(method="lm",se=FALSE))

# ANCOVA
# Model with interaction
ancout <- lm(y~treat+x+treat:x,data=fitdata)
summary(ancout)
confint(ancout)

# ANOVA using Type II SS 
Anova(ancout)

# Model without interaction
ancout_noint <- lm(y~treat+x,data=fitdata)
summary(ancout_noint)
confint(ancout_noint)

# ANOVA using Type II SS
Anova(ancout_noint)

# Calculate least squares means
emmeans(ancout_noint,~treat)

# Obtain multiple comparison object
compout <- glht(ancout_noint,linfct=mcp(treat="Tukey"))

# Tukey method - controls the EER
summary(compout)
confint(compout)
cld(compout)

# Diagnostic plots to check ANCOVA assumptions
p <- predict(ancout_noint)
r <- resid(ancout)
plot(p,r)
qqnorm(r)
