confidence_intervals.qmd

# Confidence Intervals

```{r}
#| results: "asis"
#| echo: false

source("_common.R")
status("polishing")
```

## What is a confidence interval?

When we calculate a mean from a sample, we are using it to estimate the
mean of the population. Confidence intervals are a range of values and
are a way to quantify the uncertainty in our estimate. When we report a
mean with its 95% confidence interval we give the mean plus and minus
some variation. We are saying that 95% of the time, that range will
contain the population mean.

The confidence interval is calculated from the sample mean and the
standard error of the mean. The standard error of the mean is the
standard deviation of the [sampling distribution of the
mean](logic_hyopthesis_testing.html#sampling-distribution-of-the-mean).

To understand confidence intervals we need to understand some properties
of the normal distribution.

## The normal distribution

A [distribution](ideas_about_data.html#distributions) describes the
*values* the variable can take and the *chance* of them occurring. A
distribution has a general type, given by the function, and is further
tuned by the *parameters* in the function. For the normal distribution
these parameters are the mean and the standard deviation. Every variable
that follows the normal distribution has the same bell shaped curve and
the distributions differ only in their means and/or standard deviations.
The mean determines where the centre of the distribution is, the
standard deviation determines the spread (@fig-normal-distributions).

```{r}
#| echo: false
#| label: fig-normal-distributions
#| fig-cap: "The mean determines where the centre of the distribution is, the standard deviation determines the spread. The distributions on the left have the same mean but different standard deviations. The distributions on the right have the same standard deviation but different means."
mu <- 3300
sd <- 600
mu_shift <- 900
sd_shift <- 250

a <- ggplot(data = data.frame(mass = c(mu - 3 * sd, mu + 3 * sd)),
       aes(mass)) +
  stat_function(fun = dnorm, n = 101, 
                args = list(mean = mu, sd = sd)) +
  stat_function(fun = dnorm, n = 101, 
                args = list(mean = mu, sd = sd - sd_shift)) +
  stat_function(fun = dnorm, n = 101, 
                args = list(mean = mu, sd = sd + sd_shift)) +
  scale_y_continuous(breaks = NULL, name = "",
                     expand = c(0, 0),
                     limits = c(0, 0.0016)) +
  scale_x_continuous(breaks = NULL, name = "",
                     expand = c(0, 0),
                     limits = c(1000, 6000)) + 
  annotate("text", x = mu, y = 0.0013, 
           label = "Same mean\nDifferent standard deviations") +
  theme_classic() +
  theme(axis.line.y = element_blank())


b <- ggplot(data = data.frame(mass = c(mu - 3 * sd, mu + 3 * sd)),
       aes(mass)) +
  stat_function(fun = dnorm, n = 101, 
                args = list(mean = mu - mu_shift, sd = sd)) +
  stat_function(fun = dnorm, n = 101, 
                args = list(mean = mu, sd = sd)) +
  stat_function(fun = dnorm, n = 101, 
                args = list(mean = mu + mu_shift, sd = sd)) +
  scale_y_continuous(breaks = NULL, name = "",
                     expand = c(0, 0),
                     limits = c(0, 0.0016)) +
  scale_x_continuous(breaks = NULL, name = "",
                     expand = c(0, 0),
                     limits = c(1000, 6000)) + 
  annotate("text", x = mu, y = 0.0013, 
           label = "Same standard deviation\nDifferent means") +
  theme_classic() +
  theme(axis.line.y = element_blank())

a + b + plot_layout(ncol = 2)


```

Whilst normal distributions vary in the location on the horizontal axis
and their width, they all share some properties and it is these shared
properties that allow the calculation of confidence intervals with some
standard formulae. The properties are that a fix percentage of values
lie between a given number of standard deviations. For example, 68.2%
values lie between plus and minus one standard deviation from the mean
and 95% values lie between $\pm$ 1.96 standard deviations. Another way
of saying this is that there is a 95% chance that a randomly selected
value will lie between $\pm$ 1.96 standard deviations from the mean.
This is illustrated in @fig-norm-dist-means-properties.

```{r}
#| echo: false
#| label: fig-norm-dist-means-properties
#| fig-cap: "Normal distributions share some properties regardless of the mean and standard deviation. 68% of the values are within 1 standard deviation of the mean and 95% are within 1.96 standard deviations."


# functions for shading under a normal distribution curve
# 1 sd
dnorm_limit_1 <- function(x) {
  y <- dnorm(x)
  y[x < -1 |  x > 1] <- NA
  return(y)
}
# 1.96 sd
dnorm_limit_95 <- function(x) {
  y <- dnorm(x)
  y[x < -1.96 |  x > 1.96] <- NA
  return(y)
}

a <- ggplot(data = data.frame(x = c(-3, 3)),
       aes(x)) +
  stat_function(fun = dnorm_limit_1,
                geom = "area",
                fill = pal2[3]) +
  stat_function(fun = dnorm, n = 101) +
  annotate("text", x = 0, y = 0.2, 
           label = "68%",
           size = 6) +
  scale_y_continuous(breaks = NULL, name = "",
                     expand = c(0, 0)) +
  scale_x_continuous(breaks = c(-1, 0, 1), name = "",
                     expand = c(0, 0),
                     limits = c(-3, 3)) +
  theme_classic() +
  theme(axis.line.y = element_blank(),
        axis.ticks = element_blank())

b <- ggplot(data = data.frame(x = c(-3, 3)),
       aes(x)) +
  stat_function(fun = dnorm_limit_95,
                geom = "area",
                fill = pal2[3]) +
  stat_function(fun = dnorm, n = 101) +
  annotate("text", x = 0, y = 0.2, 
           label = "95%",
           size = 6) +
  scale_y_continuous(breaks = NULL, name = "",
                     expand = c(0, 0)) +
  scale_x_continuous(breaks = c(-1.96, 0, 1.96), name = "",
                     expand = c(0, 0),
                     limits = c(-3, 3)) +
  theme_classic() +
  theme(axis.line.y = element_blank(),
        axis.ticks = element_blank())

a + b + plot_layout(ncol = 1)

```

R has some useful functions associated with distributions, including the
normal distribution.

### Distributions: the R functions

For any distribution, R has four functions:

-   the density function, which gives the height of the function at a
    given value.
-   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 random number generating function

The functions are named with a letter `d`, `p`, `q` or `r` preceding the
distribution name. @tbl-rdistrib shows these four functions for the
normal, binomial, Poisson and t distributions.

| Distribution | Density    | Distribution | Quantile   | Random number generating |
|---------------|---------------|---------------|---------------|---------------|
| Normal       | `dnorm()`  | `pnorm()`    | `qnorm()`  | `rnorm()`                |
| Binomial     | `dbinom()` | `pbinom()`   | `qbinom()` | `rbinom()`               |
| Poisson      | `dpois()`  | `ppois()`    | `qpois()`  | `rpois()`                |
| t            | `dt()`     | `pt()`       | `qt()`     | `rt()`                   |

: R functions that provide values for some example distributions
{#tbl-rdistrib}

Searching for the manual with `?normal` or any one of the functions
(`?pnorm`) will bring up a single help page for all four associated
functions.

::: {#fig-manual-page}
<iframe style="height:400px;width:750px;border:2px" src="https://stat.ethz.ch/R-manual/R-devel/library/stats/html/Normal.html" title="Manual page">

</iframe>

The R Manual Page for the normal distribution which shows the four
functions associated functions.
:::

The functions which are of most use to us are `pnorm()` and `qnorm()`
and these are illustrated in @fig-pnorm.

```{r}
#| echo: false
#| label: fig-pnorm
#| fig-height: 8
#| fig-cap: "The `pnorm()` function calculates the probability that a value is less than or equal to a given value."


# functions for shading under a normal distribution curve

# 1.96 sd
dnorm_limit <- function(x) {
  y <- dnorm(x)
  y[x > 1.96] <- NA
  return(y)
}

a <- ggplot(data = data.frame(x = c(-3, 3)),
       aes(x)) +
  stat_function(fun = dnorm_limit,
                geom = "area",
                fill = pal2[3]) +
  stat_function(fun = dnorm, n = 101) +
  annotate("text", x = 0, y = 0.45, 
           label = "pnorm(x)",
           size = 6) +
  annotate("text", x = 1.96, y = 0.3, 
           label = "You give x",
           size = 5) +
  annotate("text", x = -2.3, y = 0.3, 
           label = "pnorm returns\nthe probability\n of getting ≤ x",
           size = 5) +
  annotate("segment", x = 1.96, y = 0.27, 
           xend = 1.96, yend = 0.08, 
           linewidth = 2,
           colour = pal4[1],
           arrow = arrow(type = "open", 
                         length = unit(0.05, "npc"))) +
  annotate("segment", x = -0.8, y = 0.05, 
           xend = -2.3, yend = 0.2, 
           linewidth = 2,
           colour = pal4[2],
           arrow = arrow(type = "open", 
                         length = unit(0.05, "npc"))) +
  scale_y_continuous(breaks = NULL, name = "",
                     expand = c(0, 0),
                     limits = c(0, 0.5)) +
  scale_x_continuous(breaks = c(1.96), 
                     labels = c("𝑥"),
                     name = "",
                     expand = c(0, 0),
                     limits = c(-3, 3)) +
  theme_classic() +
  theme(axis.line.y = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_text(size = 25))

b <- ggplot(data = data.frame(x = c(-3, 3)),
       aes(x)) +
  stat_function(fun = dnorm_limit,
                geom = "area",
                fill = pal2[3]) +
  stat_function(fun = dnorm, n = 101) +
  annotate("text", x = 0, y = 0.45, 
           label = "qnorm(P)",
           size = 6) +
  annotate("text", x = 1.96, y = 0.3, 
           label = "qnorm returns x",
           size = 5) +
  annotate("text", x = -2.3, y = 0.3, 
           label = "You give\nprobability",
           size = 5) +
  annotate("segment", xend = 1.96, yend = 0.27, 
           x = 1.96, y = 0.08, 
           linewidth = 2,
           colour = pal4[2],
           arrow = arrow(type = "open", 
                         length = unit(0.05, "npc"))) +
  annotate("segment", xend = -0.8, yend = 0.05, 
           x = -2.3, y = 0.2, 
           linewidth = 2,
           colour = pal4[1],
           arrow = arrow(type = "open", 
                         length = unit(0.05, "npc"))) +
  scale_y_continuous(breaks = NULL, name = "",
                     expand = c(0, 0),
                     limits = c(0, 0.5)) +
  scale_x_continuous(breaks = c(1.96), 
                     labels = c("𝑥"),
                     name = "",
                     expand = c(0, 0),
                     limits = c(-3, 3)) +
  theme_classic() +
  theme(axis.line.y = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_text(size = 25))

a + b + plot_layout(ncol = 1)

```

## Confidence intervals on large samples

$$
\bar{x} \pm 1.96 \times s.e.
$$ {#eq-ci-large-sample}

95% of confidence intervals calculated in this way will contain the true
population mean.

Do you have to remember the value of 1.96? Not if you have R!

```{r}
qnorm(0.975)
```

Notice that it is `qnorm(0.975)` and not `qnorm(0.95)` for a 95%
confidence interval. This is because the functions are defined as giving
the area under to the curve to the left of the value given. If we gave
0.95, we would get the value that put 0.05 in one tail. We want 0.025 in
each tail, so we need to use 0.975 in `qnorm()`.

TO-DO pic

## Confidence intervals on small samples

The calculation of confidence intervals on small samples is very similar
but we use the *t*-distribution rather than the normal distribution. The
formula is:

$$
\bar{x} \pm t_{[d.f.]} \times s.e.
$$ {#eq-ci-small-sample}

The *t*-distibution is a modified version of the normal distribution and
we use it because the sampling distribution of the mean is not quite
normal when the sample size is small. The *t*-distribution has an
additional parameter called the *degrees of freedom* which is the sample
size minus one ($n -1$). Like the normal distribution, the
*t*-distribution has a mean of zero and is symmetrical. However, The
*t*-distribution has fatter tails than the normal distribution and this
means that the probability of getting a value in the tails is higher
than for the normal distribution. The degrees of freedom determine how
much fatter the tails are. The smaller the sample size, the fatter the
tails. As the sample size increases, the *t*-distribution becomes more
and more like the normal distribution.

### What are degrees of freedom?

Degrees of freedom (usually abbreivated *d.f.*) describe how much 
independent information is available to estimate variability when 
calculating probabilities. In a population, *all* values are known, 
so no parameters need estimating. In a sample, we estimate parameters 
such as the mean. In estimating a mean, only $n - 1$ values can vary
independently because one value is constrained by the sample mean. 
In hypothesis tests and confidence intervals, *d.f.* adjust for 
uncertainty when generalizing from a sample to a population, affecting 
probability distributions like the *t*-distribution. More *d.f.* mean 
better probability estimates, while fewer *d.f.* reflect greater 
uncertainty in estimating population parameters.


## 🎬 Your turn!

If you want to code along you will need to start a new [RStudio
project](workflow_rstudio.html#rstudio-projects), add a `data-raw`
folder and open a new script. You will also need to load the
**`tidyverse`** package [@tidyverse].

## Large samples

A team of biomedical researchers is studying the concentration of
Creatine Kinase (CK) in the blood of patients with muscle disorders.
They collect a large random sample (n = 100) of blood samples from
patients and measure the CK concentration in UL^-1^ (units per litre).

The goal is to estimate the average CK concentration in this patient
population and calculate a 95% confidence interval for the mean
concentration. The data are in [ck_concen.csv](data-raw/ck_concen.csv) 

### Import

```{r}
creatinekinase <- read_csv("data-raw/ck_concen.csv")
```

### Calculate sample statistics

Calculate the mean, standard deviation, sample size and standard error
and assign as variables:

```{r}
# mean
m <- mean(creatinekinase$ck_conc)

# standard deviation
sd <- sd(creatinekinase$ck_conc)

# sample size (needed for the se)
n <- length(creatinekinase$ck_conc)

# standard error
se <- sd / sqrt(n)
```

### Calcuate *C.I.*

To calculate the 95% confidence interval we need to look up the quantile
(multiplier) using `qnorm()`:

```{r}
q <- qnorm(0.975)
```

Now we can use it in our confidence interval calculation:

```{r}
lcl <- m - q * se
ucl <- m + q * se
```

I used the names `lcl` and `ucl` to stand for "lower confidence limit"
and "upper confidence limit" respectively.

Print the values:

```{r}
lcl
ucl
```

This means we are 95% confident the population mean lies between
`r round(lcl,2)` UL^-1^ and `r round(ucl,2)` UL^-1^. The amount we have
added/subtracted from the mean (`q * se`) is `r round(q * se, 2)` thus
we sometimes see this written as `r round(m, 2)` $\pm$
`r round(q * se, 2)` mm.

### Report

The mean creatine kinase concentration in patients with muscle disorders
is `r round(m, 2)` UL^-1^, 95% C.I. \[`r round(lcl,2)`,
`r round(ucl,2)`\]

## Small samples

## Summary

## Summary of Confidence Intervals

1.  A confidence interval gives a range of plausible values for a 
    population mean from a sample and is calculated from a sample
    mean and standard error. 

2.  They are possible because all normal distributions have the same
    properties. 

3.  95% Confidence Intervals for Large Samples:

    -   Formula:\
        $$
        \bar{x} \pm 1.96 \times s.e.
        $$
    -   95% of confidence intervals computed this way will contain the
        true population mean.

4.  95% Confidence Intervals for Small Samples:

    -   Uses the *t*-distribution instead of the normal distribution.\
    -   Formula:\
        $$
        \bar{x} \pm t_{[d.f.]} \times s.e.
        $$
    -   The *t*-distribution has fatter tails and varies based on
        degrees of freedom (*n-1*), approaching the normal distribution
        as *n* increases.