---
title: "Normal distributions with R"
# author: "Mark"
header-includes:
   - \usepackage{color, colortbl}
   - \definecolor{LightGray}{gray}{0.9}
   - \definecolor{Gray}{gray}{0.75}
date: "9/8/2017"
output: pdf_document
---

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

```{r, echo = F, fig.pos="h", fig.align='center',fig.cap="\\label{fig:standard_normal}The standard normal"}
dx = 0.01
x = seq(-4,4, by=dx)
y = dnorm(x)
plot(x,y, type='l', ylim=c(0,0.6), yaxs='i', xlim=c(-3,3), xlab="", ylab="", yaxt='n', lwd=3)
x <- seq(-4,-0.5, by=dx)
y <- dnorm(x)
polygon(c(x,rev(x)), c(rep(0, length(y)), rev(y)), col='lightgray')
text(-1,0.1,"A = ?")
axis(side=1, at=c(-0.5))
#points(-0.5,dnorm(-0.5), pch=21,  bg="blue")
#text(-1.1, dnorm(-0.5)+0.02, "(-0.5,0.352)")
```

A fundamental question in statistics is - how do we compute and interpret area computations like the one shown above? In this little lab, we'll learn exactly how to do so with R.

## The `pnorm` command

Suppose we want to compute the shaded area shown in figure \ref{fig:standard_normal} above - that is, the area under the standard normal curve and to the left of the line of the value $Z=-0.5$. It's a simple matter of executing the following command:

```{r, eval = F}
pnorm(-0.5)
```

Give it a try! Of course, the result should agree with the result obtained from a table. A portion of the our normal distribution table is shown in table 1 below with the relevant row and column highlighted.

\begin{table}[!htbp]
\setlength{\belowcaptionskip}{15pt plus 3pt minus 2pt}
\caption{A few relevant $Z$ values}
\begin{center}{\small
\begin{tabular}{| c | rrrr>{\columncolor{LightGray}}r | c}
  & 0.04 &  0.03 &  0.02 &  0.01 &  0.00 & $Z$  \\
  \hline
  \hline
  \vdots\ \vdots & \scriptsize{0.2296} & \scriptsize{0.2327} & \scriptsize{0.2358} & \scriptsize{0.2389} & \scriptsize{0.2420} & $-0.7$ \\
  \vdots\ \vdots & \scriptsize{0.2611} & \scriptsize{0.2643} & \scriptsize{0.2676} & \scriptsize{0.2709} & \scriptsize{0.2743} & $-0.6$ \\
  \rowcolor{LightGray}
  \vdots\ \vdots & \scriptsize{0.2946} & \scriptsize{0.2981} & \scriptsize{0.3015} & \scriptsize{0.3050} & \cellcolor{Gray}\scriptsize{0.3085} & $-0.5$ \\
  \vdots\ \vdots & \scriptsize{0.3300} & \scriptsize{0.3336} & \scriptsize{0.3372} & \scriptsize{0.3409} & \scriptsize{0.3446} & $-0.4$ \\
  \vdots\ \vdots & \scriptsize{0.3669} & \scriptsize{0.3707} & \scriptsize{0.3745} & \scriptsize{0.3783} & \scriptsize{0.3821} & $-0.3$ \\
    \hline
\end{tabular}}
\end{center}
\end{table}


## A basic example

Let's motivate our further exploration with a simple example: *Suppose the height of third graders at Isaac Dickson is normally distributed with a mean of 50.2 inches and a standard deviation of 2.66 inches. If there are 87 third graders, about how many do we expect to be taller than 54 inches?*

Note that this problem is phrased very much like a quiz or exam problem and we can approach it as such. However, when we get to where we need to compute the final answer, we'll use the R `pnorm` command, rather than a table.

*Solution (through the standard normal)*: We first compute the $z$-score as usual:
$$Z = \frac{54-50.2}{2.66} = 1.428571.$$
Of course, that computation is super easy to type into R as `(54-50.2)/2.66`. The point behind the $z$-score is that we can now compute the proportion of the kids taller than 54 inches as the area under the standard normal curve and to the *right* of $z = 1.428571$. This is the *white* area under the curve in the top plot of figure \ref{fig:two_normal_curves} shown below. Note that the `pnorm` command computes the so-called *cumulative* probability. Thus, it can compute the shaded gray region shown in the figure. To compute the white region, we'll need to subtract this from one. Taking this into account, our answer can be computed as:

