# Chocolate Data - 18-3-ANOVA

### Function Definitions

# updateToFactor - Converts the numeric column to a factor
# Data frame column dfc and list of replacement values v passed to function
# Returns the factor form of the data column

updateToFactor<-function(dfc,v)
{
  # dfc is currently a numeric vector
  # We want to generate a character vector, so create a new NULL vector
  # that can have the appropriate elements added to it
  
  newValues<-NULL
  
  # Loop through each element of dfc and add the appropriate value
  # from v to the end of newValues to get a vector of characters that is
  # the same length as the original dfc column
  
  for(i in dfc)
  {
    newValues<-c(newValues,v[i])
  }

  # Convert newValues to a factor
  
  newValues<-factor(newValues)
  
  # Return newValues
  
  return(newValues)
}

# updateToFactor - Converts the numeric column to a factor
# Data frame column dfc and list of replacement values v passed to function
# Returns the factor form of the data column

updateCodeToFactor<-function(dfc,v)
{
  # dfc is currently a numeric vector
  # We want to generate a character vector, so create a new NULL vector
  # that can have the appropriate elements added to it
  
  newValues<-NULL
  
  # Loop through each element of dfc and add the appropriate value
  # from v to the end of newValues to get a vector of characters that is
  # the same length as the original dfc column
  
  # This is not the most elegant way of doing this, but it's the easiest
  # The original Excel file gives an incorrect set of codes - 
  # HA should be HCA and P should be PT
  
  for(i in dfc)
  {
    if (i=="A") newValues<-c(newValues,v[1])
    if (i=="D") newValues<-c(newValues,v[2])
    if (i=="HCA") newValues<-c(newValues,v[3])
    if (i=="N") newValues<-c(newValues,v[4])
    if (i=="OT") newValues<-c(newValues,v[5])
    if (i=="PT") newValues<-c(newValues,v[6])
    if (i=="O") newValues<-c(newValues,v[7])
    if (i=="Leftover") newValues<-c(newValues,v[8])
  }
  
  # Convert newValues to a factor
  
  newValues<-factor(newValues)
  
  # Return newValues
  
  return(newValues)
}

# timeToString - converts number of minutes

# generateSummary - Counts how many of each chocolate were eaten where
# Creates a data frame to record the sweets eaten in each location
# Data frame df passed to function
# Returns the data frame

# This needs a data frame that contains 8 rows and 3 columns
# Column 1 - each hospital name twice
# Column 2 - each chocolate brand for each hospital
# Column 3 - the number of chocolates of that brand eaten at that hospital

generateSummary<-function(df)
{

  # Create a vector hospitals containing the hospital names
  
  hospitals<-c("Bedford AAU", "Bedford AAU",
               "Bedford Pilgrim","Bedford Pilgrim",
               "Swindon Meldon","Swindon Meldon",
               "UHNS Haem","UHNS Haem")
  
  # Create a vector chocolates containing the chocolate brands
  
  chocolates<-c("Roses","Quality Street","Roses","Quality Street",
                "Roses","Quality Street","Roses","Quality Street")
  
  # Create a vector numberEaten containing zeros (to be updated)
  
  numberEaten<-c(0,0,0,0,0,0,0,0)
  
  # Create a vector boxOpened
  
  boxOpened<-c(NA,NA,NA,NA,NA,NA,NA,NA)
  
  # Create a vector finalSweet
  
  finalSweet<-c(NA,NA,NA,NA,NA,NA,NA,NA)
  
  
  # Bind the five vectors into a data frame called summary
  summary<-data.frame(Hospital=hospitals,
                      Chocolates=chocolates,
                      Eaten=numberEaten,
                      BoxOpened=boxOpened,
                      FinalSweet=finalSweet)
  
  # Loop through each record in the df data frame
  
  for (i in seq(1:length(df$Location)))
  {
    # For each row of the summary table...
    
    for(j in 1:8)
    {
      # ... check to see which row of the summary table it refers to
      if (df[i,1]==summary[j,1] && df[i,2]==summary[j,2])
      {
        # If the chocolate has been eaten, increment the chocolates eaten counter
        # then transfer the time eaten to the summary file (this will be overwritten
        # by chocolates eaten after it).
        
        if (df[i,8]=="Eaten") 
        {
          summary[j,3]<-summary[j,3]+1
          summary[j,5]<-df[i,7]
        }
        
        # If this is the first chocolate of the box, copy box opened time to the summary
        if (df[i,3]==1) summary[j,4]<-df[i,5]
      }
    }
  }
  
  # Return the completed summary data frame to the main script
  
  return(summary)
}

