An interactive textbook

Weird Statistics Questions

Twenty-one questions that students rarely ask out loud, with simulations you can break with your own hands.

Hai Liang · School of Journalism and Communication, The Chinese University of Hong Kong · hailiang@cuhk.edu.hk

Preface

Most statistics courses teach procedures: which test to run, which function to call, and which p-value to report. However, students often get stuck on questions that the procedures do not answer. Why do researchers report the mean and the standard deviation, rather than something else? Why can a sample of 1,000 people say anything about a city of 7.5 million? Why can adding one more control variable make an estimate worse?

I call these weird questions, not because they are strange, but because they sit between the formulas and are rarely discussed. This book is organized around 21 of them. Each chapter starts with a playful scenario, explains the idea in words and in equations, and then lets you run the simulation yourself. In other words, you do not have to take my word for anything: you can generate the data, violate the assumptions, and watch what happens.

How each chapter works

The hook: a real-life scenario The weird question Play: an interactive simulation Common traps Takeaways

Every chapter also includes an R Corner, which contains R code that reproduces the analysis, and a short Check yourself quiz. The R Corner is written in R Markdown style: each piece of code is followed by its actual output, which was produced in R 4.6.1, so you can compare your results with mine line by line. The same code and results are collected in a companion file, Weird-Statistics-R-Companion.Rmd, and its knitted version, Weird-Statistics-R-Companion.md. The simulations in the browser use a fixed random seed, so your numbers will be the same each time you reload. Words with a dotted underline, such as standard error, show a definition when you hover over them.

How to cite this book

If you use this book in teaching or research, please cite it as follows (APA, 7th edition):

Liang, H. (2026). Weird statistics questions: An interactive textbook (Version 1.0) [Interactive online textbook]. School of Journalism and Communication, The Chinese University of Hong Kong.

BibTeX:

@book{liang2026weird,
  author    = {Liang, Hai},
  title     = {Weird Statistics Questions: An Interactive Textbook},
  year      = {2026},
  edition   = {Version 1.0},
  publisher = {School of Journalism and Communication, The Chinese University of Hong Kong},
  note      = {Interactive online textbook},
  url       = {}
}

To cite a single chapter, use the citation at the end of that chapter.

Declaration of AI assistance

This book is based on my teaching notes, Weird Statistics Questions (Version 2), which were written in R Markdown. I used an AI assistant, Claude (Anthropic), working inside the Cursor code editor, to turn these notes into the present interactive book. More specifically, the assistant (1) drafted the chapter text from my notes, following my instructions on structure and writing style; (2) created the teaching scenarios and programmed the interactive simulations in JavaScript; (3) checked the statistical claims in the original notes and proposed the corrections that appear in the boxes titled "Note on the original notes"; and (4) ran every R Corner in R 4.6.1, so that the results shown are the actual output of the code.

I designed the book, decided which corrections to adopt, and reviewed the content. I take full responsibility for the book, including any remaining errors. The scenarios and simulated data are fictional and serve teaching purposes only, whereas the real data sets (e.g., sleep, npk, mtcars, warpbreaks, and the UCLA admissions and absence data) are publicly available.

Part I · Describing data · Chapter 1

The Mean, the Standard Deviation, and Their Limits

The hook

The bubble-tea queue

Imagine you run a bubble-tea shop near campus. Every afternoon a queue forms outside the door, and every cup gets a scoop of toppings. Your manager asks two simple questions: how long does a typical customer wait, and how many customers arrive per minute? You could answer with the mean. However, a customer who waited 40 minutes behind a group order of 30 drinks will not be comforted by the average. This chapter examines what the mean and the standard deviation can and cannot tell you.

The weird question

How can we represent and summarize a variable? Why do we calculate the mean and the standard deviation? What if the variable is skewed or discrete?

A distribution is the data, summarized

The most faithful way to represent a variable is its distribution, usually drawn as a histogram, because a histogram stays close to the raw data. Summary statistics compress the histogram into a few numbers. The mean ($\mu$) describes the central tendency, i.e., the typical value, and the standard deviation ($\sigma$) describes the dispersion, i.e., how far values typically spread around the mean. Quantiles, such as the median and the quartiles, describe locations within the distribution.

Three notes are important. First, the mean is only one measure of central tendency. The median and the mode are two others, and the appropriate choice depends on the research question and the shape of the distribution. Second, the standard deviation is defined around the mean:

$$s=\sqrt{\frac{1}{n-1}\sum_{i=1}^{n}\left(x_i-\bar{x}\right)^2}$$

Therefore, if the mean is not meaningful, the standard deviation is not meaningful either. Third, if a mathematical formula describes the distribution, the parameters of that formula are the most precise summary.

When the formula knows best: parameters

A Normal distribution (the bell curve) is written as $X\sim\mathsf{Norm}(\mu,\sigma)$. Its two parameters are exactly the mean and the standard deviation, which explains why these two numbers became the default summary of data. A Poisson distribution, which describes counts such as the number of customers arriving per minute, is written as $X\sim\mathsf{Pois}(\lambda)$:

$$P(X=k)=\frac{\lambda^{k}e^{-\lambda}}{k!}$$

Here $\lambda$, the rate parameter, is the total number of events divided by the number of units. For example, if 500 customers arrive in 200 minutes, $\lambda=2.5$ customers per minute. A remarkable property of the Poisson distribution is that its mean and variance are equal: $E[X]=\mathrm{Var}(X)=\lambda$. You can verify this property in the explorer below.

Play · Distribution explorer

Distribution

Try this: with Poisson and $\lambda=2.5$, compare the variance with the mean. Then move $\lambda$ and check that the two stay close.
Switch to Pareto with $\alpha=1.2$ and press New sample several times. Watch how much the standard deviation jumps.
Move the Chebyshev $K$ slider. The share of data within $K$ standard deviations never falls below the Chebyshev minimum, whatever the shape.

Think about it

Many textbooks advise that a skewed distribution should be summarized by the median. In the explorer, a Poisson variable with $\lambda=2.5$ has a median of 2. Is the median more informative than the mean here?

My answer

Not here. Every Poisson distribution with $\lambda$ between about 1.7 and 2.7 has a median of 2, so the median cannot tell these distributions apart. The mean, however, estimates $\lambda$, which is the only parameter of the distribution. Once I know $\lambda$, I can recover the entire distribution from the formula. Therefore, for a Poisson variable, the mean is a complete summary, whereas the median throws information away. Skewness alone is not a reason to abandon the mean; the more important question is whether the mean corresponds to a meaningful parameter.

Comparing spread across groups. Comparing two standard deviations is straightforward when the variables share a scale and have similar means. When the means differ greatly, for example, the body weights of mice and elephants, the coefficient of variation, $CV=\sigma/\mu$, which means the standard deviation expressed as a fraction of the mean, is more informative. Note that the CV only makes sense for ratio variables with a true zero, such as weight or income.

Categories have no mean

For categorical variables, the mean, and thus the standard deviation, are not well defined. If a survey codes districts or ethnic groups as numbers, software will happily report that the average ethnic group in Hong Kong is 3.1, but the number has no interpretation. The distribution of a categorical variable is described by its category probabilities, $P(X=k)=p_k$.

To see why this matters, suppose two shops code their toppings as 8 = pearls, 9 = pudding, 10 = grass jelly, 11 = aloe, and 12 = red bean. Their customers have very different tastes, but the summary statistics of the codes can look almost identical.

Play · The twin shops

Shop A is fixed. Use the sliders to change the topping sales at Shop B (the shares are rescaled to sum to 100%). Your challenge: find a menu for Shop B with the same mean and standard deviation as Shop A, but a different distribution.

The best summary for a categorical variable is the frequency table, i.e., the proportion of observations in each category. In the twin-shop example, the frequency tables reveal at once how differently the customers of the two shops behave, even when the mean and the standard deviation of the codes are identical.

When the mean does not exist

Some continuous distributions have no finite mean or no finite variance. The Pareto distribution, which describes wealth, city sizes, and the number of followers on social media, has a tail index $\alpha$. If $\alpha\le 1$, the mean is infinite; if $\alpha\le 2$, the variance is infinite. In practice, this means that the sample mean never settles down: one new observation, such as a billionaire walking into a bar, can move it dramatically. (See this tutorial on fat-tailed distributions.)

Play · A billionaire walks into a bar

Track

Try this: set $\alpha=3$. The running mean quickly settles near the dashed theoretical line.
Set $\alpha=1.5$ and track the standard deviation. The mean settles, but the standard deviation keeps jumping, because the variance is infinite.
Set $\alpha=0.9$. Now even the mean never settles.

Expected value: why the mean still matters

If many distributions lack a meaningful mean, why are the mean and variance so popular? The answer lies in the expected value, which means the probability-weighted average of all possible values:

$$E[X]=\sum_{k} x_k\,P(X=x_k)\quad\text{(discrete)},\qquad E[X]=\int_{-\infty}^{\infty}x\,f(x)\,dx\quad\text{(continuous)}$$

The expected value is a property of the population, or more precisely of the process that generates the data, whereas the mean is a property of a sample. In a random sample, the sample mean is an unbiased estimate of $E[X]$. In a biased sample, the two can differ systematically. For example, the expected topping code at Shop B (with the reset menu) is $8(0.10)+9(0.25)+10(0.30)+11(0.10)+12(0.25)=10.15$. A random sample of 1,000 cups gives a mean close to, but not exactly, 10.15.

How many numbers are needed to recover a distribution? It depends on the distribution. A Poisson distribution needs one number ($\lambda$), a Normal distribution needs two ($\mu$ and $\sigma$), and a categorical variable with $k$ categories needs $k-1$ probabilities. A binary variable coded 0 and 1 needs only the proportion of 1s, $p$, because $E[X]=p$ and $\mathrm{Var}(X)=E[(X-E[X])^2]=p(1-p)$. Therefore, to describe data precisely, I first ask which distribution generated them, and then I estimate its parameters.

The mean and the variance are the first two members of a family of summaries called moments. The first moment is the expected value; the second central moment is the variance, $\mathrm{Var}(X)=E[(X-E[X])^2]$; the third standardized moment is skewness; and the fourth standardized moment is kurtosis. Even when $E[X]$ has no direct substantive interpretation, it remains essential for linking a sample to a population, as Chapter 3 shows.

A guarantee for any shape: Chebyshev's inequality

The mean and the standard deviation still say something useful about skewed distributions. Chebyshev's inequality states that, for any distribution with a finite variance and any $K>1$, at least $1-1/K^2$ of the values lie within $K$ standard deviations of the mean:

$$P\left(|X-\mu|<K\sigma\right)\ \ge\ 1-\frac{1}{K^2}$$

$K$	Chebyshev minimum (any shape)	Normal distribution
2	75.0%	95.4%
3	88.9%	99.7%
4	93.75%	99.99%

Chebyshev's bound is weaker than the familiar Normal figures because it holds for every shape. You can check it in the distribution explorer above by moving the $K$ slider.

Common traps

Averaging category codes, such as a mean topping of 10.0 or a mean ethnic group of 3.1.
Switching to the median automatically for skewed data, even when the mean corresponds to a parameter, such as the Poisson $\lambda$.
Reporting a mean and a standard deviation for fat-tailed variables, such as followers or wealth, without checking whether they are stable.

R Corner: code and results

set.seed(12345)
x_norm = rnorm(1000, 10, 2)    # Normal: mean = 10, sd = 2
x_pois = rpois(1000, 2.5)      # Poisson: lambda = 2.5
hist(x_norm); hist(x_pois)
c(mean(x_norm), mean(x_pois))
c(var(x_norm), var(x_pois))    # Poisson: the variance is close to the mean
median(x_pois)                 # 2: the same for any lambda from about 1.7 to 2.7

# The twin shops: different distributions, similar means and SDs
x_cat1 = sample(8:12, 1000, replace = TRUE, prob = c(0.25, 0.10, 0.30, 0.25, 0.10))
x_cat2 = sample(8:12, 1000, replace = TRUE, prob = c(0.10, 0.25, 0.30, 0.10, 0.25))
c(mean(x_cat1), mean(x_cat2)); c(sd(x_cat1), sd(x_cat2))
prop.table(table(x_cat1)); prop.table(table(x_cat2))   # the better summary

# Expected value from the true probabilities
sum(8:12 * c(0.10, 0.25, 0.30, 0.10, 0.25))   # 10.15

# A fat tail: Pareto with alpha = 1.5 and x_m = 1
x_par = 1 / runif(20000)^(1 / 1.5)
plot(cumsum(x_par) / seq_along(x_par), type = "l", ylab = "Running mean")

# Chebyshev: share of values within K = 2 SDs of the mean
mean(abs(x_pois - mean(x_pois)) < 2 * sd(x_pois))   # at least 0.75

set.seed(12345)
x_norm = rnorm(1000, 10, 2)    # Normal: mean = 10, sd = 2
x_pois = rpois(1000, 2.5)      # Poisson: lambda = 2.5
hist(x_norm); hist(x_pois)

c(mean(x_norm), mean(x_pois))

#> [1] 10.0924  2.4390

c(var(x_norm), var(x_pois))    # Poisson: the variance is close to the mean

#> [1] 3.989987 2.480760

median(x_pois)                 # 2: the same for any lambda from about 1.7 to 2.7

#> [1] 2

# The twin shops: different distributions, similar means and SDs
x_cat1 = sample(8:12, 1000, replace = TRUE, prob = c(0.25, 0.10, 0.30, 0.25, 0.10))
x_cat2 = sample(8:12, 1000, replace = TRUE, prob = c(0.10, 0.25, 0.30, 0.10, 0.25))
c(mean(x_cat1), mean(x_cat2)); c(sd(x_cat1), sd(x_cat2))

#> [1]  9.857 10.195

#> [1] 1.303937 1.309462

prop.table(table(x_cat1)); prop.table(table(x_cat2))   # the better summary

#> x_cat1
#>     8     9    10    11    12 
#> 0.246 0.092 0.322 0.239 0.101

#> x_cat2
#>     8     9    10    11    12 
#> 0.088 0.252 0.299 0.099 0.262

# Expected value from the true probabilities
sum(8:12 * c(0.10, 0.25, 0.30, 0.10, 0.25))   # 10.15

#> [1] 10.15

# A fat tail: Pareto with alpha = 1.5 and x_m = 1
x_par = 1 / runif(20000)^(1 / 1.5)
plot(cumsum(x_par) / seq_along(x_par), type = "l", ylab = "Running mean")

# Chebyshev: share of values within K = 2 SDs of the mean
mean(abs(x_pois - mean(x_pois)) < 2 * sd(x_pois))   # at least 0.75

#> [1] 0.969

Takeaways

A distribution is the data; the mean and the standard deviation are compressions of it.
When the generating distribution is known, its parameters are the best summary: $\lambda$ for Poisson, $\mu$ and $\sigma$ for Normal.
Category codes do not have meaningful means. Report frequency tables instead.
Fat tails ($\alpha\le 2$) make the standard deviation, and sometimes even the mean, unstable.
For any shape, at least 75% of values lie within two standard deviations of the mean.

Part I · Describing data · Chapter 2

Many Correlations, and Association Between Categories

The hook

Scrolling before bed

A class survey finds that students who scroll short videos longer before bed also sleep less. One classmate reports Pearson's $r$, another reports Spearman's $\rho$, and a third reports Kendall's $\tau$. The three numbers differ, and nobody made a mistake. In other words, why does statistics need more than one correlation coefficient, and what should a researcher do when the variables are categories, such as being a cat person or a dog person?

The weird question

We use correlation coefficients (why plural?) to quantify the relationship between two continuous variables. How can we measure relationships involving categorical or rank variables?

Pearson's $r$: linear association

Pearson's correlation coefficient measures the linear association between two continuous variables:

$$r=\frac{\mathrm{Cov}(X,Y)}{\sigma_X\,\sigma_Y},\qquad \mathrm{Cov}(X,Y)=E\big[(X-E[X])(Y-E[Y])\big]=E[XY]-E[X]E[Y]$$

Consider $X\sim\mathsf{Norm}(10,2)$ and $Y=0.5X+\varepsilon$, where $\varepsilon\sim\mathsf{Norm}(0,1)$. Many students expect $r=0.5$, but $r$ is about .71. Why? The covariance is $0.5\times 4=2$, and the variance of $Y$ is $0.25\times 4+1=2$. Therefore, $r=2/(2\sqrt{2})\approx .71$ and $r^2=.50$. In other words, the slope 0.5 is not the correlation; $r^2=.50$ means that half of the variance in $Y$ is shared with $X$.

Pearson's $r$ does not assume normality, although it assumes finite variances. Normality matters for the significance test. When the two variables are bivariate Normal, $t=r\sqrt{(n-2)/(1-r^2)}$ follows a t distribution with $n-2$ degrees of freedom. With small samples from non-normal data, this p-value can be misleading; Chapter 6 shows how bootstrapping avoids the assumption.

Spearman's $\rho$: monotonic association

Spearman's $\rho$ is simply Pearson's $r$ computed on ranks. It measures monotonic association, i.e., whether $Y$ tends to increase (or decrease) as $X$ increases, regardless of whether the relationship follows a straight line. Because ranks ignore distances, $\rho$ is robust to outliers and appropriate for ordinal data. When variables are highly skewed, researchers often log-transform them; ranking has a similar effect without requiring a choice of transformation.

Kendall's $\tau$: a probability

Kendall's $\tau$ works with pairs of observations. A pair is concordant if the observation with the larger $X$ also has the larger $Y$, and discordant otherwise. With $C$ concordant and $D$ discordant pairs among $n$ observations:

$$\tau=\frac{C-D}{n(n-1)/2}$$

Therefore, $\tau$ is the probability that a random pair is concordant minus the probability that it is discordant. Because it is a probability rather than a proportion of variance, $\tau$ is not on the same scale as $r$ or $\rho$. For bivariate Normal data, $\tau=\frac{2}{\pi}\arcsin(r)$, so $r=.71$ corresponds to $\tau\approx .50$. A smaller $\tau$ does not indicate a weaker relationship; it indicates a different metric. Spearman's $\rho$ is usually larger than $\tau$ for the same data, which sometimes leads researchers to the incorrect conclusion that $\rho$ is "better."

Play · The correlation lab

Start from:

Drag any point. Click an empty spot to add a point; double-click a point to remove it.

Try this: choose One outlier. Pearson's $r$ is inflated by a single point, whereas $\rho$ and $\tau$ are much less affected. Drag the outlier back into the cloud.
Choose Curved but monotonic: $\rho$ and $\tau$ are close to 1, but $r$ is not.
Choose Inverted U: all three coefficients are close to 0, although $Y$ depends strongly on $X$.

Pearson versus Spearman when the curve bends

If $Y$ is a monotonic but non-linear function of $X$, $\rho$ stays at 1 while $r$ declines. The widget below correlates $x\sim\mathsf{Norm}(\text{mean},2)$ with $x^n$ for $n=1,\dots,10$.

Play · Correlating $x$ with $x^n$

Try this: with a mean of 10, all $x$ values are positive, so every power is monotonic: $\rho=1$ throughout, while $r$ slowly declines.
Move the mean to 0. Now even powers ($x^2, x^4,\dots$) are U-shaped, and both coefficients collapse to about 0 for even $n$.

Cosine similarity. Text analysts often use cosine similarity between two vectors $A$ and $B$:

$$\mathrm{CosSim}(A,B)=\frac{\sum_i A_iB_i}{\sqrt{\sum_i A_i^2}\,\sqrt{\sum_i B_i^2}}$$

Cosine similarity is Pearson's $r$ without centering. When both variables are centered at zero, the two are identical. When all values are positive, as with word counts, cosine similarity is high even for unrelated variables. You can verify this with the No relation preset in the correlation lab.

Association between two categorical variables

Suppose I survey 170 residents in two cities, B and T, and ask whether they are cat people or dog people. For two categorical variables, the $\chi^2$ statistic compares the observed counts $O$ with the counts $E$ expected under independence:

$$\chi^2=\sum\frac{(O-E)^2}{E}$$

However, $\chi^2$ is not normalized: it grows with the sample size, so it measures evidence rather than strength. In a $k\times l$ table, its maximum is $N(\min(k,l)-1)$, which leads to Cramér's $V$:

$$V=\sqrt{\frac{\chi^2}{N\,(\min(k,l)-1)}}$$

