library(tidyverse)
library(MASS) 
library(car) 
library(agricolae)

read_csv("datasets/reduced.csv") -> reduced

##### P1 Developing hypotheses #####

#Create your 6 null and alternate hypotheses
HO1 <- "sample null hypothesis 1,"
HA1 <- "sample alternate hypothesis 1,"

HO2 <- "sample null hypothesis 2."
HA2 <- "sample alternate hypothesis 2."


print(paste("Our null hypotheses are that", HO1, "as well as", HO2, "Our alternate hypotheses are that", HA1, "as well as", HA2))


##### P2 Exploring the Data #####

#Summarize the data


#Create and save your boxplots
ggplot(reduced) +
  geom_boxplot() +
  aes(x = predictor, y = response) -> figure_1

ggsave("figure_1.png", figure_1, height = 4, width = 7)


#Make and save your interaction plot
png("interaction_plot.png", width = 7, height = 4, units = "in", res = 100)
interaction.plot(x.factor = reduced$Treatment, 
                 trace.factor = reduced$Gradient, 
                 response = reduced$EarlyGrowth, 
                 fun = mean, type = "b", leg.bty = "o", fixed = TRUE,
                 trace.label = "Legend label", xlab = "axis label", 
                 ylab = "axis label")
dev.off()




##### P3 Conducting Normality Tests #####

#Perform the Shapiro-Wilk test
reduced %>%
  group_by(Treatment, Gradient) %>%
  summarise_all(.funs = funs(statistic = 
                               shapiro.test(.)$statistic, 
                             p.value = shapiro.test(.)$p.value)) -> SWtest

print(SWtest)

#Create a model and make a QQ plot
lm(response ~ Predictor_1 * Predictor_2, data = reduced) -> model

qqnorm(model$residuals); qqline(model$residuals)


##### P4 Box-Cox Transformation #####

#Perform the Box-Cox Transformation on the model and find the best lambda
boxcox(model) -> bc

lmda <- bc$x[which.max(bc$y)] 

#Create a new model
lm(((EarlyGrowth^lmda-1)/lmda) ~ Predictor_1 * Predictor_2, data = reduced) -> normal_model

#Make a new dataset
reduced %>%
  mutate(BC_Growth = ((EarlyGrowth^lmda-1)/lmda)) -> new

#Perform Part 3 again with the transformed data



##### P5 Testing for Homogeneity #####

#Conduct the Levene Test
leveneTest(response ~ Predictor_1 * Predictor_2, data = new)

#Extract the P-value from the test and answer whether it is significant or not



##### P6 Conducting the 2-Way ANOVA #####

#Conduct the Two-Way ANOVA
lm(response ~ Predictor_1 * Predictor_2, data = new) -> new_model

#print the results of this ANOVA as well as a summary of it
aov(new_model)


##### P7 Interpreting the ANOVA #####

#Should you reject or fail to reject your hypotheses?



##### P8 Conducting a 1-Way ANOVA #####

#Conduct a One-Way ANOVA
lm(BC_Growth ~ Gradient, data = new) -> gradient_model

#Print the results and a summary of your ANOVA



##### P9 Removing data to make it work #####

#Conduct a Tukey HSD Test
HSD.test(gradient_model, "Gradient", alpha = 0.05, group = TRUE, main = "Please help", unbalanced=TRUE, console=TRUE)

#Choose what data you will be removing from the dataset

#Perform the necessary steps to explore the effectiveness of your choice