# convertToMinutes - Takes a "HH:MM" string and converts
# it to an integer number of minutes (O0:00 is 0 minutes)
# Parameters passed: data frame row x, integer column
# Returns: an integer representing minutes since midnight 

convertToMinutes<-function(x, column)
{
  
  # Extract the string to be converted
  convertable<-as.vector(x[column])

  # Some strings are 4 characters, some are 5
  # The number of characters to the left of the 
  # colon is 1 if the number is 4, 2 if it is 5
  # Use the function nchar(string) to get the 
  # the string length
  
  stringlength<-nchar(convertable)
  
  # If stringlength is 4
  # hours is character 1 * 60
  # minutes is characters 3,4
  # minuteTime is hours+minutes
  # Use the function substr(x,start,stop)
  # to get the substrings, then as.numeric()
  
  if(stringlength==4)
  {
    hours<-as.numeric(substr(convertable,1,1))*60
    minutes<-as.numeric(substr(convertable,3,4))
    minuteTime<-hours+minutes
  }

  # If stringlength is 5
  # hours is characters 1,2 * 60
  # minutes is characters 4,5
  # minuteTime is hours+minutes
  # Use the function substr(x,start,stop)
  # to get the substrings, then as.numeric()
  
  if(stringlength==5)
  {
    hours<-as.numeric(substr(convertable,1,2))*60
    minutes<-as.numeric(substr(convertable,4,5))
    minuteTime<-hours+minutes
  }
  
  # Return the value of minuteTime
  
  return(minuteTime)
}

### calcMode - calculates the modal value of an input vector
# Parameters passed: vector v
# Returns: a vector containing the modal value or values

calcMode<-function(v)
{
  # Get all unique values within v
  uniqVals<-unique(v)
  
  # Identify the index of the most common value
  modalVals<-which.max(uniqVals)
  
  # Return the value of v that occurs at index modelVals
  return(v[modalVals])
}

### generateDescStats - generates the descriptive statistics
# Parameters passed: vector v
# Returns: a vector containing the descriptive statistics values

generateDescStats<-function(v)
{
  # Create vector to be returned at the end
  
  vDescStats<-NULL
  
  # Calculating basic descriptive statistics 
  
  vDescStats<-c(vDescStats, mean(v, na.rm=TRUE))
  vDescStats<-c(vDescStats, median(v, na.rm=TRUE))
  
  # Calculating the modal value
  # This requires the function calcMode defined above
  
  vDescStats<-c(vDescStats, calcMode(v))
  
  # Calculating more basic descriptive statistics   
  vDescStats<-c(vDescStats, max(v, na.rm=TRUE))
  vDescStats<-c(vDescStats, min(v, na.rm=TRUE))
  vDescStats<-c(vDescStats, length(v))
  
  # Calculate range
  vRange<-range(v, na.rm=TRUE)
  vDescStats<-c(vDescStats, vRange[2]-vRange[1])

  # Percentiles and quartiles

  # Calculate the 10th percentiles
  vDescStats<-c(vDescStats, as.numeric(quantile(v, c(0.1), na.rm=TRUE)))
  
  # Calculate the quartiles for v
  qChoc<-quantile(v, na.rm=TRUE)
  
  # Extract the 1st [2] and 3rd [4] quartiles
  
  vDescStats<-c(vDescStats, as.numeric(qChoc[2]))
  vDescStats<-c(vDescStats, as.numeric(qChoc[4]))
  
  # Calculate the 90th percentiles
  vDescStats<-c(vDescStats, as.numeric(quantile(v, c(0.9), na.rm=TRUE)))
  
  # Calculate the interquartile range as a number
  vDescStats<-c(vDescStats, as.numeric(qChoc[4]-qChoc[2]))

  # Displaying skewness and kurtosis
  # Requires the e1071 package
  
  vDescStats<-c(vDescStats, skewness(v, na.rm=TRUE))
  vDescStats<-c(vDescStats, kurtosis(v, na.rm=TRUE))
  
  # Return vDescStats
  return(vDescStats)
}