$V$ ranges from 0 to 1. The uncertainty coefficient (Theil's $U$) takes an information-theory perspective. It is the proportion of uncertainty (entropy) in one variable that is removed by knowing the other, so it can differ by direction:

$$U(C\mid R)=\frac{H(C)-H(C\mid R)}{H(C)}$$

Note: R's chisq.test applies Yates's continuity correction to 2 × 2 tables by default, which affects the p-value. Cramér's $V$ below uses the uncorrected $\chi^2$, so packages may differ in the second decimal place. See other measures of nominal association.

Play · Cat people and dog people

	Cat person	Dog person
City B
City T

Try this: press Multiply every count by 10. $\chi^2$ becomes 10 times larger and the p-value shrinks, but $V$ and $U$ do not change.
Press Swap the two columns. The phi coefficient and the odds ratio flip direction; $V$ and $U$ stay the same.

Think about it

Is it possible for the association between two categorical variables to be negative?

My answer

For nominal variables, no. Cramér's $V$ and Theil's $U$ are never negative, because the order of categories (cat people first or dog people first) is arbitrary and should not change the answer. In a 2 × 2 table, the phi coefficient and the log odds ratio do carry a sign, but the sign flips when the categories of either variable are swapped, as the widget shows. Thus, a sign is substantively meaningful only when categories are ordered, as in ordinal variables (low, medium, high), for which measures such as Goodman and Kruskal's $\gamma$ or Kendall's $\tau_b$ are appropriate.

Common traps

Concluding that a relationship is weak because $\tau$ is smaller than $\rho$. The two are on different scales.
Reading $r\approx 0$ as "no relationship." It means no linear relationship.
Reporting $\chi^2$ as the strength of an association. Use $V$ or $U$, and report $\chi^2$ as the test.

R Corner: code and results

set.seed(12345)
rv_1 = rnorm(1000, 10, 2)
rv_2 = 0.5 * rv_1 + rnorm(1000)
cor(rv_1, rv_2); cor(rv_1, rv_2)^2        # about .71 and .50

rv_3 = rpois(1000, 2.5)
rv_4 = 0.5 * rv_3 + rnorm(1000)
cor.test(rv_3, rv_4, method = "pearson")
cor.test(rv_3, rv_4, method = "spearman", exact = FALSE)   # ties in Poisson data
cor.test(rv_3, rv_4, method = "kendall", exact = FALSE)

cor(rank(rv_1), rank(rv_2))                # Pearson on ranks = Spearman
cor(rv_1, rv_2, method = "spearman")

# x versus x^n
sapply(1:10, function(i) cor(rv_1, rv_1^i))                        # Pearson
sapply(1:10, function(i) cor(rv_1, rv_1^i, method = "spearman"))   # Spearman

# Non-monotonic: all three coefficients understate the dependence
rv_x = -0.5 * rv_1^2 + 11 * rv_1 + rnorm(1000)
cor(rv_1, rv_x, method = "pearson"); cor(rv_1, rv_x, method = "kendall")

# Cosine similarity (Pearson without centering)
sum(rv_1 * rv_x) / (sqrt(sum(rv_1^2)) * sqrt(sum(rv_x^2)))

# Two categorical variables
tbl = matrix(c(65, 45, 20, 40), nrow = 2, byrow = TRUE,
             dimnames = list(City = c("B", "T"), Pet = c("Cat", "Dog")))
chisq.test(tbl)                            # Yates-corrected for 2 x 2 tables
sqrt(chisq.test(tbl, correct = FALSE)$statistic / sum(tbl))   # Cramer's V
library(DescTools)
UncertCoef(tbl, direction = "symmetric")

set.seed(12345)
rv_1 = rnorm(1000, 10, 2)
rv_2 = 0.5 * rv_1 + rnorm(1000)
cor(rv_1, rv_2); cor(rv_1, rv_2)^2        # about .71 and .50

#> [1] 0.7174864

#> [1] 0.5147867

rv_3 = rpois(1000, 2.5)
rv_4 = 0.5 * rv_3 + rnorm(1000)
cor.test(rv_3, rv_4, method = "pearson")

#> 
#> 	Pearson's product-moment correlation
#> 
#> data:  rv_3 and rv_4
#> t = 25.737, df = 998, p-value < 2.2e-16
#> alternative hypothesis: true correlation is not equal to 0
#> 95 percent confidence interval:
#>  0.5928380 0.6674754
#> sample estimates:
#>       cor 
#> 0.6316179

cor.test(rv_3, rv_4, method = "spearman", exact = FALSE)   # ties in Poisson data

#> 
#> 	Spearman's rank correlation rho
#> 
#> data:  rv_3 and rv_4
#> S = 62485885, p-value < 2.2e-16
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#>       rho 
#> 0.6250843

cor.test(rv_3, rv_4, method = "kendall", exact = FALSE)

#> 
#> 	Kendall's rank correlation tau
#> 
#> data:  rv_3 and rv_4
#> z = 20.85, p-value < 2.2e-16
#> alternative hypothesis: true tau is not equal to 0
#> sample estimates:
#>      tau 
#> 0.478622

cor(rank(rv_1), rank(rv_2))                # Pearson on ranks = Spearman

#> [1] 0.6989149

cor(rv_1, rv_2, method = "spearman")

#> [1] 0.6989149

# x versus x^n
sapply(1:10, function(i) cor(rv_1, rv_1^i))                        # Pearson

#>  [1] 1.0000000 0.9904631 0.9650112 0.9281316 0.8836453 0.8345976 0.7833095
#>  [8] 0.7314887 0.6803623 0.6308057

sapply(1:10, function(i) cor(rv_1, rv_1^i, method = "spearman"))   # Spearman

#>  [1] 1 1 1 1 1 1 1 1 1 1

# Non-monotonic: all three coefficients understate the dependence
rv_x = -0.5 * rv_1^2 + 11 * rv_1 + rnorm(1000)
cor(rv_1, rv_x, method = "pearson"); cor(rv_1, rv_x, method = "kendall")

#> [1] 0.5200444

#> [1] 0.356012

# Cosine similarity (Pearson without centering)
sum(rv_1 * rv_x) / (sqrt(sum(rv_1^2)) * sqrt(sum(rv_x^2)))

#> [1] 0.9852853

# Two categorical variables
tbl = matrix(c(65, 45, 20, 40), nrow = 2, byrow = TRUE,
             dimnames = list(City = c("B", "T"), Pet = c("Cat", "Dog")))
chisq.test(tbl)                            # Yates-corrected for 2 x 2 tables

#> 
#> 	Pearson's Chi-squared test with Yates' continuity correction
#> 
#> data:  tbl
#> X-squared = 9.2985, df = 1, p-value = 0.002293

sqrt(chisq.test(tbl, correct = FALSE)$statistic / sum(tbl))   # Cramer's V

#> X-squared 
#>  0.246183

library(DescTools)
UncertCoef(tbl, direction = "symmetric")

#> [1] 0.04580323

Takeaways

Pearson's $r$ measures linear association, Spearman's $\rho$ measures monotonic association, and Kendall's $\tau$ is a difference between two probabilities.
The magnitudes of $\tau$ and $r$ are not comparable; $\rho$ and $r$ roughly are.
A correlation near zero does not mean that two variables are unrelated.
$\chi^2$ measures evidence; Cramér's $V$ and Theil's $U$ measure strength.
The sign of an association between nominal variables is arbitrary.

Part II · Inference · Chapter 3

Why a Sample Can Speak for a Population

The hook

One spoon of soup

A chef tastes a single spoon of soup and decides that the whole pot needs more salt. Nobody calls this reckless. However, when a poll of 1,000 residents claims to describe 7.5 million people in Hong Kong, many people are skeptical. The chef and the pollster rely on the same logic and on the same condition: the pot must be well stirred, i.e., the sample must be random.

The weird question

Why can we use a sample to infer the population, and under what conditions? What is a sampling distribution? Is it observable? Why does a large sample size work?

Samples vary, and that is fine

Suppose the population consists of 1,000 values with a mean of 15 and a standard deviation of 3. A simple random sample of 100 values gives a mean close to 15, but not exactly 15. A second sample gives a slightly different value, and a third sample differs again. The important question is not whether one sample is exactly right, but how far samples typically stray from the truth.

The sampling distribution

If I draw many samples and compute the mean of each, the distribution of these sample means is the sampling distribution of the mean. The same can be done for the standard deviation, the median, or any other statistic. Is the sampling distribution observable? Not in practice, because a researcher usually has only one sample. The sampling distribution is a thought experiment, which I can reconstruct either through theory, such as the central limit theorem, or through resampling, such as the bootstrap in Chapter 6.

The central limit theorem

The central limit theorem (CLT) states that, for a population with mean $\mu$ and finite standard deviation $\sigma$, the sampling distribution of the mean approaches a Normal distribution as the sample size $n$ grows, regardless of the shape of the population:

$$\bar{X}\ \dot\sim\ \mathsf{Norm}\!\left(\mu,\ \frac{\sigma}{\sqrt{n}}\right)$$

The standard deviation of the sampling distribution, $\sigma/\sqrt{n}$, is called the standard error (SE). Two consequences follow. First, the sampling distribution is centered on $\mu$, so a random sample is right on average. Second, the SE shrinks with $\sqrt{n}$, so a sample of 1,000 is about three times more precise than a sample of 100 ($\sqrt{10}\approx 3.16$).

Play · The sampling machine

Population Statistic Number of samples $N$

Sample with replacement Biased spoon (larger values more likely to be picked) Keep the population's x-axis

Try this: with the Normal population, the mean, and $n=100$, compare the observed SE with the theoretical $\sigma/\sqrt{n}$.
Choose the Exponential population and move $n$ from 2 to 50. The sampling distribution of the mean changes from skewed to bell-shaped.
Change $N$ from 100 to 5,000 while keeping $n$ fixed. The histogram becomes smoother, but its width (the SE) does not change.
Choose the Poisson population and the median. The sampling distribution is lumpy, not Normal.
Set $n$ to about 800 without replacement, and then tick Sample with replacement. Compare the two SEs.
Tick Biased spoon and increase $n$. The distribution gets narrower, but it stays centered at the wrong value.

Sample size $n$ versus the number of samples $N$

Two numbers are easy to confuse in simulations. The sample size $n$ is the number of observations in each sample. The number of samples $N$ is how many times I repeat the sampling. $N$ does not change the sampling distribution; it only determines how sharply I can see it, much like the resolution of a photograph. The sample size $n$ is what changes the sampling distribution: a larger $n$ makes it narrower and, for the mean, closer to Normal. A real study has $N=1$ and some $n$.

Not every sampling distribution is Normal

The CLT is a theorem about means and sums. Other statistics have their own sampling distributions. The median of a Poisson variable can only take a few values, so its sampling distribution is lumpy. The sampling distribution of the standard deviation from a skewed population is itself skewed in small samples. Many statistics become approximately Normal in large samples, but this property should be checked rather than assumed. Nevertheless, the averages of these sampling distributions tend to sit near the corresponding population values, which you can confirm in the widget.

Sampling from a finite population

When sampling without replacement from a finite population of size $P$, the standard error includes a finite population correction:

$$SE=\frac{\sigma}{\sqrt{n}}\sqrt{\frac{P-n}{P-1}}$$

When $n$ approaches $P$, the correction approaches 0: a sample that includes everyone has no sampling error. This explains why drawing 800 out of 1,000 values produces a much narrower sampling distribution than $\sigma/\sqrt{800}$ predicts. For a poll of 1,000 people in a city of 7.5 million, the correction is essentially 1. This leads to a counterintuitive fact: the precision of a poll depends on the sample size, not on the size of the population.

Why large samples work, and when they do not

A large sample reduces sampling error, i.e., random variation from sample to sample. It does not reduce bias. If the spoon only scoops from the top of the pot, a bigger spoon simply tastes the top more precisely. In 1936, the Literary Digest collected about 2.4 million straw-poll responses, mainly from car owners and telephone subscribers, and confidently predicted the wrong winner of the United States presidential election. Therefore, random sampling is the condition that makes inference possible, and the sample size only determines how precise that inference is.

Note on the original notes

An earlier version of these notes stated that all sampling distributions "should be Normal," that a larger number of samples $N$ makes them "more Normal," and that the SE approaches 0 as $n$ approaches the population size because of the CLT. This chapter corrects these points: normality applies to the mean in large samples, $N$ only affects resolution, and the shrinking SE near the population size comes from the finite population correction.

Common traps

Treating the number of simulated samples $N$ as if it were the sample size $n$.
Assuming that the sampling distribution of every statistic is Normal.
Believing that a very large sample compensates for a biased sampling frame.

R Corner: code and results

set.seed(12345)
rv_1 = rnorm(1000, 15, 3)                 # treat these 1,000 values as the population
c(mean(rv_1), sd(rv_1), median(rv_1))

s_1 = sample(rv_1, 100); c(mean(s_1), sd(s_1), median(s_1))
s_2 = sample(rv_1, 100); c(mean(s_2), sd(s_2), median(s_2))

# Repeat the sampling N times with sample size n
clt_exper = function(n, N, pop = rv_1, stat = mean, replace = FALSE) {
  replicate(N, stat(sample(pop, n, replace = replace)))
}
hist(clt_exper(80, 1000))
sd(clt_exper(80, 1000));  sd(rv_1) / sqrt(80)                 # SE versus theory

# Close to the population size: the finite population correction
sd(clt_exper(800, 1000))
sd(rv_1) / sqrt(800) * sqrt((1000 - 800) / (1000 - 1))
sd(clt_exper(800, 1000, replace = TRUE))                     # no correction

# Not every sampling distribution is Normal
rv_2 = rpois(1000, 2.5)
hist(clt_exper(100, 1000, pop = rv_2, stat = median))        # lumpy
hist(clt_exper(100, 1000, pop = rv_2, stat = mean))          # bell-shaped

set.seed(12345)
rv_1 = rnorm(1000, 15, 3)                 # treat these 1,000 values as the population
c(mean(rv_1), sd(rv_1), median(rv_1))

#> [1] 15.138594  2.996243 15.138650

s_1 = sample(rv_1, 100); c(mean(s_1), sd(s_1), median(s_1))

#> [1] 15.306339  2.906789 15.820623

s_2 = sample(rv_1, 100); c(mean(s_2), sd(s_2), median(s_2))

#> [1] 15.182909  2.887464 15.349239

# Repeat the sampling N times with sample size n
clt_exper = function(n, N, pop = rv_1, stat = mean, replace = FALSE) {
  replicate(N, stat(sample(pop, n, replace = replace)))
}
hist(clt_exper(80, 1000))

sd(clt_exper(80, 1000));  sd(rv_1) / sqrt(80)                 # SE versus theory

#> [1] 0.3176996

#> [1] 0.3349901

# Close to the population size: the finite population correction
sd(clt_exper(800, 1000))

#> [1] 0.04863719

sd(rv_1) / sqrt(800) * sqrt((1000 - 800) / (1000 - 1))

#> [1] 0.04739846

sd(clt_exper(800, 1000, replace = TRUE))                     # no correction

#> [1] 0.1025908

# Not every sampling distribution is Normal
rv_2 = rpois(1000, 2.5)
hist(clt_exper(100, 1000, pop = rv_2, stat = median))        # lumpy

hist(clt_exper(100, 1000, pop = rv_2, stat = mean))          # bell-shaped

Takeaways

A sampling distribution is a thought experiment: the distribution of a statistic across all possible samples.
By the CLT, the sample mean is approximately Normal, centered on $\mu$, with $SE=\sigma/\sqrt{n}$.
The sample size $n$ changes precision; the number of simulated samples $N$ only changes resolution.
The precision of a random sample depends on $n$, not on the size of the population (unless the sample covers a large share of it).
A large sample reduces random error, not bias.

Part II · Inference · Chapter 4

Standard Errors, Sampling Errors, and Confidence Intervals

The hook

Throwing 100 fishing nets

A fisher throws a net into a murky pond to catch one fish whose position is unknown. A wider net catches the fish more often, but it says less precisely where the fish is. A confidence interval is such a net. The population mean is the fish, and each sample throws one net. The confidence level describes how often nets thrown in this way catch the fish; it does not say whether this particular net did.

The weird question

What are the relationships between standard errors, sampling errors, confidence intervals, and confidence levels?

Four terms, one picture

All four terms describe the sampling distribution from Chapter 3, but from different angles.

Term	Meaning
Sampling error	The difference between a sample statistic and the population parameter, e.g., $\bar{x}-\mu$. It is unknown in any real study.
Standard error (SE)	The standard deviation of the sampling distribution, i.e., the typical size of the sampling error. For the mean, $SE=\sigma/\sqrt{n}$.
Margin of error	$z_{\alpha/2}\times SE$. Pollsters often call this number "the sampling error."
Confidence interval (CI)	The statistic plus or minus the margin of error.
Confidence level	$1-\alpha$: the long-run share of intervals, constructed in this way, that contain the parameter.

If the sampling distribution of the mean is Normal, 95% of sample means fall within $1.96\,SE$ of $\mu$. Turning this statement around gives the 95% confidence interval:

$$CI_{1-\alpha}=\bar{x}\pm z_{\alpha/2}\,\frac{\sigma}{\sqrt{n}},\qquad CI_{95\%}=\bar{x}\pm 1.96\,\frac{\sigma}{\sqrt{n}}$$

In practice, a researcher observes neither the sampling distribution nor $\sigma$, but only one sample, with mean $\bar{x}$ and standard deviation $s$. Replacing $\sigma$ with $s$ is not free: it adds uncertainty, which is why the critical value 1.96 should be replaced by the t distribution's $t_{n-1,\,\alpha/2}$ when $n$ is small (Chapter 5):

$$CI_{1-\alpha}=\bar{x}\pm t_{n-1,\,\alpha/2}\,\frac{s}{\sqrt{n}}$$

Play · Throwing 100 nets

Population Confidence level Critical value

Try this: set $n=5$ with the z critical value. Throw several hundred nets. The coverage falls well below 95%, because 1.96 ignores the uncertainty in $s$. Switch to t.
Change the confidence level from 95% to 99%. The nets get wider and miss less often.
Choose the skewed population with $n=5$: even the t interval under-covers. With $n=100$, the CLT rescues it.

Think about it

In one sample, the true mean falls inside the 95% CI. Is it possible that, for another random sample, the true mean is not in the interval?

My answer

Yes, and the frequency is known: about 5% of the time, which is exactly what the red nets in the widget show. The 95% describes the procedure, not any single interval. Once an interval has been calculated, it either contains $\mu$ or it does not. Therefore, the statement "there is a 95% probability that $\mu$ lies in [14.2, 15.8]" is not correct in the frequentist framework. A Bayesian credible interval does support such a statement, but it requires a prior distribution.

Common traps

Reading a 95% CI as "95% of the data lie in this interval." The interval describes the mean, not individual observations.
Using 1.96 with a small sample and an estimated standard deviation.
Treating the margin of error as the total error. It only covers random sampling error, not bias or measurement error.

R Corner: code and results

set.seed(12345)
rv_1 = rnorm(1000, 15, 3)                         # the population
clt_exper = function(n, N) replicate(N, mean(sample(rv_1, n)))

# 95% of the sample means fall within about 1.96 SE of the center
sdm = clt_exper(100, 1000)
quantile(sdm, probs = c(0.025, 0.975))
mean(sdm) + c(-1.96, 1.96) * sd(sdm)

# One sample, one interval
s_1 = sample(rv_1, 100)
m = mean(s_1); s = sd(s_1)
m + c(-1.96, 1.96) * s / sqrt(100)                # z-based interval
m + qt(c(0.025, 0.975), df = 99) * s / sqrt(100)  # t-based interval (slightly wider)
mean(rv_1)                                        # the fish

# Coverage with n = 5: z versus t
cover = function(crit) mean(replicate(5000, {
  s = sample(rv_1, 5)
  ci = mean(s) + c(-1, 1) * crit * sd(s) / sqrt(5)
  ci[1] <= mean(rv_1) && mean(rv_1) <= ci[2]
}))
c(z = cover(1.96), t = cover(qt(0.975, df = 4)))  # about .88 versus .95

set.seed(12345)
rv_1 = rnorm(1000, 15, 3)                         # the population
clt_exper = function(n, N) replicate(N, mean(sample(rv_1, n)))

# 95% of the sample means fall within about 1.96 SE of the center
sdm = clt_exper(100, 1000)
quantile(sdm, probs = c(0.025, 0.975))

#>     2.5%    97.5% 
#> 14.56538 15.64177

mean(sdm) + c(-1.96, 1.96) * sd(sdm)

#> [1] 14.58329 15.66010

# One sample, one interval
s_1 = sample(rv_1, 100)
m = mean(s_1); s = sd(s_1)
m + c(-1.96, 1.96) * s / sqrt(100)                # z-based interval

#> [1] 14.66202 15.76014

m + qt(c(0.025, 0.975), df = 99) * s / sqrt(100)  # t-based interval (slightly wider)

#> [1] 14.65523 15.76692

mean(rv_1)                                        # the fish

#> [1] 15.13859

# Coverage with n = 5: z versus t
cover = function(crit) mean(replicate(5000, {
  s = sample(rv_1, 5)
  ci = mean(s) + c(-1, 1) * crit * sd(s) / sqrt(5)
  ci[1] <= mean(rv_1) && mean(rv_1) <= ci[2]
}))
c(z = cover(1.96), t = cover(qt(0.975, df = 4)))  # about .88 versus .95

#>      z      t 
#> 0.8692 0.9514

Takeaways

The SE is the standard deviation of the sampling distribution; the margin of error is a multiple of it.
A confidence level is a property of the procedure: about 95% of 95% intervals contain the parameter.
With an estimated standard deviation and a small sample, use t rather than z.
Precision grows with $\sqrt{n}$: halving the width requires four times as many observations.

Part II · Inference · Chapter 5

t, Z, F, and $\chi^2$ Tests: One Idea, Different Distributions

The hook

Is the 30 cm pizza really 30 cm?

A pizza chain advertises 30 cm pizzas. A group of hungry students measures 12 pizzas and finds an average diameter of 29.5 cm. Is the chain cheating, or is this the ordinary variation of hand-stretched dough? Every test in this chapter answers a version of the same question: is the difference between what I observed and what I expected larger than chance alone would produce?

The weird question

What are the differences among the t-test, the Z-test, the F-test, and the $\chi^2$-test?

One recipe for every test statistic

A test statistic measures the difference between the observed data and what the null hypothesis predicts, in units of the standard error:

$$TS=\frac{\text{observed value}-\text{expected value}}{SE}$$

The larger the test statistic (in absolute value), the more unusual the data look under the null hypothesis. For a two-sided Z-test, $|Z|>1.96$ implies $p<.05$, and $|Z|>2.58$ implies $p<.01$. Why? Under the null hypothesis, $Z=(\bar{x}-\mu_0)/SE\sim\mathsf{Norm}(0,1)$, and only 5% of a standard Normal distribution lies beyond $\pm 1.96$.

The tests differ mainly in the distribution that the test statistic follows under the null hypothesis. In other words, they share one idea and use different reference distributions.

Z versus t: the price of estimating $\sigma$

A Z-test assumes that the population standard deviation $\sigma$ is known, so that $SE=\sigma/\sqrt{n}$. This is rarely the case. When $\sigma$ is estimated by $s$, the statistic $T=(\bar{x}-\mu_0)/(s/\sqrt{n})$ follows a Student t distribution with $\nu=n-1$ degrees of freedom. The t distribution has heavier tails than the Normal distribution, which reflects the extra uncertainty in $s$. As $n\to\infty$, $s$ becomes precise and $\mathsf{StudentT}(\nu)\to\mathsf{Norm}(0,1)$.

Play · t versus Normal

Play · The pizza test

The chain claims $\mu_0=30$ cm. You decide what the chain really bakes; the students only see their sample.

Alternative

What the t-test assumes

An interval or ratio measurement scale.
A simple random sample, with independent observations.
A Normal sampling distribution of the mean (not of the raw data!), which the CLT provides in large samples.
For the classic two-sample test: equal variances in the two groups.

The t-test is unreliable for small samples from strongly non-normal distributions, but it works well for large samples from non-normal distributions (see this discussion). The R Corner shows that the Z-score of the mean of a Poisson variable, which is skewed, is approximately Normal with $n=100$.

Two samples, the variance sum law, and Welch's test

To compare two independent means, the test statistic is $T=(\bar{x}_1-\bar{x}_2)/SE$. The challenge is the SE. According to the variance sum law, for independent variables, $\mathrm{Var}(X_1-X_2)=\mathrm{Var}(X_1)+\mathrm{Var}(X_2)$. Therefore, Welch's test uses

$$SE_{\text{Welch}}=\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}}$$

The classic Student test assumes equal variances and pools them:

$$s_p^2=\frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2},\qquad SE_{\text{pooled}}=s_p\sqrt{\frac{1}{n_1}+\frac{1}{n_2}}$$

Note on the original notes

An earlier version of these notes omitted the factor $\sqrt{1/n_1+1/n_2}$ from the pooled formula and concluded, from a simulation with $n_1=n_2=100$, that the pooled test overestimates significance. With equal group sizes, however, the two SEs are algebraically identical: $s_p^2(2/n)=(s_1^2+s_2^2)/n$. The pooled test only goes wrong when the group sizes and the variances differ, as the widget below shows.

Play · Welch versus pooled under the null hypothesis

Both groups have the same mean, so every significant result is a false positive. A good test rejects about 5% of the time. Each change runs 2,000 simulated studies.

Try this: give the small group the large SD (the default). The pooled test rejects far more than 5% of the time.
Give the large group the large SD. Now the pooled test becomes too conservative.
Make $n_1=n_2$. Both tests behave well, whatever the SDs.

Simulation studies (see this discussion) show that Welch's test performs as well as the pooled test when variances are equal and much better when they are not. Therefore, I recommend Welch's test by default, which is also R's default in t.test.

Beyond two means

Situation	Test	Reference distribution under $H_0$
One mean, $\sigma$ known	Z-test	$\mathsf{Norm}(0,1)$
One mean or two means, $\sigma$ estimated	t-test	$t$ with $n-1$ (or Welch) df
Two dependent measurements (before and after)	Paired t-test	$t$ with $n-1$ df, on the differences
More than two means	F-test (ANOVA)	$F$ with $k-1$ and $N-k$ df
Two categorical variables	$\chi^2$-test	$\chi^2$ with $(r-1)(c-1)$ df

The F statistic is a ratio of two variances: the variance between group means divided by the variance within groups. The $\chi^2$ statistic sums squared standardized differences between observed and expected counts. Both follow the same recipe: observed versus expected, scaled by chance variation.

Common traps

Requiring the raw data to be Normal. The assumption concerns the sampling distribution of the mean.
Testing for equal variances first and then choosing a test. Use Welch's test by default.
Using an independent-samples test for paired data, such as the same participants measured twice.

R Corner: code and results

set.seed(12345)
rv_1 = rnorm(1000, 15, 3)
rv_2 = rpois(1000, 2.5)
s_1 = sample(rv_1, 100); s_2 = sample(rv_1, 100)

t.test(s_1, mu = 14, alternative = "greater")   # one-sample t-test
t.test(s_1, s_2)                                # Welch two-sample t-test (R's default)
t.test(s_1, s_2, var.equal = TRUE)              # pooled-variance t-test

# The Z-score of the mean of a skewed variable is approximately Normal
zs = replicate(1000, {
  sp = sample(rv_2, 100)
  (mean(sp) - mean(rv_2)) / (sd(sp) / sqrt(100))
})
hist(zs); mean(abs(zs) > 1.96)                  # close to .05

# With n1 = n2, the pooled and Welch t statistics are identical
sp_1 = sample(rv_1, 100); sp_2 = sample(rv_2, 100)
c(pooled = t.test(sp_1, sp_2, var.equal = TRUE)$statistic,
  welch  = t.test(sp_1, sp_2)$statistic)

# ... but not when the smaller group has the larger variance (both means are 0)
p_pooled = replicate(2000, t.test(rnorm(10, 0, 3), rnorm(40, 0, 1), var.equal = TRUE)$p.value)
p_welch  = replicate(2000, t.test(rnorm(10, 0, 3), rnorm(40, 0, 1))$p.value)
c(pooled = mean(p_pooled < .05), welch = mean(p_welch < .05))

