How to use this sheet

This is not a fill-in-the-blanks sheet. All the code is here and all of it runs.

We assume you can already open RStudio, knit an .Rmd, read a data file, and index a data frame. If any of that is shaky, the two cheat sheets in this folder cover it, and it is a reasonable thing to ask an AI assistant to walk you through. What follows is the part that is specific to this course.

Everything below uses one seed, so your numbers will match the printed ones:

set.seed(67)

1. Random variables in R

1.1 The four-letter pattern

Every standard distribution in R comes with four functions, distinguished by their first letter. Learn the pattern once and you know all of them.

Prefix Gives you Normal Reads as
d density \(f(x)\) dnorm(x) height of the curve at x
p cdf \(F(x) = P(X \le x)\) pnorm(q) area to the left of q
q quantile \(F^{-1}(p)\) qnorm(p) the value with area p to its left
r random draws rnorm(n) n simulated observations


Note: p and q are inverses of each other.

pnorm(1.96)            # area to the left of 1.96
## [1] 0.9750021
qnorm(0.975)           # the value with 0.975 to its left
## [1] 1.959964

The distributions you will need this term, with their R names:

Distribution R name Parameters
Normal norm mean, sd: the standard deviation, not the variance
Exponential exp rate
Chi-square chisq df
Student \(t\) t df
Fisher \(F\) f df1, df2
Uniform unif min, max
Binomial binom size, prob


Two examples of the kind of thing you will do all term:

qt(0.975, df = 23)             # two-sided 5% critical value for a t with 23 df
## [1] 2.068658
1 - pf(4.2, df1 = 2, df2 = 40) # p-value of an observed F of 4.2
## [1] 0.02209493

Be careful: rnorm(n, mean, sd) takes the standard deviation. To sample from \(N(10, 4)\), whose variance is 4, you write rnorm(n, 10, 2).

1.2 Repeating an experiment

The whole idea of a sampling distribution is: the statistic is a random variable, so compute it many times and look at the spread of the answers. replicate is how you do that.

replicate(5, mean(rnorm(20, 10, 2)))   # five separate experiments
## [1]  9.644336 10.072767  9.526072  9.942828 10.039465

Each entry is the mean of its own fresh sample of 20. They differ because the sample differs, and that variation corresponds to the sampling distribution.

Three other functions worth knowing:

set.seed(67)
sample(1:10, size = 4)                    # draw without replacement
## [1] 4 1 6 7
sample(c("a","b"), 6, replace = TRUE)     # with replacement
## [1] "a" "a" "b" "a" "a" "b"
round(c(mean = mean(rnorm(100)), sd = sd(rnorm(100))), 3)
##   mean     sd 
## -0.049  0.998

set.seed makes randomness reproducible. Set it once at the top of every analysis.

2. ggplot2 for this course

Every ggplot is the same three pieces: data, an aesthetic mapping (aes, which says which variable goes on which axis), and one or more geoms added with +.

d <- data.frame(z = rnorm(500))
ggplot(d, aes(x = z)) + geom_histogram(bins = 30)

2.1 Density scale, and a curve on top

Comparing a histogram with a theoretical density only works if both are on the same vertical scale. Counts are not; densities are. This idiom appears in all four questions below:

ggplot(d, aes(x = z)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30,
                 fill = "grey80", colour = "white") +
  stat_function(fun = dnorm, args = list(mean = 0, sd = 1),
                colour = "#B3261E", linewidth = 0.8) +
  labs(x = "z", y = "density")

after_stat(density) rescales the bars so their total area is 1; stat_function draws any function of x, here the \(N(0,1)\) density.

2.2 Several panels at once

To compare the same plot under different settings, such as different sample sizes, stack all the simulated values into one data frame, with a column saying which setting each value belongs to. Then facet_wrap() draws one panel per setting, side by side. This lets you watch a distribution change as \(n\) grows, which is exactly what the central limit theorem describes.

dd <- rbind(data.frame(n = "n = 1",  m = replicate(2000, mean(rexp(1)))),
            data.frame(n = "n = 30", m = replicate(2000, mean(rexp(30)))))
ggplot(dd, aes(m)) +
  geom_histogram(bins = 30, fill = "grey80", colour = "white") +
  facet_wrap(~ n, scales = "free") + labs(x = expression(bar(X)), y = "count")

2.3 QQ plots and scatterplots with a fitted line

A histogram shows normality loosely. A QQ plot is a bit more precise. Points on the line means the sample matches the theoretical quantiles.

ggplot(d, aes(sample = z)) + stat_qq(size = 0.6) + stat_qq_line(colour = "#B3261E") +
  labs(x = "theoretical quantiles", y = "sample quantiles")

And the plot this whole course is about:

xx <- runif(60, 0, 10); yy <- 3 + 2 * xx + rnorm(60, 0, 2)
ggplot(data.frame(xx, yy), aes(xx, yy)) +
  geom_point(size = 0.9) +
  geom_smooth(method = "lm", se = FALSE, colour = "#1F5FA9", linewidth = 0.7) +
  labs(x = "x", y = "y")

Useful: labs() for axis titles, facet_wrap() for panels, theme_minimal() for a clean look, expression(bar(X)) for maths in a label.

3. Four simulations

Each question states a mathematical result and shows code that tests it. The results are Weeks 2 to 6 of this course, so you are seeing here the pictures before the proofs.

Question 1: the sampling distribution of \(\bar X\)

Let \(X_1,\dots,X_n\) be i.i.d. \(N(\mu,\sigma^2)\). Then exactly \[\bar X \sim N\!\left(\mu,\ \frac{\sigma^2}{n}\right).\]

Take \(\mu = 10\), \(\sigma^2 = 4\), \(n = 20\):