### addStat - generates mean, sd and median stats
# Parameters passed: chocolateData data frame df, integer statFun to define
# which statistic is to be generated, summary data frame summaryStats
# Returns: a vector containing the descriptive statistics values

addStat<-function(df,statFun, summaryStats)
{
  # Create temporary vector to be returned
  sumStat<-NULL

  # Create aggStat
  # This uses the switch() function to select which function should be 
  # called by the FUN= term of the aggregate() function
  aggStat<-aggregate(df$SurvivalTime, 
                      by=list(Location=df$Location, 
                              Chocolate=df$Chocolate), 
                      FUN=switch(statFun,mean,sd,median), 
                      na.rm=TRUE)

  
  # For each row in aggStat
  for (i in 1:8)
  {
    # For each row of the summary table...
    
    for(j in 1:8)
    {
      # ... check to see which row of the summary table it refers to
      if (aggStat[j,1]==summaryStats[i,1] && aggStat[j,2]==summaryStats[i,2])
      {
        # Add current value to the end of sumStat
        sumStat<-c(sumStat,aggStat[j,3])

      }
    }
  }

  
  # Return aggStat
  return(sumStat)
}





### Main Script

# Load data file into data frame ChocolateData
# We want to work with categorical data, so use stringsAsFactors=TRUE

chocolateData<-read.csv("SP18/ChocolateDataCSV.csv", header=TRUE, stringsAsFactors = FALSE)

# Update WhoAte to become a factor

# Define the values to replace the numbers in the column to be factorised

newValues<-c("Administrative","Doctor","Healthcare Assistant","Nurse",
             "Occupational Therapist","Physiotherapist","Other","Not Eaten")

# Call updateCodeToFactor with WhoAte and newValues

chocolateData$WhoAte<-updateCodeToFactor(chocolateData$WhoAte,newValues)

# Update Location to become a factor

# Define the values to replace the numbers in the column to be factorised

newValues<-c("Bedford AAU","Bedford Pilgrim","Swindon Meldon","UHNS Haem")

# Call updateToFactor with Location and newValues

chocolateData$Location<-updateToFactor(chocolateData$Location,newValues)

# Define the values to replace the numbers in the column to be factorised

newValues<-c("Roses","Quality Street")

# The $Chocolate column starts at 0, not 1, so needs to be updated
# Add 1 to every element in chocolateData$Chocolate

chocolateData$Chocolate<-chocolateData$Chocolate+1

# Call updateToFactor with $Chocolate and newValues

chocolateData$Chocolate<-updateToFactor(chocolateData$Chocolate,newValues)

# Define the values to replace the numbers in the column to be factorised

newValues<-c("Leftover","Eaten")

# The $Eaten column starts at 0, not 1, so needs to be updated
# Add 1 to every element in chocolateData$Eaten

chocolateData$Eaten<-chocolateData$Eaten+1

# Call updateToFactor with $Eaten and newValues

chocolateData$Eaten<-updateToFactor(chocolateData$Eaten,newValues)

# Generate summary

summary<-generateSummary(chocolateData)

# Display the summary data

summary

# Create a vector boxPlaced by applying convertToMinutes
# to the data frame chocolateData
# As convertToMinutes expects more than one parameter,
# we have to prefix if with function(x) in the apply()

boxPlaced<-apply(chocolateData,1,function(x) convertToMinutes(x,5))

# Replace the chocolateData$boxPlace with the new
# minutes-only data now stored in the boxPlaced vector

chocolateData$BoxPlace<-boxPlaced

# Create a vector timeOpened by applying convertToMinutes
# to the data frame chocolateData
# As convertToMinutes expects more than one parameter,
# we have to prefix if with function(x) in the apply()

timeOpened<-apply(chocolateData,1,function(x) convertToMinutes(x,6))