# More than two groups: F-test; two categorical variables: chi-squared test
summary(aov(yield ~ block, data = npk))
chisq.test(matrix(c(65, 45, 20, 40), nrow = 2, byrow = TRUE))

set.seed(12345)
rv_1 = rnorm(1000, 15, 3)
rv_2 = rpois(1000, 2.5)
s_1 = sample(rv_1, 100); s_2 = sample(rv_1, 100)

t.test(s_1, mu = 14, alternative = "greater")   # one-sample t-test

#> 
#> 	One Sample t-test
#> 
#> data:  s_1
#> t = 3.7026, df = 99, p-value = 0.0001756
#> alternative hypothesis: true mean is greater than 14
#> 95 percent confidence interval:
#>  14.60871      Inf
#> sample estimates:
#> mean of x 
#>  15.10361

t.test(s_1, s_2)                                # Welch two-sample t-test (R's default)

#> 
#> 	Welch Two Sample t-test
#> 
#> data:  s_1 and s_2
#> t = -0.35014, df = 198, p-value = 0.7266
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#>  -0.9807552  0.6849924
#> sample estimates:
#> mean of x mean of y 
#>  15.10361  15.25149

t.test(s_1, s_2, var.equal = TRUE)              # pooled-variance t-test

#> 
#> 	Two Sample t-test
#> 
#> data:  s_1 and s_2
#> t = -0.35014, df = 198, p-value = 0.7266
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#>  -0.9807552  0.6849924
#> sample estimates:
#> mean of x mean of y 
#>  15.10361  15.25149

# The Z-score of the mean of a skewed variable is approximately Normal
zs = replicate(1000, {
  sp = sample(rv_2, 100)
  (mean(sp) - mean(rv_2)) / (sd(sp) / sqrt(100))
})
hist(zs); mean(abs(zs) > 1.96)                  # close to .05

#> [1] 0.036

# With n1 = n2, the pooled and Welch t statistics are identical
sp_1 = sample(rv_1, 100); sp_2 = sample(rv_2, 100)
c(pooled = t.test(sp_1, sp_2, var.equal = TRUE)$statistic,
  welch  = t.test(sp_1, sp_2)$statistic)

#> pooled.t  welch.t 
#> 36.51843 36.51843

# ... but not when the smaller group has the larger variance (both means are 0)
p_pooled = replicate(2000, t.test(rnorm(10, 0, 3), rnorm(40, 0, 1), var.equal = TRUE)$p.value)
p_welch  = replicate(2000, t.test(rnorm(10, 0, 3), rnorm(40, 0, 1))$p.value)
c(pooled = mean(p_pooled < .05), welch = mean(p_welch < .05))

#> pooled  welch 
#> 0.2625 0.0455

# More than two groups: F-test; two categorical variables: chi-squared test
summary(aov(yield ~ block, data = npk))

#>             Df Sum Sq Mean Sq F value Pr(>F)  
#> block        5  343.3   68.66   2.318 0.0861 .
#> Residuals   18  533.1   29.62                 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

chisq.test(matrix(c(65, 45, 20, 40), nrow = 2, byrow = TRUE))

#> 
#> 	Pearson's Chi-squared test with Yates' continuity correction
#> 
#> data:  matrix(c(65, 45, 20, 40), nrow = 2, byrow = TRUE)
#> X-squared = 9.2985, df = 1, p-value = 0.002293

Takeaways

Every test statistic is (observed − expected) / SE; tests differ in their reference distributions.
Estimating $\sigma$ turns Z into t, whose heavier tails matter in small samples.
With equal group sizes, pooled and Welch statistics coincide; otherwise, prefer Welch.
F compares variances between and within groups; $\chi^2$ compares observed and expected counts.

Part II · Inference · Chapter 6

Distribution-Free Tests: Permutation, Bootstrap, and Exact Tests

The hook

The lady tasting tea

In the 1920s, Muriel Bristol, a colleague of the statistician Ronald Fisher, claimed that she could tell whether milk had been poured into the cup before or after the tea. Fisher prepared eight cups, four of each kind, in random order. How many cups must she identify correctly before her claim is believable? Fisher's answer did not require any Normal distribution. It required only counting the ways in which the cups could have been arranged.

The weird question

What are distribution-free (non-parametric) tests? We test differences between means; can we test differences between medians or variances? Can we test the association between two categorical variables?

Counting instead of assuming

There are $\binom{8}{4}=70$ ways to choose four cups out of eight. If the lady is only guessing, each choice is equally likely. Only one of the 70 choices identifies all four milk-first cups, so the probability of a perfect score by luck is $1/70\approx .014$. This is Fisher's exact test: the p-value comes from counting all possible arrangements under the null hypothesis, not from an approximate distribution.

Play · Be the lady

Four of these cups had milk poured first. Pick the four you believe are milk-first, then reveal the truth.

Select four cups.

What "distribution-free" means

A t-test is not always appropriate, and the reason is not only non-normality or unequal variances. Sometimes the mean is not the interesting quantity. Consider two restaurants whose waiting times both average 7.5 minutes. At the first restaurant, customers wait either about 5 or about 10 minutes; at the second, either about 4 or about 11 minutes. The means are the same, but the spreads differ.

For some statistics, such as the difference in means, the sampling distribution is known. For others, such as the difference in standard deviations of bimodal data or the difference in medians, it is not. "Distribution-free" refers to this: such tests do not assume a particular distribution of the test statistic. They build the null distribution from the data themselves.

The randomization test

Compute the observed statistic, e.g., $sd(x_1)-sd(x_2)$.
Pool the two groups and randomly shuffle the group labels.
Split the shuffled values into two groups and compute the statistic again.
Repeat many times; the shuffled statistics form the null distribution.
The p-value is the share of shuffled statistics at least as extreme as the observed one.

Play · The shuffling machine

Statistic

Try this: test the difference in means. The observed difference sits in the middle of the shuffled differences: no evidence of a difference.
Test the difference in standard deviations or interquartile ranges. The observed value lies far outside the shuffled values.

The bootstrap

The bootstrap resamples the observed data with replacement and recomputes the statistic on each resample. The spread of the bootstrap statistics approximates the sampling distribution, from which a standard error or a confidence interval can be read. Recall from Chapter 2 that Pearson's $r$ does not assume normality, but its classical confidence interval does. With highly skewed data, the bootstrap offers an honest alternative.

Play · Bootstrapping a correlation between skewed variables

The data are 300 pairs of highly skewed (Weibull) variables, such as the numbers of likes and shares of social-media posts.

Coefficient

For Pearson's r, the classical interval uses Fisher's z transformation, which assumes bivariate Normal data. For Spearman's ρ, the same formula is applied to the ranks.

Randomization, permutation, and bootstrap

Randomization tests shuffle the data many times and compute the statistic on each shuffle. Permutation tests do the same thing, but I reserve the label for tests that use all possible permutations, as in Fisher's tea experiment. Bootstrapping resamples with replacement and is mainly used for estimation. In other words, bootstrapping is concerned with the population from which the data came, whereas randomization tests focus on the mechanism that distributed the data between groups (see a detailed introduction).

Contingency tables. Fisher's exact test extends the tea logic to any 2 × 2 table, such as cat people and dog people in Chapter 2. Its effect size is the odds ratio. In City B, the odds of being a cat person are $65/45=1.44$; in City T, they are $20/40=0.50$. The odds ratio is $1.44/0.50=(65\times 40)/(45\times 20)=2.89$: the odds of being a cat person are 2.89 times higher in City B. (R's fisher.test reports a slightly different conditional estimate.)

Common traps

Believing that non-parametric tests have no assumptions. Randomization tests still assume exchangeability under the null hypothesis, and the bootstrap assumes that the sample represents the population.
Treating a difference in means as the only possible difference between groups.
Bootstrapping a tiny sample. The bootstrap cannot create information that the sample does not contain.

R Corner: code and results

set.seed(12345)
v_3 = c(rnorm(500, 5, 1), rnorm(500, 10, 1))   # restaurant A
v_4 = c(rnorm(500, 4, 1), rnorm(500, 11, 1))   # restaurant B
t.test(v_3, v_4)                       # the means do not differ
quantile(v_3); quantile(v_4)           # but the spreads do

# Randomization test for the difference in SDs
obs = sd(v_3) - sd(v_4)
perm = replicate(1000, {
  pv = sample(c(v_3, v_4))
  sd(pv[1:1000]) - sd(pv[1001:2000])
})
mean(abs(perm) >= abs(obs))            # randomization p-value
var.test(v_3, v_4)                     # parametric alternative (assumes normality)

# Bootstrap the difference in SDs
library(boot)
data = data.frame(v_3, v_4)
sd_diff = function(d, i) sd(d$v_3[i]) - sd(d$v_4[i])
boot_sd = boot(data, sd_diff, R = 1000)
boot.ci(boot_sd, type = c("norm", "perc", "bca"))

# Fisher's exact test: the lady tasting tea (3 of 4 cups correct), and cats versus dogs
fisher.test(matrix(c(3, 1, 1, 3), nrow = 2), alternative = "greater")   # p = 17/70
tbl = matrix(c(65, 45, 20, 40), nrow = 2, byrow = TRUE)
fisher.test(tbl)

# Bootstrapping correlation coefficients of skewed variables
r1 = rweibull(1000, 0.5)
r2 = 0.5 * r1 + rweibull(1000, 0.5)
d = data.frame(r1, r2)
cor.test(r1, r2)$conf.int              # classical interval (assumes bivariate Normal)
boot_r = boot(d, function(d, i) cor(d$r1[i], d$r2[i]), R = 1000)
boot.ci(boot_r, type = "perc")

set.seed(12345)
v_3 = c(rnorm(500, 5, 1), rnorm(500, 10, 1))   # restaurant A
v_4 = c(rnorm(500, 4, 1), rnorm(500, 11, 1))   # restaurant B
t.test(v_3, v_4)                       # the means do not differ

#> 
#> 	Welch Two Sample t-test
#> 
#> data:  v_3 and v_4
#> t = 0.53567, df = 1824.8, p-value = 0.5922
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#>  -0.2038310  0.3570123
#> sample estimates:
#> mean of x mean of y 
#>  7.546198  7.469608

quantile(v_3); quantile(v_4)           # but the spreads do

#>        0%       25%       50%       75%      100% 
#>  2.418145  5.102542  7.480331  9.997535 13.330733

#>         0%        25%        50%        75%       100% 
#>  0.6653246  3.9730550  7.5814484 10.9798811 13.7609997

# Randomization test for the difference in SDs
obs = sd(v_3) - sd(v_4)
perm = replicate(1000, {
  pv = sample(c(v_3, v_4))
  sd(pv[1:1000]) - sd(pv[1001:2000])
})
mean(abs(perm) >= abs(obs))            # randomization p-value

#> [1] 0

var.test(v_3, v_4)                     # parametric alternative (assumes normality)

#> 
#> 	F test to compare two variances
#> 
#> data:  v_3 and v_4
#> F = 0.52892, num df = 999, denom df = 999, p-value < 2.2e-16
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#>  0.4671982 0.5988058
#> sample estimates:
#> ratio of variances 
#>          0.5289244

# Bootstrap the difference in SDs
library(boot)
data = data.frame(v_3, v_4)
sd_diff = function(d, i) sd(d$v_3[i]) - sd(d$v_4[i])
boot_sd = boot(data, sd_diff, R = 1000)
boot.ci(boot_sd, type = c("norm", "perc", "bca"))

#> BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
#> Based on 1000 bootstrap replicates
#> 
#> CALL : 
#> boot.ci(boot.out = boot_sd, type = c("norm", "perc", "bca"))
#> 
#> Intervals : 
#> Level      Normal             Percentile            BCa          
#> 95%   (-1.0848, -0.9149 )   (-1.0798, -0.9104 )   (-1.0820, -0.9125 )  
#> Calculations and Intervals on Original Scale

# Fisher's exact test: the lady tasting tea (3 of 4 cups correct), and cats versus dogs
fisher.test(matrix(c(3, 1, 1, 3), nrow = 2), alternative = "greater")   # p = 17/70

#> 
#> 	Fisher's Exact Test for Count Data
#> 
#> data:  matrix(c(3, 1, 1, 3), nrow = 2)
#> p-value = 0.2429
#> alternative hypothesis: true odds ratio is greater than 1
#> 95 percent confidence interval:
#>  0.3135693       Inf
#> sample estimates:
#> odds ratio 
#>   6.408309

tbl = matrix(c(65, 45, 20, 40), nrow = 2, byrow = TRUE)
fisher.test(tbl)

#> 
#> 	Fisher's Exact Test for Count Data
#> 
#> data:  tbl
#> p-value = 0.002162
#> alternative hypothesis: true odds ratio is not equal to 1
#> 95 percent confidence interval:
#>  1.427537 5.912307
#> sample estimates:
#> odds ratio 
#>   2.870314

# Bootstrapping correlation coefficients of skewed variables
r1 = rweibull(1000, 0.5)
r2 = 0.5 * r1 + rweibull(1000, 0.5)
d = data.frame(r1, r2)
cor.test(r1, r2)$conf.int              # classical interval (assumes bivariate Normal)

#> [1] 0.4071862 0.5052858
#> attr(,"conf.level")
#> [1] 0.95

boot_r = boot(d, function(d, i) cor(d$r1[i], d$r2[i]), R = 1000)
boot.ci(boot_r, type = "perc")

#> BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
#> Based on 1000 bootstrap replicates
#> 
#> CALL : 
#> boot.ci(boot.out = boot_r, type = "perc")
#> 
#> Intervals : 
#> Level     Percentile     
#> 95%   ( 0.3118,  0.6317 )  
#> Calculations and Intervals on Original Scale

Takeaways

Distribution-free tests build the null distribution of the test statistic from the data.
Randomization tests can compare any statistic: medians, standard deviations, or quantiles.
The bootstrap estimates sampling variability by resampling with replacement.
Fisher's exact test counts arrangements; its effect size is the odds ratio.

Part II · Inference · Chapter 7

Equivalence Tests: Can the Null Hypothesis Be Accepted?

The hook

Is the new milk-tea recipe the same?

A tea chain wants to replace its expensive tea leaves with a cheaper blend. The marketing team asks customers to rate both recipes and finds no significant difference. The manager is delighted: "So customers cannot tell the difference!" However, a non-significant result is not evidence of sameness. Perhaps only 20 customers took part.

The weird question

Is it possible to test equivalence statistically, i.e., to accept rather than merely fail to reject the null hypothesis?

Failing to reject is not accepting

In conventional hypothesis testing, the null hypothesis states that there is no difference, $H_0:\mu_1=\mu_2$, and the alternative states that there is one, $H_1:\mu_1\ne\mu_2$. When $p>.05$, the test never "accepts" $H_0$; it only fails to reject it. A failure to reject can occur because of:

a lack of statistical power, e.g., a small sample;
high variability in the data; or
a true null effect.

Only the last reason implies equivalence, and the t-test cannot distinguish it from the other two.

Turning the hypotheses around

An equivalence test first defines a margin $\delta$, which means the smallest difference that would matter in practice. It then reverses the hypotheses:

$$H_0:\ |\mu_1-\mu_2|\ \ge\ \delta\qquad\text{versus}\qquad H_1:\ |\mu_1-\mu_2|\ <\ \delta$$

The most common procedure is TOST (two one-sided tests). It tests whether the difference is significantly greater than $-\delta$ and significantly smaller than $+\delta$. Equivalently, with $\alpha=.05$, equivalence is declared when the 90% confidence interval of the difference lies entirely within $(-\delta,\delta)$. Combining a t-test and TOST yields four possible verdicts.

Play · The recipe taste test

Ratings have a standard deviation of 15 points. The equivalence margin is the smallest rating difference that the chain cares about.

Try this: keep the default (true difference 1, 20 customers per group, margin 5) and run 1,000 tests. Most results are inconclusive: neither different nor equivalent.
Increase the sample to 300 customers per group. Equivalence becomes the most common verdict.
Set the true difference to 2, the margin to 5, and the sample to 2,000. Many tests are both significant and equivalent: a real difference that is too small to matter.

In the example from my original notes, two groups of 20 were drawn from populations with means of 100 and 101 (SD = 15). The t-test was not significant, and TOST with $\delta=5$ did not establish equivalence either. The result is inconclusive, because the sample is too small to decide. In general, equivalence tests require larger samples than conventional tests, because they must show that the difference is precisely small.

Choosing $\delta$. The margin should be justified substantively, e.g., the smallest rating difference that would change a business decision, and it must be fixed before the data are collected. Choosing $\delta$ after seeing the data makes any result "equivalent."

Common traps

Reporting "no difference" based on $p>.05$.
Choosing the equivalence margin after looking at the results.
Forgetting that a difference can be both statistically significant and practically negligible.

R Corner: code and results

set.seed(123)
group1 = rnorm(20, mean = 100, sd = 15)
group2 = rnorm(20, mean = 101, sd = 15)
t.test(group1, group2)                        # not significant

library(TOSTER)
eqv = t_TOST(x = group1, y = group2, eqb = 5, alpha = 0.05)
eqv                                           # not equivalent either: inconclusive
plot(eqv)

# By hand: equivalent if the 90% CI lies within (-5, 5)
t.test(group1, group2, conf.level = 0.90)$conf.int

# The same design with 1,000 customers per group
big1 = rnorm(1000, 100, 15); big2 = rnorm(1000, 101, 15)
t_TOST(x = big1, y = big2, eqb = 5)

set.seed(123)
group1 = rnorm(20, mean = 100, sd = 15)
group2 = rnorm(20, mean = 101, sd = 15)
t.test(group1, group2)                        # not significant

#> 
#> 	Welch Two Sample t-test
#> 
#> data:  group1 and group2
#> t = 0.44145, df = 37.082, p-value = 0.6615
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#>  -6.795731 10.582160
#> sample estimates:
#> mean of x mean of y 
#>  102.1244  100.2311

library(TOSTER)
eqv = t_TOST(x = group1, y = group2, eqb = 5, alpha = 0.05)
eqv                                           # not equivalent either: inconclusive

#> 
#> Welch Two Sample t-test
#> 
#> The equivalence test was non-significant, t(37.08) = -0.72, p = 0.24
#> The null hypothesis test was non-significant, t(37.08) = 0.441, p = 0.66
#> NHST: don't reject null significance hypothesis that the effect is equal to zero 
#> TOST: don't reject null equivalence hypothesis
#> 
#> TOST Results 
#>                  t    df p.value
#> t-test      0.4414 37.08   0.661
#> TOST Lower  1.6073 37.08   0.058
#> TOST Upper -0.7244 37.08   0.237
#> 
#> Effect Sizes 
#>                Estimate     SE              C.I. Conf. Level
#> Raw              1.8932 4.2886 [-5.3417, 9.1281]         0.9
#> Hedges's g(av)   0.1368 0.3248  [-0.3744, 0.646]         0.9
#> Note: SMD confidence intervals are an approximation. See vignette("SMD_calcs").

plot(eqv)

# By hand: equivalent if the 90% CI lies within (-5, 5)
t.test(group1, group2, conf.level = 0.90)$conf.int

#> [1] -5.341698  9.128127
#> attr(,"conf.level")
#> [1] 0.9

# The same design with 1,000 customers per group
big1 = rnorm(1000, 100, 15); big2 = rnorm(1000, 101, 15)
t_TOST(x = big1, y = big2, eqb = 5)

#> 
#> Welch Two Sample t-test
#> 
#> The equivalence test was significant, t(1997.9) = 5.69, p < 0.01
#> The null hypothesis test was non-significant, t(1997.9) = -1.75, p = 0.08
#> NHST: don't reject null significance hypothesis that the effect is equal to zero 
#> TOST: reject null equivalence hypothesis
#> 
#> TOST Results 
#>                 t   df p.value
#> t-test     -1.749 1998   0.081
#> TOST Lower  5.695 1998 < 0.001
#> TOST Upper -9.192 1998 < 0.001
#> 
#> Effect Sizes 
#>                Estimate      SE               C.I. Conf. Level
#> Raw            -1.17464 0.67175 [-2.2801, -0.0692]         0.9
#> Hedges's g(av) -0.07817 0.04476 [-0.1517, -0.0046]         0.9
#> Note: SMD confidence intervals are an approximation. See vignette("SMD_calcs").

Takeaways

A non-significant result means "not shown to differ," not "shown to be the same."
Equivalence tests reverse the hypotheses, using a margin $\delta$ fixed in advance.
TOST declares equivalence when the 90% CI lies within $(-\delta,\delta)$.
Results can be different, equivalent, both, or neither.

Part III · Regression (OLS) · Chapter 8

Can Regression Replace the Tests?

The hook

One model, many tests

Statistics courses usually introduce the t-test, ANOVA, and the $\chi^2$-test in separate chapters, as if they were separate tools in a toolbox. Students then memorize a flowchart: two groups, use a t-test; three groups, use ANOVA; two categorical variables, use $\chi^2$. In fact, most of these tests are special cases of one family of models. Once you see the regression behind them, the flowchart becomes much shorter.

The weird question

Can regression models (why plural?) replace all the tests in Part II?

Why "models," plural

Ordinary least squares (OLS) regression is only one member of a large family. Logistic regression, Poisson regression, and the other generalized linear models in Part IV are also regression models, and they differ in the distribution that generates the dependent variable. In this chapter, the dependent variable is continuous and the independent variable is categorical, which is exactly the setting of the tests in Chapter 5.

The t-test is a regression with one dummy variable

Code the groups as a dummy variable $D$ (0 for drug 1, 1 for drug 2) and fit

$$y_i=\beta_0+\beta_1 D_i+\varepsilon_i$$

Then $\beta_0$ is the mean of group 1, and $\beta_1$ is the difference between the two group means. The t statistic for $\beta_1$ is identical to the t statistic of the pooled (equal-variance) t-test. R's built-in sleep data record the extra hours of sleep of 10 patients under two drugs: the means are 0.75 and 2.33 hours, so $\hat\beta_1=1.58$.

Play · Same answer, two machines

Try this: move the slider. The regression slope always equals the difference in means, and the two t values stay identical.
Welch's test is a different matter, because lm assumes one common error variance. With 10 patients in each group, Welch's t happens to equal the pooled t (Chapter 5), but its degrees of freedom, and therefore its p value, differ from those of lm.

A note on the sleep data

The sleep data are actually paired: the same 10 patients took both drugs. A paired t-test is therefore more appropriate, and it, too, is a regression: lm(extra ~ group + ID), which adds a dummy variable for each patient, reproduces the paired t statistic exactly (see the R Corner).

ANOVA is a regression with several dummy variables

With $k$ groups, the regression uses $k-1$ dummy variables. The overall F-test of the regression, reported at the bottom of summary(lm(...)), is exactly the ANOVA F-test, because both compare the variance between group means with the variance within groups.

Play · ANOVA versus lm on the npk field trial

R's npk data record pea yields in six blocks of land. Shift the yields in block 3 and compare the two F-tests.

Non-normal outcomes and categorical outcomes

