---
title: "Calculating summary statistics, probabilities and confidence intervals"
author: "Emma Rand"
output:
  html_document:
    depth: 3
    highlight: pygments
    theme: lumen
    toc: yes
    toc_float:
      collapsed: no
      smooth_scroll: no
  pdf_document:
    toc: yes
  word_document: default
---
![](../../pics/17C.png)

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

# Introduction

## Aims

You will practice calculating and choosing summary statistics, quantiles and confidence intervals.
 

## Learning Outcomes 

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

* Explain the properties of 'normal distributions' and their use in statistics (MLO 1 and 2)
* Define, select and calculate with R probabilities, quantiles and confidence intervals (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

Calculating summary statistics, probabilities and confidence intervals: [pdf](../lects/web04Describing normal distributions and Confidence Intervals.pdf) (recommended) / [pptx](../lects/web04Describing normal distributions and Confidence Intervals.pptx)

![Artwork by Allison Horst](../../pics/horst/not_normal.png)




## 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 previously for your 17C Data Analysis work.

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

# Exercises

## Distributions: the R functions

For any distribution, two very useful quantities can be calculated:   
* the Distribution Function, which gives the **probability** that a variable takes a particular value or less.  
* the Quantile function which is the inverse of the Distribution function, *i.e.*, it returns the value ('quantile') for a given probability.

The functions are names with a letter `p` or `q` preceding the distribution name. Below are some examples:

```{r echo=FALSE}
library(knitr)
library(kableExtra)
stattable <- read.table("stattable.txt", 
                        sep = "/", 
                        header = TRUE, 
                        row.names = 1)
kable_styling(kable(stattable, "html") )

```

## Calculating Probabilities for single value

### Using `pnorm()`

Look up `pnorm()` in the manual using `?pnorm`

You give it a values for which you want a probability and by default it gives you the probability of getting that value or less from a normal distribution with a mean of 0 and a standard devation of 1. If you want the probability of a vlaue from a different normal distribution you need to set the mean and standard deviation appropriately.

For example, I.Q. in the U.K. population is normally distributed with a mean of 100 and a standard deviation of 15. We can use `pnorm()` to calculate probabilities associated with having a particular range of IQs.

```{r echo=FALSE, fig.width=4,fig.height=4}
# figure to show IQ
# only for illustration    
m <- 100
sd <- 15
IQ <- seq(40, 160, 1)
curve(dnorm(x, 100, 15),
      xlim = c(40, 160), bty = "n",
      axes = F, xlab = "IQ", ylab = "") 
axis(1, pos = 0)
    
```


We can use the values of mean = 100 and standard deviation = 15 in `pnorm()` to work out the probability of having an I.Q. of _115 or less_. 

![R](../../pics/R.png) First, create variables for the parameter values - this is considered good practice.
```{r}
# create variables for the parameter values
m <- 100
sd <- 15
```

![R](../../pics/R.png) Now pass those variables to the `pnorm()` function along with the value for which we want a probability:

```{r}
# Now use pnorm()
pnorm(115, m, sd)
```

![R](../../pics/R.png) Look at the manual page. Because the default is lower.tail = TRUE, we get the probability we want, P[IQ < 115]

```{r echo=FALSE, fig.width=4,fig.height=4}
# fig to illustrate IQ < 115
# only for illustration
IQ <- seq(40, 160, 1)
cord.x <- c(40, seq(40, 115, 1), 115) 
cord.y <- c(0, dnorm(seq(40, 115, 1), 100, 15), 0) 
curve(dnorm(x, 100, 15),
      xlim = c(40, 160),
      bty = "n",axes = F,
      xlab = "IQ", ylab = "") 
polygon(cord.x, cord.y, col = "gray")
axis(1, pos = 0)
```
I recommend sketching the distribution and shading the area you want to work out what arguments you want to give the function.

![R](../../pics/R.png) Determine  the probability of having an IQ of 115 OR MORE? Do a sketch first.


```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
pnorm(115, m, sd, lower.tail = FALSE)
```    

![R](../../pics/R.png) Determine the probability of having an IQ between 85 and 115? Do a sketch first.
```{r echo=FALSE, fig.width=4,fig.height=4}
# fig to illustrate IQ between 85 and 115
# You sketch should look like this:
IQ <- seq(40, 160, 1)
cord.x <- c(85, seq(85, 115, 1), 115) 
cord.y <- c(0, dnorm(seq(85, 115, 1), 100, 15), 0) 
curve(dnorm(x, 100, 15),
      xlim = c(40, 160),
      bty = "n",axes = F,
      xlab = "IQ", ylab = "") 
polygon(cord.x, cord.y, col = "gray")
axis(1, pos = 0)
```
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
pnorm(115, m, sd) - pnorm(85, m, sd)
```   

![Q](../../pics/Q.png)Is this what you expect?
<!-- #---THINKING ANSWER--- -->
<!-- # You should expect this if you remember from the lecture that about 68% of the values are between - 1 s.d. and + 1 s.d and that 15 is the s.d. in this case -->

![R](../../pics/R.png) What is 1.96 * the standard deviation
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
v <- 1.96 * 15
```   
![R](../../pics/R.png) What is the probability of having an IQ between -1.96 standard deviations and +1.96 standard deviations? Is this what you expect?
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
pnorm(m + v, m, sd) - pnorm(m - v, m, sd)
# 95% because 95% of values lies between +1.96 s.d. and -1.96 s.d.

```   


### Using qnorm()
We can use `qnorm()` to find the IQ associated with a particular probability.

We will again use the values of mean = 100 and standard deviation = 15 in `qnorm()` to work out what I.Q. value 0.2 (20%) of people fall below. Make sure you relate the manual information to the command.

![R](../../pics/R.png) To find the I.Q. value that 20% people fall below:
```{r}
 qnorm(0.2, m, sd)
```

```{r echo=FALSE, fig.width=4,fig.height=4}
# fig to illustrate IQ value that 20% people are below

IQ <- seq(40, 160, 1)
cord.x <- c(0, seq(0, 87, 1), 87) 
cord.y <- c(0, dnorm(seq(0, 87, 1), 100, 15), 0) 
curve(dnorm(x, 100, 15),
      xlim = c(40, 160),
      ylim = c(-0.006, 0.03),
      bty = "n",axes = F,
      xlab = "IQ", ylab = "") 
polygon(cord.x, cord.y, col = "gray")
text(65, 0.01, "p = 0.2")
arrows(65, 0.009, 75, 0.004, length = .15)
text(87, -0.006, "?", col = "red")
arrows(87, -0.004,87, 0, length = .15, col = "red")
axis(1, pos = 0, labels = FALSE)


```

20% people have an IQ less than `r round(qnorm(0.2, m, sd),1) `

![R](../../pics/R.png) What I.Q. value are 0.025 (2.5%) of people below?
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
qnorm(0.025, m, sd)
``` 

![R](../../pics/R.png) In what range do 99% of the population fall? Note that 99% means 1% (0.01) in both tails so 0.5% (0.005) in each tail. The figure may help you.
    
```{r echo=FALSE, fig.width=6,fig.height=4.5}
# the figure for 99% 
# only for illustration
z <- qnorm(0.995)
cord.x <- c(-z, seq(-z, z, 0.01), z) 
cord.y <- c(0, dnorm(seq(-z, z, 0.01)), 0) 
curve(dnorm(x, 0, 1),
      xlim = c(-3, 3),
      ylim = c(-0.15, 0.45),
      bty = "n", 
      axes = F, xlab = "IQ",
      ylab = "") 
polygon(cord.x, cord.y, col = "gray")
arrows(-2.7, 0.1, -2.7, 0.01, length = .15)
text(-2.6, 0.11, "p = 0.005")
arrows(2.7, 0.1, 2.7, 0.01, length = .15)
text(2.6, 0.11, "p = 0.005")
text(0, 0.13, "p = 0.99")
text(z, -0.06, "?", col = "red")
arrows(z, -0.04,z, 0, length = .15, col = "red")
text(-z, -0.06, "?", col = "red")
arrows(-z, -0.04,-z, 0, length = .15, col = "red")
axis(1, pos = 0, labels = FALSE)
```
    
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
qnorm(0.005, m, sd)
qnorm(0.995, m, sd)
# between 61.4 and 138.6
```


## Calculating Probabilities for samples

The only difference in using `pnorm()` and `qnorm()` for samples is in what we give as the sd argument. Since we are now thinking about the distribution of the sample means, we need to use the standard error.

We used mean = 100 and standard deviation = 15 in `pnorm()` to work out the probability of an _individual_ having an I.Q. of _115 or less_. 

We can use a similar approach to find the probability of getting a _sample_ of n = 5 having a mean I.Q. of 115 or less
The only difference is that we use the standard error instead of the standard deviation. 

![R](../../pics/R.png) First, calculate the standard error:

```{r}
n <- 5
se <- 15 / sqrt(n)
```
![R](../../pics/R.png) Now the probability of getting a sample mean of 115 or less from that distribution:

```{r}
pnorm(115, m, se)
```


There's a `r round(pnorm(115, m, se), 4) ` probability that a sample of 5 people will have a mean of 115 or less. Thus there is a probability of just `r round(pnorm(115, m, se, lower.tail = FALSE), 4)` that a sample of n = 5 will have a mean above 115. This is quite unlikely and we might suspect this group was not sampled from the general population.

![R](../../pics/R.png) What is the probability of sample of size 10 having a mean of 105 or more?
    
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
n <- 10
se <- 15 / sqrt(n)
pnorm(105, m, se, lower.tail = FALSE)
# it is about 0.146
``` 

## Confidence intervals (large samples)

The data in [beewing.txt](../data/beewing.txt)  are left wing widths of 100  honey bees (mm). The confidence interval for large samples is given by: $\bar{x} \pm 1.96 \times s.e.$)

