---
title: "Goodness of Fit and Contingency chi-squared tests"
author: "Emma Rand"
output:
  html_document:
    depth: 3
    highlight: pygments
    theme: lumen
    toc: yes
    toc_float:
      collapsed: no
      smooth_scroll: no
  word_document: default
  pdf_document:
    toc: yes
bibliography: "../refs/references.bib"
---
![](../../pics/17C.png)

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


# Introduction

## Aims

To explain how to analyse data that are counts falling into mutually exclusive categories using two types of chi-squared test and the difference between these two tests.

## Learning Outcomes 

By actively following the lecture and practical and carrying out the independent study the successful student will be able to:

* recognise when to use chi-squared Goodness of Fit and Contingency tests (MLO 2)
* be able to carry out, interpret and report scientifically both types in R (MLO 3 and 4)


## Philosophy

Workshops are not a test. It is expected that you often don't know how to start, make a lot of mistakes and need help.
Do not be put off and don't let what you can not do interfere with what you can do. You will benefit from collaborating with others and/or discussing your results.

The lectures and the workshops are [closely integrated](../BIO00017C.html) and it is expected that you are familar with the lecture content before the workshop. You need not understand every detail as the workshop should build and consolidate your understanding. You may wish to refer to the slides as you work through the workshop schedule.

## Slides

Goodness of Fit and Contingency chi-squared tests: [pdf](../lects/web03Chi-squared tests.pdf) (recommended) / [pptx](../lects/web03Chi-squared tests.pptx)

## Getting started

![W](../../pics/W.png) Start RStudio from the [Start menu](../../pics/startmenu.png).

![R](../../pics/R.png) In RStudio, [set your working directory](../../pics/menusetwd.png) to the folder you created last week for your 17C Data Analysis work.

![R](../../pics/R.png) Make a [new script](../../pics/newscript.png) file called practical3.R to carry out the rest of the work. 


# Exercises

## Inductions 

In a local maternity hospital, the total numbers of births induced on each day of the week over a six week period were recorded a

|Day  |No.inductions  | 
|---------------|---------- | 
|Monday  | 43  |  
|Tuesday  | 36  | 
|Wednesday  | 35  |  
|Thursday   | 38  | 
|Friday | 48 |
|Saturday | 26 |
|Sunday  | 24  |  
|Total  | 250  |

### Inductions - coding the ${\chi}^2$

We can use a chi-squared test to ascertain whether there is a pattern in these data that might suggest that surgeons are more reluctant to perform inductions on some days than on others.

![R](../../pics/R.png)  Make a vector `obs` that holds the number of inductions on each day.
```{r echo=FALSE}
#---CODING ANSWER---
obs <- c(43, 36, 35, 38, 48, 26, 24) 
```
![R](../../pics/R.png) Examine the 'structure' of the `obs` object using `str()`
```{r echo=FALSE, results="hide"}
#---CODING ANSWER---
str(obs)
# it's a vector containing integers
```
![Q](../../pics/Q.png) What is your null hypothesis and what type of test is required?

<!-- One obvious null hypothesis is that the proportion of inductions is the same on each day (as one might expect if medical considerations were of overriding concern) - we know what the proportions should be so we have a Goodness of Fit test. -->

We can carry out a Goodness of Fit chi-squared test on these data by coding the how the test works ourselves. 
If the null hypothesis is true, 1/7th of inductions would be expected to occur each day *i.e.,* 1/7th of 250

![R](../../pics/R.png) Assign the total number of inductions to a variable called `totalinductions` 
```{r}
# the total number of inductions
totalinductions <- sum(obs)

```
![R](../../pics/R.png) Now calculate the expected values: 
```{r}
# the expected values
# note that R mulitples every element in the c() by totalinductions
exp <- c(1/7, 1/7, 1/7, 1/7, 1/7, 1/7, 1/7) * totalinductions

```
The formula for chi-squared is: ${\chi}^2 = \Sigma \frac{(O - E)^2}{E}$

![R](../../pics/R.png) We can write this in R using: 
```{r}
# calculate the chisquared
chisqd <- sum(((obs - exp)^2) / exp)

```

We have  ${\chi}^2 =$ `r chisqd`

The degrees of freedom for a Goodness of fit chi-squared test are the number of categories minus one.

![R](../../pics/R.png) Calculate and assign the *d.f.* to a variable called `df`
```{r echo=FALSE, results="hide"}
#---CODING ANSWER---
df <- length(obs) - 1
# note: you might have used df = 6 but using the length of the 
# observation data is more reproducible. You can use the same code
#  for any goodness of fit test and only have to update the obs vector (the data).
```

![R](../../pics/R.png) The p value can be calculated with:
```{r}
pchisq(chisqd, df = df, lower.tail = FALSE)

```


