dat<-read.csv("AD Study nesting scores.csv") #set working directory to where ever you have the .csv file saved

##### Package installation #########

install.packages("irr") #inter-rater reliability
install.packages("psych") #descriptives
install.packages("tidyr") #used for formatting wide to long data, if needed
install.packages("ez") #Mauchly's results
install.packages("lsmeans") #Bonferroni posthocs
install.packages("ggplot2") #graphs
install.packages("ggsignif") #significance bars for ggplot2

####### IRR (ICC) calculation ########

# 3, 4, 5 = square
# 6, 7, 8 = paper
# 9, 10, 11 = bedding
# 12, 13, 14 = twist

square<-dat[c(3,4,5)]
paper<-dat[c(6,7,8)]
bedding<-dat[c(9,10,11)]
twist<-dat[c(12,13,14)]

library(irr)

print(icc(square, model="twoway", type="agreement", unit="average")) #0.91
print(icc(paper, model="twoway", type="agreement", unit="average")) #0.94 # paper has the strongest agreement
print(icc(bedding, model="twoway", type="agreement", unit="average")) #0.87
print(icc(twist, model="twoway", type="agreement", unit="average")) #0.87

######## Wide to long formatting #######

dat$Animal_ID<-factor(dat$Animal_ID)
library(tidyr)
View(dat)
dat2<-dat[c(1,2,15,16,17,18)]
dat_long<-gather(dat2, material, avg_score, Avg_Square:Avg_Twist)
View(dat_long)

####### Mixed ANOVA model #######

# within subjects variable - the nesting material
# between subjects variable - genotype
# dependent variable -- average of scores

## check Mauchly's test of Sphericity - major assumption for within-designs/variables

library(ez)

ezANOVA(data=dat_long, dv=.(avg_score), wid=.(Animal_ID), within=.(material),
        between=.(Genotype), detailed=F, type=2) ## select type 2 if not expecting an interaction

## did not violate Mauchly's - woohoo! But do have unequal sample sizes. ANOVA is pretty robust and can often handle this

# set main effects of genotype and material and genotype*material interaction
# below is a litle redundant but it's the ANOVA built into CRAN; does not give you Mauchly's results

options(contrasts=c("constr.sum", "contr.poly"))
aov_nesting<-aov(avg_score ~ Genotype + material + Genotype*material + Error(Animal_ID/material), 
                data=dat_long)
design<-factor(c("Shredded paper", "Square", "Bedding", "Twist"))
summary(aov_nesting)

## to get partial eta, divide SS-effect by (SS-effect+SS-error)

29.68/(29.68+18.02) # partial eta for genotype
69.06/(69.06+22.83) # partial eta for material


## get some Bonferroni-adjusted posthocs

library(lsmeans)

lsmeans(aov_nesting, pairwise ~ material, adjust="Bonferroni")

####### Descriptives ######

library(psych)

describe.by(dat_long, dat_long$Genotype) ## average score for each genotype across all materials
describe.by(dat2, dat2$material) ## average score for each material across genotype
describe.by(dat2, dat2$Genotype)## average score for each material, by genotype


####### graphs ########

library(ggplot2)
library(ggsignif)

Genotype = c("Wildtype", "APOE e4", "Wildtype", "APOE e4", "Wildtype", "APOE e4", "Wildtype", "APOE e4")
Material = c("Shredded paper","Shredded paper","Square","Square","Bedding","Bedding","Twist", "Twist")
Mean = c(4.77,3.52, 2.57, 1.39, 2.73, 1.73, 2.63, 1.30)
se = c(0.16, 0.24, 0.36, 0.06, 0.25, 0.26, 0.29, 0.05) 
df = data.frame(Genotype, Material, Mean, se)
df$Genotype<- factor(df$Genotype, levels=(c("Wildtype", "APOE e4")))
df$Material<- factor(df$Material, levels=(c("Shredded paper","Square", "Bedding", "Twist")))
str(df)

nestgraph<-ggplot(df, aes(x=Material, y = Mean, fill=Genotype))+
  geom_bar(position=position_dodge(), stat="identity")+
  scale_fill_manual(values=c("#CCCCCC", "#666666", "#CCCCCC", "#666666", "#CCCCCC", "#666666", "#CCCCCC", "#666666"))+
  geom_errorbar(aes(ymin=Mean-se, ymax=Mean+se), 
                width=0.2,
                position=position_dodge(0.9))+
  geom_segment(aes(x = 2, y = 3.98, xend = 4, yend = 3.98))+  #this line of code makes the bar that extends across the 3-nonsig materials
  scale_y_continuous(name="Average score\n", expand=c(0,0))+
  geom_text(aes(x=1, y=5.9, label="Stretch It"), vjust=-1)+
  scale_x_discrete(name="\nProvided materials")+
  ggtitle("          Averaged scores of nests constructed by wildtype or APOE e4 mice\n")+ #added chunk of space to center title in jpeg file
  theme(plot.title=element_text(lineheight=0.8, face="bold", hjust=0.5))+
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        panel.background=element_blank(), axis.line=element_line(color="black"))+

  geom_signif(data=df,
              aes(xmin=0.75, xmax=1.25, annotations="*", y_position=5.25),
              textsize = 5, vjust = 0.05, tip_length = c(0.04, 0.2),
              manual=TRUE) +
  geom_signif(data=df,
               aes(xmin=1.75, xmax=2.25, annotations="*", y_position=3.30),
               textsize = 5, vjust = 0.05, tip_length = c(0.04, 0.3),
               manual=TRUE) +
  geom_signif(data=df,
              aes(xmin=2.75, xmax=3.25, annotations="*", y_position=3.30),
              textsize = 5, vjust = 0.05, tip_length = c(0.04, 0.2),
              manual=TRUE)+
  geom_signif(data=df,
              aes(xmin=3.75, xmax=4.25, annotations="*", y_position=3.30),
              textsize = 5, vjust = 0.05, tip_length = c(0.04, 0.3),
              manual=TRUE)+
  geom_signif(data=df,
              aes(xmin=1, xmax=3, annotations="**", y_position=5.7),
              textsize = 5, vjust = -0.0000025, tip_length = c(0.04, 0.35),
              manual=TRUE)
 
print(nestgraph)

ggsave(filename="nesting3.jpeg", plot=nestgraph, width=7.5, height=6.75) #this will save the graph to where you have your working directory set