Where 1.96 is the quantile for 95% confidence.

You may need to refer to previous practicals to remind yourself how to carry out some of the following steps.

![W](../../pics/W.png) Save a copy of the file. I saved mine to my 'data' directory

![R](../../pics/R.png) Read in the data and check the structure of the resulting dataframe
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
bee <- read.table("../data/beewing.txt", header = FALSE)
str(bee)
``` 
![R](../../pics/R.png) Rename the column to 'wing'
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
names(bee) <- "wing"
str(bee)
``` 
![R](../../pics/R.png) Calculate and assign to variables: the mean, standard deviation and standard error
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
m <- mean(bee$wing)
sd <- sd(bee$wing)
n <- length(bee$wing)
se <- sd / sqrt(n)
``` 
![R](../../pics/R.png) To calculate the 95% confidence interval we need to look up quantile (multiplier) using `qnorm()`
```{r warning=FALSE,message=FALSE}
q <- qnorm(0.975)
```

![R](../../pics/R.png) Now we can use it in our confidence interval calculation
```{r warning=FALSE,message=FALSE}
(lcl <- m - q * se)
(ucl <- m + q * se)
``` 
![R](../../pics/R.png) Between what values would you be 99% confident of the population mean being?
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
q <- qnorm(0.995)
lcl <- m - q * se
ucl <- m + q * se
#---THINKING ANSWER---
# the mean is 4.55 +/- 0.077. We are 95% certain the populations mean lies between 4.627 and 4.473
``` 