For two groups of Poisson counts, t.test and lm give nearly the same answer, because both rely on the CLT for the difference in means (Welch's test differs slightly in its df). For two categorical variables, the regression is less familiar. A log-linear model treats the cell counts of a contingency table as a Poisson outcome:

$$\log(\mu_{ij})=\lambda+\lambda^{\text{city}}_i+\lambda^{\text{pet}}_j+\lambda^{\text{city}\times\text{pet}}_{ij}$$

The interaction term captures the dependence between the two variables, and its coefficient is exactly the log odds ratio. For the cat-and-dog table, $\log(2.89)=1.06$. Its sign depends on which categories serve as references. In other words, the log-linear model is a regression version of the $\chi^2$-test. Equivalently, with the data disaggregated to one row per person, a logistic regression of pet preference on city yields the same log odds ratio.

Common traps

Forgetting that a factor coefficient is a difference from the reference group, not a group mean.
Analyzing paired data as independent groups.
Believing that regression requires continuous predictors.

R Corner: code and results

head(sleep)
t.test(extra ~ group, data = sleep, var.equal = TRUE)
summary(lm(extra ~ group, data = sleep))        # same t and p; slope = 2.33 - 0.75

# The sleep data are paired (the same 10 patients took both drugs)
t.test(sleep$extra[sleep$group == 2], sleep$extra[sleep$group == 1], paired = TRUE)
summary(lm(extra ~ group + ID, data = sleep))   # same t as the paired test

# More than two groups: ANOVA = lm
aggregate(yield ~ block, npk, mean)
summary(aov(yield ~ block, data = npk))
summary(lm(yield ~ block, data = npk))          # same F and p at the bottom

# Non-normal outcomes
set.seed(12345)
g1 = rpois(500, 2.5); g2 = rpois(500, 3.5)
d = data.frame(dv = c(g1, g2), groups = rep(c(1, 0), each = 500))
t.test(dv ~ groups, data = d)
summary(lm(dv ~ groups, data = d))

# Contingency tables: a log-linear (Poisson) model
ctbl = data.frame(freq = c(65, 45, 20, 40), city = c("B", "B", "T", "T"),
                  pet = c("Cat", "Dog", "Cat", "Dog"))
summary(glm(freq ~ city * pet, family = poisson, data = ctbl))
log((65 * 40) / (45 * 20))                      # the interaction equals the log odds ratio

# ... or a logistic regression on the disaggregated data (one row per person)
people = ctbl[rep(1:4, ctbl$freq), c("city", "pet")]
coef(glm(I(pet == "Dog") ~ city, family = binomial, data = people))

head(sleep)

#>   extra group ID
#> 1   0.7     1  1
#> 2  -1.6     1  2
#> 3  -0.2     1  3
#> 4  -1.2     1  4
#> 5  -0.1     1  5
#> 6   3.4     1  6

t.test(extra ~ group, data = sleep, var.equal = TRUE)

#> 
#> 	Two Sample t-test
#> 
#> data:  extra by group
#> t = -1.8608, df = 18, p-value = 0.07919
#> alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
#> 95 percent confidence interval:
#>  -3.363874  0.203874
#> sample estimates:
#> mean in group 1 mean in group 2 
#>            0.75            2.33

summary(lm(extra ~ group, data = sleep))        # same t and p; slope = 2.33 - 0.75

#> 
#> Call:
#> lm(formula = extra ~ group, data = sleep)
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#> -2.430 -1.305 -0.580  1.455  3.170 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)  
#> (Intercept)   0.7500     0.6004   1.249   0.2276  
#> group2        1.5800     0.8491   1.861   0.0792 .
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1.899 on 18 degrees of freedom
#> Multiple R-squared:  0.1613,	Adjusted R-squared:  0.1147 
#> F-statistic: 3.463 on 1 and 18 DF,  p-value: 0.07919

# The sleep data are paired (the same 10 patients took both drugs)
t.test(sleep$extra[sleep$group == 2], sleep$extra[sleep$group == 1], paired = TRUE)

#> 
#> 	Paired t-test
#> 
#> data:  sleep$extra[sleep$group == 2] and sleep$extra[sleep$group == 1]
#> t = 4.0621, df = 9, p-value = 0.002833
#> alternative hypothesis: true mean difference is not equal to 0
#> 95 percent confidence interval:
#>  0.7001142 2.4598858
#> sample estimates:
#> mean difference 
#>            1.58

summary(lm(extra ~ group + ID, data = sleep))   # same t as the paired test

#> 
#> Call:
#> lm(formula = extra ~ group + ID, data = sleep)
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#> -1.510 -0.215  0.000  0.215  1.510 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)   
#> (Intercept)   0.5100     0.6450   0.791  0.44946   
#> group2        1.5800     0.3890   4.062  0.00283 **
#> ID2          -1.7000     0.8697  -1.955  0.08235 . 
#> ID3          -0.8500     0.8697  -0.977  0.35395   
#> ID4          -1.8500     0.8697  -2.127  0.06232 . 
#> ID5          -1.4000     0.8697  -1.610  0.14193   
#> ID6           2.6000     0.8697   2.989  0.01522 * 
#> ID7           3.3000     0.8697   3.794  0.00425 **
#> ID8          -0.1000     0.8697  -0.115  0.91099   
#> ID9           1.0000     0.8697   1.150  0.27987   
#> ID10          1.4000     0.8697   1.610  0.14193   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.8697 on 9 degrees of freedom
#> Multiple R-squared:  0.912,	Adjusted R-squared:  0.8142 
#> F-statistic: 9.328 on 10 and 9 DF,  p-value: 0.001254

# More than two groups: ANOVA = lm
aggregate(yield ~ block, npk, mean)

#>   block  yield
#> 1     1 54.025
#> 2     2 57.450
#> 3     3 60.775
#> 4     4 50.125
#> 5     5 50.525
#> 6     6 56.350

summary(aov(yield ~ block, data = npk))

#>             Df Sum Sq Mean Sq F value Pr(>F)  
#> block        5  343.3   68.66   2.318 0.0861 .
#> Residuals   18  533.1   29.62                 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(lm(yield ~ block, data = npk))          # same F and p at the bottom

#> 
#> Call:
#> lm(formula = yield ~ block, data = npk)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -7.2250 -3.4937 -0.5375  2.1062 11.8750 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   54.025      2.721  19.855 1.09e-13 ***
#> block2         3.425      3.848   0.890   0.3852    
#> block3         6.750      3.848   1.754   0.0964 .  
#> block4        -3.900      3.848  -1.013   0.3243    
#> block5        -3.500      3.848  -0.910   0.3751    
#> block6         2.325      3.848   0.604   0.5532    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 5.442 on 18 degrees of freedom
#> Multiple R-squared:  0.3917,	Adjusted R-squared:  0.2228 
#> F-statistic: 2.318 on 5 and 18 DF,  p-value: 0.08607

# Non-normal outcomes
set.seed(12345)
g1 = rpois(500, 2.5); g2 = rpois(500, 3.5)
d = data.frame(dv = c(g1, g2), groups = rep(c(1, 0), each = 500))
t.test(dv ~ groups, data = d)

#> 
#> 	Welch Two Sample t-test
#> 
#> data:  dv by groups
#> t = 7.7038, df = 989.92, p-value = 3.198e-14
#> alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
#> 95 percent confidence interval:
#>  0.6081435 1.0238565
#> sample estimates:
#> mean in group 0 mean in group 1 
#>           3.456           2.640

summary(lm(dv ~ groups, data = d))

#> 
#> Call:
#> lm(formula = dv ~ groups, data = d)
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#> -3.456 -1.456 -0.456  1.360  5.360 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   3.4560     0.0749  46.143  < 2e-16 ***
#> groups       -0.8160     0.1059  -7.704 3.18e-14 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1.675 on 998 degrees of freedom
#> Multiple R-squared:  0.05613,	Adjusted R-squared:  0.05518 
#> F-statistic: 59.35 on 1 and 998 DF,  p-value: 3.176e-14

# Contingency tables: a log-linear (Poisson) model
ctbl = data.frame(freq = c(65, 45, 20, 40), city = c("B", "B", "T", "T"),
                  pet = c("Cat", "Dog", "Cat", "Dog"))
summary(glm(freq ~ city * pet, family = poisson, data = ctbl))

#> 
#> Call:
#> glm(formula = freq ~ city * pet, family = poisson, data = ctbl)
#> 
#> Coefficients:
#>              Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)    4.1744     0.1240  33.655  < 2e-16 ***
#> cityT         -1.1787     0.2557  -4.609 4.04e-06 ***
#> petDog        -0.3677     0.1939  -1.896  0.05793 .  
#> cityT:petDog   1.0609     0.3356   3.161  0.00157 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for poisson family taken to be 1)
#> 
#>     Null deviance:  2.5378e+01  on 3  degrees of freedom
#> Residual deviance: -1.3767e-14  on 0  degrees of freedom
#> AIC: 30.036
#> 
#> Number of Fisher Scoring iterations: 3

log((65 * 40) / (45 * 20))                      # the interaction equals the log odds ratio

#> [1] 1.060872

# ... or a logistic regression on the disaggregated data (one row per person)
people = ctbl[rep(1:4, ctbl$freq), c("city", "pet")]
coef(glm(I(pet == "Dog") ~ city, family = binomial, data = people))

#> (Intercept)       cityT 
#>  -0.3677248   1.0608720

Takeaways

The pooled t-test is OLS with one dummy variable; ANOVA is OLS with $k-1$ dummy variables.
The paired t-test is OLS with an additional dummy variable for each participant.
The $\chi^2$-test has a regression counterpart, the log-linear model, whose interaction is the log odds ratio.
Learning regression well replaces a flowchart of tests.

Part III · Regression (OLS) · Chapter 9

Degrees of Freedom and Over-Fitting

The hook

Memorizing last year's exam

A student memorizes the answers to last year's exam word for word. On a practice test made of last year's questions, the student scores 100%. On this year's exam, the student fails. The student learned the noise, i.e., the particular questions, rather than the signal, i.e., the subject. Statistical models can make exactly the same mistake.

The weird question

What is a degree of freedom? Is it possible to fit a regression model with $R^2=1$? What is over-fitting, and why should we avoid it?

Degrees of freedom

Degrees of freedom ($df$) are the number of independent pieces of information that remain after estimating the parameters. In a regression, $df=N-P$, where $N$ is the sample size and $P$ is the number of estimated parameters. The intercept counts: a model with two predictors and 100 observations has $df=100-3=97$, not 98. Degrees of freedom indicate how much independent information goes into the estimate of the error variance. More $df$ means more precise estimates and more powerful tests.

Three ways to reach $R^2=1$

A just-identified model. With 3 observations and 3 parameters, $df=0$. The regression line passes through every point, $R^2=1$, and no standard error can be computed. In R, lm reports NA for the standard errors. Such a system has a unique solution; with $df<0$, it has infinitely many.
A deterministic model. If the data contain no random error, e.g., $y=2x_1+1.5x_2+1.1x_3$ exactly, then the true model has $R^2=1$. This is a perfect model, not an over-fitted one.
An over-fitted model. A model with too many parameters for the data extracts some of the residual variation (the noise) as if it were structure. It fits the sample perfectly and new data poorly.

In my original notes, I created a factor with 50 levels that was unrelated to any predictor. Regressing $y$ on this factor in the first 50 observations gave $R^2=1$. Predicting the other 50 observations with the same model gave a squared correlation close to 0. The widget below shows the same phenomenon with polynomials.

Play · The memorizing model

Blue dots: training data (last year's exam). Grey dots: new data (this year's exam). Green dashed curve: the true relationship. Red curve: a polynomial fitted to the training data.

Try this: raise the degree one step at a time. The training $R^2$ always rises, but the test $R^2$ peaks and then collapses.
Set the degree to $n-1$. Now $df=0$, and the curve passes through every training point.
Increase the training sample size. Over-fitting starts later, because each parameter is backed by more information.

Therefore, a high $R^2$ on the data used to fit a model is weak evidence of a good model. Cross-validation, i.e., evaluating a model on data it has not seen, is the standard remedy, and adjusted $R^2$ offers a crude penalty for the number of parameters.

Common traps

Forgetting that the intercept uses a degree of freedom.
Celebrating a high training $R^2$ from a model with many predictors and few observations.
Confusing a deterministic model with an over-fitted one.

R Corner: code and results

set.seed(12345)
x1 = rnorm(100, 5, 1); x2 = rnorm(100, 1, 1); x3 = rnorm(100, 2, 2)
y = 2 * x1 + 1.5 * x2 + 1.1 * x3 + rnorm(100)
data = data.frame(y, x1, x2, x3)
summary(lm(y ~ x1 + x2, data = data))                     # df = 100 - 3: the intercept counts

summary(lm(y ~ x1 + x2, data = data[sample(100, 3), ]))   # df = 0: R^2 = 1, no SEs
summary(lm(y ~ x1 + x2, data = data[sample(100, 4), ]))   # df = 1

summary(lm(y ~ x1 + x2 + x3, data))$r.squared             # the true model with noise
y2 = 2 * x1 + 1.5 * x2 + 1.1 * x3                         # no noise: a deterministic model
summary(lm(y2 ~ x1 + x2 + x3))$r.squared                  # R^2 = 1 (R warns about a perfect fit)

# Over-fitting: a factor with 50 levels that is unrelated to y
data$v1 = factor(c(1:50, 1:50))
model1 = lm(y ~ v1, data[1:50, ])
summary(model1)$r.squared                                 # 1 on the training data
y_predict = predict(model1, data[51:100, ])
cor(y[51:100], y_predict)^2                               # close to 0 on new data

# Polynomial over-fitting: training R^2 versus out-of-sample R^2
f = function(x) sin(2.5 * x) + 0.5 * x
train = data.frame(x = runif(12, -1, 1)); train$y = f(train$x) + rnorm(12, 0, 0.3)
test  = data.frame(x = runif(200, -0.95, 0.95)); test$y = f(test$x) + rnorm(200, 0, 0.3)
oos = function(m) 1 - sum((test$y - predict(m, test))^2) / sum((test$y - mean(test$y))^2)
for (k in c(1, 3, 6, 9, 11)) {
  m = lm(y ~ poly(x, k), data = train)
  cat("degree", k, " train R2 =", round(summary(m)$r.squared, 3), " test R2 =", round(oos(m), 3), "\n")
}

set.seed(12345)
x1 = rnorm(100, 5, 1); x2 = rnorm(100, 1, 1); x3 = rnorm(100, 2, 2)
y = 2 * x1 + 1.5 * x2 + 1.1 * x3 + rnorm(100)
data = data.frame(y, x1, x2, x3)
summary(lm(y ~ x1 + x2, data = data))                     # df = 100 - 3: the intercept counts

#> 
#> Call:
#> lm(formula = y ~ x1 + x2, data = data)
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#> -5.209 -1.423  0.253  1.122  6.299 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   2.1535     0.9793   2.199   0.0302 *  
#> x1            2.0541     0.1825  11.256  < 2e-16 ***
#> x2            1.3818     0.2012   6.869  6.2e-10 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 2.013 on 97 degrees of freedom
#> Multiple R-squared:  0.6645,	Adjusted R-squared:  0.6575 
#> F-statistic: 96.04 on 2 and 97 DF,  p-value: < 2.2e-16

summary(lm(y ~ x1 + x2, data = data[sample(100, 3), ]))   # df = 0: R^2 = 1, no SEs

#> 
#> Call:
#> lm(formula = y ~ x1 + x2, data = data[sample(100, 3), ])
#> 
#> Residuals:
#> ALL 3 residuals are 0: no residual degrees of freedom!
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)
#> (Intercept)  16.6964        NaN     NaN      NaN
#> x1           -0.2948        NaN     NaN      NaN
#> x2            0.7595        NaN     NaN      NaN
#> 
#> Residual standard error: NaN on 0 degrees of freedom
#> Multiple R-squared:      1,	Adjusted R-squared:    NaN 
#> F-statistic:   NaN on 2 and 0 DF,  p-value: NA

summary(lm(y ~ x1 + x2, data = data[sample(100, 4), ]))   # df = 1

#> 
#> Call:
#> lm(formula = y ~ x1 + x2, data = data[sample(100, 4), ])
#> 
#> Residuals:
#>       47       41       86        2 
#>  0.29329  1.27410  0.07655 -1.64394 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)
#> (Intercept)  17.1542     7.9452   2.159    0.276
#> x1           -0.5013     1.4370  -0.349    0.786
#> x2            4.0842     2.8720   1.422    0.390
#> 
#> Residual standard error: 2.102 on 1 degrees of freedom
#> Multiple R-squared:  0.7896,	Adjusted R-squared:  0.3688 
#> F-statistic: 1.876 on 2 and 1 DF,  p-value: 0.4587

summary(lm(y ~ x1 + x2 + x3, data))$r.squared             # the true model with noise

#> [1] 0.9293098

y2 = 2 * x1 + 1.5 * x2 + 1.1 * x3                         # no noise: a deterministic model
summary(lm(y2 ~ x1 + x2 + x3))$r.squared                  # R^2 = 1 (R warns about a perfect fit)

#> [1] 1

# Over-fitting: a factor with 50 levels that is unrelated to y
data$v1 = factor(c(1:50, 1:50))
model1 = lm(y ~ v1, data[1:50, ])
summary(model1)$r.squared                                 # 1 on the training data

#> [1] 1

y_predict = predict(model1, data[51:100, ])
cor(y[51:100], y_predict)^2                               # close to 0 on new data

#> [1] 0.002394023

# Polynomial over-fitting: training R^2 versus out-of-sample R^2
f = function(x) sin(2.5 * x) + 0.5 * x
train = data.frame(x = runif(12, -1, 1)); train$y = f(train$x) + rnorm(12, 0, 0.3)
test  = data.frame(x = runif(200, -0.95, 0.95)); test$y = f(test$x) + rnorm(200, 0, 0.3)
oos = function(m) 1 - sum((test$y - predict(m, test))^2) / sum((test$y - mean(test$y))^2)
for (k in c(1, 3, 6, 9, 11)) {
  m = lm(y ~ poly(x, k), data = train)
  cat("degree", k, " train R2 =", round(summary(m)$r.squared, 3), " test R2 =", round(oos(m), 3), "\n")
}

#> degree 1  train R2 = 0.938  test R2 = 0.869 
#> degree 3  train R2 = 0.964  test R2 = 0.913 
#> degree 6  train R2 = 0.979  test R2 = 0.582 
#> degree 9  train R2 = 0.988  test R2 = -811.102 
#> degree 11  train R2 = 1  test R2 = -350144.5

Takeaways

$df=N-P$, and the intercept is one of the $P$ parameters.
$R^2=1$ can mean a just-identified model, a deterministic model, or an over-fitted model.
Over-fitting means fitting noise; it shows up as a gap between training and test performance.
Judge predictive models on data they have not seen.

Part III · Regression (OLS) · Chapter 10

Interpreting and Comparing Regression Coefficients

The hook

Pricing a Hong Kong flat

An estate agent models flat prices (in HK$ million) using the building's age, the flat's size, whether it has a sea view, and its district: the New Territories, Kowloon, or Hong Kong Island. The client asks three questions. How much is a sea view worth? Is a flat on Hong Kong Island worth more than one in Kowloon? Does size matter more than age? Each question requires a different reading of the regression output. (The numbers are invented for teaching.)

The weird question

How do we interpret regression coefficients in terms of direction, magnitude, significance, and form? How do we interpret categorical predictors? How do we compare coefficients?

Four things to read from a coefficient

Direction: the sign. Older buildings are cheaper; larger flats are more expensive.
Magnitude: the expected change in $y$ for a one-unit change in $x$, holding the other predictors constant.
Significance: whether the estimate is distinguishable from 0, given its SE.
Form: whether the relationship is linear, or whether a curve or an interaction (Chapter 14) is needed.

Binary and categorical predictors

The coefficient of a binary predictor, such as sea view, is the difference in $y$ between flats with and without a sea view, holding the other predictors constant. A categorical predictor with $k$ levels enters as $k-1$ dummy variables, and each coefficient compares one level with the reference level. With the New Territories as the reference, the coefficient of Hong Kong Island is the difference between Hong Kong Island and the New Territories, not between Hong Kong Island and the other two districts combined. Changing the reference changes the coefficients, but not the model's predictions.

Note also that the raw group means differ from the coefficients, because the raw means do not hold age, size, and view constant.

Comparing coefficients

Standardized coefficients ($\beta$), which are the coefficients obtained after converting all variables to z-scores, put predictors on a common scale. However, a larger $|\beta|$ is not yet evidence of a larger effect; that claim needs a significance test, which requires the SE of the difference. By the full variance sum law:

$$\mathrm{Var}(c_1b_1+c_2b_2)=c_1^2\,\mathrm{Var}(b_1)+c_2^2\,\mathrm{Var}(b_2)+2c_1c_2\,\mathrm{Cov}(b_1,b_2)$$

The variances and covariances come from the variance-covariance matrix of the coefficients (vcov in R). For a difference $b_1-b_2$, $c_1=1$ and $c_2=-1$, so the covariance term is subtracted. To test whether the Hong Kong Island premium is twice the Kowloon premium, $c_1=1$ and $c_2=-2$, so the covariance term is $-4\,\mathrm{Cov}$.

Note on the original notes

An earlier version of these notes added the covariance term ($+4\,\mathrm{Cov}$) when testing $b_{x_4=3}=2b_{x_4=2}$. The correct sign is negative. The comparison of $|b_{x_2}|$ and $|b_{x_1}|$ used $+2\,\mathrm{Cov}$, which happens to be correct because $b_{x_1}$ is negative, so $|b_{x_2}|-|b_{x_1}|=b_{x_2}+b_{x_1}$.

Play · The flat-price lab

Reference district Standardize all continuous variables

Test a linear combination: $b_A - k\,b_B = 0$

A k B

Try this: set A = Hong Kong Island, k = 2, and B = Kowloon. The true values are 8 and 4, so the test should not be significant.
Change the reference district to Kowloon. The district coefficients change, but the other coefficients do not.
Tick Standardize. The coefficients and SEs change, but the t values of the continuous predictors stay the same.

Does standardization change the variance-covariance matrix? Yes. The coefficients and their variances are rescaled, but the t statistics are not. The implication is that significance does not depend on units, whereas comparisons of standardized coefficients must use the variance-covariance matrix of the standardized model.

Common traps

Reading a dummy coefficient as a comparison with all other categories.
Concluding that one effect is larger than another without testing the difference.
Ignoring the covariance term when testing combinations of coefficients.

R Corner: code and results

set.seed(12345)
x1 = rnorm(1000, 10, 2)                       # building age
x2 = rnorm(1000, 7, 3) + 0.5 * x1             # size (hundreds of square feet)
x3 = rbinom(1000, 1, 0.3)                     # sea view
x4 = sample(1:3, 1000, replace = TRUE, prob = c(0.2, 0.5, 0.3))   # 1 = NT, 2 = Kowloon, 3 = HK Island
y = -2.2 * x1 + 2.2 * x2 + 3.3 * x3 + 4 * (x4 == 2) + 8 * (x4 == 3) + rnorm(1000)
data = data.frame(y, x1, x2, x3, x4)

summary(m <- lm(y ~ x1 + x2 + x3 + factor(x4), data = data))
aggregate(y ~ x4, data, mean)                 # raw means are not the coefficients
summary(lm(y ~ x1 + x2 + x3 + relevel(factor(x4), ref = "2"), data = data))

library(lm.beta)
lm.beta(m)                                    # standardized coefficients

# Var(c1*b1 + c2*b2) = c1^2 Var(b1) + c2^2 Var(b2) + 2 c1 c2 Cov(b1, b2)
vc = vcov(m)
test_combo = function(cvec, value = 0) {
  est = sum(cvec * coef(m)) - value
  se = sqrt(drop(t(cvec) %*% vc %*% cvec))
  c(estimate = est, se = se, z = est / se, p = 2 * pnorm(-abs(est / se)))
}
cv = function(a) { v = setNames(rep(0, length(coef(m))), names(coef(m))); v[names(a)] = a; v }
test_combo(cv(c(x2 = 1, x1 = 1)))                           # |b_x2| - |b_x1|, since b_x1 < 0
test_combo(cv(c(x1 = 1)), value = -2.2)                     # is b_x1 different from -2.2?
test_combo(cv(c("factor(x4)3" = 1, "factor(x4)2" = -2)))    # is b_HKIsland = 2 * b_Kowloon?

# Standardizing changes the variance-covariance matrix, but not the t values
z = function(v) (v - mean(v)) / sd(v)
m1 = lm(y ~ x1 + x2, data = data)
m2 = lm(z(y) ~ z(x1) + z(x2), data = data)
vcov(m1); vcov(m2)
summary(m1)$coef[-1, "t value"]; summary(m2)$coef[-1, "t value"]

set.seed(12345)
x1 = rnorm(1000, 10, 2)                       # building age
x2 = rnorm(1000, 7, 3) + 0.5 * x1             # size (hundreds of square feet)
x3 = rbinom(1000, 1, 0.3)                     # sea view
x4 = sample(1:3, 1000, replace = TRUE, prob = c(0.2, 0.5, 0.3))   # 1 = NT, 2 = Kowloon, 3 = HK Island
y = -2.2 * x1 + 2.2 * x2 + 3.3 * x3 + 4 * (x4 == 2) + 8 * (x4 == 3) + rnorm(1000)
data = data.frame(y, x1, x2, x3, x4)

summary(m <- lm(y ~ x1 + x2 + x3 + factor(x4), data = data))