![Q](../../pics/Q.png) What do you conclude about the pattern of birth inductions?

<!-- so it's a little over 0.05 and thus (just!), not significant. We have no reason to reject the null hypothesis that clinical needs dictate when inductions occur. -->


We could have decided to test a different null hypothesis. Instead of testing whether we had 1/7 of the inductions each day (7 categories) we could test whether we had 2/7 and 5/7 in two categories, 'weekends' and 'weekdays'.

![Q](../../pics/R.png) You can pool your data like this: 

```{r}
# Combining the weekdays and the weekend days
# first combine the obs
obs2 <- c(sum(obs[1:5]), sum(obs[6:7])) 

# then calculate the expected values
exp2 <- c(5/7, 2/7) * totalinductions

```
![Q](../../pics/R.png) Carry out this Goodness of fit test.
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
# calculate the chisquared
chisqd <- sum(((obs2 - exp2)^2) / exp2)
df <- length(obs2) - 1
pchisq(chisqd, df = df, lower.tail = FALSE)
# Signifcant
```

![Q](../../pics/Q.png) Do your conclusions alter?

<!-- p = 0.0027 - reject the null hypothesis, fewer inductions are being performed during the weekend (look at the peaks on Fridays and Mondays!).  Note that now we have divided the data into two categories - weekend and weekdays - so the chi-squared value is associated with just 1 d.f.  It appears that clinical need may not be uppermost in surgeon's minds!
This is why it is important to determine what you are going to test before you collect your data. Testing in different ways until you find significance is called p-hacking and is not legitimate.-->

### Inductions - alternative method (inbuilt test)

We have carried out a Goodness of Fit chi-squared test by coding the how the test works ourselves.
R also has an inbuilt function.
![Q](../../pics/R.png) The inbuilt function `chisq.test` is used like this: 

```{r}
chisq.test(obs)
```

This looks much easier but it is only that straightforward when the expected ratio numbers are all equal, *i.e.*, 1:1:1:1:1:1:1 

For other cases, such as we had for on combining the weekdays and the weekend days,  then we have to give the expected ratios.
![Q](../../pics/R.png) The inbuilt function `chisq.test` when expected values are not equal in each category: 

```{r}
exp2 <- c(5/7, 2/7)
chisq.test(obs2, p = exp2)
```

### Inductions - alternative data format

The data we were given was already tabulated however, data often comes to us with an observation in each row. In this case that would mean a column of 250 values which indicated which day and induction was performed. Data formated this way is in [inductions.txt](../data/inductions.txt).

![W](../../pics/W.png) You should have made a [folder](../../pics/folderdata.png) in your working directory called 'data' last week. If you have not, look it up [here](http://www-users.york.ac.uk/~er13/17C%20-%202018/pracs/02TestingDataTypesReadingInData.html#import_data_2)

![W](../../pics/W.png) Save a copy of [inductions.txt](../data/inductions.txt) to the 'data' folder

![R](../../pics/R.png) Read the data in to R. Remember, you need to use the relative path to the file in the `read.table()` command and decide whether you need header to be TRUE or FALSE:

```{r echo=FALSE}
#---CODING ANSWER---
#I have my files organised differently so I use a different 'path' to my file
inductions <- read.table("../data/inductions.txt", header = T)
```

![R](../../pics/R.png) look at the structure of the dataframe using `str()`
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
str(inductions) 
```

![R](../../pics/R.png) We can tabulate the data and assign it using the `table()` command:
```{r}
indtab <- table(inductions$day)
```

The `table()` outputs the days in alphbetical order which is a common defualt in R. This does not matter sinces days are just arbitrary categories for the chi-squared test.

![R](../../pics/R.png) Now carry out the 7 category goodness of fit like this:
```{r}
chisq.test(indtab)
```

### Illustrating

![R](../../pics/R.png) We could create a figure using ggplot like this:
```{r}
# load tidyverse for ggplot2 - you only have to load once in a session
# tidyverse gives us ggplot and several other packages. tidyverse is known as a metapackage (it's several packages)
library(tidyverse)

ggplot(data = inductions, aes(x = day)) +
  geom_bar() +
  theme_classic()

# themes change several things at once about the appearance. 
# Classic is useful one but there are others you can play
# with. Type theme_ and use the tab key to see the options
```

Remember - the ggplot geom_bar() will count the number of each day for us which is great! `r emo::ji("party")`

What is not so great is the alphbetical ordering `r emo::ji("annoyed")`

```{r}
 levels(inductions$day)
```


![R](../../pics/R.png) Reordering factor levels is something we quite often want to do.

```{r}
inductions$day <- fct_relevel(inductions$day,
                              "mon", 
                              "tues",
                              "wed",
                              "thurs",
                              "fri",
                              "sat",
                              "sun")

levels(inductions$day)
```

![R](../../pics/R.png) Redraw the figure.