## Confidence intervals (small samples)

The confidence interval for small samples is given by: $\bar{x} \pm \sf t_{[d.f]} \times s.e.$

The fatty acid Docosahexaenoic acid (DHA) is a major component of membrane phospholipids in nerve cells and deficiency leads to many behavioural and functional deficits. The cross sectional area of neurons in the CA 1 region of the hippocampus of normal rats is 155 $\mu m^2$. A DHA deficient diet was fed to 8 animals and the cross sectional area (csa) of neurons is given in [neuron.txt](../data/neuron.txt)

![W](../../pics/W.png) Save a copy of the file. I saved mine to my 'data' directory

![R](../../pics/R.png) Read in the data and check the structure of the resulting dataframe
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
neur <- read.table("../data/neuron.txt", header = TRUE)
str(neur)
``` 

![R](../../pics/R.png) Assign the mean to m
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
m <- mean(neur$csa)

``` 

![R](../../pics/R.png) Calculate and assign the standard error to se
```{r echo=FALSE,results="hide"}
#---CODING ANSWER---
sd <- sd(neur$csa)
n <- length(neur$csa)
se <- sd / sqrt(n)
``` 

To work out the confidence interval for our sample mean we need to use the _t_ distribution because it is a small sample. This means we need to determine the degrees of freedom (the number in the sample minus one). 