# Replace the chocolateData$TimeOpen with the new
# minutes-only data now stored in the timeOpened vector

chocolateData$TimeOpen<-timeOpened

# Create a vector timeEaten by applying convertToMinutes
# to the data frame chocolateData
# As convertToMinutes expects more than one parameter,
# we have to prefix if with function(x) in the apply()

timeEaten<-apply(chocolateData,1,function(x) convertToMinutes(x,7))

# Replace the chocolateData$boxPlace with the new
# minutes-only data now stored in the timeEaten vector

chocolateData$TimeEat<-timeEaten

# Survival time = $TimeEat-$BoxPlace

survival<-chocolateData$TimeEat-chocolateData$BoxPlace

# Bind survival time column to data frame

chocolateData<-cbind(chocolateData,SurvivalTime=survival)

# Update Leftover times
# Loop through each row of chocolateData and if Eaten="Leftover" then set SurvivalTime=NA

for (i in seq(length(chocolateData$Location)))
{
  if (chocolateData[i,8]=="Leftover") chocolateData[i,9]<-NA
}




# Create labels for descStats data frame

descStatsLabels<-c("Mean", "Median", "Mode", "Max", "Min", "Count", "Range", 
                   "10th percentile", "1st quartile", "3rd quartile", "90th percentile", 
                   "Interquartile range", "Skewness", "Kurtosis")

# Calculate descriptive statistics for $SurvivalTime

descStats<-generateDescStats(chocolateData$SurvivalTime)

# Create data frame from labels and $SurvivalTime data

descStats<-data.frame(Labels=descStatsLabels,SurvivalTime=descStatsST)
descStats






# Aggregate Data by chocolate type and who ate the chocolate

whoAteWhat<-aggregate(chocolateData$SurvivalTime, 
                      by=list(Chocolate=chocolateData$Chocolate,
                              Who=chocolateData$WhoAte), 
                      FUN=mean, 
                      na.rm=TRUE)
whoAteWhat

# Add columns to the summary table for mean, sd, median

sumMeans<-addStat(chocolateData,1,summary)
sumSD<-addStat(chocolateData,2,summary)
sumMedian<-addStat(chocolateData,3,summary)

summary<-cbind(summary,Means=sumMeans,SD=sumSD,Medians=sumMedian)








### Statistical Analyses

# t-test: Survival time for different types of chocolates

# Extract chocolateData$SurvivalTime for each type of chocolate

Roses<-chocolateData$SurvivalTime[chocolateData$Chocolate=="Roses"]
QStreet<-chocolateData$SurvivalTime[chocolateData$Chocolate=="Quality Street"]

# Is each data sample normal? Use Shapiro test to find out

shapiro.test(Roses)
shapiro.test(QStreet)

# Both of these show as p<0.05, so they are not normally distributed
# Therefore create boxplots to see where the outliers are 

boxplot(Roses, main="Roses")
boxplot(QStreet, main="QStreet")

# Roses outliers are above 200, so filter the data sample to remove them
# NB: This may take several iterations to remove all outliers

Roses<-Roses[Roses<200]
boxplot(Roses, main="Roses")

# QStreet outliers are above 200, so filter the data sample to remove them
# NB: This may take several iterations to remove all outliers

QStreet<-QStreet[QStreet<200]
boxplot(QStreet, main="QStreet")

# Perform t-test for independent samples

# Now normal, so can do the t-test

# HO: There is no difference between the means, i.e. both chocolates have the same survival time
# H1: There is a difference between the means, i.e. one type lasts longer than the other

t.test(Roses, QStreet)

# The calculated p-value for this is 0.7967. This is more than 0.05,
# so we cannot reject the null hypothesis at the p=0.05 level.
# We therefore accept the null hypothesis that the type of chocolate 
# makes no difference to the mean survival time.



# Difference between chocolates eaten by doctors and those eaten by nurses

# t-test: Survival time for different types of chocolates

# Extract chocolateData$SurvivalTime for each type of chocolate

Doctors<-chocolateData$SurvivalTime[chocolateData$WhoAte=="Doctor"]
Nurses<-chocolateData$SurvivalTime[chocolateData$WhoAte=="Nurse"]