#> 
#> Call:
#> lm(formula = y ~ x1 + x2 + x3 + factor(x4), data = data)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -3.5214 -0.6106 -0.0408  0.6337  2.4490 
#> 
#> Coefficients:
#>              Estimate Std. Error  t value Pr(>|t|)    
#> (Intercept)  0.070260   0.182046    0.386      0.7    
#> x1          -2.213362   0.016082 -137.632   <2e-16 ***
#> x2           2.202264   0.009968  220.938   <2e-16 ***
#> x3           3.329405   0.066496   50.070   <2e-16 ***
#> factor(x4)2  3.991623   0.080790   49.407   <2e-16 ***
#> factor(x4)3  8.160564   0.087449   93.318   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.9514 on 994 degrees of freedom
#> Multiple R-squared:  0.9848,	Adjusted R-squared:  0.9847 
#> F-statistic: 1.286e+04 on 5 and 994 DF,  p-value: < 2.2e-16

aggregate(y ~ x4, data, mean)                 # raw means are not the coefficients

#>   x4         y
#> 1  1  4.859102
#> 2  2  9.347056
#> 3  3 12.745883

summary(lm(y ~ x1 + x2 + x3 + relevel(factor(x4), ref = "2"), data = data))

#> 
#> Call:
#> lm(formula = y ~ x1 + x2 + x3 + relevel(factor(x4), ref = "2"), 
#>     data = data)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -3.5214 -0.6106 -0.0408  0.6337  2.4490 
#> 
#> Coefficients:
#>                                  Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)                      4.061883   0.171422   23.70   <2e-16 ***
#> x1                              -2.213362   0.016082 -137.63   <2e-16 ***
#> x2                               2.202264   0.009968  220.94   <2e-16 ***
#> x3                               3.329405   0.066496   50.07   <2e-16 ***
#> relevel(factor(x4), ref = "2")1 -3.991623   0.080790  -49.41   <2e-16 ***
#> relevel(factor(x4), ref = "2")3  4.168941   0.069078   60.35   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.9514 on 994 degrees of freedom
#> Multiple R-squared:  0.9848,	Adjusted R-squared:  0.9847 
#> F-statistic: 1.286e+04 on 5 and 994 DF,  p-value: < 2.2e-16

library(lm.beta)
lm.beta(m)                                    # standardized coefficients

#> 
#> Call:
#> lm(formula = y ~ x1 + x2 + x3 + factor(x4), data = data)
#> 
#> Standardized Coefficients::
#> (Intercept)          x1          x2          x3 factor(x4)2 factor(x4)3 
#>          NA  -0.5748724   0.9226980   0.1963364   0.2596385   0.4901162

# Var(c1*b1 + c2*b2) = c1^2 Var(b1) + c2^2 Var(b2) + 2 c1 c2 Cov(b1, b2)
vc = vcov(m)
test_combo = function(cvec, value = 0) {
  est = sum(cvec * coef(m)) - value
  se = sqrt(drop(t(cvec) %*% vc %*% cvec))
  c(estimate = est, se = se, z = est / se, p = 2 * pnorm(-abs(est / se)))
}
cv = function(a) { v = setNames(rep(0, length(coef(m))), names(coef(m))); v[names(a)] = a; v }
test_combo(cv(c(x2 = 1, x1 = 1)))                           # |b_x2| - |b_x1|, since b_x1 < 0

#>    estimate          se           z           p 
#> -0.01109777  0.01571636 -0.70612833  0.48010833

test_combo(cv(c(x1 = 1)), value = -2.2)                     # is b_x1 different from -2.2?

#>    estimate          se           z           p 
#> -0.01336202  0.01608172 -0.83088276  0.40603986

test_combo(cv(c("factor(x4)3" = 1, "factor(x4)2" = -2)))    # is b_HKIsland = 2 * b_Kowloon?

#>  estimate        se         z         p 
#> 0.1773176 0.1222716 1.4501945 0.1470043

# Standardizing changes the variance-covariance matrix, but not the t values
z = function(v) (v - mean(v)) / sd(v)
m1 = lm(y ~ x1 + x2, data = data)
m2 = lm(z(y) ~ z(x1) + z(x2), data = data)
vcov(m1); vcov(m2)

#>              (Intercept)            x1            x2
#> (Intercept)  0.348172375 -0.0241544583 -0.0077906226
#> x1          -0.024154458  0.0032096708 -0.0006891507
#> x2          -0.007790623 -0.0006891507  0.0012334403

#>               (Intercept)         z(x1)         z(x2)
#> (Intercept)  1.903551e-04  9.297389e-20 -3.086428e-20
#> z(x1)        9.297389e-20  2.165202e-04 -7.499345e-05
#> z(x2)       -3.086428e-20 -7.499345e-05  2.165202e-04

summary(m1)$coef[-1, "t value"]; summary(m2)$coef[-1, "t value"]

#>        x1        x2 
#> -40.04186  62.13996

#>     z(x1)     z(x2) 
#> -40.04186  62.13996

Takeaways

Read each coefficient for direction, magnitude, significance, and form, holding other predictors constant.
Dummy coefficients compare each category with the reference category.
Test differences between coefficients with the variance-covariance matrix, minding the sign of the covariance term.
Standardization rescales coefficients and their SEs, but not their t values.

Part III · Regression (OLS) · Chapter 11

Why the I.I.D. Assumption Matters

The hook

Marbles from one bag

Suppose I want to know whether a new teaching method improves exam scores. I compare 10 classrooms that use the new method with 10 that use the old one, with 30 students in each. That looks like 600 observations. However, students in the same classroom share a teacher, a timetable, and a noisy air conditioner. In terms of information, I may have closer to 20 observations than 600.

The weird question

Why is the I.I.D. assumption essential?

A regression has a random part and a systematic part

A regression model can be written as a combination of two parts:

$$Y_i\sim\mathsf{Norm}(\mu_i,\sigma)\quad\text{(random part)},\qquad \mu_i=\beta_0+\beta_1X_i\quad\text{(systematic part)}$$

I.I.D. refers to observations that are independent and identically distributed. Think of drawing marbles from a bag:

Independent: each draw does not depend in any way on other draws. Repeated measures of the same person, or students from the same classroom, violate this condition.
Identically distributed: every draw comes from the same bag, i.e., the same distribution. The shape of the distribution does not matter, despite claims that it must be Normal; what matters is that all observations share it, including the same variance (Chapter 16 relaxes this).

What goes wrong without independence

The SE formulas of OLS assume that every observation brings new information. When observations are clustered, they partly repeat each other, so the SEs are too small, and tests reject too often. The degree of clustering is measured by the intraclass correlation (ICC), which means the share of the variance that lies between clusters. The variance inflation is approximately the design effect:

$$\text{DEFF}=1+(m-1)\,\text{ICC}$$

where $m$ is the cluster size. With 30 students per classroom and an ICC of .2, $\text{DEFF}=6.8$: 600 students carry about as much information as 88 independent students.

Play · Classrooms that look like a big sample

The teaching method has no effect at all. Each change runs 2,000 simulated studies and counts how often each analysis wrongly finds a significant effect.

Try this: set the ICC to 0. Both analyses reject about 5% of the time.
Raise the ICC to .1 and the class size to 50. The student-level analysis finds "effects" in a large share of studies.

Remedies include analyzing cluster means, using cluster-robust standard errors, and fitting multilevel (mixed-effects) models. The common principle is that the analysis must respect the unit at which the treatment varies.

Common traps

Counting repeated measures or clustered observations as independent.
Believing that I.I.D. requires Normal data.
Checking only the coefficients: clustering usually leaves estimates unbiased but makes SEs wrong.

R Corner: code and results

set.seed(12345)
one_study = function(K = 10, m = 30, icc = 0.2) {
  class = rep(1:(2 * K), each = m)
  treat = rep(rep(0:1, each = K), each = m)          # no true treatment effect
  u = rnorm(2 * K, 0, sqrt(icc))[class]              # shared classroom effect
  y = u + rnorm(2 * K * m, 0, sqrt(1 - icc))
  naive = summary(lm(y ~ treat))$coef["treat", "Pr(>|t|)"]
  means = aggregate(y ~ class + treat, FUN = mean)
  cluster = summary(lm(y ~ treat, data = means))$coef["treat", "Pr(>|t|)"]
  c(naive = naive, cluster = cluster)
}
p = replicate(1000, one_study())
rowMeans(p < .05)                                    # naive: far above .05; cluster means: about .05
1 + (30 - 1) * 0.2                                   # the design effect

set.seed(12345)
one_study = function(K = 10, m = 30, icc = 0.2) {
  class = rep(1:(2 * K), each = m)
  treat = rep(rep(0:1, each = K), each = m)          # no true treatment effect
  u = rnorm(2 * K, 0, sqrt(icc))[class]              # shared classroom effect
  y = u + rnorm(2 * K * m, 0, sqrt(1 - icc))
  naive = summary(lm(y ~ treat))$coef["treat", "Pr(>|t|)"]
  means = aggregate(y ~ class + treat, FUN = mean)
  cluster = summary(lm(y ~ treat, data = means))$coef["treat", "Pr(>|t|)"]
  c(naive = naive, cluster = cluster)
}
p = replicate(1000, one_study())
rowMeans(p < .05)                                    # naive: far above .05; cluster means: about .05

#>   naive cluster 
#>    0.45    0.05

1 + (30 - 1) * 0.2                                   # the design effect

#> [1] 6.8

Takeaways

I.I.D. means independent draws from the same distribution; it does not mean Normal.
Dependence makes OLS standard errors too small and tests too eager.
The design effect $1+(m-1)\,\text{ICC}$ shows how much information clustering removes.
Analyze at the level at which the treatment varies, or use models that account for clustering.

Part III · Regression (OLS) · Chapter 12

Skewed Outcomes, Normal Residuals

The hook

Incomes and noise

Household incomes are famously right-skewed, and many students log-transform any skewed outcome before running a regression, because "OLS requires a Normal dependent variable." It does not. A skewed outcome can arise entirely from skewed predictors, while the random part of the model is perfectly Normal.

The weird question

What if the residuals are normally distributed, but the dependent variable is not? Is this possible?

The assumption is about the random part

Recall the model from Chapter 11: $Y_i\sim\mathsf{Norm}(\mu_i,\sigma)$ with $\mu_i=X_i\beta$. The Normal assumption concerns the distribution of $Y$ given the predictors, i.e., the errors. Because $\mu_i$ varies across observations, the marginal distribution of $Y$ can have any shape. For example, if $y=7x_1+2x_2+\varepsilon$, where $x_1$ and $x_2$ are strongly skewed and $\varepsilon$ is Normal, then $y$ is skewed, the residuals are Normal, and OLS recovers the coefficients 7 and 2.

Even the Normality of the errors matters mainly for small-sample inference. With large samples, the CLT makes the coefficients approximately Normal whatever the error distribution, provided that the errors are I.I.D. with finite variance.

Play · Where does the skew come from?

Predictors Errors

Try this: with skewed predictors and Normal errors, $y$ is heavily skewed, but the residuals are bell-shaped, and the estimates sit close to 7 and 2.
Switch to Normal predictors and skewed errors. Now the residuals are skewed, yet with 1,000 observations the estimates are still close to the truth.

Therefore, a log transformation should be chosen for substantive reasons, e.g., because effects are multiplicative (a 10% raise rather than a HK$1,000 raise), not merely because a histogram of $y$ is skewed. A transformation changes the model, not just the shape of a histogram.

Common traps

Testing the normality of $y$ instead of examining the residuals.
Log-transforming every skewed outcome by habit.
Forgetting that the CLT protects inference about coefficients in large samples.

R Corner: code and results

set.seed(12345)
x1 = rexp(1000, 1.5)
x2 = rweibull(1000, 0.5)
y = 7 * x1 + 2 * x2 + rnorm(1000)
hist(y)                       # strongly skewed
m = lm(y ~ x1 + x2)
summary(m)                    # coefficients close to 7 and 2
hist(resid(m))                # bell-shaped
qqnorm(resid(m)); qqline(resid(m))

set.seed(12345)
x1 = rexp(1000, 1.5)
x2 = rweibull(1000, 0.5)
y = 7 * x1 + 2 * x2 + rnorm(1000)
hist(y)                       # strongly skewed

m = lm(y ~ x1 + x2)
summary(m)                    # coefficients close to 7 and 2

#> 
#> Call:
#> lm(formula = y ~ x1 + x2)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2.7345 -0.6871  0.0325  0.6707  3.3957 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 0.009148   0.047044   0.194    0.846    
#> x1          7.004881   0.049612 141.192   <2e-16 ***
#> x2          1.982469   0.007676 258.259   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1.001 on 997 degrees of freedom
#> Multiple R-squared:  0.9889,	Adjusted R-squared:  0.9888 
#> F-statistic: 4.422e+04 on 2 and 997 DF,  p-value: < 2.2e-16

hist(resid(m))                # bell-shaped

qqnorm(resid(m)); qqline(resid(m))

Takeaways

OLS assumes Normal errors, not a Normal dependent variable.
A skewed outcome can be produced by skewed predictors with Normal noise.
Check residuals, and remember that large samples make the error distribution less important.
Transform $y$ only when the transformed model makes substantive sense.

Part III · Regression (OLS) · Chapter 13

Effect Sizes and Unique Contributions

The hook

Who deserves credit for the group project?

Two students worked on a group project, and their contributions overlapped heavily: both wrote parts of the literature review. If I ask how much the grade would drop without each student, the two answers do not add up to the total grade. Regression faces the same problem when predictors are correlated: the variance that they explain jointly belongs to neither of them alone.

The weird question

Why is effect size important? Is it true that larger coefficients indicate larger effect sizes? How can we measure the unique effect of a predictor without the confounding effects of other predictors, and can stepwise regression solve the problem?

Effect sizes

Significance indicates whether an effect is distinguishable from zero; effect size indicates how large it is. With a large sample, trivial effects become significant. There are many effect size measures. A raw coefficient depends on units, so a larger coefficient does not imply a larger effect. The standardized $\beta$ can serve as an effect size, but it is hard to interpret. $R^2$, the share of variance explained, is the most intuitive, but it describes the whole model.

Unique contributions: $sr^2$

The squared semi-partial correlation ($sr^2$) of a predictor is its unique contribution to the variance of $y$, i.e., the increase in $R^2$ when that predictor is added to a model that already contains all the others. A tempting shortcut is to add predictors one at a time and record each increase in $R^2$. However, with correlated predictors, the answer depends on the order of entry: the first predictor entered takes the credit for the shared variance. The correct "one-less" approach is:

Fit the full model and record $R^2_{\text{full}}$.
Remove one predictor $x_i$ and record $R^2_{(i)}$.
Compute $sr^2_i=R^2_{\text{full}}-R^2_{(i)}$.

The sum of the $sr^2$ values usually falls short of $R^2_{\text{full}}$, because the shared variance is credited to nobody. However, the sum can also exceed $R^2_{\text{full}}$. This happens under suppression, e.g., when two predictors are positively correlated but have opposite effects, so that each one hides the other. In my original flat-price data, age and size are positively correlated, but age lowers the price and size raises it, which is why the $sr^2$ values summed to more than the total $R^2$.

Play · Dividing the credit

$y=b_1x_1+2x_2+1\cdot x_3+\varepsilon$, where $x_3$ is unrelated to $x_1$ and $x_2$.

Try this: with a correlation of .6 and $b_1=2$, the increase in $R^2$ for $x_1$ is much larger when it enters first than when it enters last.
Set $b_1=-2$ with a positive correlation. The shared part turns negative: suppression.
Set the correlation to 0. Order no longer matters, and the $sr^2$ values add up to $R^2$.

Can stepwise regression solve the problem? No. Stepwise procedures add or drop predictors based on significance, so their results depend on the order of entry and on chance correlations in the sample. They answer a question about prediction, not about the unique effect of each predictor.

Common traps

Comparing raw coefficients measured in different units.
Reporting order-dependent $\Delta R^2$ values as unique contributions.
Expecting unique contributions to add up to the total $R^2$.

R Corner: code and results

set.seed(12345)
x1 = rnorm(1000, 10, 2)
x2 = rnorm(1000, 7, 3) + 0.5 * x1
x3 = rbinom(1000, 1, 0.3)
x4 = sample(1:3, 1000, replace = TRUE, prob = c(0.2, 0.5, 0.3))
y = -2.2 * x1 + 2.2 * x2 + 3.3 * x3 + 4 * (x4 == 2) + 8 * (x4 == 3) + rnorm(1000)
data = data.frame(y, x1, x2, x3, x4)
r2 = function(f) summary(lm(f, data))$r.squared

# Hierarchical entry: the answer depends on the order
base = r2(y ~ x3 + factor(x4))
c(x1_entered_first = r2(y ~ x1 + x3 + factor(x4)) - base,
  x1_entered_last  = r2(y ~ x1 + x2 + x3 + factor(x4)) - r2(y ~ x2 + x3 + factor(x4)))

# The one-less approach: squared semi-partial correlations
full = r2(y ~ x1 + x2 + x3 + factor(x4))
sr2 = c(x1 = full - r2(y ~ x2 + x3 + factor(x4)),
        x2 = full - r2(y ~ x1 + x3 + factor(x4)),
        x3 = full - r2(y ~ x1 + x2 + factor(x4)),
        x4 = full - r2(y ~ x1 + x2 + x3))
sr2; sum(sr2); full                    # sum(sr2) > R^2: suppression
cor(data$x1, data$x2)                  # positively correlated, opposite-signed effects

set.seed(12345)
x1 = rnorm(1000, 10, 2)
x2 = rnorm(1000, 7, 3) + 0.5 * x1
x3 = rbinom(1000, 1, 0.3)
x4 = sample(1:3, 1000, replace = TRUE, prob = c(0.2, 0.5, 0.3))
y = -2.2 * x1 + 2.2 * x2 + 3.3 * x3 + 4 * (x4 == 2) + 8 * (x4 == 3) + rnorm(1000)
data = data.frame(y, x1, x2, x3, x4)
r2 = function(f) summary(lm(f, data))$r.squared

# Hierarchical entry: the answer depends on the order
base = r2(y ~ x3 + factor(x4))
c(x1_entered_first = r2(y ~ x1 + x3 + factor(x4)) - base,
  x1_entered_last  = r2(y ~ x1 + x2 + x3 + factor(x4)) - r2(y ~ x2 + x3 + factor(x4)))

#> x1_entered_first  x1_entered_last 
#>        0.0650897        0.2901754

# The one-less approach: squared semi-partial correlations
full = r2(y ~ x1 + x2 + x3 + factor(x4))
sr2 = c(x1 = full - r2(y ~ x2 + x3 + factor(x4)),
        x2 = full - r2(y ~ x1 + x3 + factor(x4)),
        x3 = full - r2(y ~ x1 + x2 + factor(x4)),
        x4 = full - r2(y ~ x1 + x2 + x3))
sr2; sum(sr2); full                    # sum(sr2) > R^2: suppression

#>         x1         x2         x3         x4 
#> 0.29017540 0.74775853 0.03840323 0.13782565

#> [1] 1.214163

#> [1] 0.9847733

cor(data$x1, data$x2)                  # positively correlated, opposite-signed effects

#> [1] 0.3463577

Takeaways

Effect size, not significance, indicates whether an effect matters.
$sr^2$ is the drop in $R^2$ when one predictor is removed from the full model.
Shared variance belongs to no single predictor; under suppression, the "shared" part is negative.
Stepwise regression does not identify unique effects.

Part III · Regression (OLS) · Chapter 14

Main Effects in Models with Interactions

The hook

Coffee and night owls

Does coffee improve exam performance? Suppose that each cup helps morning people a lot but helps night owls only a little. A model with an interaction term reports a "coefficient of coffee," but this number is the effect of coffee for morning people only. So what is the main effect of coffee?

The weird question

Is it possible to estimate the main effects from a regression model with interaction terms?

Coefficients in an interaction model are conditional

Consider the model

$$y=b_0+b_1\,\text{coffee}+b_2\,\text{owl}+b_3\,(\text{coffee}\times\text{owl})+\varepsilon$$

For morning people ($\text{owl}=0$), the slope of coffee is $b_1$. For night owls ($\text{owl}=1$), it is $b_1+b_3$. Therefore, $b_1$ is not the main effect of coffee; it is the conditional effect when the moderator equals 0. If 0 is not a meaningful value of the moderator (e.g., an age of 0), $b_1$ may not describe anyone at all.

The average marginal effect

A marginal effect is the change in the prediction for an infinitesimal change in the focal variable. The average marginal effect (AME) averages it over the sample, which yields one number per predictor that can be considered an adjusted regression coefficient. In this linear model,

$$\text{AME}_{\text{coffee}}=b_1+b_3\,\overline{\text{owl}}$$

If the moderator is centered before forming the product, $b_1$ equals the AME directly. A model without the interaction term estimates a similar average effect. In other words, the main effect lies between the two conditional slopes, weighted by the shares of the groups.

Play · Coffee for larks and owls

Center the moderator (owl − mean)

Try this: with the defaults, the coefficient of coffee (about 4) describes morning people only, while the AME is about 2.8.
Tick Center the moderator. The coefficient of coffee becomes the AME, but the interaction coefficient does not change.
Increase the share of night owls. The AME moves toward the night owls' slope.

In nonlinear models, such as logistic regression, the relationship between coefficients and marginal effects is more complicated, as Chapter 17 shows. The R function margins computes AMEs for both kinds of models.

Common traps

Reporting $b_1$ from an interaction model as the "main effect."
Interpreting a conditional effect at a moderator value of 0 that nobody has.
Concluding that a lower-order term is unimportant because it is non-significant in the interaction model.

R Corner: code and results

set.seed(12345)
coffee = rnorm(1000, 3, 1)                       # cups per day
owl = rbinom(1000, 1, 0.4)                       # 1 = night owl
score = 60 + 4 * coffee + 3 * owl - 3 * coffee * owl + rnorm(1000, 0, 4)
d = data.frame(score, coffee, owl)

summary(full <- lm(score ~ coffee * owl, data = d))   # b_coffee: the slope for morning people
coef(lm(score ~ coffee, data = d[d$owl == 0, ]))["coffee"]
coef(lm(score ~ coffee, data = d[d$owl == 1, ]))["coffee"]

library(margins)
summary(margins(full))                           # AME of coffee: about 4 - 3 * 0.4 = 2.8
coef(lm(score ~ coffee + owl, data = d))["coffee"]   # model without the interaction

d$owl_c = d$owl - mean(d$owl)                    # centering the moderator
coef(lm(score ~ coffee * owl_c, data = d))["coffee"] # now b_coffee equals the AME

set.seed(12345)
coffee = rnorm(1000, 3, 1)                       # cups per day
owl = rbinom(1000, 1, 0.4)                       # 1 = night owl
score = 60 + 4 * coffee + 3 * owl - 3 * coffee * owl + rnorm(1000, 0, 4)
d = data.frame(score, coffee, owl)

summary(full <- lm(score ~ coffee * owl, data = d))   # b_coffee: the slope for morning people

#> 
#> Call:
#> lm(formula = score ~ coffee * owl, data = d)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -12.7018  -2.6481   0.0226   2.6941  11.2118 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  59.8858     0.5248 114.118  < 2e-16 ***
#> coffee        3.9342     0.1641  23.968  < 2e-16 ***
#> owl           3.2608     0.8063   4.044 5.66e-05 ***
#> coffee:owl   -2.9374     0.2512 -11.692  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.923 on 996 degrees of freedom
#> Multiple R-squared:  0.5258,	Adjusted R-squared:  0.5244 
#> F-statistic: 368.2 on 3 and 996 DF,  p-value: < 2.2e-16

coef(lm(score ~ coffee, data = d[d$owl == 0, ]))["coffee"]

#> coffee 
#> 3.9342

coef(lm(score ~ coffee, data = d[d$owl == 1, ]))["coffee"]

#>    coffee 
#> 0.9967873

library(margins)
summary(margins(full))                           # AME of coffee: about 4 - 3 * 0.4 = 2.8

#>  factor     AME     SE        z      p   lower   upper
#>  coffee  2.7622 0.1245  22.1925 0.0000  2.5182  3.0061
#>     owl -5.6871 0.2533 -22.4541 0.0000 -6.1836 -5.1907

coef(lm(score ~ coffee + owl, data = d))["coffee"]   # model without the interaction

#>   coffee 
#> 2.680273

d$owl_c = d$owl - mean(d$owl)                    # centering the moderator
coef(lm(score ~ coffee * owl_c, data = d))["coffee"] # now b_coffee equals the AME

#>   coffee 
#> 2.762172

Takeaways