![R](../../pics/R.png) We can assign this to a variable using:
```{r}
df <- length(neur$csa) - 1; df
```

![R](../../pics/R.png) The _t_ value is found by:
```{r}
t <- qt(0.975, df = df); t
```

![R](../../pics/R.png) And the confidence interval by:
```{r}
round(m + t * se, 2)
round(m - t * se, 2)
```

![Q](../../pics/Q.png) Given the upper and lower confidence values for the estimate of the population mean, what do you think about the effect of the DHA deficient diet?
<!--
#---THINKING ANSWER---
the 95% confidence interval from 132.75 $\mu m^2$ to 151.95 $\mu m^2$ doesn't include 155 $\mu m^2$ so we can conclude that DHA deficiency significantly lowers csa
--> 

# Independent study

You need to carry out this work before the next practical.

## Reading

[An introduction to the normal distribution](https://tinystats.github.io/teacups-giraffes-and-statistics/02_bellCurve.html) from the Teacups, giraffes and statistics book (the whole online book is listed in the Additional Resources folder on the VLE).

## Confidence interval

Adiponectin is exclusively secreted from adipose tissue and modulates a number of metabolic processes. Nicotinic acid can affect adiponectin secretion. 3T3-L1 adipocytes were treated with nicotinic acid or with a control treatment and adiponectin concentration (pg/mL) measured. The data are in [adipocytes.txt](../data/adipocytes.txt). Each row represents an independent sample of adipocytes and the first column gives the concentration adiponectin and the second column indicates whether they were treated with nicotinic acid or not. Estimate the mean Adiponectin concentration in each group - this means calculate the sample mean and construct a confidence interval around it for each group.
Hint: you will find 'tapply()' which was in a previous practical useful.
```{r echo=FALSE,results="hide"}
#---CODING AND THINKING ANSWER---
adip <- read.table("../data/adipocytes.txt", header = T)
str(adip)
library(tidyverse)
adip_summary <- adip %>% 
  group_by(treatment) %>% 
  summarise(mean = mean(adiponectin),
            sd = sd(adiponectin),
            n = length(adiponectin),
            se = sd/sqrt(n))
#since there are only 15 values in each group we need to use the t value to consruct our ci
t <- qt(0.975, df = 14)
#And the confidence interval are:
round(adip_summary$mean + t * adip_summary$se, 2)
round(adip_summary$mean - t * adip_summary$se, 2)
#we conclude we're 95% certain the mean for the conrol group is between 4.73 and 6.36 and the mean for the nicotinic group is between 6.52 and 8.50
#More usually we might put is like this:
round(t * adip_summary$se, 2)
# the mean for the control group is 5.55 +/- 0.82 and that for the nicotinic group is 7.51 +/- 0.99

```

## Probability of misdiagnosis

Healthy people have Thyroid Stimulating Hormone (TSH) levels of (mean $\pm$ s.d) 3 $\pm$ 2.8 units per mL of blood, those experiencing hypothyroidism have elevated TSH. Individuals with TSH 6 units per mL or higher are treated for hypothyroidism. What is the probability of being misdiagnosed with hypothyroidism?
```{r echo=FALSE,results="hide"}
   #---CODING AND THINKING ANSWER---
#This question is asking what is the probability of having a value of 6 or more if you come from the healthy people distribution
m <- 3
sd <- 2.8
# Now use pnorm()
pnorm(6, m, sd, lower.tail = F)
#it is 0.142
``` 

## Geting a bit ove`R`whelmed?

Try [Chapter 3 Getting started with R](https://learningstatisticswithr.com/book/introR.html#gettingR) of Danielle Navarro's book. 


# The Code files

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

[Rmd file](04CalculatingSummaryStatisticsProbabilitiesAndCI.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/04CalculatingSummaryStatisticsProbabilitiesAndCI.R)
This is plain script (.R) version of the practical


# Script example
This is an example of a well formatted analysis script 

[Script example](../scripts/exampleprac4.R)

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

[Goodness of Fit and Contingency chi-squared tests](03ChiSquaredTests.html)

* 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)

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