## Blood group and ulcers

A human geneticist is interested in whether there is an association between ABO blood group phenotypes and the incidence of peptic ulcers.  He found that in a sample of 477 blood group O people 65 had peptic ulcers whereas in a sample of 387 blood group A people 31 had peptic ulcers.

![Q](../../pics/Q.png) Draw a table of these data (on a piece of paper).

<!-- 2 x 2 contingency table. -->

<!-- 	            Blood group	 Total -->
<!-- 		            O	    A	         -->
<!-- Ulcers	Yes	   65 	 31     	96 -->
<!--       	No  	412 	356 	   768 -->
<!-- Total		    477  	387	     864 -->


![Q](../../pics/Q.png) What is your null hypothesis and what type of test is required?

<!-- Here we are scoring two characters on each individual and we want to see if there is an association between the characters i.e. we want to perform a test of association. Our null hypothesis is that there is no difference in the proportion of people with peptic ulcers in the two blood-group phenotypes  so we have a contingency test -->

### Blood group and ulcers - inbuilt test

![R](../../pics/R.png) Make a vector `obs` that holds the observed numbers. For the moment, don't worry about what order they are in: 
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
obs <- c(65, 412, 31, 356)
```

For a contingency chi squared test, the inbuilt chi-squared test can be used in a straightforward way. However, we need to structure our data as a 2 x 2 table rather than as a 1 x 4 vector. A 2 x 2 table can be created with the `matrix()` function. We can also name the rows and columns which helps us interpret the results.

![R](../../pics/R.png) To create a list containing two elements which are vectors for the two groups in each variable we do:

```{r}
# list of two elments
# the two variables are whether someone has an ulcer or not and whether they are blood group O or A
vars <- list(ulcer = c("yes","no"), blood = c("O", "A"))
vars

```
![R](../../pics/R.png) Now we can create the matrix from our vector of numbers `obs` and use our list `vars` to give the column and row names:
```{r}
ulcers <- matrix(obs, nrow = 2, dimnames = vars)
```

![R](../../pics/R.png) Check the content of `ulcers` and recreate if the numbers are not in the correct place (*i.e.*, do not match your table) 

![R](../../pics/R.png) Run a contingency chi-squared with:

```{r}
chisq.test(ulcers, correct = FALSE)
# you should look up the command in the manual to see what correct = FALSE does
```
To help us discover what is the direction of any deviation from the null hypothesis it is helpful to see what the expected values were. These are accessible in the `$expected` variable in the output value of the `chisq.test()` method (See the manual!).

![R](../../pics/R.png) View the expected values with:

```{r}
chisq.test(ulcers, correct = FALSE)$expected
```

![Q](../../pics/Q.png) What do you conclude about the association between ABO blood group and peptic ulcers?

<!-- People with group O seem to get more ulcers: 65 observed is higher than the expected 53. People with group A have fewer than expected: 31 observed is lower than the expected 43. -->

<!-- For a contingency table such as this d.f. = (r - 1)(c - 1) = 1 i.e. (number of rows - 1)( number of columns - 1) -->

### Blood group and ulcers- alternative data format.

The data we were given was already tabulated. There are raw data in [blood_ulcers.txt](../data/blood_ulcers.txt). Examine this file and compare its content to [inductions.txt](../data/inductions.txt)

![W](../../pics/W.png) Save a copy of [blood_ulcers.txt](../data/blood_ulcers.txt) to the 'data' folder

![R](../../pics/R.png) Read the data in to R. Remember, you need to use the relative path to the file in the `read.table()` command and decide whether you need header to be TRUE or FALSE:

```{r echo=FALSE}
#---CODING ANSWER---
#I have my files organised differently so I use a different 'path' to my file
blood_ulcers <- read.table("../data/blood_ulcers.txt", header = T)
```

![R](../../pics/R.png) look at the structure of the dataframe using `str()`
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
str(blood_ulcers) 
```

Notice that we have two variables this time.

![R](../../pics/R.png) We can tabulate the data and assign it using the `table()` command:
```{r}
ulctab <- table(blood_ulcers$blood, blood_ulcers$ulcer)
# examine the result
ulctab
```

We need to give both variables to cross tabulate.

![R](../../pics/R.png) Now carry out the contingency chi-squared like this:
```{r}
chisq.test(ulctab, correct = FALSE)
```



## Genetics

In a particular species of ground beetle, red and yellow elytral colours are controlled by alleles at one locus, with red (R) dominant to yellow (r). Spotted and plain elytral patterns are controlled by alleles at another locus, with spotted elytra (S) dominant to plain (s).  You set up a cross between two indivuals that are heterozygous at both loci (RrSs). If the two loci are unlinked you would expect a phenotype ratio of 9 : 3 : 3 : 1 ratio of red/spotted : red/plain : yellow/spotted : yellow plain, respectively.  The numbers you actually observe for these four phenotypes are: 79, 17, 19 and 20 (total = 135).