# Is each data sample normal? Use Shapiro test to find out

shapiro.test(Doctors)
shapiro.test(Nurses)

# Both of these show as p<0.05, so they are not normally distributed
# Therefore create boxplots to see where the outliers are 

boxplot(Doctors, main="Doctors")
boxplot(Nurses, main="Nurses")

# Doctors outliers are above 200, so filter the data sample to remove them
# NB: This may take several iterations to remove all outliers

Doctors<-Doctors[Doctors<200]
boxplot(Doctors, main="Doctors")

# Nurses outliers are above 300, so filter the data sample to remove them
# NB: This may take several iterations to remove all outliers

Nurses<-Nurses[Nurses<300]
boxplot(Nurses, main="Nurses")

# Perform t-test for independent samples

# Now normal, so can do the t-test

# HO: There is no difference between the means, i.e. doctors and nurses eat chocolates at similar times
# H1: There is a difference between the means, i.e. doctors and nurses eat chocolates at different times

t.test(Doctors,Nurses)

# The calculated p-value for this is 0.7624. This is more than 0.05,
# so we cannot reject the null hypothesis at the p=0.05 level.
# We therefore accept the null hypothesis that doctors and nurses eat chocolates
# at about the same times



# ANOVA

# Differences between doctors, nurses, and administrators

# H0: There is no difference between doctors, nurses, and administrators
# H1: There is a significant difference between at least two of the groups

# Filter the dataSet to create a new dataframe with just the survival time and the person who ate the chocolate

anovaData<-chocolateData[,c(4,9)]
anovaData

# Filter the new dataset to extract only chocolates eaten by doctors, nurses and administrators
anovaData<-anovaData[anovaData$WhoAte %in% c("Doctor","Nurse","Administrative"),]
anovaData

# Carry out ANOVA test

# Compute the analysis of variance
resaov <- aov(anovaData$SurvivalTime ~ anovaData$WhoAte, data = anovaData)

# Summary of the analysis
summary(resaov)

#                    Df    Sum Sq   Mean Sq   F value   Pr(>F)  
# anovaData$WhoAte     2    45814     22907     3.382   0.0378 *
# Residuals          103   697674      6774                 
#  ---
#  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# There is a significant difference between groups at the p=0.05 level
# We can reject the null hypothesis and accept the alternative hypothesis
# It is therefore necessary to perform post-hoc tests to find out where these differences occur

# Perform Tukey Pair-wise post-hoc analysis

TukeyHSD(resaov)

#Tukey multiple comparisons of means
#95% family-wise confidence level

#Fit: aov(formula = anovaData$SurvivalTime ~ anovaData$WhoAte, data = anovaData)

#$`anovaData$WhoAte`
#                       diff        lwr       upr     p adj
#Doctor-Administrative  58.65367   4.006131 113.30122 0.0323376
#Nurse-Administrative   24.73349 -23.998754  73.46574 0.4517825
#Nurse-Doctor          -33.92018 -78.978504  11.13815 0.1779281

# The Tukey test shows that there is a significant difference at the p=0.05 level
# between doctors and administrators, but not between other pairs of groups



# Explore the differences between mean chocolate survival time at the four hospitals

# Extract the columns for chocolate survival time and location. There is no need to filter the
# rows as we want all elements

anovaData<-chocolateData[,c(1,9)]
anovaData

# Carry out ANOVA test

# H0: There is no difference between the survival times at any hospitals
# H1: There is a significant difference between at least two of the hospitals

# Compute the analysis of variance
resaov <- aov(anovaData$SurvivalTime ~ anovaData$Location, data = anovaData)

# Summary of the analysis
summary(resaov)

# There is a significant difference between groups at the p=0.05 level
# We can reject the null hypothesis and accept the alternative hypothesis
# It is therefore necessary to perform post-hoc tests to find out where these differences occur

TukeyHSD(resaov)

#Significant differences exist between the following pairs of hospitals:
# Bedford Pilgrim and Bedford AAU
# UHNS Haem-Bedford AAU 
# Swindon Meldon-Bedford Pilgrim
# UHNS Haem-Swindon Meldon