# Tdubius_bait_color_new.R
# Two-way ANOVA for T. dubius counts (Reeve et al. 2009)

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

# Read in data set
Tdata <- read.table(header=T,colClasses=c("factor","factor","numeric"),text="
bait color Tdubius
FRT	B	18
FRT	W	12
FRT	W	15
FRT	W	7
FRT	W	4
IDT	B	2
IDT	B	1
IDT	B	4
IDT	W	2
IDT	W	1
IST	B	2
IST	B	2
IST	B	10
IST	B	7
IST	W	4
")

# Apply transformations here
Tdata <- transform(Tdata,y=log10(Tdubius+1))

# Print data
Tdata

# Graphics using ggplot2
ggplot(Tdata,aes(bait,y,group=color,color=color))+
geom_jitter(size=2,position=position_jitter(width=0.1))+
stat_summary(fun.y="mean",geom="line",size=1)+
stat_summary(fun.data="mean_cl_normal",geom="errorbar",width=0.2,size=1)

# MODEL WITH INTERACTION - USE THIS OUTPUT IF INTERACTION SIGNIFICANT
aovout <- aov(y~bait+color+bait:color,data=Tdata) 

# ANOVA using Type I SS
anova(aovout)

# ANOVA using Type II SS
Anova(aovout)

# Calculate least squares means
emmeans(aovout,~bait:color)

# Tests for simple effects
testInteractions(aovout,fixed="color",across="bait",adjustment="none")

# MODEL WITHOUT INTERACTION - USE THIS OUTPUT IF INTERACTION NS
aovout_noint <- aov(y~bait+color,data=Tdata)  # Model without interaction

# Type I SS
anova(aovout_noint)

# Type II SS
Anova(aovout_noint)

# Calculate least squares means
emmeans(aovout_noint,~bait)
emmeans(aovout_noint,~color)

# All pairwise comparisons - Tukey method
compout <- glht(aovout_noint,linfct=mcp(bait="Tukey"))
summary(compout)
confint(compout)
cld(compout)

# Diagnostic plots to check ANOVA assumptions
p <- predict(aovout)
r <- resid(aovout)
plot(p,r)
qqnorm(r)
