#########################################################
# 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() )


#########################################################
# Preparing the Explanatory Data for blocking
#########################################################
#Getting the dplyr library to join the ExperimentData with the Responses
library(dplyr)

#Sort the data.frame by InitialCholesterol
Explanatory_ByCholesterol <- arrange(Explanatory, InitialCholesterol)

#Adding a variable to data.frame to identify the block
Explanatory_ByCholesterol <- mutate(Explanatory_ByCholesterol, 'Block' = rep(1:12,each=2))



#########################################################
# Doing a single iteration of a simulation
#########################################################
#Make Treatment Assignments (within each Block)
Treatment <- rep(c("A","B"),times=12)

#Getting one random assignment
RandomAssignments <- Explanatory_ByCholesterol %>% group_by(Block) %>% sample_n(2,replace=F) %>% select(PersonID,Block)


ExperimentData <- data.frame(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 with blocking and get the ANOVA table
ANOVA_Output <- aov(Improvement ~ Treatment, data=ExperimentData)
summary(ANOVA_Output)


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

MyBlockSimulation <- function(){
  
  #Make Treatment Assignments (within each Block)
  Treatment <- rep(c("A","B"),times=12)
  
  #Getting one random assignment
  RandomAssignments <- Explanatory_ByCholesterol %>% group_by(Block) %>% sample_n(2,replace=F) %>% select(PersonID,Block)
  
  
  ExperimentData <- data.frame(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 + Block, data=ExperimentData)
  SSE_Value <- summary(ANOVA_Output)[[1]][3,2]
  
  return(SSE_Value)  
}


#########################################################
# Doing the simulation and plotting outcomes
#########################################################

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


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

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

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


#########################################################
# Comparing the CRD Design to the Randomzied Block Design
#########################################################
SimulationOutcomes<-mutate(SimulationOutcomes, 'Type'= "CRD")
BlockSimulationOutcomes<-mutate(BlockSimulationOutcomes, 'Type' = "Block")
AllSimulationOutcomes <- bind_rows(SimulationOutcomes, BlockSimulationOutcomes)

#ggplot to compare
ggplot(AllSimulationOutcomes, aes(x=SSE_Values)) +
  geom_dotplot(binwidth=100) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank()) +
  scale_x_continuous(limits=c(0,8000)) +
  facet_wrap(~ Type, ncol=1)


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