These notes are meant to serve as an introduction to some intermediate concepts in mathematics that are useful for quantitative scientists. They are not meant to replace undergraduate or graduate training in calculus, linear algebra, analysis, and other key math disciplines. However, the hope is that this document will serve as a useful review for those who have had prior exposure to these topics, and a useful introduction for those who have not.
Rather than focus on developing technical proficiency (e.g., the ability to compute complex derivatives using the chain rule), we focus here on ideas, conceptual framing, and interpreting notation. The goal is to give you some skills you can use so that, when you encounter an equation, you can slow down, take it apart, and understand what it is saying.
One goal of using math in a quantitative field like epidemiology is that it can serve a dual function: First, mathematical notation is a compression device for scientific reasoning. Consider how a scientific question gets translated into mathematics:
Each step strips away context and gains precision, making the vague question “does treatment reduce mortality?” much more specific, rigorously defined, quantifiable, and (ideally) testable with data.
However, simultaneously, this act of translation from broad substantive meaning to precise and specific mathematical form, strips away important context that matters for what we do. Therefore, another goal of using math in a field like epidemiology is to understand exactly what goes missing when we abstract elements of our work into a more precise and specific form (Naimi 2026).
Finally, understanding math is a skill required by much of the modern literature in machine learning, data science, and causal inference (Naimi et al. 2026). These topic areas have advanced tremendously in the past decades, making them important for what we do as scientists. Developing a tool-kit to be able to understand this literature is therefore very useful.
In these notes, we cover four domains: algebra, focusing on rules for manipulating symbols, and notational conventions (functions, sums, products, indicators); calculus, which is the mathematics of change (derivatives) and accumulation (integrals), central to key concepts and tools in probability and statistics (rates, hazards, expected values, and optimization); linear algebra, which is central to how regression actually happens on your computer; finally, we practice mathematical literacy: reading real equations from the statistics, causal inference, and machine learning literatures.
A note on how to use these notes: The examples here were constructed to be deliberately easy to facilitate developing a general conceptual understanding. Work through the small numerical examples with pencil and paper or in R. The R code (folded by default; click “Show R code”) is there to show you how each idea looks in software.
Algebra is the grammar of math and methods development. Deriving an estimator or evaluating the properties of an estimator, or identifying an estimand (and som much more) requires algebra. This section reviews the grammar, rules, and notation.
A variable is a symbol that stands for a quantity that can take different values: age, blood pressure, exposure status. A constant is a symbol whose value is fixed, either a specific number (like \(\pi\)) or a quantity treated as fixed in a given problem.
Three notational conventions1 appear frequently in statistics and epidemiology:
We will also occasionally use basic set notation: curly braces list members of a set, as in \(A \in \{0, 1\}\), which reads “\(A\) takes values in the set containing 0 and 1” (i.e., \(A\) is binary).
Throughout these notes we will return to a small cohort of five patients:
cohort <- data.frame(
id = 1:5,
age = c(45, 62, 37, 54, 41),
smoker = c(0, 1, 0, 1, 1),
sbp = c(128, 142, 119, 131, 125)
)
kable(cohort, col.names = c("ID", "Age (years)", "Smoker", "SBP (mmHg)"))| ID | Age (years) | Smoker | SBP (mmHg) |
|---|---|---|---|
| 1 | 45 | 0 | 128 |
| 2 | 62 | 1 | 142 |
| 3 | 37 | 0 | 119 |
| 4 | 54 | 1 | 131 |
| 5 | 41 | 1 | 125 |
We use five rows (observations) so we can check calculations by hand. Whenever notation feels abstract in what follows, consider coming back to this table and asking: what does this symbol mean for these five people?
Nearly every equation in the methods literature is built from these conventions. When you see \(\sum_{i=1}^n (Y_i - \hat{\beta}_0 - \hat{\beta}_1 X_i)^2\), you should be able to parse it instantly as: “for each person \(i\), take their outcome, subtract their model-predicted value, square it, and add up across all \(n\) people.”
An exponent is repeated multiplication: \(2^3 = 2 \times 2 \times 2 = 8\), read “two to the third power” or “two cubed.” Three rules cover almost every manipulation you will see:
\[ a^m \times a^n = a^{m+n}, \qquad (a^m)^n = a^{mn}, \qquad a^{-n} = \frac{1}{a^n}. \]
The first rule says multiplying powers adds exponents; the second says raising a power to a power multiplies them; the third defines negative exponents as reciprocals. Fractional exponents define roots: \(a^{1/2} = \sqrt{a}\), the number that, multiplied by itself, gives \(a\). So \(\sqrt{9} = 3\) because \(3^2 = 9\).
These show up a lot in epidemiology:
Exponential growth. Early in an outbreak, cases may double every fixed interval. If cases double every 3 days, then after \(t\) days there are \(c_0 \times 2^{t/3}\) cases, where \(c_0\) is the starting count. After 30 days that multiplier is \(2^{10} = 1{,}024\)—a thousand-fold increase in a month. Exponents are how small, constant multiplicative changes compound into enormous totals, which is why epidemic curves, cumulative survival, and compound risks all involve them.
Square roots and precision. The standard error of a sample mean is \(\sigma / \sqrt{n}\): the population standard deviation divided by the square root of the sample size. The square root is why precision is expensive: to cut your standard error in half, you must quadruple your sample size, because \(\sqrt{4n} = 2\sqrt{n}\).
Flexible regression modeling. Polynomial and fractional polynomial terms are often used in regression modeling to avoid assuming perfectly linear relationships between an independent variable (on the right hand side of the regression equation) and the dependent variable (Greenland 1995). For example, \(E(Y \mid X) = \beta_0 + \beta_1 X\) assumes that as \(X\) increases, the conditional mean of \(Y\) changes linearly. In contrast, \(E(Y \mid X) = \beta_0 + \beta_1 X + \beta_2 X^2\) or \(E(Y \mid X) = \beta_0 + \beta_1 X^{1/2} + \beta_2 X + \beta_3 X^{3/2} + \beta_4 X^2\) both allow this relationship to take different shapes.
One number deserves special mention: \(e \approx 2.718\), Euler’s number. It arises naturally when growth compounds continuously rather than in discrete steps, and it is the base used almost universally in statistical models. You will constantly see \(e^x\), also written \(\exp(x)\) (“the exponential of \(x\)”). For now, treat \(\exp(x)\) as a specific, well-behaved function your software knows; its special properties will become clear when we reach calculus.
Exponents answer “what do I get after \(t\) periods of doubling?”. Logarithms answer the reverse: “how many doublings does it take to reach this value?” Formally, \(\log_b(y)\) is the exponent to which the base \(b\) must be raised to produce \(y\):
\[ \log_2(8) = 3 \quad \text{because} \quad 2^3 = 8; \qquad \log_{10}(1000) = 3 \quad \text{because} \quad 10^3 = 1000. \]
The natural logarithm, written \(\ln(y)\) or, in nearly all statistical writing, simply \(\log(y)\), uses base \(e\). In these notes and in most of statistics, \(\log\) will mean the natural log: \(\log_e\).
Strogatz (Strogatz 2019) describes logarithms as one of the greatest labor-saving devices in the history of calculation, and the reason is one property:
\[ \log(a \times b) = \log(a) + \log(b), \qquad \log(a^k) = k \log(a). \]
Logarithms convert multiplication into addition. Before computers, this is what made hard arithmetic feasible (it is the principle behind the slide rule). In modern statistics, the same property is why so much modeling happens “on the log scale”:
\[ \operatorname{logit}(p) = \log\!\left(\frac{p}{1-p}\right), \]
live on the entire real line \((-\infty, \infty)\). Logistic regression models the logit precisely so that no combination of covariates can ever imply a probability below 0 or above 1. Figure 1 shows the logit transformation: the middle of the probability scale is stretched only gently, while the extremes near 0 and 1 are stretched enormously, so they never lie outside of the (0,1) bounds.
p <- seq(0.005, 0.995, by = 0.001)
d <- data.frame(p = p, logodds = log(p / (1 - p)))
ggplot(d, aes(p, logodds)) +
geom_hline(yintercept = 0, colour = ref_grey, linewidth = 0.3) +
geom_line(colour = pal["blue"], linewidth = 0.9) +
annotate("point", x = 0.5, y = 0, size = 2, colour = pal["orange"]) +
annotate("text", x = 0.56, y = -0.55, label = "p = 0.5 → log-odds = 0",
hjust = 0, size = 3.4, colour = "grey20") +
labs(x = "Probability p", y = "Log-odds: log(p / (1 - p))")
Consider: a risk of \(p = 0.2\) gives odds of \(0.2/0.8 = 0.25\) and a log-odds of \(\log(0.25) \approx -1.39\). When a logistic regression reports a coefficient of, say, \(0.7\) for smoking, it is saying the log-odds of the outcome differ by \(0.7\) between (e.g.) smokers and non-smokers. Equivalently, the odds ratio is \(e^{0.7} \approx 2.0\), because exponentiating undoes the log.
“Log” in a statistics paper almost never means base 10. When you see \(\log\) in a model, a likelihood, or R output, read it as the natural log, base \(e\). (In R, log() is the natural log; base 10 is log10().) Mixing up bases will not change which models fit best, but it will change how you back-transform coefficients.
Solving an equation means finding the value of an unknown that makes the equation true, and the method is always the same idea: undo the operations, in reverse order, doing the same thing to both sides. To solve \(3x + 6 = 21\), subtract 6 from both sides (\(3x = 15\)), then divide both sides by 3 (\(x = 5\)).
This matters in epidemiology because many quantities are defined implicitly by a model, and you have to invert the model to solve for them (using data). Alternatively, we sometimes want to extract information out of a model that is not immediately present, and requires some transformation and solving. For example: logistic regression gives you a linear model for the log-odds,
\[ \log\!\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1 a. \]
But you might want the risk \(p\) for a person with exposure \(a\). Solve for \(p\): exponentiate both sides to undo the log, giving \(p/(1-p) = e^{\beta_0 + \beta_1 a}\); then a few algebraic moves (multiply both sides by \(1-p\), collect the \(p\) terms, divide) yield
\[ p = \frac{e^{\beta_0 + \beta_1 a}}{1 + e^{\beta_0 + \beta_1 a}} = \operatorname{expit}(\beta_0 + \beta_1 a). \]
That final function is called the expit (or inverse-logit), and it is basically the logit equation solved for \(p\). Every time software converts your logistic model into predicted risks, it is doing this algebra.
A second example that illustrates a key methods tool in statistics. We can define the mean by expressing it as a solution to an equation: find the value \(\mu\) that makes
\[ \sum_{i=1}^{n} (x_i - \mu) = 0. \]
In words: the value around which the data’s deviations exactly cancel. A little algebra (distribute the sum, move terms across the equals sign, divide by \(n\)) shows the solution is \(\mu = \frac{1}{n}\sum_i x_i\), the sample mean. In effect, this means that we can estimate the mean by finding a value of \(mu\) that solves the above equation. This “estimator as the solution to an equation set to zero” is the foundation of estimating equations, a general framework in which nearly every estimator you will see (means, regression coefficients, standardized risks) can be characterized by. In fact, M estimators are useful in that they can be used to solve for things like the g formula, and simultaneously provide standard errors as a bonus (Zivich et al. 2026; Ross et al. 2024; Stefanski and Boos 2002), obviating the need to use tools like the bootstrap.
A function is a rule that turns each input into exactly one output. We write \(f(x)\), read “\(f\) of \(x\),” for the output of rule \(f\) applied to input \(x\). It helps to picture a machine:
flowchart LR A["Input<br/>age = 50"] --> B["Rule f<br/>risk as a function of age"] --> C["Output<br/>f(50) = 0.12"]
Two words describe a function’s scope. The domain is the set of allowed inputs; the range is the set of possible outputs:
Matching domains and ranges is what link functions do in the context of generalized linear models, and the mechanism is composition: applying one function to the output of another, written \(f(g(x))\), read “\(f\) of \(g\) of \(x\).” In a logistic model, the risk for a person with exposure \(a\) is
\[ p(a) = \operatorname{expit}(\underbrace{\beta_0 + \beta_1 a}_{\text{inner: linear}}), \]
a composition of a linear function (which can output any real number) with the expit (which maps any real number into \((0,1)\)). Generalized linear models are, at heart, a catalog of useful compositions. Composition will return when we meet the chain rule (Section 2.2.1), which tells us how change propagates through composed functions.
Functions can also take several inputs at once, as in \(f(a, w)\): risk as a function of both treatment and a covariate.
The simplest functions are built from powers of \(x\):
x <- seq(-2, 2, length.out = 200)
d <- rbind(
data.frame(x = x, y = 1 + 0.8 * x, f = "Linear: 1 + 0.8x"),
data.frame(x = x, y = 1 + 0.8 * x - 0.9 * x^2, f = "Quadratic: adds x\u00b2"),
data.frame(x = x, y = 0.5 * x^3 - x, f = "Cubic: adds x\u00b3")
)
d$f <- factor(d$f, levels = unique(d$f))
ggplot(d, aes(x, y)) +
geom_line(colour = pal["blue"], linewidth = 0.9) +
facet_wrap(~ f, scales = "free_y") +
labs(x = "x", y = "f(x)")
Polynomials matter for at least two reasons: First, they allow model flexibility: splines, fractional polynomials, and many machine learning basis expansions are polynomials. Second (covered more in Section 2.5) polynomials are “easy” functions that mathematicians use to approximate hard ones. Much of statistical theory works because complicated functions behave, locally, like quadratic polynomial functions.
The capital Greek sigma, \(\sum\), is an instruction to add:
\[ \sum_{i=1}^{n} x_i = x_1 + x_2 + \cdots + x_n. \]
\(\sum_{i=1}^{n} x_i\) reads: “the sum, over observations \(i\) from 1 to \(n\), of \(x_i\)”, that is, add up the values of \(x_i\) for every observation from the first to the last. The expression below the sigma names the index and where it starts; the expression above says where it stops; the expression to the right is what gets added.
Using our five-patient cohort, the sample mean of systolic blood pressure is
\[ \bar{x} = \frac{1}{n}\sum_{i=1}^{n} x_i = \frac{128 + 142 + 119 + 131 + 125}{5} = \frac{645}{5} = 129 \text{ mmHg}, \]
where \(\bar{x}\) (“x-bar”) is standard notation for a sample mean. A slight extension covers weighted means, which appear whenever observations should not count equally (survey weights, inverse-probability weights, standardization weights):
\[ \bar{x}_w = \frac{\sum_{i=1}^{n} w_i x_i}{\sum_{i=1}^{n} w_i}. \]
Each person’s value is multiplied by their weight \(w_i\); the denominator rescales so the weights behave like counts. If patient 2 in our cohort represented 200 people in the source population while the others represented 100 each, their blood pressure would count twice as heavily in \(\bar{x}_w\).
The capital pi, \(\prod\), is the multiplicative sibling: an instruction to multiply.
\[ \prod_{i=1}^{n} p_i = p_1 \times p_2 \times \cdots \times p_n. \]
Products arise naturally with independent probabilities, because independent probabilities multiply. If three independent patients each have probability \(0.9\) of remaining event-free, the probability all three remain event-free is \(\prod_{i=1}^3 0.9 = 0.9^3 = 0.729\). This is exactly the structure of a likelihood: the probability of observing your whole dataset is the product, over people, of each person’s probability contribution. Note how this here is related ot the section on logs, which turn products into sums (Section 1.3),
\[ \log\!\left(\prod_{i=1}^{n} p_i\right) = \sum_{i=1}^{n} \log(p_i), \]
statisticians maximize the log-likelihood rather than the likelihood because sums are easier to work with than products (and better behaved on a computer, where multiplying thousands of small probabilities can “underflow” to zero; i.e., can become so small that the computer is not able to distinguish it from true zero).
sbp <- cohort$sbp
c(sum = sum(sbp), mean = mean(sbp)) # sigma notation in R sum mean
645 129 p <- c(0.9, 0.9, 0.9)
c(product = prod(p), exp_sum_log = exp(sum(log(p)))) # products via logs product exp_sum_log
0.729 0.729 An indicator function converts a condition into a 1 if the condition holds, and 0 if it does not. Notation varies across papers, such as \(\mathbf{1}[\cdot]\), \(\mathbb{1}(\cdot)\), \(I(\cdot)\):
\[ \mathbf{1}[\text{smoker}_i] = \begin{cases} 1 & \text{if person } i \text{ smokes} \\ 0 & \text{otherwise.} \end{cases} \]
Indicators are the bridge between categories and arithmetic. Once “is a smoker” is coded as 0/1, we can add it, average it, and put it in a regression. In our cohort, the smoking indicator takes values \((0, 1, 0, 1, 1)\), and its mean is
\[ \frac{1}{n}\sum_{i=1}^{n} \mathbf{1}[\text{smoker}_i] = \frac{0+1+0+1+1}{5} = 0.6, \]
which is just the proportion of smokers. This is a fact worth framing: the average of an indicator is a proportion, and (in its population version) the expected value of an indicator is a probability (Hardt and Recht 2022). This is one reason why it’s useful to use dummy varaible coding, since using other levels such as \((1,2)\) can change the meaning.
Prevalence, risk, and attack rates are all means of indicators. That one observation lets every tool built for means, such as regression, standardization, and standard errors, apply immediately to probabilities.
Indicators are everywhere once you look:
mean(cohort$smoker) # mean of an indicator = proportion[1] 0.6sum(cohort$smoker == 1) # 1[condition] via a logical test[1] 3mean(cohort$sbp[cohort$smoker == 1]) # mean SBP among smokers only[1] 132.66671. Write out \(\sum_{i=1}^{3} x_i^2\) for \(x = (2, 1, 3)\) and compute it.
\(x_1^2 + x_2^2 + x_3^2 = 4 + 1 + 9 = 14\). Note the order of operations: square first, then sum. \(\left(\sum x_i\right)^2 = 36\) is a different quantity (a distinction that matters constantly in variance formulas).
2. Your logistic model gives log-odds \(\beta_0 + \beta_1 a = -2.2\) for an unexposed person. What is their predicted risk, and which function did you use?
Apply the expit: \(p = e^{-2.2}/(1 + e^{-2.2}) \approx 0.111/1.111 \approx 0.10\). The expit is the logit solved for \(p\)—an exercise in undoing operations.
3. In words, what does \(\mathbf{1}[\text{age}_i \geq 65]\) equal for a 70-year-old? What does the mean of this indicator across a cohort estimate?
It equals 1 for the 70-year-old (the condition holds). Its mean across the cohort is the proportion aged 65 or older—an example of “average of an indicator = proportion.”
4. Why do statisticians work with \(\sum_i \log(p_i)\) rather than \(\prod_i p_i\) when fitting models by maximum likelihood?
The log converts the product into a sum without changing where the maximum is (the log is an increasing function, so whatever maximizes the likelihood maximizes the log-likelihood). Sums are easier to differentiate and numerically stable, whereas products of many small probabilities underflow to zero on a computer.
Looking ahead: every regression model, likelihood, and causal formula in this course is assembled from the parts in this section, including functions composed with links, sums over people, products over probability contributions, and indicators selecting subgroups.
Next we add the mathematics of change and accumulation.
Calculus is fundamentally about change and accumulation over time, which is central to what we do in epidemiology. A rate is a change per unit time; a hazard is an instantaneous rate; cumulative incidence accumulates risk over a follow-up period; person-time accumulates observation. Moreover, the statistical machinery we use every day (maximum likelihood, least squares, machine learning tools like gradient boosting) works by using derivatives to find optimal parameter values.
Strogatz (Strogatz 2019) summarizes the strategy of calculus as the “Infinity Principle”: to solve a hard problem about change or accumulation, chop it into infinitely many infinitesimally small pieces, solve the easy piece-level problem, and reassemble. Derivatives are the chopping (what happens in an instant?); integrals reassemble (what accumulates over the whole interval?).
A limit asks: as the input approaches some value, what does the output approach? We write
\[ \lim_{x \to a} f(x) = L, \]
read “the limit of \(f\) of \(x\), as \(x\) approaches \(a\), equals \(L\).” The function need never actually reach \(L\); what matters is where it is heading. The simplest example: \(1/n\) gets arbitrarily close to 0 as \(n\) grows, so \(\lim_{n \to \infty} 1/n = 0\), even though \(1/n\) never actually equals 0.
A function is continuous if it has no jumps or holes—informally, you can draw it without lifting your pen. Risk as a function of age is usually modeled as continuous; a step function (like “risk jumps at the 65th birthday”) is not.
Limits present themselves in at least three areas related to epidemiology:
A hazard is a rate, not a probability, and it is not bounded by 1. A hazard of 2 per person-year is perfectly meaningful (events are occurring fast); a probability of 2 is nonsense. Keeping “limit of a rate” (hazard) distinct from “accumulated probability” (risk) will help you in survival analysis.
The derivative of a function is its instantaneous rate of change, which can be thought of as the slope of a tangent line to a curve at a single point. Between two points we can always compute an average rate of change (the slope of the connecting line, called a secant): change in output over change in input, “rise over run.” The derivative is what happens to that average as the two points slide together:
\[ f'(x) = \lim_{h \to 0} \frac{f(x + h) - f(x)}{h}. \]
\(f'(x)\) reads “\(f\)-prime of \(x\)”: the derivative of \(f\) at input \(x\). The equivalent notation \(\frac{df}{dx}\) or \(\frac{dy}{dx}\) reads “the derivative of \(f\) (or \(y\)) with respect to \(x\)” and emphasizes what is changing. Inside the limit, \(h\) is the distance between the two points: \([f(x+h) - f(x)]/h\) is the slope of a secant over an interval of width \(h\), and the limit slides the interval’s width to zero.
Figure 4 shows the idea on an epidemiologic curve: risk of disease as a function of age. The slope of the secant from age 40 to age 60 is the average increase in risk per year over that span; the tangent line at age 50 shows the instantaneous increase per year at exactly age 50.
expit <- function(x) exp(x) / (1 + exp(x))
risk <- function(a) expit(-6 + 0.08 * a)
a <- seq(30, 75, length.out = 300)
d <- data.frame(a = a, r = risk(a))
p50 <- risk(50); slope50 <- 0.08 * p50 * (1 - p50)
sec_slope <- (risk(60) - risk(40)) / 20
ggplot(d, aes(a, r)) +
geom_line(colour = pal["blue"], linewidth = 0.9) +
annotate("segment", x = 38, xend = 62, y = risk(40) + (38 - 40) * sec_slope,
yend = risk(40) + (62 - 40) * sec_slope,
colour = ref_grey, linetype = "dashed", linewidth = 0.6) +
annotate("segment", x = 42, xend = 58, y = p50 + (42 - 50) * slope50,
yend = p50 + (58 - 50) * slope50,
colour = pal["orange"], linewidth = 0.8) +
annotate("point", x = c(40, 60), y = c(risk(40), risk(60)), colour = ref_grey, size = 2) +
annotate("point", x = 50, y = p50, colour = pal["orange"], size = 2) +
annotate("text", x = 61.5, y = risk(60) - 0.01,
label = "secant: average\nchange, 40 to 60",
hjust = 0, size = 3.2, colour = "grey35") +
annotate("text", x = 40, y = 0.30, label = "tangent: instantaneous\nchange at 50",
hjust = 0, size = 3.2, colour = pal["orange"]) +
coord_cartesian(xlim = c(30, 75)) +
labs(x = "Age (years)", y = "Risk of disease")
If the risk curve at age 50 has derivative \(0.008\) per year, then at that age, risk is climbing by about 0.8 percentage points per year of age. Another derivative commonly encountered is in survival analysis. Let \(H(t)\) denote the cumulative hazard, which is the total accumulated event intensity up to some arbitrary time \(t\). The hazard function is its derivative:
\[ h(t) = \frac{d}{dt} H(t). \]
The hazard is to the cumulative hazard what speed is to distance traveled: one is the instantaneous rate, the other is the running total. Because the cumulative hazard is also mathematically related to the survival, or cumulative distribution function, we can define the hazard as a specific type of derivative of the survival curve.
You do not really need to be fluent in differentiation rules and how to apply them in a given setting. In fact, there are software tools available to do these computations. But two facts are worth knowing because they appear in derivations constantly: the derivative of \(x^2\) is \(2x\) (so, e.g., at \(x = 3\) the parabola climbs at slope 6), and the exponential function is its own derivative: \(\frac{d}{dx} e^x = e^x\). That second fact is the reason \(e\) is the universal base in modeling: growth whose rate is proportional to its current size is exponential growth, exactly.
f <- function(x) x^2
h <- 1e-6
(f(3 + h) - f(3)) / h # numerical slope at x = 3; compare to 2x = 6[1] 6.000001Real models are built by composing functions (Section 1.5), so it’s sometimes important to know how change occurs through a composition: in other words, how to take a derivative fo a function of a function. The chain rule is often used to do this: if \(y = f(g(x))\), then
\[ \frac{dy}{dx} = f'(g(x)) \times g'(x). \]
In words: rates of change multiply along a chain. If \(y\) changes 3 times as fast as \(u\), and \(u\) changes 2 times as fast as \(x\), then \(y\) changes 6 times as fast as \(x\).
A logistic regression model is a composition: risk is \(p(x) = \operatorname{expit}(\beta_0 + \beta_1 x)\). The chain rule (using the elegant fact that the expit’s derivative is \(p(1-p)\)) gives the marginal effect of \(x\) on the probability scale:
\[ \frac{dp}{dx} = \beta_1 \, p(x)\,\bigl(1 - p(x)\bigr). \]
This formula explains something important: a constant effect on the log-odds scale is not a constant effect on the risk scale. With \(\beta_1 = 0.08\) per year of age, a person at risk \(p = 0.10\) experiences a marginal effect of \(0.08 \times 0.10 \times 0.90 \approx 0.007\), while a person at \(p = 0.5\) experiences \(0.08 \times 0.25 = 0.02\), nearly three times larger. Note: same coefficient, but different slopes. This is because the chain rule passes \(\beta_1\) through the S-shaped link. (This is precisely the tangent slope drawn in Figure 4.)
The product rule is the companion fact for products of functions: if \(y = u(x) \times v(x)\), then \(y' = u'v + uv'\). This is the derivative of the first times the second, plus the first times the derivative of the second. You will rarely apply it by hand, but you will see it invoked in derivations whenever a formula is a product of moving parts (for instance, survival curves multiplied by densities; weights multiplied by outcomes). It comes into play in the differentiation steps in estimating-equation derivations (Zivich et al. 2026).
Derivatives play a huge role in statistics, particularly with optimization. Optimization in statistics involves choosing parameter values that make some criterion as good as possible. Parameter estimation is optimization:
Derivatives can be used to find these optima. At the bottom of a valley or the top of a hill, a curve is momentarily flat: its slope is zero. This gives the first derivative test (in optimization language, the first-order condition): candidate optima are the points where
\[ f'(\theta) = 0. \]
In other words, the parameter value where the loss function (least squares, maximum likelihood, or some other function) is zero.
The second derivative test distinguishes valley from hill: if \(f''(\theta) > 0\) (the slope is increasing—the curve bends upward), the point is a minimum; if \(f''(\theta) < 0\), a maximum. The second derivative measures curvature, and plays a key role in variance estimation (when possible, as we’ll see in a moment).
Let’s make this concrete with our cohort. Suppose we summarize the five blood pressures with a single number \(c\), and we score any candidate \(c\) by its total squared error, \(L(c) = \sum_{i=1}^{5}(x_i - c)^2\). Figure 5 plots this loss function. Its derivative is \(L'(c) = -2\sum_i (x_i - c)\); setting that to zero gives \(\sum_i (x_i - c) = 0\)—the equation from Section 1.4, whose solution is the sample mean, \(c = 129\). And \(L''(c) = 2n > 0\) confirms it is a minimum. The sample mean is the least-squares summary of a set of numbers, and we have just derived that fact with the two derivative tests.
Lc <- function(c) sapply(c, function(cc) sum((cohort$sbp - cc)^2))
cc <- seq(115, 143, length.out = 300)
d <- data.frame(c = cc, L = Lc(cc))
ggplot(d, aes(c, L)) +
geom_line(colour = pal["blue"], linewidth = 0.9) +
annotate("segment", x = 125, xend = 133, y = Lc(129), yend = Lc(129),
colour = pal["orange"], linewidth = 0.7) +
annotate("point", x = 129, y = Lc(129), colour = pal["orange"], size = 2.4) +
annotate("text", x = 129, y = Lc(129) + 150,
label = "slope = 0 at the minimum:\nc = mean = 129",
size = 3.4, colour = "grey20") +
labs(x = "Candidate summary value c (mmHg)",
y = "Loss: sum of squared errors")
Two further connections make this section central to everything downstream:
Optimization conditions are estimating equations. “Set the derivative of the log-likelihood to zero and solve” produces an equation of exactly the form we met in Section 1.4: an estimator characterized as the value that zeroes an equation (here called the score equation). This is why M-estimators (“M” for maximization) and Z-estimators (“Z” for zero) are two names for the same framework (Zivich et al. 2026).
Curvature is information. The second derivative does more than classify optima. it quantifies how sharply the criterion deteriorates as you move away from the best value. A steeply curved log-likelihood means the data strongly distinguish the best estimate from nearby values (i.e., small uncertainty). A shallow curve means many parameter values fit nearly as well (large uncertainty). This slope-and-curvature reasoning is exactly how Zivich et al. explain why variance estimators involve inverting a derivative matrix: a steeper curvature leads to a smaller variance (Zivich et al. 2026). Standard errors, confidence intervals, and the “information matrix” all flow from the curvature of a criterion function. We will meet the multivariable version (the Hessian) shortly.
Integration is the other half of calculus: accumulation. While the derivative asks “how fast is this changing in a particular instant?”, the integral asks “how much has piled up in total?”
The notation is
\[ \int_a^b f(t) \, dt, \]
read “the integral of \(f\) of \(t\), \(dt\), from \(a\) to \(b\).” The elongated S symbol \(\int\) is an S because it stands for “sum.” The expression instructs us as follows: over the interval from \(a\) to \(b\), chop the axis into tiny slivers of width \(dt\), multiply each sliver’s width by the height of the function there (\(f(t) \times dt\), the area of a thin rectangle), and add all those slivers up. Geometrically, the result is the area under the curve between \(a\) and \(b\). Figure 6 shows the construction: rectangles approximate the area, and as they get thinner the approximation gets better and better. The integral is the limiting value of this approximation.
hz <- function(t) 0.05 + 0.01 * t
rects <- function(width) {
t0 <- seq(0, 10 - width, by = width)
data.frame(xmin = t0, xmax = t0 + width, ymin = 0, ymax = hz(t0),
panel = paste0("Rectangle width = ", width))
}
rr <- rbind(rects(2), rects(0.5))
rr$panel <- factor(rr$panel, levels = c("Rectangle width = 2", "Rectangle width = 0.5"))
curve_d <- do.call(rbind, lapply(levels(rr$panel), function(pn)
data.frame(t = seq(0, 10, length.out = 200), h = hz(seq(0, 10, length.out = 200)), panel = pn)))
curve_d$panel <- factor(curve_d$panel, levels = levels(rr$panel))
ggplot() +
geom_rect(data = rr, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
fill = pal["blue"], alpha = 0.35, colour = "white", linewidth = 0.4) +
geom_line(data = curve_d, aes(t, h), colour = pal["blue"], linewidth = 0.9) +
facet_wrap(~ panel) +
labs(x = "Time t (years)", y = "Hazard h(t)")
Accumulation is central to epidemiology. Consider what our field routinely adds up in small pieces over time:
An integral is a sum in the limit. Whenever you can say “the total is built from many small contributions; a little sliver of time, times the intensity during that sliver, all added up,” you are describing an integral. If the contributions come in discrete chunks (people, strata, days), the integral is just a sum. The \(\int\) symbol generalizes \(\sum\) to quantities that vary continuously.
Integration is fundamental to probability (Naimi et al. 2026).
First, let \(P(A)\) denotes the probability of event \(A\). The vertical bar denotes conditioning: \(P(Y = 1 \mid A = 1)\) reads “the probability that \(Y\) equals 1, given that \(A\) equals 1”: the risk of the outcome within the subgroup of exposed people. Conditioning is the mathematical act of restricting attention to a subgroup, which is why it is the natural language for stratification, adjustment, and prediction.
The expected value (or expectation, or mean) of a random variable \(X\), written \(E(X)\), is its probability-weighted average. For a discrete variable, that is a sum of values times probabilities; for a continuous variable, the sum becomes an integral of values times density:
\[ E(X) = \sum_x x \, P(X = x) \qquad \text{or} \qquad E(X) = \int x \, f(x) \, dx, \]
where \(f(x)\) is the probability density function, playing the role of the weights.
Notice these are the same idea (a weighted average) written for discrete versus continuous variables. A conditional expectation \(E(Y \mid X = x)\) is the average of \(Y\) within the subgroup with \(X = x\); if \(Y\) is a 0/1 outcome, \(E(Y \mid X = x) = P(Y = 1 \mid X = x)\), because the mean of an indicator is a probability (Section 1.8).
Because writing separate formulas for discrete and continuous variables can be tedious, technical papers often use a single unified notation:
\[ E(X) = \int x \, d\mathbb{P}(x), \]
read “the integral of \(x\) with respect to the probability distribution \(\mathbb{P}\).” The differential (the “\(d\)-something” at the end of an integral) states what is supplying the weights. Three variants cover most of what you will encounter (Naimi et al. 2026):
| Notation | Read as | Weights come from | Typical setting |
|---|---|---|---|
| \(\int g(x)\,dx\) | “integral of \(g\), \(dx\)” | equal weight per unit of \(x\) (length) | ordinary calculus; areas |
| \(\int g(x)\,f(x)\,dx\) | “…weighted by the density \(f\)” | a probability density | continuous random variables |
| \(\int g(x)\,d\mathbb{P}(x)\) | “…with respect to the distribution of \(X\)” | the distribution, whatever its type | statistics and causal inference papers |
When \(X\) is categorical, \(\int g(x)\, d\mathbb{P}(x)\) collapses to a weighted sum: \(\sum_x g(x) P(X = x)\).
With that notation in hand, we can state one of the most useful identities in statistics, the law of iterated expectation: for any two random variables,
\[ E(Y) \;=\; E\bigl[\,E(Y \mid X)\,\bigr] \;=\; \int E(Y \mid X = x)\, d\mathbb{P}(x). \]
In words: to get the overall mean of \(Y\), take the mean of \(Y\) within each stratum of \(X\), then average those stratum-specific means, weighting by how common each stratum is (Naimi et al. 2026; Wasserman 2004). The middle expression is the iterative form (an expectation of an expectation; the inner one computed within strata of \(X\), the outer one averaging over \(X\)); the right-hand expression is the non-iterative form, which makes the weighting explicit.
For epidemiologists, this is effectively standardization in its mathematical form. You can verify this in our cohort, with \(Y\) = blood pressure and \(X\) = smoking. Smokers (three of five; weight 0.6) have mean SBP \((142 + 131 + 125)/3 = 132.67\); non-smokers (weight 0.4) have mean \((128+119)/2 = 123.5\). The weighted average of the stratum means is
\[ 132.67 \times 0.6 + 123.5 \times 0.4 = 79.6 + 49.4 = 129, \]
exactly the overall mean from Section 1.7. Splitting a mean into strata and reassembling it with probability weights always returns the original mean (the law of iterated expectation).
w <- mean(cohort$smoker) # P(smoker) = 0.6
m1 <- mean(cohort$sbp[cohort$smoker == 1]) # E(SBP | smoker)
m0 <- mean(cohort$sbp[cohort$smoker == 0]) # E(SBP | non-smoker)
c(iterated = m1 * w + m0 * (1 - w), overall = mean(cohort$sbp))iterated overall
129 129 Note that, by replacing “average within strata of \(X\)” with “average within strata of confounders, at a fixed exposure level”, this law of iterated expectation becomes the g-formula, the central identification tool of causal inference (Robins 1986; Naimi et al. 2026). The important of the law of iterated expectations cannot be understated for the work that we do.
Everything so far has involved functions of one input. But an outcome usually depends on many things at once (age, treatment, blood pressure, smoking, etc) and a likelihood depends on many parameters at once. Multivariable calculus extends “slope” to functions with several inputs. The extension is the machinery of regression coefficient estimation, gradient descent, and variance estimation.
For a function of several inputs, we might ask about change in one input while holding the others fixed. That is a partial derivative, written with the “curly d”:
\[ \frac{\partial f}{\partial x} \]
\(\frac{\partial f}{\partial x}\) reads “the partial derivative of \(f\) with respect to \(x\)”: the rate at which \(f\) changes as \(x\) alone changes, with every other input frozen at its current value. Computationally, you treat the other variables as constants and differentiate as usual.
Consider a linear model for mean blood pressure given age and smoking:
\[ E(\text{SBP} \mid \text{age}, \text{smoke}) = \beta_0 + \beta_1\,\text{age} + \beta_2\,\text{smoke}. \]
Then \(\partial E / \partial \text{age} = \beta_1\): the model says mean SBP rises by \(\beta_1\) mmHg per year of age among people with the same smoking status. A regression coefficient is a partial derivative, and that is the precise meaning of the phrase “holding other variables constant” that accompanies regression interpretations we are often introduced to.
This viewpoint also explains why a coefficient changes when you add a covariate to the model. The partial derivative is a property of the whole regression function. Adding a covariate changes which function you are differentiating. In a model with age alone, \(\partial E/\partial \text{age}\) describes how mean SBP differs across ages, with smokers and non-smokers mixed together as they happen to co-occur with age. In a model with age and smoking, it describes the age slope within smoking groups. If smoking is associated with age, these are genuinely different quantities (crude and adjusted associations).
The mathematics is telling you that “the association between age and the outcome” is not one thing; it depends on what is held fixed.
We can collect all the partial derivatives into vectors and matrices, and these become important objects in the context of regression.
The gradient of a function \(f\) with several inputs, written \(\nabla f\) (the symbol is “nabla,” usually read “the gradient of \(f\)”), is the vector containing all its first partial derivatives (Hardt and Recht 2022). For a two-parameter function \(f(\beta_0, \beta_1)\):
\[ \nabla f = \begin{bmatrix} \dfrac{\partial f}{\partial \beta_0} \\[1.2ex] \dfrac{\partial f}{\partial \beta_1} \end{bmatrix}. \]
The gradient generalizes “slope” to a surface: at any point, it points in the direction of steepest ascent, and its length says how steep the surface is. Two immediate uses:
Optimization in many dimensions. The first-order condition of Section 2.2.2 becomes: at an optimum, every partial derivative is zero at once—\(\nabla f = \mathbf{0}\). Maximum likelihood in a regression with 10 coefficients means solving a system of 10 equations (the score equations) stating that the gradient of the log-likelihood is zero. And when no closed-form solution exists, iterative algorithms can be used to “walk” up or down the “hill” created by these functions in 10 dimensional space. Gradient descent, the engine behind a lot of machine learning, repeatedly nudges the parameters a small step in the direction opposite the gradient of the loss. Figure 7 shows the landscape being navigated for a two-parameter least-squares problem on our cohort: the elliptical contours are level sets of the loss, and the minimum sits where the surface is flat.
grid <- expand.grid(b0 = seq(40, 140, length.out = 120),
b1 = seq(-0.2, 1.85, length.out = 120))
grid$ssq <- with(grid, sapply(seq_len(nrow(grid)), function(k)
sum((cohort$sbp - (b0[k] + b1[k] * cohort$age))^2)))
fit <- lm(sbp ~ age, data = cohort)
ggplot(grid, aes(b0, b1, z = ssq)) +
geom_contour(breaks = c(80, 200, 600, 2000, 6000, 20000),
colour = pal["blue"], linewidth = 0.5) +
annotate("point", x = coef(fit)[1], y = coef(fit)[2],
colour = pal["orange"], size = 2.4) +
annotate("label", x = coef(fit)[1], y = coef(fit)[2] + 0.18,
label = "minimum: gradient = 0", size = 3.4, colour = "grey20",
fill = "white", label.size = 0) +
labs(x = expression("Intercept " * beta[0]),
y = expression("Slope " * beta[1]))
Curvature and uncertainty in many dimensions. The multivariable analog of the second derivative is the Hessian matrix, which is the square matrix collecting all second partial derivatives, written \(\nabla^2 f\). It describes the curvature of the surface in every direction (remember, direction here is potentially multi-dimensional). Its statistical role mirrors the one-dimensional idea of the second derivative: the curvature of the log-likelihood around its maximum determines how precisely the data pin down the parameters, and the variance matrix of the estimates is obtained by inverting a curvature matrix: sharp curvature, small variance; flat curvature, large variance.
In the estimating-equations framework, the “bread” matrix of the sandwich variance estimator is precisely an expected matrix of partial derivatives, \(B_\theta = E[-\nabla_\theta \psi(O_i; \theta)]\), and its inversion is what converts curvature into uncertainty (Zivich et al. 2026; Mansournia et al. 2021).
You may notice we introduced partial derivatives before assembling them into gradients and Hessians, since the gradient is defined as a vector of partials. Keep both in your head as one idea at two zoom levels: the partial derivative is the single-coordinate question (“how does \(f\) respond to this one input?”); the gradient and Hessian are about exactly this issue, but in all directions at once.
The final calculus idea is also commonly used in statistical theory: complicated functions can be approximated by polynomials, and near a point of interest, a short polynomial is usually enough.
This is the basic idea behind Taylor’s theorem. Near a chosen point \(a\), a smooth function satisfies
\[ f(x) \;\approx\; \underbrace{f(a)}_{\text{level}} \;+\; \underbrace{f'(a)\,(x - a)}_{\text{linear term}} \;+\; \underbrace{\tfrac{1}{2} f''(a)\,(x - a)^2}_{\text{quadratic term}} \;+\; \cdots \]
Each term is a correction built from a higher derivative at \(a\): start at the function’s value, tilt by its slope, bend by its curvature, and so on. Extending the sum forever gives an infinit series. Historically, this is the trick that let Newton and his successors tame otherwise impossible functions (Strogatz 2019). Cutting the sum short gives an approximation whose accuracy is very close to \(a\), and decays as you move away. Figure 8 shows the linear and quadratic approximations of \(\log(x)\) around \(x = 1\) hugging the curve near 1 and drifting away from it further out.
x <- seq(0.25, 2.6, length.out = 300)
d <- rbind(
data.frame(x = x, y = log(x), f = "log(x)"),
data.frame(x = x, y = x - 1, f = "linear"),
data.frame(x = x, y = (x - 1) - (x - 1)^2 / 2, f = "quadratic")
)
cols <- c("log(x)" = unname(pal["blue"]), "linear" = unname(pal["orange"]),
"quadratic" = unname(pal["green"]))
lty <- c("log(x)" = "solid", "linear" = "dashed", "quadratic" = "dotted")
lab <- data.frame(x = c(2.65, 2.65, 2.65), y = c(log(2.6), 1.6, 0.32),
f = c("log(x)", "linear", "quadratic"))
ggplot(d, aes(x, y, colour = f, linetype = f)) +
geom_vline(xintercept = 1, colour = ref_grey, linewidth = 0.3) +
geom_line(linewidth = 0.8) +
geom_text(data = lab, aes(label = f), hjust = 0, size = 3.3, show.legend = FALSE) +
scale_colour_manual(values = cols, guide = "none") +
scale_linetype_manual(values = lty, guide = "none") +
coord_cartesian(xlim = c(0.25, 3.05)) +
labs(x = "x", y = "f(x)")
Two pillars of everyday statistical practice are effectively employing Taylor expansions behind the scenes.
The central limit theorem (CLT) says: the sample mean of \(n\) independent observations is approximately normally distributed when \(n\) is large, regardless of the shape of the underlying variable’s distribution. The CLT is obtained with the approximation:
\[ \bar{X} \;\overset{\cdot}{\sim}\; \mathcal{N}\!\left(\mu, \; \frac{\sigma^2}{n}\right), \]
read “\(\bar{X}\) is approximately distributed as normal with mean \(\mu\) and variance \(\sigma^2/n\).” Note that, even if the distribution of \(X\) is itself skewed, lopsided, a rare binary event, a complex zero-truncated function: average enough of them and the average’s distribution is bell-shaped.
set.seed(750)
sim <- do.call(rbind, lapply(c(1, 5, 30), function(n) {
data.frame(n = paste0("Means of n = ", n),
xbar = replicate(4000, mean(rexp(n, rate = 1/5))))
}))
sim$n <- factor(sim$n, levels = paste0("Means of n = ", c(1, 5, 30)))
ggplot(sim, aes(xbar)) +
geom_histogram(bins = 45, fill = pal["blue"], colour = "white", linewidth = 0.2) +
geom_vline(xintercept = 5, colour = pal["orange"], linewidth = 0.5) +
facet_wrap(~ n, scales = "free") +
labs(x = "Sample mean length of stay (days)", y = "Count")
The classical proofs of the CLT work by expanding a transform of the distribution of \(\bar X\) in a Taylor series and observing which terms survive averaging: the mean term and the variance term can play an important role, but everything beyond shrinks faster as \(n\) grows. Only the first two moments (mean and variance) persist, and a distribution characterized by only a mean and a variance, with quadratic structure in its logarithm, is exactly the normal curve (note \(e^{-z^2/2}\): the exponential of a quadratic).
The same logic explains a companion fact: that, near its maximum, a smooth log-likelihood is approximately quadratic (its own Taylor expansion, with the linear term vanishing because the gradient is zero there), so the likelihood itself is approximately a normal curve in the parameter. Note that this is why regression coefficient estimates are assumed to be approximately normal, enabling us to use “estimate \(\pm\) 1.96 standard errors” as the workhorse behind constructing confidence intervals (i.e., the Wald, or normal interval, equation).
The practical consequence is that we rely on the CLT and a quadratic Taylor expansion whenever we construct 95% Wald confidence intervals.
Suppose you have a standard error for one quantity, but you report a function of that quantity. For instance, we estimate log risk ratio as coefficients in a regression model. But we exponentiate those log risk ratios to report risk ratios. Can we simply take the standard error from the regression model for the log risk ratio, and report it as the stanard error of the risk ratio? In other words, how does uncertainty travel through a transformation \(g(\cdot)\)?
Taylor’s linear term is the answer. Near the truth \(\theta\), approximate the transformation by its tangent line, \(g(\hat\theta) \approx g(\theta) + g'(\theta)(\hat\theta - \theta)\). A linear function of a random quantity just rescales its variability by the slope, so
\[ \operatorname{Var}\bigl[g(\hat{\theta})\bigr] \;\approx\; \bigl[g'(\theta)\bigr]^2 \, \operatorname{Var}(\hat{\theta}). \]
This is the delta method: variances propagate through transformations via the squared slope of the transformation. It is a two-line Taylor argument, and it is the standard machinery for standard errors of ratios, differences of transformed coefficients, and marginal effects.
A worked example you will recognize from every epidemiologic software output. Suppose \(\widehat{RR} = 2.0\) with a standard error of \(0.2\) on the log scale. Build the interval where the normal approximation is good—the log scale, then transform:
\[ \log(2.0) \pm 1.96 \times 0.2 = 0.693 \pm 0.392 = (0.301,\; 1.085) \;\;\xrightarrow{\;\exp\;}\;\; (1.35,\; 2.96). \]
This is why confidence intervals for ratio measures are asymmetric around the point estimate: they were symmetric on the log scale, where the statistics behave, and the exponential stretched the upper arm more than the lower one. Logs (Section 1.3), the CLT (Section 2.5.1), and the delta method combine in that one routine calculation.
In contrast, we don’t usually do this (simply because it’s more complicated and, in some respects, less trustworty), but if we wanted the standard error of the risk ratio (i.e., after exponentiating the log risk ratio), we could deploy the delta method directly:
Here the transformation is \(g(\theta) = e^{\theta}\), applied to \(\theta = \log RR\). Its slope is \(g'(\theta) = e^{\theta}\) (recall from Section 2.2 that the exponential is its own derivative) and evaluated at our estimate \(\hat\theta = \log(2.0) = 0.693\), that slope is \(e^{0.693} = 2.0\), which is the risk ratio itself. The delta method then gives us:
\[ \operatorname{Var}\bigl(\widehat{RR}\bigr) \;\approx\; \bigl[e^{\hat\theta}\bigr]^2 \operatorname{Var}\bigl(\log \widehat{RR}\bigr) = (2.0)^2 \times (0.2)^2 = 0.16, \]
so \(\operatorname{SE}\bigl(\widehat{RR}\bigr) \approx \sqrt{0.16} = 0.4\). Because the exponential’s slope at the estimate is the risk ratio, the delta method reduces to:
\[ \operatorname{SE}\bigl(\widehat{RR}\bigr) \;\approx\; \widehat{RR} \times \operatorname{SE}\bigl(\log \widehat{RR}\bigr) = 2.0 \times 0.2 = 0.4. \]
This standard error is perfectly legitimate to report. But what is less trustworthy is the symmetric interval built from it: \(2.0 \pm 1.96 \times 0.4 = (1.22,\; 2.78)\), which differs from the log-scale interval \((1.35,\; 2.96)\) we obtained above.
This is because the sampling distribution of \(\widehat{RR}\) is skewed (bounded by \([0, \infty]\)), so the normal approximation behind the Wald recipe works better on the log scale, where the distribution is closer to symmetric. That (not a defect in the delta method) is why ratio-measure intervals are routinely built on the log scale and then exponentiated. With a larger standard error, the direct symmetric interval can even dip below zero, an impossible value for a risk ratio.
1. If the slope (derivative) of a risk-versus-age curve is positive at age 50, what does that tell you? What would a zero slope at age 50 mean?
A positive derivative at 50 means risk is increasing with age at that point: people slightly older than 50 have higher risk than people slightly younger. A zero slope would mean risk is locally flat at 50—possibly a plateau, a peak, or a trough (the second derivative would distinguish these).
2. What epidemiologic quantity could be viewed as an accumulation (integral) over time? Name two.
Several: person-time (accumulated follow-up), cumulative hazard (accumulated event intensity), cumulative incidence (accumulated risk), pack-years (accumulated exposure), restricted mean survival time (accumulated survival probability). Each is an “area under a curve.”
3. In the logistic model of Section 2.2.1 with \(\beta_1 = 0.08\), why is the marginal effect of age on risk larger for someone at \(p = 0.5\) than for someone at \(p = 0.1\)?
The chain rule gives \(dp/dx = \beta_1 p (1-p)\). The factor \(p(1-p)\) is largest at \(p = 0.5\) (\(0.25\)) and small near the extremes (\(0.09\) at \(p = 0.1\)). The S-shaped expit is steepest in the middle, so the same log-odds effect translates into a bigger risk change there.
4. State the law of iterated expectation in words, using “stratum” somewhere in your answer.
The overall mean of an outcome equals the weighted average of its stratum-specific means, with each stratum weighted by its probability. Compute the mean within strata of \(X\), then average across strata according to how common each stratum is.
Looking ahead: derivatives gave us optimization and marginal effects; integrals gave us expectation and standardization; Taylor gave us the normal approximation and the delta method. What remains is the arithmetic that lets computers do all of this for thousands of observations and dozens of parameters at once: matrices.
Linear algebra, which involves the algebraic manipulation of matrices, is a workhorse of data science and regression. Most datasets are rectangular in structure (e.g., rows for people, columns for variables). Linear algebra is the mathematics of such rectangles. Formulas like \(Y = X\beta + \epsilon\) compress a model for thousands of people into five symbols (you cannot read or write in the methods literature without these conventions). Computationally, when R fits your regression model to your data, the actual work is matrix arithmetic.
This section follows the treatment of matrix operations in Zivich et al. (Zivich et al. 2026), with supporting material from open matrix algebra texts (Hartman 2011; Deisenroth et al. 2020).
A vector is an ordered list of numbers. Order matters: in our cohort, patient 2’s covariates form the row vector
\[ x_2 = \begin{bmatrix} 1 & 62 & 1 \end{bmatrix}, \]
with a fixed convention for what each position means, here, a leading 1 (explained momentarily), then age, then smoking status. A patient is a row vector: one person, summarized as an ordered list of measurements. Written vertically, a list is a column vector; a single variable across all patients (everyone’s age, say) is naturally a column.
A matrix is a rectangular array of numbers, or equivalently, a stack of row vectors or a bundle of column vectors. Stacking all five patients’ covariate rows gives the design matrix:
\[ X = \begin{bmatrix} 1 & 45 & 0 \\ 1 & 62 & 1 \\ 1 & 37 & 0 \\ 1 & 54 & 1 \\ 1 & 41 & 1 \end{bmatrix}. \]
Rows are observations; columns are variables. The dimensions of a matrix are quoted as rows \(\times\) columns: \(X\) is a \(5 \times 3\) matrix (“five by three”), and in general a dataset with \(n\) people and \(p\) variables yields an \(n \times p\) design matrix.
Here, the column of 1s in the leftmost column is the “variable” multiplying the intercept (in this case, a constant input every person shares). This intercept column lets the intercept be handled by the same arithmetic as every other coefficient. Individual entries in a matrix can be identified using double subscripts: \(x_{ij}\) is the entry in row \(i\), column \(j\), so \(x_{22} = 62\) is patient 2’s age.
One operator belongs in this section: the transpose, written \(X^\top\) (or \(X^T\); read “X-transpose”), which flips a matrix over its diagonal so that the rows become the columns and the columns become the rows (Zivich et al. 2026). Our \(5 \times 3\) design matrix has a \(3 \times 5\) transpose: each variable’s column becomes a row. Transposition is important because it plays a role in expressions like \(X^\top X\), coming shortly, which is how “summarize across people” gets written in matrix language.
X <- matrix(
c(1, 45, 0,
1, 62, 1,
1, 37, 0,
1, 54, 1,
1, 41, 1),
nrow = 5, byrow = TRUE
)
X [,1] [,2] [,3]
[1,] 1 45 0
[2,] 1 62 1
[3,] 1 37 0
[4,] 1 54 1
[5,] 1 41 1dim(X) # rows, columns: 5 by 3[1] 5 3t(X) # the transpose: 3 by 5 [,1] [,2] [,3] [,4] [,5]
[1,] 1 1 1 1 1
[2,] 45 62 37 54 41
[3,] 0 1 0 1 1When a methods paper writes “let \(X\) denote the \(n \times p\) matrix of covariates,” it is talking about your data frame, minus the outcome column, plus a column of 1s. Keeping the picture “row = person, column = variable” in mind converts much apparently abstract matrix notation into statements about datasets.
Matrix multiplication is the operation that makes the notation pay off, and it is worth learning as a structured way to combine information: every entry of the result is a weighted sum: each row of the first matrix “meets” each column of the second, multiplying element-by-element and adding up.
The rule, for the \(2 \times 2\) case (Zivich et al. 2026):
\[ A \cdot B = \begin{bmatrix} a_1 & a_2 \\ a_3 & a_4 \end{bmatrix} \begin{bmatrix} b_1 & b_2 \\ b_3 & b_4 \end{bmatrix} = \begin{bmatrix} a_1 b_1 + a_2 b_3 & a_1 b_2 + a_2 b_4 \\ a_3 b_1 + a_4 b_3 & a_3 b_2 + a_4 b_4 \end{bmatrix}. \]
Trace one entry with your finger: the top-left result pairs the first row of \(A\) with the first column of \(B\). Two consequences of the row-meets-column rule (Zivich et al. 2026):
Now, take our design matrix \(X\) (\(5 \times 3\)) and a column vector of coefficients \(\beta = (\beta_0, \beta_1, \beta_2)^\top\), say \(\beta = (100,\ 0.5,\ 8)^\top\), meaning: baseline 100 mmHg, plus 0.5 mmHg per year of age, plus 8 mmHg for smokers.
The product \(X\beta\) is a \(5 \times 1\) vector, and its \(i\)-th entry is row \(i\) of \(X\) (patient \(i\)’s covariates) dotted with the coefficients:
\[ (X\beta)_i = \beta_0 \cdot 1 + \beta_1 \cdot \text{age}_i + \beta_2 \cdot \text{smoke}_i. \]
For patient 1: \(100 + 0.5(45) + 8(0) = 122.5\). For patient 2: \(100 + 0.5(62) + 8(1) = 139\). One matrix multiplication computes every patient’s model-predicted value simultaneously, which is the same weighted sum of covariates, applied to each row at once. That is why the linear regression model for all \(n\) people is written as
\[ Y = X\beta + \epsilon, \]
where \(Y\) is the \(n \times 1\) vector of outcomes, \(X\) the \(n \times p\) design matrix, \(\beta\) the \(p \times 1\) vector of coefficients, and \(\epsilon\) (“epsilon”) the \(n \times 1\) vector of errors (each person’s deviation of their observed outcome from their model value). Five symbols, \(n\) equations. Whether \(n\) is five or five hundred thousand, the notation and computation is identical, which is exactly why matrix notation is the language of regression at scale.
beta <- c(100, 0.5, 8)
drop(X %*% beta) # %*% is matrix multiplication in R; one product, five predictions[1] 122.5 139.0 118.5 135.0 128.5One special matrix is the identity matrix \(I\), with 1s on the main diagonal and 0s elsewhere. It is the matrix analog of the number 1: multiplying by it changes nothing (\(IX = X\)). It matters because it defines what an inverse must accomplish, in Section 3.4.
Matrix multiplication is not elementwise. \(AB\) does not multiply corresponding entries (R’s A * B does that, a different operation from A %*% B). Every entry of a matrix product is a weighted sum across a whole row and column. That is what makes it powerful, but can trip up intuition transferred from ordinary arithmetic.
A square matrix is symmetric if it equals its own transpose, \(A = A^\top\): the entry in row \(i\), column \(j\) equals the entry in row \(j\), column \(i\) (a mirror image across the diagonal). Symmetric matrices arise in statistics for one dominant reason: relationships between pairs of variables are unordered. The covariance of age with blood pressure is the covariance of blood pressure with age, so a table of pairwise covariances is automatically symmetric. In our cohort, the covariance matrix of (age, SBP) is
\[ \Sigma = \begin{bmatrix} 102.7 & 83.8 \\ 83.8 & 72.5 \end{bmatrix}, \]
with variances on the diagonal and the (twice-appearing) covariance on the off-diagonal. Correlation matrices, and the variance-covariance matrices of regression coefficient estimates printed by your software, have the same symmetric structure.
A symmetric matrix \(M\) is positive definite if, for every nonzero vector \(z\),
\[ z^\top M z > 0, \]
and positive semidefinite if the product is \(\geq 0\) (Hardt and Recht 2022). The quantity \(z^\top M z\) (a single number, called a quadratic form) looks opaque until you learn what it computes when \(M\) is a covariance matrix: \(z^\top \Sigma z\) is the variance of the weighted combination \(z_1 X_1 + z_2 X_2 + \cdots\). With \(z = (1, 1)^\top\) and our \(\Sigma\), the arithmetic gives \(102.7 + 72.5 + 2(83.8) = 342.8\), which is the variance of (age + SBP). Since a variance can never be negative, every legitimate covariance matrix must be positive semidefinite. “Positive definite” is thus the matrix generalization of “a positive number”: a matrix that behaves, in every direction \(z\), the way a positive variance should.
S <- cov(cohort[, c("age", "sbp")])
S # symmetric: S equals t(S) age sbp
age 102.70 83.75
sbp 83.75 72.50z <- c(1, 1)
c(quad_form = drop(t(z) %*% S %*% z),
var_sum = var(cohort$age + cohort$sbp)) # same numberquad_form var_sum
342.7 342.7 Positive definiteness is important in three places:
Matrices have no division, but you can sometimes invert them. The inverse of a square matrix \(A\), written \(A^{-1}\) (“A-inverse”), is defined by what it accomplishes: it is the matrix that undoes \(A\), in the sense that
\[ A \, A^{-1} = I, \]
where \(I\) is the identity matrix (Zivich et al. 2026). Just as multiplying by \(1/5\) undoes multiplying by 5, multiplying by \(A^{-1}\) undoes multiplying by \(A\). This is how matrix equations get solved: if \(A b = c\), then \(b = A^{-1} c\).
For a \(2 \times 2\) matrix there is an explicit formula (Zivich et al. 2026):
\[ A^{-1} = \begin{bmatrix} a_1 & a_2 \\ a_3 & a_4 \end{bmatrix}^{-1} = \frac{1}{a_1 a_4 - a_2 a_3} \begin{bmatrix} a_4 & -a_2 \\ -a_3 & a_1 \end{bmatrix}. \]
Look at the scalar out front: its denominator, \(a_1 a_4 - a_2 a_3\), is the determinant of \(A\), written \(\det(A)\) or \(|A|\), a single number computable from any square matrix (geometrically, the factor by which the matrix scales areas or volumes). The formula announces the determinant’s central role: if \(\det(A) = 0\), the inverse requires dividing by zero and does not exist. Such a matrix is called singular. A zero determinant is the matrix version of “you cannot divide by zero,” and, as we’ll see, it is the mathematical signature of redundant information.
A quick numerical check, with \(A = \begin{bmatrix} 2 & 1 \\ 1 & 3 \end{bmatrix}\): the determinant is \(2(3) - 1(1) = 5\), so \(A^{-1} = \tfrac{1}{5}\begin{bmatrix} 3 & -1 \\ -1 & 2 \end{bmatrix} = \begin{bmatrix} 0.6 & -0.2 \\ -0.2 & 0.4 \end{bmatrix}\), and multiplying \(A A^{-1}\) indeed returns the identity:
A <- matrix(c(2, 1, 1, 3), nrow = 2, byrow = TRUE)
det(A)[1] 5solve(A) # the inverse [,1] [,2]
[1,] 0.6 -0.2
[2,] -0.2 0.4round(A %*% solve(A), 10) # A times its inverse: the identity [,1] [,2]
[1,] 1 0
[2,] 0 1Where inverses star epidemiology and statistics:
Ordinary least squares in closed form. The coefficients that minimize \(\sum_i (y_i - x_i^\top \beta)^2\) are
\[ \hat{\beta} = (X^\top X)^{-1} X^\top y, \]
read aloud: “X-transpose-X, inverse, times X-transpose-y.” In fact, this is the ordinary least squares estimator written out in matrix form. From right to left: \(X^\top y\) summarizes how the covariates co-move with the outcome; \(X^\top X\) (\(p \times p\), symmetric matrix) summarizes how the covariates co-move with each other; and the inverse divides out that among-covariate overlap, so each coefficient reflects its variable’s own contribution. This formula is also the algebra behind an elegant geometric fact: the fitted values \(X\hat\beta\) are the projection of the outcome vector onto the space of all outcomes the covariates can express, which is the closest the model can get to the data. That geometry explains a familiar phenomenon: adding a covariate enlarges the space being projected onto, so all the coefficients can shift, not just the new one [only when the new covariate is unrelated (orthogonal, in the language of linear algebr) to the others do the old coefficients stay put].
The sandwich variance. The robust variance estimator that accompanies estimating equations is
\[ V_\theta = B_\theta^{-1} \, M_\theta \, \bigl(B_\theta^{-1}\bigr)^\top, \]
the “bread-meat-bread” sandwich: \(M_\theta\) (the meat) captures the variability of the data’s contributions, and \(B_\theta^{-1}\)—the inverted curvature matrix from Section 2.4.2—converts that variability into uncertainty about the parameters, exactly as promised by the “steeper curvature, smaller variance” principle (Zivich et al. 2026; Mansournia et al. 2021). You can now read every symbol in that formula: a transpose, two inverses, and two matrix products, each doing a job you’ve seen.
(A practical footnote: modern software rarely computes \(A^{-1}\) explicitly—it solves the associated equations by more stable routes. The inverse remains how we humans write and reason about the operation.)
When does the inverse fail, and what does failure mean in an actual analysis? The concepts that answer these questions are rank and null space.
The rank of a matrix is the number of genuinely independent columns it contains. Columns carrying information not already expressible as a weighted combination of other columns. A matrix whose columns are all independent has full rank; if some column duplicates information, the matrix is rank-deficient, its determinant is 0, and \((X^\top X)^{-1}\) does not exist.
It’s easy to see this in an example. Suppose we code sex with two indicator columns, an indicator for male and one for female, and we keep the intercept’s column of 1s:
\[ X = \begin{bmatrix} 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 1 & 0 \end{bmatrix} \qquad \text{male} + \text{female} = \text{intercept column, in every row.} \]
The third column is exactly the first minus the second: perfect multicollinearity. Only two of the three columns carry independent information (rank 2), \(X^\top X\) is singular, and the regression cannot be fit as written. The same trap arises whenever a full set of category indicators (all race categories, all sites, all months) enters a model with an intercept (the indicators sum to the intercept column, which is why software drops one category as the reference).
The null space of \(X\) is the set of coefficient vectors \(v\) for which \(Xv = \mathbf{0}\): recipes for combining the columns into nothing. Full-rank matrices have only the trivial recipe (\(v = \mathbf{0}\)); a rank-deficient matrix has real ones—here, \(v = (1, -1, -1)^\top\), since intercept minus male minus female is zero for every person.
Statistically, this is important: if \(Xv = \mathbf{0}\), then \(X(\beta + cv) = X\beta\) for any constant \(c\)—infinitely many different coefficient vectors produce identical predictions for every individual, so the data cannot distinguish among them. There is no unique \(\hat\beta\); the parameters are not identifiable from the design.
Rank deficiency is a statement that your question, as parameterized, has no unique answer. R responds by reporting NA for the redundant coefficient:
male <- c(1, 0, 1, 0, 1)
female <- 1 - male
coef(lm(sbp ~ male + female, data = cohort)) # 'female' is NA: rank-deficient design(Intercept) male female
136.5 -12.5 NA Near-deficiency matters too: when columns are almost linearly dependent (age and year of birth, say), the determinant is near zero, the inverse is numerically volatile, and coefficient variances balloon. This is the practical consequence of multicollinearity.
The last piece of vocabulary measures the size of a vector. A norm, written with double bars \(\lVert v \rVert\), converts a vector into a single nonnegative number summarizing its magnitude.
Two norms dominate applied statistics (but there are others), and it is worth knowing both because they behave differently in an important way. For \(v = (v_1, \ldots, v_p)\):
\[ \lVert v \rVert_2 = \sqrt{\textstyle\sum_j v_j^2} \qquad \text{(the } L_2 \text{, or Euclidean, norm)} \]
\[ \lVert v \rVert_1 = \textstyle\sum_j |v_j| \qquad \text{(the } L_1 \text{ norm)}. \]
For \(v = (3, 4)\): \(\lVert v \rVert_2 = \sqrt{9 + 16} = 5\) (basically, a straight-line length, or the Pythagorean theorem generalized), while \(\lVert v \rVert_1 = 3 + 4 = 7\) (for instance, total blocks walked on a city grid, hence the commonly used name “Manhattan norm”).
The distance between two points is the norm of their difference, \(\lVert u - v \rVert\), which is how the question: “how similar are these two patients’ covariate profiles?” becomes a computable number (this is the quantity that motivates nearest-neighbor prediction, clustering algorithms such as K means, and some matching algorithms).
One practical consequence should be noted: variables measured on wildly different scales should be standardized before computing distances, lest age in years be dwarfed by blood pressure in mmHg.
Norms also unify things you have already seen. Least squares estimation (Section 2.2.2) is, in fact, norm-minimization: \(\sum_i (y_i - x_i^\top\beta)^2 = \lVert y - X\beta \rVert_2^2\), so “minimize the squared error” reads “make the residual vector as short as possible, in the \(L_2\) sense.”
The modern reason norms deserve their own subsection, though, is penalization. High-dimensional models, or models with hundreds or thousands of covariates (some EHR-based, environmental, genomic, or metabolomic data) will overfit if fit by plain least squares or maximum likelihood. One remedy is to add a penalty on the size of the coefficient vector:
\[ \hat\beta = \arg\min_\beta \; \Bigl\{ \lVert y - X\beta \rVert_2^2 \;+\; \lambda \lVert \beta \rVert \Bigr\}, \]
One reading of this equation is: the estimate \(\hat{\beta}\) is the vector of \(\beta\) values that makes the model’s predictions \(X\beta\) as close as possible to the observed outcomes \(y\) (minimizing the squared \(L_2\) distance between the two) while also keeping the coefficient vector itself “short”, so to speak, as measured by its norm (whichever we happen to choose, more below).
The two terms pull in opposite directions: fitting the data more closely generally requires larger coefficients, but the \(\arg\min\) will settle at a compromise where further gains in fit (as measured by \(lVert y - X\beta \rVert_2^2\)) are no longer worth the growth in coefficient size (as measured by \(\lambda \lVert \beta \rVert\)).
Here \(\lambda\) (“lambda”) controls the penalty’s strength, and the choice of norm for the penalty defines the method: penalizing \(\lVert \beta \rVert_2^2\) gives ridge regression, which shrinks all coefficients smoothly toward zero; penalizing \(\lVert \beta \rVert_1\) gives the lasso, which shrinks and selects by setting some coefficients exactly to zero.
The geometric reason is visible in Figure 10: the \(L_1\) ball has corners on the axes, and corners are where solutions land, so lasso solutions sit at points where some coordinates are exactly zero. Meanwhile the smooth \(L_2\) ball has no corners, so ridge shrinks without zeroing. This is a classic figure in machine learning (e.g., Figure 3.11 in Hastie et al. (2009)), and it can be a little cryptic to understand. However, we will cover this in a later section.
th <- seq(0, 2 * pi, length.out = 200)
l2 <- data.frame(x = cos(th), y = sin(th))
l1 <- data.frame(x = c(1, 0, -1, 0, 1), y = c(0, 1, 0, -1, 0))
ggplot() +
geom_hline(yintercept = 0, colour = ref_grey, linewidth = 0.3) +
geom_vline(xintercept = 0, colour = ref_grey, linewidth = 0.3) +
geom_path(data = l2, aes(x, y), colour = pal["blue"], linewidth = 0.9) +
geom_path(data = l1, aes(x, y), colour = pal["orange"], linewidth = 0.9) +
annotate("text", x = 0.82, y = 0.82, label = "L2: circle (ridge)",
colour = pal["blue"], size = 3.5, hjust = 0) +
annotate("text", x = 0.74, y = -0.78, label = "L1: diamond (lasso)",
colour = pal["orange"], size = 3.5, hjust = 0) +
annotate("point", x = c(1, 0), y = c(0, 1), colour = pal["orange"], size = 2) +
coord_fixed(xlim = c(-1.35, 1.6), ylim = c(-1.25, 1.25)) +
labs(x = expression(v[1]), y = expression(v[2]))
1. In the design matrix \(X\) of Section 3.1, what does one row represent? One column? The entry \(x_{42}\)?
A row is one patient’s covariate values (in a fixed order); a column is one variable’s values across all patients; \(x_{42}\) is the entry in row 4, column 2—patient 4’s age, 54.
2. Why might matrix notation be preferable when fitting the same regression model to 100,000 observations?
Because \(Y = X\beta + \epsilon\) and \(\hat\beta = (X^\top X)^{-1}X^\top y\) are the same expressions whether \(n\) is 5 or 100,000: one line of notation (and one set of matrix routines) replaces 100,000 person-level equations. The notation scales because matrix multiplication applies the same weighted sum to every row simultaneously.
3. A colleague’s regression includes an intercept plus indicator variables for all four study sites, and the software returns an NA coefficient. Diagnose the problem in the language of this section.
The four site indicators sum to the intercept’s column of 1s, so the design matrix is rank-deficient (perfect multicollinearity): a nonzero vector in the null space means infinitely many coefficient vectors give identical predictions, \(X^\top X\) is singular, and no unique solution exists. Drop one site as the reference category (or drop the intercept).
4. Your software reports that the estimated covariance matrix of \((\hat\beta_1, \hat\beta_2)\) is symmetric. Why must it be? And what would it mean, informally, if it were not positive semidefinite?
Symmetric because the covariance of \(\hat\beta_1\) with \(\hat\beta_2\) is the same quantity as the covariance of \(\hat\beta_2\) with \(\hat\beta_1\). If it were not positive semidefinite, some weighted combination of the estimates would have negative variance (\(z^\top \Sigma z < 0\)), which is an impossibility, signaling a defective estimate of uncertainty.
The purpose of this final section is practical. Methods papers and, sometimes, applied papers, communicate their central ideas in equations. The notation in most epidemiology, biostatistics, and machine learning papers draws on a modest shared vocabulary, nearly all of which you have now seen. What remains is to consolidate the conventions, practice on real equations, and leave you with a repeatable strategy for formulas you have never seen before.
The table below collects the conventions that cover most equations you will see:
| Symbol | How to read it | Notes |
|---|---|---|
| \(X\) vs. \(x\) | “big X” (random variable) vs. “little x” (a realized value) | \(P(X = x)\) links them |
| \(x_i\) | “x sub i” | value for observation \(i\); subscripts index, superscripts often mean powers—context decides |
| \(\bar{x}\) | “x-bar” | sample mean; but in longitudinal causal papers, \(\bar{a}_t\) is a history \((a_0, \ldots, a_t)\), context will decide |
| \(\hat{\theta}\) | “theta-hat” | an estimate of parameter \(\theta\) computed from data |
| \(\beta, \mu, \sigma, \lambda, \theta\) | Greek letters | unknown parameters (coefficients, means, spreads, penalties) |
| \(Y^a\) | “Y under a” | potential outcome: the outcome if exposure were set to \(a\) |
| \(P(A)\), \(P(Y{=}1 \mid A{=}1)\) | “probability of…”, “given” | the bar \(\mid\) restricts to a subgroup |
| \(E(Y)\), \(E(Y \mid X)\) | “expected value of…” | population mean; conditional mean within strata of \(X\) |
| \(\operatorname{Var}, \operatorname{Cov}\) | variance, covariance | spread; co-movement of pairs |
| \(\sim\) | “is distributed as” | \(Y \sim \mathcal{N}(\mu, \sigma^2)\): normal with that mean and variance |
| \(\mathbf{1}[\cdot]\), \(I(\cdot)\) | “the indicator that…” | 1 if true, 0 if false |
| \(\sum, \prod\) | “sum over…”, “product over…” | check the index and its range |
| \(\int \cdots \, d\mathbb{P}(x)\) | “averaged over the distribution of X” | a weighted average in disguise (Table 2) |
| \(f'(x)\), \(\frac{d}{dx}\), \(\partial\), \(\nabla\) | derivatives | slope; partial slope; vector of partial slopes |
| \(X^\top\), \(A^{-1}\) | “transpose”, “inverse” | flip rows/columns; the undoing matrix |
| \(\lVert v \rVert_1, \lVert v \rVert_2\) | “L-one/L-two norm” | vector size: sum of absolute values; Euclidean length |
| \(\arg\max_\theta, \arg\min_\theta\) | “the value of theta that maximizes/minimizes” | returns the location of the optimum, not its height |
| \(\in\), \(\{\cdot\}\) | “is in”; a set | \(A \in \{0,1\}\): A is binary |
| \(\perp\!\!\!\perp\) | “is independent of” | \(Y^a \perp\!\!\!\perp A \mid W\): exchangeability statements |
| i.i.d. | “independent and identically distributed” | the standard sampling assumption |
Two warnings about conventions. First, they are conventions, not laws: some authors use \(f\) for densities and others for generic functions; \(\bar{a}\) means “mean” in one literature and “history” in another. Good papers define their notation early, but sometimes you do just have to figure it out based on context. Second, notation overloading is normal, not a sign you have misunderstood.
\[ \operatorname{logit}\bigl\{ P(Y = 1 \mid A, W) \bigr\} = \beta_0 + \beta_1 A + \beta_2 W \]
Left side, inside out: \(P(Y=1 \mid A, W)\) is the risk of the outcome among people with exposure level \(A\) and covariate level \(W\); the logit maps that risk from \((0,1)\) onto the whole real line (Section 1.3). The right side is a linear function of exposure and covariate. So the model asserts: the log-odds of the outcome is linear in \(A\) and \(W\). Everything you know about the pieces now activates: \(\beta_1\) is a partial derivative on the log-odds scale (\(\partial/\partial A\), holding \(W\) fixed; Section 2.4.1), so \(e^{\beta_1}\) is a conditional odds ratio (Section 1.3); predicted risks come from inverting the logit with the expit (Section 1.4); the marginal effect on the risk scale requires the chain rule and varies across people (Section 2.2.1); and fitting the model means maximizing a log-likelihood, an optimization problem whose first-order condition sets a gradient to zero (Section 2.4.2), computed with the design matrix \(X = [\,1 \;\; A \;\; W\,]\) (Section 3.1).
The g-formula (Robins 1986; Naimi et al. 2017; Hernan and Robins 2025) identifies the mean of a potential outcome \(E(Y^a)\), or the average outcome if everyone’s exposure were set to \(a\). From observed data, we have:
\[ E(Y^a) \;=\; \sum_w E(Y \mid A = a, W = w)\, P(W = w). \]
Read it aloud: “the mean of \(Y\) that would be observed under exposure level \(a\) equals the sum, over confounder strata \(w\), of the observed mean outcome among people with exposure \(a\) in stratum \(w\), weighted by the probability of being in stratum \(w\).” Structurally, you already know this object: it is a weighted average of stratum-specific means (the non-iterative form of the law of iterated expectation), with one edit: the inner expectation conditions on \(A = a\), but the weights \(P(W=w)\) come from the whole population, not from the exposed subgroup. That edit is standardization: it asks what the average outcome would be if the \(A = a\) group’s stratum-specific risks applied to everyone.
The law of iterated expectation is a mathematical identity, true for any random variables. The g-formula is a causal claim: it equates the identity’s right-hand side with the counterfactual quantity \(E(Y^a)\), and that equation holds only under causal identification assumptions, including causal consistency (the observed outcome equals the potential outcome under the exposure actually received), conditional exchangeability (within strata of \(W\), the exposed and unexposed are comparable, no unmeasured confounding), and positivity (every stratum contains people at each exposure level). This is a critical diffrence between the standard LIE, and the LIE used to derive the g formula (Naimi et al. 2026).
Let’s compute it, using the data that Sato and Matsuyama used to illustrate standardization (Sato and Matsuyama 2003), as analyzed in the primer this section follows (Naimi et al. 2026): 4,901 women with breast cancer, where \(A\) is tamoxifen treatment, \(Y\) is cancer recurrence, and \(W\) indicates positive lymph node metastasis at surgery, a cause of both treatment choice and recurrence, i.e., a confounder:
flowchart LR W((W)) --> A((A)) W((W)) --> Y((Y)) A((A)) --> Y((Y))
tam <- data.frame(
W = c(0,0,0,0,1,1,1,1),
A = c(0,0,1,1,0,0,1,1),
Y = c(0,1,0,1,0,1,0,1),
N = c(1421, 171, 1238, 96, 507, 253, 847, 368)
)
kable(tam)| W | A | Y | N |
|---|---|---|---|
| 0 | 0 | 0 | 1421 |
| 0 | 0 | 1 | 171 |
| 0 | 1 | 0 | 1238 |
| 0 | 1 | 1 | 96 |
| 1 | 0 | 0 | 507 |
| 1 | 0 | 1 | 253 |
| 1 | 1 | 0 | 847 |
| 1 | 1 | 1 | 368 |
The formula needs two kinds of ingredients, all computable from the table. The stratum-specific mean outcomes among the treated (recall: for binary \(Y\), a mean is a risk):
\[ \hat{E}(Y \mid A=1, W=1) = \frac{368}{847 + 368} = 0.303, \qquad \hat{E}(Y \mid A=1, W=0) = \frac{96}{1238 + 96} = 0.072, \]
and the confounder distribution in the whole cohort: \(\hat{P}(W=1) = (507+253+847+368)/4901 = 0.403\), so \(\hat{P}(W=0) = 0.597\). Plug in:
\[ \hat{E}(Y^{a=1}) = 0.303 \times 0.403 + 0.072 \times 0.597 \approx 0.165. \]
Repeating with the untreated means (\(0.333\) and \(0.107\)) gives \(\hat{E}(Y^{a=0}) \approx 0.198\), and the estimated causal risk difference is \(0.165 - 0.198 \approx -0.03\): under the identification assumptions, tamoxifen lowers recurrence risk by about 3 percentage points. R confirms the hand arithmetic:
risk <- function(a, w) with(tam, N[A==a & W==w & Y==1] / sum(N[A==a & W==w]))
pW1 <- with(tam, sum(N[W==1]) / sum(N))
EY1 <- risk(1,1) * pW1 + risk(1,0) * (1 - pW1)
EY0 <- risk(0,1) * pW1 + risk(0,0) * (1 - pW1)
round(c(EY_a1 = EY1, EY_a0 = EY0, risk_difference = EY1 - EY0), 3) EY_a1 EY_a0 risk_difference
0.165 0.198 -0.033 One more literacy point: the same quantity can be written in the iterative form \(E(Y^a) = E[\,E(Y \mid A=a, W)\,]\), which are nested expectations, computed in software by predicting each person’s outcome under \(A = a\) and then averaging the predictions. The two forms are algebraically equivalent expressions of one estimand, differing in computational strategy, but not in meaning (Naimi et al. 2026).
Survival papers routinely display the object that Cox regression maximizes (Cox 1972):
\[ L(\beta) \;=\; \prod_{i \,:\, \delta_i = 1} \; \frac{\exp(x_i^\top \beta)}{\displaystyle\sum_{j \,\in\, R(t_i)} \exp(x_j^\top \beta)}. \]
Let’s dissect this piece by piece:
So the equation reads: the partial likelihood is the product, over observed events, of the probability that the person who failed was the one who did, among all those at risk at that moment.
It’s not easy to see, but something important is absent: the baseline hazard. It’s not easy to see this because the \(x\) vector in \(x_i^\top \beta\) DOES NOT include a 1 at the beginning for the intercept. This is because the underlying baseline event rate is shared by everyone, and is thus canceled out of every numerator and denominator (it multiplies both). That cancellation is the entire genius of the method: \(\beta\) (the vector of coefficients for each covariate, without the interecept) is estimated without ever modeling how the hazard evolves over time, which is why the likelihood is called “partial,” and this was Cox’s famous discovery that made his regression model possible (Cox 1975).
From here the workflow is the standard one: take logs to convert the product to a sum (Section 1.3), differentiate to get a score (Section 2.4.2), set it to zero, and let the curvature at the maximum supply standard errors (Section 2.2.2).
Finally, an equation in the context of the machine learning literature, namely, the LASSO (Hardt and Recht 2022; Deisenroth et al. 2020):
\[ \hat{\beta} \;=\; \arg\min_{\beta} \; \Biggl\{ \; \frac{1}{n}\sum_{i=1}^{n} \bigl(y_i - x_i^\top \beta\bigr)^2 \;+\; \lambda \sum_{j=1}^{p} |\beta_j| \; \Biggr\}. \]
Read: “\(\hat\beta\) is the coefficient vector that minimizes the average squared prediction error plus \(\lambda\) times the \(L_1\) norm of the coefficients.” Every component is now familiar: \(\arg\min\) (Table 3) announces an optimization (Section 2.2.2); the first term is the least-squares loss, an average of squared residuals built from linear predictors (Section 3.2); the second is the \(L_1\) penalty whose corner geometry sets some coefficients exactly to zero (Section 3.6); and \(\lambda\) is a tuning parameter trading fit against simplicity.
Papers will show you formulas these notes did not. Here is a strategy/checklist to run, slowly, whenever an equation is challenging your interpretation:
mean() an expectation. It also doubles as a check on step 4.A demonstration, on an estimator we have not covered, the inverse-probability-weighting (IPW) estimator of a counterfactual mean:
\[ \hat{E}(Y^a) \;=\; \frac{1}{n} \sum_{i=1}^{n} \frac{\mathbf{1}[A_i = a] \; Y_i}{\hat{P}(A_i = a \mid W_i)}. \]
Running the checklist: the symbols are data (\(A_i, Y_i, W_i\)), one estimated ingredient (\(\hat P(A_i{=}a \mid W_i)\), each person’s estimated probability of the exposure they are being evaluated at, given their covariates, i.e., the propensity score), an index \(i\), and operators. Inside out: the indicator switches off everyone whose exposure is not \(a\); surviving people contribute their outcome, divided by their probability of having the exposure they in fact had; the sum-and-divide-by-\(n\) is an average over the whole sample. Skeleton: a weighted average, with weights \(1/\hat P\). And the small case delivers the intuition: a person in a covariate stratum where exposure \(a\) was unlikely (\(\hat P = 0.1\)) counts as 10 copies of themselves, reconstructing, from the exposed alone, the outcomes of the full population as if everyone had received \(a\). The formula has become a sentence: upweight the observed \(a\)-group to stand in for everybody, then average. (You can even see why positivity matters: a \(\hat{P}\) of zero would put a zero in a denominator.)
1. In words: what does \(P(Y = 1 \mid A = 1)\) mean? How does it differ from \(P(Y=1)\) and from \(E(Y^{a=1})\)?
\(P(Y=1 \mid A=1)\) is the risk of the outcome among the exposed subgroup (an observed, conditional quantity). \(P(Y=1)\) is the risk in the whole population. \(E(Y^{a=1})\) is the counterfactual risk if everyone were exposed, generally different from the first quantity unless conditions like exchangeability hold; connecting the two is exactly what the g-formula does.
2. In the Cox partial likelihood, what is the risk set \(R(t_i)\), and why do censored individuals still matter to the estimate?
\(R(t_i)\) is the set of people still event-free and under observation just before event time \(t_i\). Censored people never contribute a numerator factor (no observed event), but until their censoring time they appear in the risk-set denominators, where they represent survival experience. Their presence changes every event’s conditional probability.
3. Apply steps 1–2 and 6 of the checklist to \(\hat\theta = \arg\max_\theta \sum_{i=1}^n \log f(y_i; \theta)\).
Symbols: \(y_i\) data; \(\theta\) parameter; \(i\) index; \(f\) a density function; \(\log\), \(\sum\), \(\arg\max\) operators. Inside out: person \(i\)’s probability contribution, logged; summed over people (a log-likelihood, i.e., a logged product); then the \(\theta\) that maximizes the sum. Skeleton: an \(\arg\max\) of a criterion. This is maximum likelihood estimation.
Note these are conventions, meaning these notational tools can be used for other purposes↩︎