In an interaction model, a lower-order coefficient is the effect when the other variable equals 0.
The average marginal effect summarizes the effect across the sample.
Centering the moderator turns $b_1$ into the AME in linear models.
Report conditional slopes for meaningful values of the moderator.

Part IV · Generalized linear models · Chapter 15

Generalized Linear Models

The hook

Admissions, broken yarn, and skipped classes

A graduate school records whether each applicant was admitted (yes or no). A textile mill counts how many times its yarn breaks on each loom. A high school counts the days each student was absent, and a household survey records the share of income spent on food. None of these outcomes looks like a bell curve: one is binary, two are counts, and one is a proportion. Should all four be analyzed with OLS anyway?

The weird question

How do we deal with non-normally distributed dependent variables? How do we interpret the coefficients?

A random part and a systematic part

Chapter 11 described OLS as a data-generating process: each $y_i$ is drawn from a Normal distribution whose mean varies systematically with the predictors. A generalized linear model (GLM) keeps this structure but lets the random part come from another distribution. More specifically, a GLM has three components:

a distribution for $y_i$, e.g., Bernoulli, Poisson, negative binomial, or Beta;
a linear predictor $\eta_i=\beta_0+\beta_1x_i$, which is the familiar regression equation;
a link function, which connects the linear predictor to the parameter of the distribution and keeps that parameter in its allowed range.

Therefore, the choice of model follows from the process that generates the dependent variable, not from the shape of a histogram alone.

Outcome	Distribution	Link	$e^{\beta_1}$ means
Binary (admitted or not)	Bernoulli	logit	odds ratio
Count (yarn breaks)	Poisson	log	rate ratio
Over-dispersed count (days absent)	negative binomial	log	rate ratio
Proportion (food share)	Beta	logit	ratio of the odds $\mu/(1-\mu)$

Logistic regression

$$Y_i\sim\mathsf{Bern}(\pi_i),\qquad \mathrm{logit}(\pi_i)=\log\frac{\pi_i}{1-\pi_i}=\beta_0+\beta_1X_i$$

The Bernoulli distribution has a single parameter $\pi_i$, the probability of the event. The link is the logit, i.e., the log of the odds (not the log of the probability). $\beta_0$ is the log odds of the event when $X_i=0$, and $e^{\beta_0}$ is the odds. $\beta_1=\log(\text{odds}_{x+1})-\log(\text{odds}_x)$ is the change in the log odds for a one-unit increase in $x$, so $e^{\beta_1}=\text{odds}_{x+1}/\text{odds}_x$ is an odds ratio.

Poisson and negative binomial regression

$$Y_i\sim\mathsf{Pois}(\lambda_i),\qquad \log(\lambda_i)=\beta_0+\beta_1X_i$$

$\beta_0$ is the log of the expected count when $X_i=0$, and $e^{\beta_1}=\lambda_{x+1}/\lambda_x$ is a rate ratio: the expected count is multiplied by $e^{\beta_1}$ for each one-unit increase in $x$. The Poisson distribution forces the variance to equal the mean. Real counts often vary more than that, which is called over-dispersion. The negative binomial distribution adds a parameter, so that $\mathrm{Var}(Y)=\mu+\mu^2/\theta$. As $\theta\to\infty$, the negative binomial distribution approaches the Poisson. The coefficients are interpreted in the same way as in Poisson regression.

Play · One line, three scales

The linear predictor $\eta=\beta_0+\beta_1x$ is always a straight line. The link function bends it into a probability, a count, or a mean. The dots are 150 simulated observations.

Family

Try this: with the logit link, compare the slope of the curve at $x=-2$, 0, and 2. The odds ratio is constant, but the change in probability is not.
Switch to Poisson. A constant rate ratio means that each step in $x$ multiplies the count, so the curve grows faster and faster.

Reading the admissions model