![R](../../pics/R.png)  Make a vector `obs` that holds the number of phenotype .
```{r echo=FALSE}
#---CODING ANSWER---
obs <- c(79, 17, 19, 20) 
```
![R](../../pics/R.png) Examine the 'structure' of the `obs` object using `str()`
```{r echo=FALSE, results="hide"}
#---CODING ANSWER---
str(obs)
# it's a vector containing integers
```
![Q](../../pics/Q.png) What is your null hypothesis and what type of test is required?
<!-- Here we expected values dictated by Mendel's Laws - we would expect a 9:3:3:1 ratio in our total of 135 beetles. This is Goodness of fit because we are testing the fit of our data to a known ratio - our expectations are based on a particular theory that determines 9:3:3:1. We are not making setting the expected values based on the observed ration. -->

![Q](../../pics/R.png) Carry out the appropriate chi-squared test on these data either with or withour using inbuilt test, whatever you prefer.

```{r echo=FALSE, results="hide"}
# using the inbuilt test
# give the expected ratio
exp <- c(9/16, 3/16, 3/16, 1/16)
chisq.test(obs, p = exp)

# we can examine the expected values
chisq.test(obs, p = exp)$expected
```

![Q](../../pics/Q.png) What do you conclude about results?
<!-- Reject  X-squared = 20.272, df = 3, p-value = 0.000149 -->
<!-- There is an excess of red/spotted and yellow/plain and a deficiency of red/plain and yellow/spotted.  Considering the character states of the original parents, these deviations are in a direction consistent with the two loci being linked on the same chromosome.  -->

# Independent study

## Slug colour forms - tabulated

The slug *Arion ater* has three major colour forms, black, chocolate brown and red.  Sampling in population X revealed 27 black, 17 brown and 9 red individuals, whereas in population Y the corresponding numbers were 39, 10 and 21.
Create an appropriate data structure and run the appropriate chi-squared test

```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
##################################################################################
#                            Slugs - tabulated                                   #
##################################################################################
#  name the rows and columns
vars <- list(pop = c("x","y"), colour = c("black", "brown", "red"))
slugs <- matrix(c(27, 39, 17, 10, 9, 21), nrow = 2, dimnames = vars)
# gives me
#    colour
# pop black brown red
#   x    27    17   9
#   y    39    10  21

# for a contingency test like this one where the expected values are
# derived from the data, the inbuilt chi squared test works well
#
chisq.test(slugs)
# X-squared = 6.5726, df = 2, p-value = 0.03739
# p < 0.05 so we reject the null hypothesis i.e., the proportions of the colour 
# forms are significantly different in the two populations 
# (mostly as a result of differences in the brown and red classes - 
# look at the differences between observed and expected values for 
# the three colour forms in the table above).
chisq.test(slugs)$expected
#       colour
# pop    black    brown      red
#   x 28.43902 11.63415 12.92683
#   y 37.56098 15.36585 17.07317

```

## Slug colour forms - raw data
The raw, untabulated data are in [slugs.txt](../data/slugs.txt). Perform the test on these data.


```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
##################################################################################
#                            Slugs - raw                                         #
##################################################################################
# import the data
slugs <- read.table("../data/slugs.txt", header = T)

# put it into a table
slugtab <- table(slugs$colour, slugs$pop)

# carry out the test
chisq.test(slugtab)
```
# The Code files

These contain answers and code even though they do not appear on the webpage itself.

[Rmd file](03ChiSquaredTests.Rmd)
The Rmd file is the file I use to compile the practical. Rmd stands for R markdown allow R code and ordinary text to be inter weaved to produce well-formatted reports including webpages.

[Plain script file](../scripts/03ChiSquaredTests.R)
This is plain script (.R) version of the practical



# Objectives from previous sessions

[Introduction to module and RStudio](01IntroductionToModuleAndRStudio.html)

* to explain why we need statistical tests and the logic of hypothesis testing (MLO 1)
* use the R command line as a calculator and to assign variables (MLO 3)
* create and use the basic data types in R (MLO 3)
* find their way around the RStudio windows (MLO 3)
* create, use and save a script file to run r commands (MLO 3)
* search and understand manual pages (MLO 3)

[Testing, Data types and reading in data](02TestingDataTypesReadingInData.html)

* to able to explain what response and explanatory variables are, distinguish between data types and describe how these impact choice of test (MLO 1 and 2)
* demonstrate the process of hypothesis testing with an example and evaluate potential inferences (MLO 1 and 2)
* read in data in to RStudio, create simple summaries and plots using manual pages where necessary (MLO 3)
* create neat reports in Word which include text and figures (MLO 4)

![](../../pics/17Cend.png)