```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

#  This is an accompaniment to the lecture 5.
## It addresses tests for difference between multiple groups and between nominal variables 
##  It includes all the examples used in the lecture

### Preliminaries
```{r}
library(dplyr) 
library(semTools)
library(userfriendlyscience) #install the library userfriendlyscience and load it using the library command. It has a really nice one-way anova function that providesnice summary output
library(psych)
library(stats)
library(FSA)
```


```{r}
#We are using a .dat file (survey.dat) created from the SPSS file survey.sav  taken from SPSS Survival Manual 6th Edition Julie Pallant
#http://spss.allenandunwin.com.s3-website-ap-southeast-2.amazonaws.com/data-files.html#.Wb0vvnWP-po
#Results on a survey on well being
#We need to load the file so that we can use it in R.
sdata <- read.table("survey.dat")
#Setting the column names to be that used in the dataset
colnames(sdata) <- tolower(colnames(sdata))
#Question
#Is there a difference in optimism scores for young, middle-aged and old participants?
#One independent variable with three or more levels (age category)
#One continuous variable (optimism scores)
skew<-semTools::skew(sdata$toptim)
kurt<-semTools::kurtosis(sdata$toptim)
skew[1]/skew[2] #We have an issue with skew
kurt[1]/kurt[2]
#Not including plots in this example - you would normally
#looking at skew and kurtosis
#looking at standardised scores, we can see we have an outliers 
summary(scale(sdata$toptim)) #look at the max and min values to get the range, min value outside +/-3.29
head(sort(scale(sdata$toptim))) #need to check how many, we can see that it meets our criteria

#Use dplyr filter to filter out the rows with a standardised value outsude the range
outliers <- sdata %>% 
  filter(scale(toptim) >3.29 | scale(toptim) < -3.29)

#count them using dplyr summarize
outliers %>%
summarize(count=n())
#Meets our criteria
  
#Summary statistics by group (we know we can use mean and sd)
group_by(sdata, sdata$agegp3) %>% 
summarise( 
count = n(), 
mean = mean(toptim, na.rm = TRUE), 
sd = sd(toptim, na.rm = TRUE) )


#Check for homogeneity of variance
#We use Bartlett's test 

bartlett.test(sdata$toptim, sdata$agegp3)
#Can be argued that the variances are homogeneous if the p-value > 0.05

```
###ANOVA
```{r}
#Basic Approach
# Compute the analysis of variance
res.aov <- aov(sdata$toptim ~ sdata$agegp3, data = sdata)
# Summary of the analysis
summary(res.aov)
#Tukey pairwise comparison
TukeyHSD(res.aov)

#Alternative
#run User friendly science one-way anova test using the correct post-hoc test Tukey in our case
#Use Games-Howell for unequal variances
one.way <- userfriendlyscience::oneway(sdata$agegp3, y = sdata$toptim, posthoc = 'Tukey') 
 
#printout a summary of the anova 
one.way 


```

### There are repeated measures ANOVA available in SPSS which follow the same pattern of input 


### Kruskal-Wallis (Non-Parametric)
```{r}
#We are using youthcohort.dat (Taken from Quantitative Data Analysis in Education, Paul Connolly)
#created from the .sav version available at http://cw.routledge.com/textbooks/9780415372985/resources/datasets.asp


#We need to load the file so that we can use it in R.
ydata <- read.table("youthcohort.dat")
#Setting the column names to be that used in the dataset
colnames(sdata) <- tolower(colnames(ydata))
#Question:
#"Are there any differences between pupils of different ethnicity in England and Wales in relation to the grades they achieved in GCSE Maths?"
#Variables: 'ethsfr' ethnicity, 'gradmath' grade achieved in maths

#We know we are going to use non-parametric tests since the variable gradmath is ordinal
#Get some descriptive statistics on the maths grade for each ethnicity

psych::describeBy(as.numeric(ydata$gradmath),factor(ydata$ethsfr) )

#run a kruskal wallis test
stats::kruskal.test(gradmath~ethsfr,data=ydata)

#Post hoc testing

#Need library FSA to run the post-hoc tests

tmp<-FSA::dunnTest(x=as.numeric(ydata$gradmath), g=factor(ydata$ethsfr), method="bonferroni") 
#We need the as.numeric to make sure our gradmath variable is treated as a numeric variable and we are using the numerical categories for ethsfr so we use the factors
print(tmp, dunn.test.results = TRUE)
#in the output X=gradmath, g=ethsfr and the test statistic is Z with significance shown underneath

```


### Friedmans ANOVA (non-parametric) for repeated measures
```{r}
#We are using diet.dat (Taken from Discovering Statistics with R, Andy Field)
#created from the .sav version available at https://studysites.uk.sagepub.com/dsur/study/default.htm
#adapted to work as a dataframe on load

#Does the 'andikins' diet work?
#Variables
#Outcome: weight (Kg)
#Time since beginning the diet
#At Baseline (start)
#At 1 Month
#At 2 Months

#read in the data
diet <- read.table("diet.dat")
#treat our data as a matrix
Responses <- na.omit(with(diet, cbind(month1, month2, start)))
cat("\nMedians:\n")#display the heading Medians
#print out the medians of each of our variables
print(apply(Responses, 2, median)) 
#conduct the friedman test
friedman.test(Responses)

#Conduct post-hoc tests as needed
```
##Comparing Nominal Variables
###Contingency table
```{r}
#Using the bullying.dat dataset (Paul Connolly)
#created from the .sav version available at http://cw.routledge.com/textbooks/9780415372985/resources/datasets.asp
#Contains the results of national school level survey, subset of this survey which deals with bullying.
#Question:
#"Is there a difference between boys and girls in terms of whether they have been bullied at school?"
#H0:  There is no difference
#HA: There is a difference
#Variables: 
#'rsex' {1-Male, 2-Female, 99-not answered}; 
#'ubullsch' - Q26. Have you yourself ever been bullied in school? {1- Yes, 2 -No, 99-not answered}

#import the library gmodels and then load it 
library(gmodels)
bully<-read.table("bullying.dat")
#Use the Crosstable function
#CrossTable(predictor, outcome, fisher = TRUE, chisq = TRUE, expected = TRUE)
gmodels::CrossTable(bully$rsex, bully$ubullsch, fisher = TRUE, chisq = TRUE, expected = TRUE, sresid = TRUE, format = "SPSS")

#more simplistic way of doing Chi-Square

#Create your contingency table
mytable<-xtabs(~ubullsch+rsex, data=bully)

ctest<-chisq.test(mytable, correct=TRUE)#chi square test
#correct=TRUE to get Yates correction needed for 2x2 table

ctest#will give you the details of the test statistic and p-value
ctest$expected#expected frequencies
ctest$observed#observed frequencies
ctest$p.value



```
###For repeated measures for categorical variables use McNemars test mcnemars.test
