#########################################################
# Reading in the data
#########################################################
#Reading in the data containing the experimental units
Explanatory <- read.csv( file.choose() )

#Reading in the data containing the possible responses
Response <- read.csv( file.choose() )

#########################################################
# Doing a single iteration of a simulation
#########################################################
#Get a sequence for the treatments
Treatment <- c("A","A","A","A","A","A","A","A","A","A","A","A","B","B","B","B","B","B","B","B","B","B","B","B")

#Getting one random assignment
RandomAssignments <- sample(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),replace=FALSE)

ExperimentData <- data.frame('PersonID' = RandomAssignments, 'Treatment'= Treatment)

#Getting the dplyr library to join the ExperimentData with the Responses
library(dplyr)

#Doing the join
ExperimentData <- left_join(ExperimentData,Response,
                            by = c('PersonID' = 'PersonID', 
                                   'Treatment' = 'Treatment')
                            
                            
                            )

#Get ANOVA analysis and get the ANOVA table
ANOVA_Output <- aov(Improvement ~ Treatment, data=ExperimentData)
summary(ANOVA_Output)


#########################################################
# Creating a function for the simulation study
#########################################################

MySimulation <- function(){
  
  #Setting up the treatment
  Treatment <- c("A","A","A","A","A","A","A","A","A","A","A","A","B","B","B","B","B","B","B","B","B","B","B","B")
  
  #Getting the random assignment
  RandomAssignments <- sample(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),replace=FALSE)
  
  #Creating the ExperimentData
  ExperimentData <- data.frame('PersonID' = RandomAssignments, 'Treatment'= Treatment)
  
  #Create a jOin to obtain the responses
  ExperimentData <- left_join(ExperimentData,Response,
                              by = c('PersonID' = 'PersonID', 
                                     'Treatment' = 'Treatment')
                              )
  ANOVA_Output <- aov(Improvement ~ Treatment, data=ExperimentData)
  SSE_Value <- summary(ANOVA_Output)[[1]][2,2]
  
  return(SSE_Value)  
}


#########################################################
# Doing the simulation
#########################################################

#Running 100 iteration of my simulation
SimulationOutcomes <- data.frame('SSE_Values' = replicate(100,MySimulation() ) )


#Doing some plotting and getting some summaries
library(ggplot2)

#Creating a dotplot of the SSE_Values
ggplot(SimulationOutcomes, aes(SSE_Values)) + geom_dotplot()

#Getting the mean and standard deviation of the simulated values
summarize(SimulationOutcomes, 'Mean' =  mean(SSE_Values))
summarize(SimulationOutcomes, 'Standard Deviation' =  sd(SSE_Values))

########################################################
# End of code
########################################################