#####################
# This R script can be used to evaluate the number of collected petals.
# You only need to change the working directory to the path you extracted
# the data to. Further, three additional packages are required:
# ez (performs the actual ANOVAs)
# schoRsch (allows convenient formatting of the ANOVA tables)
# plyr (a toolset for renaming data frame structures)
# You can install these packages via the shell (internet connection required):
# install.packages("ez", dependencies = TRUE)
# install.packages("schoRsch", dependencies = TRUE)
# install.packages("plyr", dependencies = TRUE)
#
# The offsets are aggregated to obtain a binary factor (visual offsets to the
# left or the rigth), you can just outcomment the respective lines to work
# with the four-level offset factor.
#####################

library(ez)
library(schoRsch)
library(plyr)

# enter the working directory here
setwd("")

# obtain the data...
data <- read.table("LiS_Petal_Data.dat", header = TRUE)

# these lines rename two factor levels of the offset factor to obtain a binary offset factor,
# just outcomment or remove them to work with the four-level factor
data$OffsetType <- revalue(data$OffsetType, c("PositiveXPositiveZ"="PositiveX"));
data$OffsetType <- revalue(data$OffsetType, c("NegativeXPositiveZ"="NegativeX"));

# aggregate data and remove empty cells
data <- aggregate(Petals ~ ParticipantID*Block*OffsetType, data, mean)
data <- droplevels(data);

summary(data)

data$ParticipantID <- factor(data$ParticipantID)

#### some descriptive statistics

st.err <- function(x) {
  sd(x)/sqrt(length(x))
}
aggregate(Petals ~ Block, data, mean)
aggregate(Petals ~ Block, data, st.err)
aggregate(Petals ~ Block, data, sd)

bp = boxplot(Petals ~ Block, col=c("white", "lightgray"), data)
t(bp$stats)

aggregate(Petals ~ OffsetType, data, mean)
aggregate(Petals ~ OffsetType, data, st.err)
aggregate(Petals ~ OffsetType, data, sd)

bp = boxplot(Petals ~ OffsetType, col=c("white", "lightgray"), data)
t(bp$stats)

aggregate(Petals ~ ParticipantID, data, mean)
aggregate(Petals ~ ParticipantID, data, st.err)
aggregate(Petals ~ ParticipantID, data, sd)

#### the actual ANOVA analysis
anova_rt <- ezANOVA(data=data, dv = .(Petals), wid = .(ParticipantID), within = .(Block,OffsetType), detailed = TRUE)
anova_out(anova_rt)

#### post-hoc tests, not necessary for the two-level offset factor
tests = pairwise.t.test(data$Petals, data$OffsetType, p.adj = "holm", pool.sd = FALSE, paired = TRUE)
tests$p.value