The admissions data (from UCLA's statistical consulting group) contain 400 applicants, their GRE score, their GPA, and the prestige rank of their undergraduate institution (1 = highest, 4 = lowest). Because odds ratios are not intuitive for most readers, I prefer to report predicted probabilities and average marginal effects, which are both on the probability scale.

Play · The admissions calculator

The model is fitted in your browser by iteratively reweighted least squares, the same algorithm that R's glm uses. Move the sliders to describe an applicant.

Try this: raise the GPA from 3.0 to 4.0. The odds of admission are multiplied by the same odds ratio for every rank, but the gain in probability is largest where the probability is closest to 50%.

Play · Poisson versus negative binomial

Both distributions have the same mean $\mu$. The bars show the negative binomial distribution; the red dots show the Poisson distribution.

Try this: with $\mu=6$ and $\theta=1$, the negative binomial distribution piles up many zeros and a long right tail, like the absence data. Increase $\theta$ until the two distributions coincide.

Beta regression for proportions

$$Y_i\sim\mathsf{Beta}(\mu_i,\phi),\qquad \mathrm{logit}(\mu_i)=\beta_0+\beta_1X_i$$

The Beta distribution lives between 0 and 1 (the end points are excluded), so it suits proportions such as the share of income spent on food. It also uses a logit link. The precision parameter $\phi$ means that the higher $\phi$, the smaller the variance for a given mean $\mu$.

Note on the original notes

Three corrections. First, average marginal effects are differences in probabilities, so the rank 2 effect means that the admission probability is about 16 percentage points lower than for rank 1, not "15.7% lower." Likewise, one more household member is associated with a food share that is about 2.4 percentage points higher, holding income constant; this is an association, not a causal effect. Second, the absence data store the program type as a labelled number, so it must be converted to a factor before it enters the model. Third, R's $\theta$ is not the dispersion parameter itself but its inverse: Stata and SAS report $\alpha=1/\theta$.

Common traps

Reading a logit coefficient as a change in probability.
Reporting an average marginal effect as a percentage change instead of a percentage-point change.
Fitting a Poisson model to over-dispersed counts, which makes the standard errors too small.

R Corner: code and results

# Logistic regression: graduate admissions
mydata = read.csv("https://stats.oarc.ucla.edu/stat/data/binary.csv")
head(mydata)
mylogit = glm(admit ~ gre + gpa + factor(rank), data = mydata, family = binomial)
summary(mylogit)
exp(coef(mylogit))                                 # odds ratios

library(jtools)
summ(mylogit)                                      # includes pseudo-R^2
library(margins)
summary(margins(mylogit, variables = "rank"))      # AMEs: differences in probability
library(ggeffects)
ggpredict(mylogit, terms = "rank")                 # predicted probabilities

# Poisson regression: yarn breaks
mypois = glm(breaks ~ wool + tension, data = warpbreaks, family = poisson)
summary(mypois)
exp(coef(mypois))                                  # rate ratios
sum(residuals(mypois, type = "pearson")^2) / df.residual(mypois)   # well above 1: over-dispersion

# Negative binomial regression: days absent
library(MASS)
library(haven)
dat = read_stata("https://stats.oarc.ucla.edu/stat/stata/dae/nb_data.dta")
dat$prog = factor(dat$prog, levels = 1:3, labels = c("General", "Academic", "Vocational"))
hist(dat$daysabs)
c(mean = mean(dat$daysabs), variance = var(dat$daysabs))
summary(mynegbin <- glm.nb(daysabs ~ math + prog, data = dat))
1 / mynegbin$theta                                 # the dispersion parameter reported by Stata

# Beta regression: share of income spent on food
library(betareg)
data("FoodExpenditure", package = "betareg")
mybeta = betareg(I(food / income) ~ income + persons, data = FoodExpenditure)
summary(mybeta)
summary(margins(mybeta))                           # AMEs in proportion units

# Logistic regression: graduate admissions
mydata = read.csv("https://stats.oarc.ucla.edu/stat/data/binary.csv")
head(mydata)

#>   admit gre  gpa rank
#> 1     0 380 3.61    3
#> 2     1 660 3.67    3
#> 3     1 800 4.00    1
#> 4     1 640 3.19    4
#> 5     0 520 2.93    4
#> 6     1 760 3.00    2

mylogit = glm(admit ~ gre + gpa + factor(rank), data = mydata, family = binomial)
summary(mylogit)

#> 
#> Call:
#> glm(formula = admit ~ gre + gpa + factor(rank), family = binomial, 
#>     data = mydata)
#> 
#> Coefficients:
#>                Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)   -3.989979   1.139951  -3.500 0.000465 ***
#> gre            0.002264   0.001094   2.070 0.038465 *  
#> gpa            0.804038   0.331819   2.423 0.015388 *  
#> factor(rank)2 -0.675443   0.316490  -2.134 0.032829 *  
#> factor(rank)3 -1.340204   0.345306  -3.881 0.000104 ***
#> factor(rank)4 -1.551464   0.417832  -3.713 0.000205 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 499.98  on 399  degrees of freedom
#> Residual deviance: 458.52  on 394  degrees of freedom
#> AIC: 470.52
#> 
#> Number of Fisher Scoring iterations: 4

exp(coef(mylogit))                                 # odds ratios

#>   (Intercept)           gre           gpa factor(rank)2 factor(rank)3 
#>     0.0185001     1.0022670     2.2345448     0.5089310     0.2617923 
#> factor(rank)4 
#>     0.2119375

library(jtools)
summ(mylogit)                                      # includes pseudo-R^2

#> MODEL INFO:
#> Observations: 400
#> Dependent Variable: admit
#> Type: Generalized linear model
#>   Family: binomial 
#>   Link function: logit 
#> 
#> MODEL FIT:
#> χ²(5) = 41.46, p = 0.00
#> Pseudo-R² (Cragg-Uhler) = 0.14
#> Pseudo-R² (McFadden) = 0.08
#> AIC = 470.52, BIC = 494.47 
#> 
#> Standard errors:MLE
#> --------------------------------------------------
#>                        Est.   S.E.   z val.      p
#> ------------------- ------- ------ -------- ------
#> (Intercept)           -3.99   1.14    -3.50   0.00
#> gre                    0.00   0.00     2.07   0.04
#> gpa                    0.80   0.33     2.42   0.02
#> factor(rank)2         -0.68   0.32    -2.13   0.03
#> factor(rank)3         -1.34   0.35    -3.88   0.00
#> factor(rank)4         -1.55   0.42    -3.71   0.00
#> --------------------------------------------------

library(margins)
summary(margins(mylogit, variables = "rank"))      # AMEs: differences in probability

#>  factor     AME     SE       z      p   lower   upper
#>   rank2 -0.1566 0.0736 -2.1272 0.0334 -0.3009 -0.0123
#>   rank3 -0.2872 0.0733 -3.9170 0.0001 -0.4309 -0.1435
#>   rank4 -0.3212 0.0802 -4.0052 0.0001 -0.4784 -0.1640

library(ggeffects)
ggpredict(mylogit, terms = "rank")                 # predicted probabilities

#> # Predicted probabilities of admit
#> 
#> rank | Predicted |     95% CI
#> -----------------------------
#>    1 |      0.51 | 0.38, 0.64
#>    2 |      0.35 | 0.27, 0.43
#>    3 |      0.22 | 0.15, 0.30
#>    4 |      0.18 | 0.11, 0.30
#> 
#> Adjusted for:
#> * gre = 580.00
#> * gpa =   3.39

# Poisson regression: yarn breaks
mypois = glm(breaks ~ wool + tension, data = warpbreaks, family = poisson)
summary(mypois)

#> 
#> Call:
#> glm(formula = breaks ~ wool + tension, family = poisson, data = warpbreaks)
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  3.69196    0.04541  81.302  < 2e-16 ***
#> woolB       -0.20599    0.05157  -3.994 6.49e-05 ***
#> tensionM    -0.32132    0.06027  -5.332 9.73e-08 ***
#> tensionH    -0.51849    0.06396  -8.107 5.21e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for poisson family taken to be 1)
#> 
#>     Null deviance: 297.37  on 53  degrees of freedom
#> Residual deviance: 210.39  on 50  degrees of freedom
#> AIC: 493.06
#> 
#> Number of Fisher Scoring iterations: 4

exp(coef(mypois))                                  # rate ratios

#> (Intercept)       woolB    tensionM    tensionH 
#>  40.1235380   0.8138425   0.7251908   0.5954198

sum(residuals(mypois, type = "pearson")^2) / df.residual(mypois)   # well above 1: over-dispersion

#> [1] 4.261522

# Negative binomial regression: days absent
library(MASS)
library(haven)
dat = read_stata("https://stats.oarc.ucla.edu/stat/stata/dae/nb_data.dta")
dat$prog = factor(dat$prog, levels = 1:3, labels = c("General", "Academic", "Vocational"))
hist(dat$daysabs)

c(mean = mean(dat$daysabs), variance = var(dat$daysabs))

#>      mean  variance 
#>  5.955414 49.518773

summary(mynegbin <- glm.nb(daysabs ~ math + prog, data = dat))

#> 
#> Call:
#> glm.nb(formula = daysabs ~ math + prog, data = dat, init.theta = 1.032713156, 
#>     link = log)
#> 
#> Coefficients:
#>                 Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)     2.615265   0.197460  13.245  < 2e-16 ***
#> math           -0.005993   0.002505  -2.392   0.0167 *  
#> progAcademic   -0.440760   0.182610  -2.414   0.0158 *  
#> progVocational -1.278651   0.200720  -6.370 1.89e-10 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for Negative Binomial(1.0327) family taken to be 1)
#> 
#>     Null deviance: 427.54  on 313  degrees of freedom
#> Residual deviance: 358.52  on 310  degrees of freedom
#> AIC: 1741.3
#> 
#> Number of Fisher Scoring iterations: 1
#> 
#> 
#>               Theta:  1.033 
#>           Std. Err.:  0.106 
#> 
#>  2 x log-likelihood:  -1731.258

1 / mynegbin$theta                                 # the dispersion parameter reported by Stata

#> [1] 0.9683231

# Beta regression: share of income spent on food
library(betareg)
data("FoodExpenditure", package = "betareg")
mybeta = betareg(I(food / income) ~ income + persons, data = FoodExpenditure)
summary(mybeta)

#> 
#> Call:
#> betareg(formula = I(food/income) ~ income + persons, data = FoodExpenditure)
#> 
#> Quantile residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2.5328 -0.4599  0.1698  0.6416  1.7733 
#> 
#> Coefficients (mean model with logit link):
#>              Estimate Std. Error z value Pr(>|z|)    
#> (Intercept) -0.622548   0.223854  -2.781 0.005418 ** 
#> income      -0.012299   0.003036  -4.052 5.09e-05 ***
#> persons      0.118462   0.035341   3.352 0.000802 ***
#> 
#> Phi coefficients (precision model with identity link):
#>       Estimate Std. Error z value Pr(>|z|)    
#> (phi)    35.61       8.08   4.407 1.05e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
#> 
#> Type of estimator: ML (maximum likelihood)
#> Log-likelihood: 45.33 on 4 Df
#> Pseudo R-squared: 0.3878
#> Number of iterations: 28 (BFGS) + 4 (Fisher scoring)

summary(margins(mybeta))                           # AMEs in proportion units

#>   factor     AME     SE       z      p   lower   upper
#>   income -0.0025 0.0006 -4.1266 0.0000 -0.0037 -0.0013
#>  persons  0.0239 0.0070  3.3945 0.0007  0.0101  0.0377

Takeaways

A GLM combines a distribution for $y$, a linear predictor, and a link function.
Exponentiated coefficients are odds ratios (logit) or rate ratios (log).
Report predicted probabilities or average marginal effects, in percentage points.
Check for over-dispersion before trusting a Poisson model.

Part IV · Generalized linear models · Chapter 16

Modeling the Variance

The hook

Salaries that fan in

Among young graduates, salaries vary enormously: some join investment banks, others start their own companies, and others take part-time jobs. Among employees close to retirement, salaries are more predictable. An OLS regression of salary on age describes how the average salary changes with age, but it says nothing about this change in spread. Can a regression model the spread itself?

The weird question

Can we model the variance of the dependent variable instead of the mean? For example, the variance of the salaries of older individuals is smaller than that of younger individuals, given equal salary levels.

A second regression equation, for $\sigma$

OLS assumes that every observation has the same standard deviation $\sigma$. Nothing prevents me from giving $\sigma$ its own systematic part:

$$Y_i\sim\mathsf{Norm}(\mu_i,\sigma_i),\qquad \mu_i=\beta_0+\beta_1X_i,\qquad \log(\sigma_i)=\gamma_0+\gamma_1X_i$$

The log link keeps $\sigma_i$ positive. $\gamma_1$ describes how the spread changes with $x$: $e^{\gamma_1}$ is the factor by which the standard deviation is multiplied for each one-unit increase in $x$. Such a model is called a location-scale model, which means that predictors may shift both the center (location) and the spread (scale) of the distribution. In R, the gamlss package (Generalized Additive Models for Location, Scale and Shape) fits it with two formulas, one for the mean and one for $\sigma$. The package also allows distributions other than the Normal.

Play · Salaries by age

Simulated monthly salaries (thousand HK$): $\mu=20+0.5\,\text{age}$ and $\log\sigma=\gamma_0+\gamma_1\,\text{age}$. The green band is the true $\mu\pm2\sigma$; the dashed red lines are the OLS band, which assumes a constant $\sigma$.

Try this: with $\gamma_1=-0.02$, the OLS band is too narrow for young employees and too wide for older ones, even though the OLS slope is fine.
Set $\gamma_1=0$. The two bands coincide, and the model reduces to OLS.

The widget also reports a quick two-step estimate: regress $\log(e_i^2)$, the log of the squared OLS residuals, on age. Because $\log(e_i^2)\approx 2\log\sigma_i+\text{noise}$, half of the slope estimates $\gamma_1$. A proper location-scale model estimates both equations simultaneously by maximum likelihood, which is what gamlss does.

Think about it

Why might a researcher care about the variance rather than the mean?

My answer

Many theories are about inequality, risk, or consensus rather than about averages. For example, a policy may leave the average income unchanged but compress the income distribution, and a news event may leave average opinion unchanged but polarize it. In these cases, the substantive effect is in $\sigma$, and a model of the mean alone would conclude that nothing happened.

Note on the original notes

The original equation modeled $\sigma_i$ as a linear function of the predictors and kept the mean fixed. A linear equation can predict a negative standard deviation, so gamlss (and this chapter) uses a log link for $\sigma$. The coefficients of the $\sigma$ equation are therefore on the log scale.

Common traps

Treating unequal spread only as a nuisance to be "fixed" with robust SEs, when it may be the finding.
Reading the $\sigma$ coefficients on the original scale instead of the log scale.

R Corner: code and results

library(gamlss)
head(mtcars)
summary(lm(mpg ~ hp + wt, data = mtcars))

# Only the standard deviation depends on the predictors (constant mean)
myglass1 = gamlss(mpg ~ 1, sigma.formula = ~ hp + wt, family = NO(),
                  data = mtcars, control = gamlss.control(trace = FALSE))
# Both the mean and the standard deviation depend on the predictors
myglass2 = gamlss(mpg ~ hp + wt, sigma.formula = ~ hp + wt, family = NO(),
                  data = mtcars, control = gamlss.control(trace = FALSE))
summary(myglass1)
summary(myglass2)                    # sigma coefficients are on the log scale

# Salaries whose spread shrinks with age
set.seed(12345)
age = runif(1000, 25, 65)
salary = 20 + 0.5 * age + rnorm(1000, 0, exp(2.5 - 0.02 * age))
d = data.frame(salary, age)
m_sal = gamlss(salary ~ age, sigma.formula = ~ age, family = NO(), data = d,
               control = gamlss.control(trace = FALSE))
summary(m_sal)                       # sigma: intercept near 2.5, age near -0.02

# A Bayesian alternative (slow; run at your own risk)
# library(brms)
# myBayes = brm(bf(mpg ~ hp + wt, sigma ~ hp + wt), data = mtcars, family = gaussian)
# summary(myBayes)

library(gamlss)
head(mtcars)

#>                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
#> Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
#> Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
#> Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
#> Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
#> Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
#> Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

summary(lm(mpg ~ hp + wt, data = mtcars))

#> 
#> Call:
#> lm(formula = mpg ~ hp + wt, data = mtcars)
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#> -3.941 -1.600 -0.182  1.050  5.854 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 37.22727    1.59879  23.285  < 2e-16 ***
#> hp          -0.03177    0.00903  -3.519  0.00145 ** 
#> wt          -3.87783    0.63273  -6.129 1.12e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 2.593 on 29 degrees of freedom
#> Multiple R-squared:  0.8268,	Adjusted R-squared:  0.8148 
#> F-statistic: 69.21 on 2 and 29 DF,  p-value: 9.109e-12

# Only the standard deviation depends on the predictors (constant mean)
myglass1 = gamlss(mpg ~ 1, sigma.formula = ~ hp + wt, family = NO(),
                  data = mtcars, control = gamlss.control(trace = FALSE))
# Both the mean and the standard deviation depend on the predictors
myglass2 = gamlss(mpg ~ hp + wt, sigma.formula = ~ hp + wt, family = NO(),
                  data = mtcars, control = gamlss.control(trace = FALSE))
summary(myglass1)

#> ******************************************************************
#> Family:  c("NO", "Normal") 
#> 
#> Call:  gamlss(formula = mpg ~ 1, sigma.formula = ~hp + wt,  
#>     family = NO(), data = mtcars, control = gamlss.control(trace = FALSE)) 
#> 
#> Fitting method: RS() 
#> 
#> ------------------------------------------------------------------
#> Mu link function:  identity
#> Mu Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  15.0982     0.3599   41.95   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------------------------------
#> Sigma link function:  log
#> Sigma Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  3.356603   0.355505   9.442 3.37e-10 ***
#> hp          -0.013336   0.003307  -4.033 0.000385 ***
#> wt           0.060041   0.146335   0.410 0.684713    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------------------------------
#> No. of observations in the fit:  32 
#> Degrees of Freedom for the fit:  4
#>       Residual Deg. of Freedom:  28 
#>                       at cycle:  13 
#>  
#> Global Deviance:     192.8026 
#>             AIC:     200.8026 
#>             SBC:     206.6656 
#> ******************************************************************

summary(myglass2)                    # sigma coefficients are on the log scale

#> ******************************************************************
#> Family:  c("NO", "Normal") 
#> 
#> Call:  gamlss(formula = mpg ~ hp + wt, sigma.formula = ~hp +  
#>     wt, family = NO(), data = mtcars, control = gamlss.control(trace = FALSE)) 
#> 
#> 
#> Fitting method: RS() 
#> 
#> ------------------------------------------------------------------
#> Mu link function:  identity
#> Mu Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 35.97612    2.47576  14.531 5.40e-14 ***
#> hp          -0.02782    0.00713  -3.901 0.000604 ***
#> wt          -3.69296    0.79155  -4.666 8.13e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------------------------------
#> Sigma link function:  log
#> Sigma Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)  
#> (Intercept)  1.224105   0.545910   2.242   0.0337 *
#> hp          -0.003228   0.002775  -1.163   0.2553  
#> wt           0.039559   0.224216   0.176   0.8613  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------------------------------
#> No. of observations in the fit:  32 
#> Degrees of Freedom for the fit:  6
#>       Residual Deg. of Freedom:  26 
#>                       at cycle:  8 
#>  
#> Global Deviance:     146.998 
#>             AIC:     158.998 
#>             SBC:     167.7924 
#> ******************************************************************

# Salaries whose spread shrinks with age
set.seed(12345)
age = runif(1000, 25, 65)
salary = 20 + 0.5 * age + rnorm(1000, 0, exp(2.5 - 0.02 * age))
d = data.frame(salary, age)
m_sal = gamlss(salary ~ age, sigma.formula = ~ age, family = NO(), data = d,
               control = gamlss.control(trace = FALSE))
summary(m_sal)                       # sigma: intercept near 2.5, age near -0.02

#> ******************************************************************
#> Family:  c("NO", "Normal") 
#> 
#> Call:  gamlss(formula = salary ~ age, sigma.formula = ~age, family = NO(),  
#>     data = d, control = gamlss.control(trace = FALSE)) 
#> 
#> Fitting method: RS() 
#> 
#> ------------------------------------------------------------------
#> Mu link function:  identity
#> Mu Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 20.05989    0.73034   27.47   <2e-16 ***
#> age          0.49716    0.01411   35.23   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------------------------------
#> Sigma link function:  log
#> Sigma Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  2.556955   0.092362   27.68   <2e-16 ***
#> age         -0.021121   0.001967  -10.74   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ------------------------------------------------------------------
#> No. of observations in the fit:  1000 
#> Degrees of Freedom for the fit:  4
#>       Residual Deg. of Freedom:  996 
#>                       at cycle:  3 
#>  
#> Global Deviance:     6027.161 
#>             AIC:     6035.161 
#>             SBC:     6054.792 
#> ******************************************************************

# A Bayesian alternative (slow; run at your own risk)
# library(brms)
# myBayes = brm(bf(mpg ~ hp + wt, sigma ~ hp + wt), data = mtcars, family = gaussian)
# summary(myBayes)

Takeaways

A location-scale model gives the standard deviation its own regression equation.
Use a log link for $\sigma$; $e^{\gamma_1}$ multiplies the standard deviation.
Differences in spread can be substantive findings, not just violations of assumptions.

Part IV · Generalized linear models · Chapter 17

Interactions in Logistic Regression

The hook

The bent probability scale

Suppose that practice ($x$) raises the chance of passing a driving test and that prior experience ($m$) does too. A logistic regression without any product term still implies that practice matters more for some learners than for others: an hour of practice changes little for someone who will almost surely pass anyway, but a lot for someone on the edge. On the probability scale, every logistic model contains a built-in interaction.

The weird question

What are the problems with using Z-statistics and their associated p-values for the coefficients of multiplicative terms when testing for statistical interactions in nonlinear models with categorical dependent variables, such as logistic regression?

Two scales, two questions

In a linear model, the coefficient of $x\times m$ is the change in the slope of $x$ for a one-unit change in $m$, and its t-test answers the question of whether the slopes differ. In a logistic model,

$$\mathrm{logit}(\pi)=\beta_0+\beta_1x+\beta_2m+\beta_3xm$$

the coefficient $\beta_3$ describes the interaction on the log-odds scale. However, most research questions are about probabilities. The marginal effect of $x$ on the probability is

$$\frac{\partial\pi}{\partial x}=(\beta_1+\beta_3m)\,\pi(1-\pi)$$

which depends on $m$ twice: through $\beta_3m$ and through $\pi(1-\pi)$, because $m$ also moves the baseline probability. Therefore, as Long and Mustillo (2021) argued, the equality of coefficients does not imply the equality of marginal effects. The z-test of $\beta_3$ tests an interaction on the log-odds scale, which may differ in size, and even in sign, from the interaction on the probability scale.

Play · The bent probability scale

$\mathrm{logit}(\pi)=0.25+\beta_1x+\beta_2m+\beta_3xm$, with $x,m\sim\mathsf{Norm}(0,1)$. Left: the marginal effect of $x$ on the probability, averaged over $x$, for each value of $m$. Right: the slope of $x$ on the log-odds scale, estimated from 1,000 simulated people, with a 95% band. The shading marks the Johnson-Neyman region, i.e., the values of $m$ for which the slope is significant.

Try this: set $\beta_3=0$ and $\beta_2=1.5$. The log-odds slope is flat across $m$ (no interaction), yet the marginal effect on the probability changes a lot with $m$.
Set $\beta_3=-0.5$ (the default). The interaction is significant on the log-odds scale, and the marginal effect changes sign within the observed range of $m$.

For this reason, I recommend three practices. First, plot predicted probabilities across values of the moderator (interact_plot or ggpredict). Second, report marginal effects at meaningful values of the moderator, and test their differences directly (margins with the at argument). Third, if the theory is about the log-odds scale, e.g., in a latent-variable interpretation, say so explicitly, and then the test of $\beta_3$ is the right one. The Johnson-Neyman interval reported by sim_slopes for a logistic model refers to the log-odds slope, not to the probability.

Note on the original notes

The original simulation set $y=1$ whenever the inverse logit of the linear predictor plus a Normal error exceeded .5. This rule generates a latent-variable model with a Normal error, which is a probit model with rescaled coefficients, rather than a logistic model with the stated coefficients. The R Corner draws $y$ with rbinom instead. The original also attributed the varying marginal effects to a "log" transformation; the transformation is the logit.

Common traps

Concluding "no interaction" from a non-significant product term in a logistic model, when the question is about probabilities.
Reporting one average marginal effect when the effect changes sign across the moderator.
Reading a Johnson-Neyman interval from a logistic model as a statement about probabilities.

R Corner: code and results

set.seed(12345)
x = rnorm(1000)
m = rnorm(1000)
y = rbinom(1000, 1, plogis(0.25 + 0.3 * x + 0.3 * m - 0.5 * x * m))
logit_fit = glm(y ~ x * m, family = binomial)
summary(logit_fit)                               # x:m tests the log-odds interaction

library(interactions)
interact_plot(logit_fit, pred = x, modx = m, interval = TRUE)   # predicted probabilities

library(margins)
summary(margins(logit_fit))                      # AMEs on the probability scale
summary(margins(logit_fit, at = list(m = c(-2.5, 0, 2.5)), variables = "x"))

library(ggeffects)
ggpredict(logit_fit, terms = "x")

sim_slopes(logit_fit, pred = x, modx = m, jnplot = TRUE)   # Johnson-Neyman (log-odds slopes)

# No product term in the data-generating process, but unequal marginal effects
y0 = rbinom(1000, 1, plogis(0.25 + 1.5 * x + 1.5 * m))
fit0 = glm(y0 ~ x * m, family = binomial)
summary(fit0)                                    # x:m is close to 0
summary(margins(fit0, at = list(m = c(-1.5, 0, 1.5)), variables = "x"))

set.seed(12345)
x = rnorm(1000)
m = rnorm(1000)
y = rbinom(1000, 1, plogis(0.25 + 0.3 * x + 0.3 * m - 0.5 * x * m))
logit_fit = glm(y ~ x * m, family = binomial)
summary(logit_fit)                               # x:m tests the log-odds interaction

#> 
#> Call:
#> glm(formula = y ~ x * m, family = binomial)
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  0.30582    0.06715   4.554 5.25e-06 ***
#> x            0.41779    0.07354   5.681 1.34e-08 ***
#> m            0.20446    0.07036   2.906  0.00366 ** 
#> x:m         -0.54440    0.08311  -6.550 5.76e-11 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 1364.9  on 999  degrees of freedom
#> Residual deviance: 1268.6  on 996  degrees of freedom
#> AIC: 1276.6
#> 
#> Number of Fisher Scoring iterations: 4

library(interactions)
interact_plot(logit_fit, pred = x, modx = m, interval = TRUE)   # predicted probabilities

library(margins)
summary(margins(logit_fit))                      # AMEs on the probability scale

#>  factor    AME     SE      z      p  lower  upper
#>       m 0.0410 0.0143 2.8653 0.0042 0.0130 0.0691
#>       x 0.0898 0.0142 6.3447 0.0000 0.0620 0.1175

summary(margins(logit_fit, at = list(m = c(-2.5, 0, 2.5)), variables = "x"))

#>  factor       m     AME     SE       z      p   lower   upper
#>       x -2.5000  0.2903 0.0178 16.3124 0.0000  0.2554  0.3252
#>       x  0.0000  0.0979 0.0160  6.1206 0.0000  0.0666  0.1293
#>       x  2.5000 -0.1809 0.0321 -5.6337 0.0000 -0.2438 -0.1179

library(ggeffects)
ggpredict(logit_fit, terms = "x")

#> # Predicted probabilities of y
#> 
#>  x | Predicted |     95% CI
#> ---------------------------
#> -3 |      0.27 | 0.19, 0.37
#> -2 |      0.36 | 0.29, 0.44
#> -1 |      0.47 | 0.42, 0.52
#>  0 |      0.57 | 0.54, 0.61
#>  1 |      0.68 | 0.63, 0.72
#>  2 |      0.76 | 0.70, 0.82
#>  3 |      0.83 | 0.76, 0.89
#>  4 |      0.88 | 0.81, 0.93
#> 
#> Adjusted for:
#> * m = -0.03

sim_slopes(logit_fit, pred = x, modx = m, jnplot = TRUE)   # Johnson-Neyman (log-odds slopes)

#> JOHNSON-NEYMAN INTERVAL
#> 
#> When m is OUTSIDE the interval [0.49, 1.17], the slope of x is p < .05.
#> 
#> Note: The range of observed values of m is [-3.33, 3.29]

#> SIMPLE SLOPES ANALYSIS
#> 
#> Slope of x when m = -1.03874572 (- 1 SD): 
#> 
#>   Est.   S.E.   z val.      p
#> ------ ------ -------- ------
#>   0.98   0.12     8.00   0.00
#> 
#> Slope of x when m = -0.03039248 (Mean): 
#> 
#>   Est.   S.E.   z val.      p
#> ------ ------ -------- ------
#>   0.43   0.07     5.87   0.00
#> 
#> Slope of x when m =  0.97796076 (+ 1 SD): 
#> 
#>    Est.   S.E.   z val.      p
#> ------- ------ -------- ------
#>   -0.11   0.10    -1.15   0.25

# No product term in the data-generating process, but unequal marginal effects
y0 = rbinom(1000, 1, plogis(0.25 + 1.5 * x + 1.5 * m))
fit0 = glm(y0 ~ x * m, family = binomial)
summary(fit0)                                    # x:m is close to 0

#> 
#> Call:
#> glm(formula = y0 ~ x * m, family = binomial)
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  0.17843    0.09188   1.942   0.0521 .  
#> x            1.51617    0.11496  13.189   <2e-16 ***
#> m            1.62720    0.11923  13.647   <2e-16 ***
#> x:m         -0.11409    0.11991  -0.951   0.3414    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 1381.11  on 999  degrees of freedom
#> Residual deviance:  860.97  on 996  degrees of freedom
#> AIC: 868.97
#> 
#> Number of Fisher Scoring iterations: 5

summary(margins(fit0, at = list(m = c(-1.5, 0, 1.5)), variables = "x"))

#>  factor       m    AME     SE       z      p  lower  upper
#>       x -1.5000 0.1815 0.0215  8.4529 0.0000 0.1394 0.2235
#>       x  0.0000 0.2658 0.0108 24.5433 0.0000 0.2446 0.2871
#>       x  1.5000 0.1106 0.0213  5.1928 0.0000 0.0689 0.1524

Takeaways

In logistic regression, the product term describes an interaction on the log-odds scale.
Marginal effects on the probability scale vary with the moderator even without a product term.
Plot predicted probabilities and report marginal effects at meaningful moderator values.

Part V · Causal inference · Chapter 18

Correlation, Causation, and the Backdoor

The hook

Ice cream and drowning

In a coastal city, days with high ice cream sales are also days with many drownings. Nobody believes that ice cream causes drowning: hot weather sends people both to ice cream stands and to the beach. Hot weather is a confounder, which means that it causes both the treatment and the outcome. Once the temperature is held constant, the association between ice cream and drowning disappears. But which variables should be held constant in less obvious cases?

The weird question

Everyone knows that correlation is not causation. When can regression coefficients be interpreted as causal effects? Is including a lagged explanatory variable a viable option?

Description versus causal inference

Regression models serve two different purposes. For description and prediction, the overall fit (e.g., $R^2$) matters. For causal inference, what matters is an unbiased estimate of one focal coefficient, the treatment effect. These goals are different. A model with a low $R^2$ can estimate a treatment effect accurately, and a model with a high $R^2$ can estimate it badly. For example, if the omitted causes of $y$ are independent of the treatment, leaving them out lowers $R^2$ but does not bias the treatment coefficient.

The backdoor criterion

Causal relationships can be drawn as a directed acyclic graph (DAG), in which arrows point from causes to effects. A backdoor path from the treatment $T$ to the outcome $Y$ is any path that starts with an arrow into $T$, such as $T\leftarrow\text{heat}\rightarrow Y$. Such paths transmit association that is not causation. Pearl's backdoor criterion states that a set of control variables $Z$ identifies the causal effect of $T$ on $Y$ if (1) no variable in $Z$ is a descendant of $T$, and (2) $Z$ blocks every backdoor path between $T$ and $Y$. One important implication is that it is not necessary to control for every cause of $Y$; it is enough to block every backdoor path once.

Play · The backdoor lab

Each scenario simulates 1,000 observations from the causal diagram. Dashed nodes are unobserved and cannot be controlled. Tick the variables to control, and compare the estimated effect of the treatment with the truth.

Scenario

Try this: in scenario 2, control for $x_1$ only, then for $x_2$ only. Either one blocks the single backdoor path, so both recover the true effect of 3.
In scenario 3, controlling for the lagged treatment fails, but using it as an instrument works. In scenarios 4 and 5, both strategies fail.

Is a lagged treatment a solution?

Suppose that the confounders are unobserved, and a researcher controls for last year's value of the treatment, $T_{t0}$. This strategy rarely works, because $T_{t0}$ does not lie on the backdoor path $T_{t1}\leftarrow x_{t1}\rightarrow Y$. However, if the confounders at the two time points are independent, $T_{t0}$ affects $Y$ only through $T_{t1}$, so it is a valid instrumental variable (IV), which means a variable that shifts the treatment but is unrelated to the outcome through any other path. Two-stage least squares (2SLS) then recovers the effect. In more realistic settings, the confounders persist over time ($x_{t0}\rightarrow x_{t1}$) or do not change at all. Then the lagged treatment is connected to $Y$ through the confounders, and both OLS and IV are biased.

Think about it

Why do researchers routinely control for demographics, such as age and gender?

My answer

Age and gender are exogenous: nothing else in the model causes them, so they can never be descendants of the treatment. Therefore, controlling for them never opens a new path. When they influence both the treatment and the outcome, they lie on backdoor paths, and controlling for them blocks these paths. When they influence only the outcome, controlling for them does not change the expected estimate but reduces the residual variance, which makes the estimate more precise. In other words, demographics are safe controls, which is not true for every variable (Chapters 20 and 21).

Note on the original notes

The original notes described $x_{t0}$ as the instrumental variable. The instrument is the lagged treatment $T_{t0}$; $x_{t0}$ is unobserved and only makes $T_{t0}$ vary. The instrument is valid only if $T_{t0}$ is unrelated to $x_{t1}$ and has no direct effect on $Y$.

Common traps

Judging a causal model by its $R^2$.
Believing that a lagged treatment automatically removes confounding.
Controlling for every available variable without drawing the causal diagram first.

R Corner: code and results

library(jtools)
set.seed(12345)

# Example 1: other causes of y independent of the treatment
x1 = rnorm(1000, 2, 1); x2 = rnorm(1000, 3, 1); Treatment = rnorm(1000, 5, 2)
y = 2 * x1 + 3 * Treatment + rnorm(1000)
summ(lm(y ~ x1 + x2 + Treatment), model.info = FALSE)
summ(lm(y ~ Treatment), model.info = FALSE)       # lower R^2, still unbiased

# Example 2: x1 -> x2 -> Treatment, x1 -> y
x1 = rnorm(1000, 2, 1)
x2 = 3.1 * x1 + rnorm(1000)
Treatment = -0.9 * x2 + rnorm(1000)
y = 2 * x1 + 3 * Treatment + rnorm(1000)
summ(lm(y ~ Treatment), model.info = FALSE)       # biased
m0 = lm(y ~ Treatment + x2)
m1 = lm(y ~ Treatment + x1)
export_summs(m0, m1)                              # both recover 3

# Example 3: lagged treatment with independent confounders
library(ivreg)
x_t0 = rnorm(1000, 1, 1.5); x_t1 = rnorm(1000, 2, 1.5)
Treatment_t0 = 1.5 * x_t0 + rnorm(1000, 3, 1.2)
Treatment_t1 = 2 * Treatment_t0 + 1.5 * x_t1 + rnorm(1000)
y = 2 * Treatment_t1 + 3 * x_t1 + rnorm(1000)    # the true effect is 2
coef(lm(y ~ Treatment_t1 + Treatment_t0))         # biased
summary(ivreg(y ~ Treatment_t1 | Treatment_t0))   # close to 2

# Example 4: persistent confounders (x_t0 -> x_t1)
x_t0 = rnorm(1000, 1, 1.5); x_t1 = 1.8 * x_t0 + rnorm(1000)
treatment_t0 = 1.5 * x_t0 + rnorm(1000, 3, 1.2)
treatment_t1 = 2 * treatment_t0 + 1.5 * x_t1 + rnorm(1000)
y = 2 * treatment_t1 + 3 * x_t1 + rnorm(1000)
coef(lm(y ~ treatment_t1 + treatment_t0))         # biased
coef(ivreg(y ~ treatment_t1 | treatment_t0))      # also biased

# Example 5: a time-invariant confounder
x = rnorm(1000, 3, 1.5)
treatment_t0 = 1.5 * x + rnorm(1000, 3, 1.2)
treatment_t1 = 2 * treatment_t0 + 1.5 * x + rnorm(1000)
y = 2 * treatment_t1 + 3 * x + rnorm(1000)
coef(lm(y ~ treatment_t1 + treatment_t0))         # biased

library(jtools)
set.seed(12345)

# Example 1: other causes of y independent of the treatment
x1 = rnorm(1000, 2, 1); x2 = rnorm(1000, 3, 1); Treatment = rnorm(1000, 5, 2)
y = 2 * x1 + 3 * Treatment + rnorm(1000)
summ(lm(y ~ x1 + x2 + Treatment), model.info = FALSE)

#> MODEL FIT:
#> F(3,996) = 13282.86, p = 0.00
#> R² = 0.98
#> Adj. R² = 0.98 
#> 
#> Standard errors:OLS
#> -----------------------------------------------
#>                     Est.   S.E.   t val.      p
#> ----------------- ------ ------ -------- ------
#> (Intercept)         0.13   0.14     0.97   0.33
#> x1                  1.97   0.03    65.32   0.00
#> x2                  0.00   0.03     0.13   0.90
#> Treatment           2.99   0.02   189.39   0.00
#> -----------------------------------------------

summ(lm(y ~ Treatment), model.info = FALSE)       # lower R^2, still unbiased

#> MODEL FIT:
#> F(1,998) = 6737.79, p = 0.00
#> R² = 0.87
#> Adj. R² = 0.87 
#> 
#> Standard errors:OLS
#> -----------------------------------------------
#>                     Est.   S.E.   t val.      p
#> ----------------- ------ ------ -------- ------
#> (Intercept)         4.24   0.19    22.09   0.00
#> Treatment           2.97   0.04    82.08   0.00
#> -----------------------------------------------

# Example 2: x1 -> x2 -> Treatment, x1 -> y
x1 = rnorm(1000, 2, 1)
x2 = 3.1 * x1 + rnorm(1000)
Treatment = -0.9 * x2 + rnorm(1000)
y = 2 * x1 + 3 * Treatment + rnorm(1000)
summ(lm(y ~ Treatment), model.info = FALSE)       # biased

#> MODEL FIT:
#> F(1,998) = 30971.55, p = 0.00
#> R² = 0.97
#> Adj. R² = 0.97 
#> 
#> Standard errors:OLS
#> -----------------------------------------------
#>                     Est.   S.E.   t val.      p
#> ----------------- ------ ------ -------- ------
#> (Intercept)         0.69   0.09     7.93   0.00
#> Treatment           2.40   0.01   175.99   0.00
#> -----------------------------------------------

m0 = lm(y ~ Treatment + x2)
m1 = lm(y ~ Treatment + x1)
export_summs(m0, m1)                              # both recover 3

#>                   ───────────────────────────────────────────
#>                                     Model 1       Model 2    
#>                                 ─────────────────────────────
#>                     (Intercept)      0.34 ***     -0.03      
#>                                     (0.08)        (0.07)     
#>                     Treatment        3.01 ***      3.00 ***  
#>                                     (0.04)        (0.02)     
#>                     x2               0.60 ***                
#>                                     (0.03)                   
#>                     x1                             2.00 ***  
#>                                                   (0.07)     
#>                                 ─────────────────────────────
#>                     N             1000          1000         
#>                     R2               0.98          0.98      
#>                   ───────────────────────────────────────────
#> *** p < 0.001; ** p < 0.01; * p < 0.05.
#> 
#> Column names: names, Model 1, Model 2

# Example 3: lagged treatment with independent confounders
library(ivreg)
x_t0 = rnorm(1000, 1, 1.5); x_t1 = rnorm(1000, 2, 1.5)
Treatment_t0 = 1.5 * x_t0 + rnorm(1000, 3, 1.2)
Treatment_t1 = 2 * Treatment_t0 + 1.5 * x_t1 + rnorm(1000)
y = 2 * Treatment_t1 + 3 * x_t1 + rnorm(1000)    # the true effect is 2
coef(lm(y ~ Treatment_t1 + Treatment_t0))         # biased

#>  (Intercept) Treatment_t1 Treatment_t0 
#>    0.9906754    3.6714621   -3.3180557

summary(ivreg(y ~ Treatment_t1 | Treatment_t0))   # close to 2

#> 
#> Call:
#> ivreg(formula = y ~ Treatment_t1 | Treatment_t0)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -13.3729  -3.1217  -0.1554   3.0333  13.8259 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   6.24878    0.37749   16.55   <2e-16 ***
#> Treatment_t1  1.99690    0.02902   68.82   <2e-16 ***
#> 
#> Diagnostic tests:
#>                  df1 df2 statistic p-value    
#> Weak instruments   1 998      4127  <2e-16 ***
#> Wu-Hausman         1 997      3360  <2e-16 ***
#> Sargan             0  NA        NA      NA    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 4.561 on 998 degrees of freedom
#> Multiple R-Squared: 0.8866,	Adjusted R-squared: 0.8865 
#> Wald test:  4736 on 1 and 998 DF,  p-value: < 2.2e-16

# Example 4: persistent confounders (x_t0 -> x_t1)
x_t0 = rnorm(1000, 1, 1.5); x_t1 = 1.8 * x_t0 + rnorm(1000)
treatment_t0 = 1.5 * x_t0 + rnorm(1000, 3, 1.2)
treatment_t1 = 2 * treatment_t0 + 1.5 * x_t1 + rnorm(1000)
y = 2 * treatment_t1 + 3 * x_t1 + rnorm(1000)
coef(lm(y ~ treatment_t1 + treatment_t0))         # biased

#>  (Intercept) treatment_t1 treatment_t0 
#>    -1.324403     3.709734    -2.945374

coef(ivreg(y ~ treatment_t1 | treatment_t0))      # also biased

#>  (Intercept) treatment_t1 
#>    -4.363371     2.839607

# Example 5: a time-invariant confounder
x = rnorm(1000, 3, 1.5)
treatment_t0 = 1.5 * x + rnorm(1000, 3, 1.2)
treatment_t1 = 2 * treatment_t0 + 1.5 * x + rnorm(1000)
y = 2 * treatment_t1 + 3 * x + rnorm(1000)
coef(lm(y ~ treatment_t1 + treatment_t0))         # biased

#>  (Intercept) treatment_t1 treatment_t0 
#>    -1.268138     3.003514    -1.242152

Takeaways

A causal estimate needs every backdoor path blocked, not a high $R^2$.
One control per backdoor path is enough; controls must not be descendants of the treatment.
A lagged treatment is a valid instrument only under strong and usually implausible conditions.

Part V · Causal inference · Chapter 19

Are Experiments Always Better?

The hook

A/B testing on a biased user pool

A learning app tests a new reminder feature. Users are randomly assigned to the new version (B) or the old version (A), so the test is a genuine randomized controlled trial (RCT). However, the feature helps women more than men, the users who volunteer for beta tests are mostly men, and some users in group B switch the reminders off. The randomization is perfect. Is the estimate?

The weird question

Are randomized controlled trials always better than observational methods for identifying causality?

What randomization buys, and what it does not

Randomization makes the treatment independent of every confounder, observed or not, so it blocks every backdoor path at once (Chapter 18). This is why RCTs are the benchmark for causal inference. However, an RCT rests on its own assumptions, and it does not answer every question:

Plausibility and ethics. Many causes cannot be manipulated, e.g., gender, or should not be, e.g., exposure to harm.
Randomization must succeed. In a small sample, the groups can still differ by chance.
Compliance. If some participants do not take the treatment they are assigned, the comparison of assigned groups, which is called the intention-to-treat (ITT) effect, is diluted.
Heterogeneous effects and external validity. If the effect differs across people, the experiment estimates the average effect in its sample. A non-representative sample estimates the wrong average for the population.

Play · The A/B test

The reminder raises the score by 2 points for women and by 1 point for men, i.e., $y=2D-1\,D\times\text{male}+\varepsilon$, where $D$ means using the feature. Half of the population is male, so the population average treatment effect (ATE) is 1.5. Assignment is always random.

Try this: set the share of men to 0.8. The experiment is still unbiased for its own sample, but it estimates 1.2, not the population effect of 1.5.
Lower the compliance rate to 0.5. The ITT estimate halves; dividing it by the compliance rate recovers the effect among those who use the feature.

The last estimate is the complier average causal effect (CACE), which is the effect among participants who take the treatment when assigned to it. It equals the ITT divided by the difference in treatment uptake between the groups, which is exactly an IV estimate with the random assignment as the instrument. Note that, in the widget, compliance is unrelated to gender; if men were also less likely to comply, the CACE would describe a population that is even more different from the target population.

In sum, RCTs are better at internal validity, i.e., estimating the effect for the people studied. Observational studies with large, representative samples can be better at external validity, i.e., generalizing to the population. The two designs complement each other.

Common traps

Generalizing an experimental effect from a convenience sample to the population.
Ignoring non-compliance, or analyzing only the participants who complied, which breaks the randomization.
Assuming that randomization justifies any question, including unethical or implausible manipulations.

R Corner: code and results

library(jtools)
set.seed(12345)
cond = sample(c(0, 1), 1000, replace = TRUE)
gender = sample(c(0, 1), 1000, replace = TRUE)          # 1 = male
y = 2 * cond - 1 * gender * cond + rnorm(1000)
pp = data.frame(y, cond, gender)

# A random experiment on a representative sample
table(pp$cond) / 1000
table(pp$gender) / 1000
summ(lm(y ~ cond, data = pp), model.info = FALSE)       # close to 1.5

# A random experiment on a biased sample: 400 men and 100 women
male = pp[pp$gender == 1, ]
female = pp[pp$gender == 0, ]
ps = rbind(male[sample(nrow(male), 400), ], female[sample(nrow(female), 100), ])
table(ps$cond) / 500
table(ps$gender) / 500
summ(lm(y ~ cond, data = ps), model.info = FALSE)       # close to 2 - 0.8 = 1.2

# Non-compliance: half of the assigned users do not use the feature
library(ivreg)
z = sample(c(0, 1), 1000, replace = TRUE)               # random assignment
comply = rbinom(1000, 1, 0.5)
d = z * comply                                          # actual use
y2 = 2 * d - 1 * gender * d + rnorm(1000)
coef(lm(y2 ~ z))["z"]                                   # ITT: about 0.75
coef(ivreg(y2 ~ d | z))["d"]                            # CACE: about 1.5

library(jtools)
set.seed(12345)
cond = sample(c(0, 1), 1000, replace = TRUE)
gender = sample(c(0, 1), 1000, replace = TRUE)          # 1 = male
y = 2 * cond - 1 * gender * cond + rnorm(1000)
pp = data.frame(y, cond, gender)

# A random experiment on a representative sample
table(pp$cond) / 1000

#> 
#>     0     1 
#> 0.516 0.484

table(pp$gender) / 1000

#> 
#>     0     1 
#> 0.508 0.492

summ(lm(y ~ cond, data = pp), model.info = FALSE)       # close to 1.5

#> MODEL FIT:
#> F(1,998) = 482.32, p = 0.00
#> R² = 0.33
#> Adj. R² = 0.33 
#> 
#> Standard errors:OLS
#> -----------------------------------------------
#>                     Est.   S.E.   t val.      p
#> ----------------- ------ ------ -------- ------
#> (Intercept)         0.00   0.05     0.10   0.92
#> cond                1.48   0.07    21.96   0.00
#> -----------------------------------------------

# A random experiment on a biased sample: 400 men and 100 women
male = pp[pp$gender == 1, ]
female = pp[pp$gender == 0, ]
ps = rbind(male[sample(nrow(male), 400), ], female[sample(nrow(female), 100), ])
table(ps$cond) / 500

#> 
#>     0     1 
#> 0.552 0.448

table(ps$gender) / 500

#> 
#>   0   1 
#> 0.2 0.8

summ(lm(y ~ cond, data = ps), model.info = FALSE)       # close to 2 - 0.8 = 1.2

#> MODEL FIT:
#> F(1,498) = 151.47, p = 0.00
#> R² = 0.23
#> Adj. R² = 0.23 
#> 
#> Standard errors:OLS
#> -----------------------------------------------
#>                     Est.   S.E.   t val.      p
#> ----------------- ------ ------ -------- ------
#> (Intercept)         0.02   0.06     0.34   0.74
#> cond                1.17   0.10    12.31   0.00
#> -----------------------------------------------

# Non-compliance: half of the assigned users do not use the feature
library(ivreg)
z = sample(c(0, 1), 1000, replace = TRUE)               # random assignment
comply = rbinom(1000, 1, 0.5)
d = z * comply                                          # actual use
y2 = 2 * d - 1 * gender * d + rnorm(1000)
coef(lm(y2 ~ z))["z"]                                   # ITT: about 0.75

#>         z 
#> 0.7360385

coef(ivreg(y2 ~ d | z))["d"]                            # CACE: about 1.5

#>        d 
#> 1.492442

Takeaways

Randomization removes confounding, but the estimate applies to the sample that was randomized.
With heterogeneous effects, a biased sample gives a biased estimate of the population effect.
Non-compliance dilutes the ITT; the CACE rescales it by the uptake.

Part V · Causal inference · Chapter 20

Controlling for Too Much

The hook

Why attractive actors seem less talented

Among famous actors, the most attractive ones often seem to be the least talented. Is there a trade-off between looks and talent? Probably not. To become famous, an actor needs talent, looks, or both. Among the famous, an actor with modest looks must therefore be very talented, and vice versa. Selecting on fame creates a negative correlation that does not exist in the population. This is Berkson's paradox, and it is a warning about control variables.

The weird question

Is it always better to control for more variables than fewer? Should we remove non-significant variables from the regression?

Colliders

A collider is a variable that is caused by two other variables, as in $\text{talent}\rightarrow\text{fame}\leftarrow\text{looks}$. A path through a collider is blocked by default. Controlling for the collider, or selecting the sample on it, opens the path and creates a spurious association. In other words, conditioning on a common effect is the mirror image of failing to condition on a common cause.

Play · Hollywood selection

Talent and looks are independent in 1,000 simulated actors. Fame = talent + looks + noise. Only the top share of actors by fame (gold) becomes famous.

Try this: select the top 10%. Among the famous, the correlation is strongly negative. Set the share to 100%, and the correlation disappears.
The last card shows the same bias without any selection: in the full sample, a regression of talent on looks that controls for fame also finds a negative coefficient.

Irrelevant and small-but-relevant variables

What about a variable that is truly unrelated to $y$? Including it does not bias the other coefficients, but it costs a degree of freedom and, if it is correlated with the other predictors, it inflates their standard errors. The inflation factor for the SE is $\sqrt{VIF}$, where the variance inflation factor $VIF_j=1/(1-R_j^2)$ and $R_j^2$ comes from regressing predictor $j$ on the other predictors. If the irrelevant variable is uncorrelated with the other predictors, the SEs hardly change.

Conversely, a variable may be relevant but have a small effect, so that it is non-significant. Removing it biases the other coefficients by its effect times its association with them. The bias is small when the effect is small, but it is not zero. Therefore, a non-significant coefficient alone is not a reason to remove a variable; the decision should be based on the causal diagram.

Play · The price of an irrelevant control

$y=2.5x_1+3.5x_2+\varepsilon$. The variable $x$ has no effect on $y$, but it is correlated with $x_1$. The sample size is 1,000.

Note on the original notes

The original notes concluded that "it is always safe to include all variables." This is not true. Including a collider (this chapter) or a mediator (Chapter 21) biases the treatment effect, and including an irrelevant variable that is correlated with the other predictors inflates their standard errors. The safe rule is to choose controls with a causal diagram, not to include or drop variables based on their significance.

Common traps

Controlling for a consequence of both the treatment and the outcome.
Studying a selected sample (the famous, the admitted, the survivors) without asking how selection happened.
Dropping a variable only because its coefficient is non-significant.

R Corner: code and results

library(jtools)
set.seed(12345)

# A collider: x1 -> x <- y
x1 = rnorm(1000); x2 = rnorm(1000)
y = 3.5 * x2 + rnorm(1000)
x = 3 * x1 + 2 * y + rnorm(1000)                   # x is a collider
m0 = lm(y ~ x1 + x2)                               # correct: x1 has no effect
m1 = lm(y ~ x1 + x2 + x)                           # wrong: x1 now "affects" y
export_summs(m0, m1)

# Berkson's paradox: selection on fame
talent = rnorm(1000); looks = rnorm(1000)
fame = talent + looks + rnorm(1000, 0, 0.5)
famous = fame > quantile(fame, 0.9)
cor(talent, looks)                                 # close to 0
cor(talent[famous], looks[famous])                 # strongly negative
coef(lm(talent ~ looks + fame))["looks"]           # controlling for the collider

# An irrelevant variable that is correlated with the predictors
x1 = rnorm(1000); x2 = rnorm(1000)
x = 2.5 * x1 + 3.5 * x2 + rnorm(1000)
y = 2.5 * x1 + 3.5 * x2 + rnorm(1000)
m0 = lm(y ~ x1 + x2 + x)
m1 = lm(y ~ x1 + x2)
export_summs(m0, m1)                               # same estimates, larger SEs in m0

# A relevant variable with a small effect
x1 = rnorm(1000, 3, 1); x2 = rnorm(1000, 5, 2)
x = 3 * x1 + 2 * x2 + 5 * rnorm(1000)
y = 2 * x1 + 3 * x2 + 0.02 * x + rnorm(1000, 6, 2.5)
m0 = lm(y ~ x1 + x2 + x)                           # x may be non-significant
m1 = lm(y ~ x1 + x2)                               # small bias: 0.02 * 3 and 0.02 * 2
export_summs(m0, m1)

library(jtools)
set.seed(12345)

# A collider: x1 -> x <- y
x1 = rnorm(1000); x2 = rnorm(1000)
y = 3.5 * x2 + rnorm(1000)
x = 3 * x1 + 2 * y + rnorm(1000)                   # x is a collider
m0 = lm(y ~ x1 + x2)                               # correct: x1 has no effect
m1 = lm(y ~ x1 + x2 + x)                           # wrong: x1 now "affects" y
export_summs(m0, m1)

#>                   ───────────────────────────────────────────
#>                                     Model 1       Model 2    
#>                                 ─────────────────────────────
#>                     (Intercept)     -0.03         -0.01      
#>                                     (0.03)        (0.01)     
#>                     x1              -0.01         -1.21 ***  
#>                                     (0.03)        (0.02)     
#>                     x2               3.49 ***      0.67 ***  
#>                                     (0.03)        (0.05)     
#>                     x                              0.40 ***  
#>                                                   (0.01)     
#>                                 ─────────────────────────────
#>                     N             1000          1000         
#>                     R2               0.93          0.99      
#>                   ───────────────────────────────────────────
#> *** p < 0.001; ** p < 0.01; * p < 0.05.
#> 
#> Column names: names, Model 1, Model 2

# Berkson's paradox: selection on fame
talent = rnorm(1000); looks = rnorm(1000)
fame = talent + looks + rnorm(1000, 0, 0.5)
famous = fame > quantile(fame, 0.9)
cor(talent, looks)                                 # close to 0

#> [1] -0.01049679

cor(talent[famous], looks[famous])                 # strongly negative

#> [1] -0.6552115

coef(lm(talent ~ looks + fame))["looks"]           # controlling for the collider

#>      looks 
#> -0.7956486

# An irrelevant variable that is correlated with the predictors
x1 = rnorm(1000); x2 = rnorm(1000)
x = 2.5 * x1 + 3.5 * x2 + rnorm(1000)
y = 2.5 * x1 + 3.5 * x2 + rnorm(1000)
m0 = lm(y ~ x1 + x2 + x)
m1 = lm(y ~ x1 + x2)
export_summs(m0, m1)                               # same estimates, larger SEs in m0

#>                   ───────────────────────────────────────────
#>                                     Model 1       Model 2    
#>                                 ─────────────────────────────
#>                     (Intercept)      0.01          0.01      
#>                                     (0.03)        (0.03)     
#>                     x1               2.45 ***      2.46 ***  
#>                                     (0.08)        (0.03)     
#>                     x2               3.43 ***      3.44 ***  
#>                                     (0.11)        (0.03)     
#>                     x                0.00                    
#>                                     (0.03)                   
#>                                 ─────────────────────────────
#>                     N             1000          1000         
#>                     R2               0.95          0.95      
#>                   ───────────────────────────────────────────
#> *** p < 0.001; ** p < 0.01; * p < 0.05.
#> 
#> Column names: names, Model 1, Model 2

# A relevant variable with a small effect
x1 = rnorm(1000, 3, 1); x2 = rnorm(1000, 5, 2)
x = 3 * x1 + 2 * x2 + 5 * rnorm(1000)
y = 2 * x1 + 3 * x2 + 0.02 * x + rnorm(1000, 6, 2.5)
m0 = lm(y ~ x1 + x2 + x)                           # x may be non-significant
m1 = lm(y ~ x1 + x2)                               # small bias: 0.02 * 3 and 0.02 * 2
export_summs(m0, m1)

#>                   ───────────────────────────────────────────
#>                                     Model 1       Model 2    
#>                                 ─────────────────────────────
#>                     (Intercept)      5.60 ***      5.61 ***  
#>                                     (0.33)        (0.33)     
#>                     x1               2.23 ***      2.17 ***  
#>                                     (0.10)        (0.08)     
#>                     x2               3.10 ***      3.06 ***  
#>                                     (0.05)        (0.04)     
#>                     x               -0.02                    
#>                                     (0.02)                   
#>                                 ─────────────────────────────
#>                     N             1000          1000         
#>                     R2               0.87          0.87      
#>                   ───────────────────────────────────────────
#> *** p < 0.001; ** p < 0.01; * p < 0.05.
#> 
#> Column names: names, Model 1, Model 2

Takeaways

Never control for a collider, and beware of samples selected on one.
Irrelevant controls do not bias estimates, but they can inflate SEs.
Choose controls with a causal diagram, not with p values.

Part V · Causal inference · Chapter 21

Mediators

The hook

Coffee, alertness, and productivity

Coffee makes office workers more productive, partly because it makes them more alert. A consultant regresses productivity on coffee and alertness and reports that the effect of coffee is small. The manager concludes that the coffee budget can be cut. What went wrong?

The weird question

Should we control for mediators to estimate the treatment effect?

Direct, indirect, and total effects

A mediator $M$ lies on a causal path from the treatment $X$ to the outcome $Y$: $X\rightarrow M\rightarrow Y$. With linear relationships, $M=aX+e_M$ and $Y=c'X+bM+e_Y$. Three effects follow:

the direct effect $c'$, which is the effect of $X$ holding $M$ constant;
the indirect effect $ab$, which is the effect transmitted through $M$;
the total effect $c=c'+ab$, which is the effect of changing $X$ and letting $M$ respond.

Controlling for a mediator therefore estimates the direct effect only. Because the mediator is a descendant of the treatment, it violates the backdoor criterion (Chapter 18). If the question is "what happens if workers drink more coffee," the answer is the total effect, and the mediator should not be controlled.

Play · The coffee path diagram

Coffee (cups) $\sim\mathsf{Norm}(5,1)$; alertness $=a\cdot\text{coffee}+3\varepsilon_M$; productivity $=c'\cdot\text{coffee}+b\cdot\text{alertness}+\sigma\varepsilon_Y$; 1,000 workers. The confidence interval of $ab$ comes from 500 bootstrap resamples.

Try this: set $c'=0$. The coffee effect disappears once alertness is controlled, yet the total effect is large: full mediation.
Set $a=0$. Alertness no longer depends on coffee, so controlling for it does not change the coffee coefficient.

Testing mediation

The traditional approach compares the coefficient of $X$ in models with and without the mediator. The difference equals $ab$ in linear models, but the approach provides neither a standard error nor a test. A better approach estimates both equations simultaneously, e.g., as a structural equation model (SEM) in lavaan, and tests $ab$ with a bootstrap confidence interval, because the sampling distribution of a product is usually skewed. Unless a particular mechanism is the research question, I estimate the treatment effect without controlling for mediators.

Think about it

Does a significant indirect effect prove that alertness is the mechanism?

My answer

No. The decomposition assumes that nothing confounds the relationship between the mediator and the outcome. Even if coffee is randomized, alertness is not. A third variable, e.g., a good night's sleep, could raise both alertness and productivity. Then $b$ is biased, and so is $ab$. Mediation analysis is thus a causal analysis of two relationships, and the second one usually rests on observational data.

Common traps

Controlling for a mediator and reporting the result as the treatment effect.
Testing mediation only by comparing coefficients across models.
Interpreting an indirect effect causally without considering mediator-outcome confounding.

R Corner: code and results

library(jtools)
set.seed(12345)
x = rnorm(1000, 5, 1)
m = 2 * x + 3 * rnorm(1000)
y = 3 * x + 5 * m + 7 * rnorm(1000, 6, 3)

# The traditional way: compare models with and without the mediator
m0 = lm(y ~ x)
m1 = lm(y ~ x + m)
export_summs(m0, m1)
coef(m0)["x"] - coef(m1)["x"]          # equals a * b in linear models

# The SEM way, with a bootstrap test of the indirect effect
library(lavaan)
data = data.frame(y, x, m)
mod = "
  m ~ a * x            # a path
  y ~ b * m            # b path
  y ~ c * x            # direct path
  ab := a * b          # indirect effect
  total := c + ab      # total effect
"
fsem = sem(mod, data = data, se = "bootstrap", bootstrap = 1000)
summary(fsem)
parameterEstimates(fsem, boot.ci.type = "perc")[7:8, ]

library(jtools)
set.seed(12345)
x = rnorm(1000, 5, 1)
m = 2 * x + 3 * rnorm(1000)
y = 3 * x + 5 * m + 7 * rnorm(1000, 6, 3)

# The traditional way: compare models with and without the mediator
m0 = lm(y ~ x)
m1 = lm(y ~ x + m)
export_summs(m0, m1)

#>                   ───────────────────────────────────────────
#>                                     Model 1       Model 2    
#>                                 ─────────────────────────────
#>                     (Intercept)     39.28 ***     42.67 ***  
#>                                     (4.09)        (3.28)     
#>                     x               13.34 ***      2.84 ***  
#>                                     (0.79)        (0.78)     
#>                     m                              4.96 ***  
#>                                                   (0.21)     
#>                                 ─────────────────────────────
#>                     N             1000          1000         
#>                     R2               0.22          0.50      
#>                   ───────────────────────────────────────────
#> *** p < 0.001; ** p < 0.01; * p < 0.05.
#> 
#> Column names: names, Model 1, Model 2

coef(m0)["x"] - coef(m1)["x"]          # equals a * b in linear models

#>        x 
#> 10.50494

# The SEM way, with a bootstrap test of the indirect effect
library(lavaan)
data = data.frame(y, x, m)
mod = "
  m ~ a * x            # a path
  y ~ b * m            # b path
  y ~ c * x            # direct path
  ab := a * b          # indirect effect
  total := c + ab      # total effect
"
fsem = sem(mod, data = data, se = "bootstrap", bootstrap = 1000)
summary(fsem)

#> lavaan 0.7-3 ended normally after 1 iteration
#> 
#>   Estimator                                         ML
#>   Optimization method                           NLMINB
#>   Number of model parameters                         5
#> 
#>   Number of observations                          1000
#> 
#> Model Test User Model:
#>                                                       
#>   Test statistic                                 0.000
#>   Degrees of freedom                                 0
#>                                                       
#>   Browne's residual (NT model-based) test             
#>   Test statistic                                 0.000
#>   Degrees of freedom                                 0
#> 
#> Parameter Estimates:
#> 
#>   Standard errors                            Bootstrap
#>   Number of requested bootstrap draws             1000
#>   Number of successful bootstrap draws             984
#> 
#> Regressions:
#>                    Estimate  Std.Err  z-value  P(>|z|)
#>   m ~                                                 
#>     x          (a)    2.117    0.099   21.287    0.000
#>   y ~                                                 
#>     m          (b)    4.961    0.217   22.902    0.000
#>     x          (c)    2.838    0.785    3.614    0.000
#> 
#> Variances:
#>                    Estimate  Std.Err  z-value  P(>|z|)
#>    .m                 9.128    0.397   22.976    0.000
#>    .y               403.210   17.364   23.222    0.000
#> 
#> Defined Parameters:
#>                    Estimate  Std.Err  z-value  P(>|z|)
#>     ab               10.505    0.640   16.408    0.000
#>     total            13.343    0.772   17.273    0.000

parameterEstimates(fsem, boot.ci.type = "perc")[7:8, ]

#>     lhs op  rhs label    est    se      z pvalue ci.lower ci.upper
#> 7    ab :=  a*b    ab 10.505 0.640 16.408      0    9.281   11.715
#> 8 total := c+ab total 13.343 0.772 17.273      0   11.867   14.901

Takeaways

Total effect = direct effect + indirect effect.
Controlling for a mediator removes part of the treatment effect.
Test indirect effects with SEM and bootstrap confidence intervals.

Situation	Test	Reference distribution under \(H_0\)
One mean, \(\sigma\) known	Z-test	\(\mathsf{Norm}(0,1)\)
One mean or two means, \(\sigma\) estimated	t-test	\(t\) with \(n-1\) (or Welch) df
Two dependent measurements (before and after)	Paired t-test	\(t\) with \(n-1\) df, on the differences
More than two means	F-test (ANOVA)	\(F\) with \(k-1\) and \(N-k\) df
Two categorical variables	\(\chi^2\)-test	\(\chi^2\) with \((r-1)(c-1)\) df

Weird Statistics Questions

Preface

How each chapter works

Contents

How to cite this book

Declaration of AI assistance

The Mean, the Standard Deviation, and Their Limits

The bubble-tea queue

A distribution is the data, summarized

When the formula knows best: parameters

Categories have no mean

When the mean does not exist

Expected value: why the mean still matters

A guarantee for any shape: Chebyshev's inequality

Many Correlations, and Association Between Categories

Scrolling before bed

Pearson's \(r\): linear association

Spearman's \(\rho\): monotonic association

Kendall's \(\tau\): a probability

Pearson versus Spearman when the curve bends

Association between two categorical variables

Why a Sample Can Speak for a Population

One spoon of soup

Samples vary, and that is fine

The sampling distribution

The central limit theorem

Sample size \(n\) versus the number of samples \(N\)

Not every sampling distribution is Normal

Sampling from a finite population

Why large samples work, and when they do not

Standard Errors, Sampling Errors, and Confidence Intervals

Throwing 100 fishing nets

Four terms, one picture

t, Z, F, and \(\chi^2\) Tests: One Idea, Different Distributions

Is the 30 cm pizza really 30 cm?

One recipe for every test statistic

Z versus t: the price of estimating \(\sigma\)

What the t-test assumes

Two samples, the variance sum law, and Welch's test

Beyond two means

Distribution-Free Tests: Permutation, Bootstrap, and Exact Tests

The lady tasting tea

Counting instead of assuming

What "distribution-free" means

The randomization test

The bootstrap

Randomization, permutation, and bootstrap

Equivalence Tests: Can the Null Hypothesis Be Accepted?

Is the new milk-tea recipe the same?

Failing to reject is not accepting

Turning the hypotheses around

Can Regression Replace the Tests?

One model, many tests

Why "models," plural

The t-test is a regression with one dummy variable

ANOVA is a regression with several dummy variables

Non-normal outcomes and categorical outcomes

Degrees of Freedom and Over-Fitting

Memorizing last year's exam

Degrees of freedom

Three ways to reach \(R^2=1\)

Interpreting and Comparing Regression Coefficients

Pricing a Hong Kong flat

Four things to read from a coefficient

Binary and categorical predictors

Comparing coefficients

Test a linear combination: \(b_A - k\,b_B = 0\)

Why the I.I.D. Assumption Matters

Marbles from one bag

A regression has a random part and a systematic part

What goes wrong without independence

Skewed Outcomes, Normal Residuals

Incomes and noise

The assumption is about the random part

Effect Sizes and Unique Contributions

Who deserves credit for the group project?

Effect sizes

Unique contributions: \(sr^2\)

Main Effects in Models with Interactions

Coffee and night owls