mu <- 10; sigma <- 2; n <- 20
means <- replicate(5000, mean(rnorm(n, mu, sigma)))

c(simulated_mean = mean(means), theoretical_mean = mu,
  simulated_sd   = sd(means),   theoretical_sd   = sigma / sqrt(n))
##   simulated_mean theoretical_mean     simulated_sd   theoretical_sd 
##       10.0046918       10.0000000        0.4522440        0.4472136
ggplot(data.frame(means), aes(means)) +
  geom_histogram(aes(y = after_stat(density)), bins = 40,
                 fill = "grey80", colour = "white") +
  stat_function(fun = dnorm, args = list(mean = mu, sd = sigma / sqrt(n)),
                colour = "#B3261E", linewidth = 0.8) +
  labs(x = expression(bar(X)), y = "density")

To discuss.

  1. The simulated standard deviation is close to \(\sigma/\sqrt{n} = 0.447\), not to \(\sigma = 2\). Where does the \(\sqrt{n}\) come from?
  2. What would change in the picture if \(n\) were 80 instead of 20? Predict first, then edit the code.
  3. Why is this result exact here, with no appeal to a limit theorem?

Question 2: the central limit theorem

Question 1 was exact because the population was normal. The CLT says \(\bar X\) is approximately normal for large \(n\) whatever the population, provided the variance is finite.

The exponential with rate 1 is strongly right-skewed, which makes it a hard test:

sizes <- c(1, 5, 30)
clt <- do.call(rbind, lapply(sizes, function(k)
  data.frame(n = factor(paste("n =", k), levels = paste("n =", sizes)),
             m = replicate(5000, mean(rexp(k, rate = 1))))))

ggplot(clt, aes(m)) +
  geom_histogram(aes(y = after_stat(density)), bins = 40,
                 fill = "grey80", colour = "white") +
  facet_wrap(~ n, scales = "free") +
  labs(x = expression(bar(X)), y = "density")

ggplot(subset(clt, n == "n = 30"), aes(sample = m)) +
  stat_qq(size = 0.5) + stat_qq_line(colour = "#B3261E") +
  labs(title = "n = 30", x = "theoretical", y = "sample")

To discuss.

  1. At \(n = 1\) the histogram is the population. Describe how the shape changes as \(n\) grows, in words.
  2. The QQ plot at \(n = 30\) bends at the right. What does that bend tell you that the histogram hides?
  3. The centre of every panel is 1. Why does the centre not move, when the shape does?

Question 3: where the chi-square comes from

Let \(X_1,\dots,X_k\) be independent \(N(0,1)\) and set \(Y_k=\sum_{i=1}^k X_i^2\). Then \(Y_k\sim\chi^2_k\), with mean \(k\) and variance \(2k\). Since \(Y_k\) is itself a sum of independent terms, the CLT applies to it: \[Z_k=\frac{Y_k-k}{\sqrt{2k}}\ \longrightarrow\ N(0,1).\]

ks <- c(1, 5, 50)
chi <- do.call(rbind, lapply(ks, function(k) {
  Y <- replicate(5000, sum(rnorm(k)^2))
  data.frame(k = factor(paste("k =", k), levels = paste("k =", ks)),
             z = (Y - k) / sqrt(2 * k)) }))

ggplot(chi, aes(z)) +
  geom_histogram(aes(y = after_stat(density)), bins = 40,
                 fill = "grey80", colour = "white") +
  stat_function(fun = dnorm, colour = "#B3261E", linewidth = 0.7) +
  facet_wrap(~ k) + coord_cartesian(xlim = c(-3, 4)) +
  labs(x = expression(Z[k]), y = "density")

To discuss.

  1. At \(k=1\), \(Z_k\) is a squared normal, shifted and scaled. Why can it never go below \(-1/\sqrt2 \approx -0.71\)?
  2. By \(k=50\) the approximation is convincing. Is that fast or slow, compared with Question 2?

Question 4: the sampling distribution of \(\hat\beta_1\)

This course is, in particular, concerned with the study of one random variable: the estimated slope. You can watch it behave before deriving anything.

Generate from a known model, fit it, and repeat:

beta0 <- 100; beta1 <- 20; sigma <- 5; n <- 25
x <- seq(1, 10, length.out = n)

one_slope <- function() {
  y <- beta0 + beta1 * x + rnorm(n, mean = 0, sd = sigma)
  coef(lm(y ~ x))[2]
}

slopes <- replicate(2000, one_slope())
Sxx <- sum((x - mean(x))^2)

c(mean_slope  = mean(slopes), true_beta1 = beta1,
  sd_slope    = sd(slopes),   theoretical = sigma / sqrt(Sxx))
##  mean_slope  true_beta1    sd_slope theoretical 
##  20.0054222  20.0000000   0.3733863   0.3698001
ggplot(data.frame(slopes), aes(slopes)) +
  geom_histogram(aes(y = after_stat(density)), bins = 40,
                 fill = "grey80", colour = "white") +
  stat_function(fun = dnorm, args = list(mean = beta1, sd = sigma / sqrt(Sxx)),
                colour = "#B3261E", linewidth = 0.8) +
  labs(x = expression(hat(beta)[1]), y = "density")

To discuss.

  1. x, beta1 and sigma never change between runs. So why is one_slope() different every time?
  2. The average of the 2,000 slopes is close to 20. State that property in the language of estimators.
  3. The spread matches \(\sigma/\sqrt{S_{xx}}\). What does that say about the \(x\) values you should choose, if you get to design the experiment?
  4. The histogram is symmetric and bell-shaped. Which assumption of the model is doing that work, and would it survive if the errors were not normal?