# TrapRCBD_clerids.R - trapping experiment data
# Randomized block mixed model and Tukey method

# Load necessary libraries
library(ggplot2)
library(lme4)
library(lmerTest)
library(multcomp)

# Read in data set
trapdata <- read.table(header=T,colClasses=c("factor","factor","numeric"),text="
block treat count
1	AP	    4
1	BLANK   0
1	FRAP	 79
1	IDAP	  7
1	ISAP	 10
2	AP	    1
2	BLANK	  0
2	FRAP  124
2	IDAP	 13
2	ISAP	 20
3	AP	    0
3	BLANK	  0
3	FRAP	 14
3	IDAP	 NA
3	ISAP	  2
4	AP	    0
4	BLANK	  0
4	FRAP   15
4	IDAP	 11
4	ISAP	  7
5	AP	    0
5	BLANK	  0
5	FRAP	 29
5	IDAP	  7
5	ISAP	  7
6	AP	    2
6	BLANK	  0
6	FRAP	 70
6	IDAP	 14
6	ISAP	 20
")

# Apply transformations here
trapdata <- transform(trapdata,y=log(count+1))

# Print data
trapdata

# Plot treatment means and confidence intervals
ggplot(trapdata,aes(treat,y))+
geom_jitter(size=2,position=position_jitter(width=0.1))+
stat_summary(fun.y="mean",geom="point",size=3)+
stat_summary(fun.data="mean_cl_normal",geom="errorbar",width=0.2,size=1)

# Plot treatments and blocks
# Remove missing values using is.na
ggplot(trapdata[!is.na(trapdata$y),],aes(treat,y,group=block,color=block))+
geom_line(size=1)

# Delete blank traps before analysis
trapdata2 <- subset(trapdata,subset=treat!="BLANK")

# Mixed model
lmerout <- lmer(y~treat+(1|block),data=trapdata2)
summary(lmerout)
anova(lmerout)

# Calculate least squares means
lsmeansLT(lmerout)

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

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