#Two way Analysis of variance for balanced designs (equal group size)
install.packages("car")
library(car)


my_data <- read.table("twanova.txt", header=T)
attach(my_data)



#Check if all fixed factors have been correctly identified as factor
str(my_data)
my_data$predator_density <- as.factor(my_data$predator_density)
my_data$timing <- as.factor(my_data$timing)



######data visualization########

if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")
library(ggpubr)

plot <- ggline(my_data, x = "timing", y = "activity", color = "predator_density",
               add = c("mean_se", "dotplot"),
               palette = c("#00AFBB", "#E7B800"))
plot



#####Two-way ANOVA with interaction#####
#Interaction = the effect of one factor on the dependent variable depends on the level of another factor
#Type does not matter in a balanced design#

model1 <- lm(activity ~ predator_density * timing, data = my_data)
Anova(model1,type="III")



####Check model assumptions####

#Assumption 1: Normal distribution of the model residuals
plot(model1,2)

aov_residuals <- residuals(object = model1)
shapiro.test(aov_residuals)


#Assumption 2: Homogeneity of variance of the groups
plot(model1,1)
leveneTest(activity ~ predator_density * timing, data = my_data)



#Post-hoc test: In case of a significant interaction effect, we investigate all separate group combinations
install.packages("lsmeans")
library("lsmeans")

install.packages("multcomp") 
library("multcomp")

install.packages("multcompView")
library("multcompView")

posthoc <- lsmeans(model1,
           pairwise ~ predator_density*timing, 
           adjust="tukey")
posthoc

cld(posthoc[[1]],  alpha=.05,  Letters=letters)