```{r, eval = F}
(1 - pnorm(1.428571))*87
```


```{r, echo=F, fig.width=7,fig.height=6,fig.cap="\\label{fig:two_normal_curves}Two normal curves for our example problem"}
par(
  mfrow=c(2,1),
  oma = c(5,4,0,0) + 0.1,
  mar = c(0,0,1,1) + 0.1
)

dx = 0.01
x = seq(-4,4, by=dx)
y = dnorm(x)
plot(x,y, type='l', ylim=c(0,0.6), yaxs='i', xlim=c(-3,3), xlab="", ylab="", yaxt='n', lwd=3, axes=F)
x <- seq(-4,1.43, by=dx)
y <- dnorm(x)
polygon(c(x,rev(x)), c(rep(0, length(y)), rev(y)), col='lightgray')
axis(side=1, at=c(-4,-1,0,1.43, 4))
axis(side=2, at=c(0,0.5))


dx = 2.66*0.001
x = seq(50.2-3*2.66,50.2+3*2.66, by=dx)
y = dnorm(x, 50.2, 2.66)
plot(x,y, type='l', yaxs='i', ylim=c(0,0.6/2.66), xlim=c(50.2-3*2.66,50.2+3*2.66), 
     xlab="", ylab="", yaxt='n', lwd=3, axes=F)
x <- seq(50.2-3*2.66,54, by=dx)
y <- dnorm(x, 50.2, 2.66)
polygon(c(x,rev(x)), c(rep(0, length(y)), rev(y)), col='lightgray')
axis(side=1, at=c(30,50.2-2.66,50, 54, 70))
axis(side=2, at=c(0,0.2))
```


*Solution (through another normal)*: Using R, we can actually skip the standardization step altogether. The reason is that the `pnorm` command accepts more arguments allowing us to specify the mean and standard deviation that we'd like to use. Thus, another command that computes the answer is:

```{r, eval = F}
(1 - pnorm(54, 50.2, 2.66))*87
```
Give it a try! The result should agree with your previous computation.

### Problem 1

Suppose that cholesterol levels of an adult can be described by a normal model with a mean of 190 mg/dL and a standard deviation of 20.

 * What percent of adults do you expect to have cholesterol levels over 200 mg/dL?
 * What percent of adults do you expect to have cholesterol levels between 170 mg/dL and 200 mg/dL?

## The `qnorm` command

The `qnorm` command is essentially an inverse of the `pnorm` command - that is, given a quantile `q`, `qnorm(q)` returns the $z$-value at that quantile. For example, can you interpret the following result?

```{r, eval = F}
qnorm(0.5)
```

You can also specify the mean `m` and standard deviation `s`. For example, we can compute the theoretical interquartile range for our Isaac Dickson heights (where the mean and standard deviation were 50.2 and 2.66) as follows:

```{r, eval=F}
qnorm(0.75, 50.2, 2.66) - qnorm(0.25, 50.2, 2.66)
```

### Problem 2

Find the theoretical interquartile range for the cholesterol levels of problem 1.

## Working with data

Let's grab some real world data and use `qqnorm` to generate a normal probability plot to see if it looks normally distributed.

```{r, eval=F}
cdc = read.csv("https://marksmath.org/data/cdc.csv")
heights = subset(cdc, gender == 'm')$height
qqnorm(heights)
```

Looks good! If we want to model this data with a normal distribution, we'll need to know the mean and standard deviation.

```{r, eval=F}
m = mean(heights)
s = sd(heights)
c(m,s)
```

We can now ask the question - what percentage of men are taller than 6 feet (according to the model)? We can answer that question with the following computation.

```{r, eval=F}
1-pnorm(72, m, s)
```

\pagebreak

We can compare this to the actual value for this data set.

```{r, eval=F}
100 - quantile(heights, 0.72)
```

### Problem 3

Find an appropriate model for the heights of women in the CDC dataset and use it to estimate the proportion of women who are between 5' 2'' and 5'6'' tall.