A mechanistic model, like constrained growth or SIR, explains a system. Each term stands for something real, so we can ask “what if?” questions about it.
An empirical model is a function chosen only because it follows the trend of the data. The book’s definition (p. 353): an empirical model “is based only on data and is used to predict, not explain, a system.” It can be very useful, but only where we have data. Outside that range, the data have nothing to say.
A load cell is a sensor that bends slightly under weight. Engineers at NIST hung known loads on one and recorded the deflection, twice at each of 20 loads. If the relationship is reliable, we can turn any deflection reading into a weight.
pontius <- read.csv("pontius.csv")
head(pontius) load_thousands deflection
1 150 0.11019
2 300 0.21956
3 450 0.32949
4 600 0.43899
5 750 0.54803
6 900 0.65694nrow(pontius)[1] 40plot(pontius$load_thousands, pontius$deflection, pch = 16, col = "steelblue",
xlab = "Load (thousands)", ylab = "Deflection",
main = "Load cell calibration (NIST Pontius data)")
That looks like a line. R’s lm() (“linear model”) finds the least-squares line: the slope and intercept that make the sum of the squared vertical distances from the points to the line as small as possible (pp. 356–357). The formula deflection ~ load_thousands reads “deflection as a function of load.”
line_fit <- lm(deflection ~ load_thousands, data = pontius)
coef(line_fit) (Intercept) load_thousands
0.0061496842 0.0007221026 The slope says each additional 1,000 units of load adds about 0.00072 to the deflection. To draw the fitted line, ask the model for predictions at many loads with predict():
plot(pontius$load_thousands, pontius$deflection, pch = 16, col = "steelblue",
xlab = "Load (thousands)", ylab = "Deflection",
main = "Least-squares line")
load_grid <- seq(0, 3000, by = 10)
lines(load_grid, predict(line_fit, data.frame(load_thousands = load_grid)),
col = "firebrick", lwd = 2)
predict(line_fit, data.frame(load_thousands = c(1000, 10000))) 1 2
0.7282523 7.2271755 A load of 1,000 is inside the range we measured (150 to 3,000), so we have good reason to trust that prediction. A load of 10,000 is more than three times the largest load tested. The line happily gives an answer, but would the real sensor stay linear, or bend, or break? The data can’t tell us. Using a model beyond the range of its data is called extrapolation.
The line looks perfect. Is it? A residual is the part of each observation the model doesn’t explain:
\[ \text{residual} = \text{observed } y - \text{predicted } y. \]
If a model has captured the trend, its residuals should look like random scatter around zero. If they show a pattern, the model is missing something.
plot(pontius$load_thousands, resid(line_fit), pch = 16, col = "steelblue",
xlab = "Load (thousands)", ylab = "Residual",
main = "Residuals from the straight line")
abline(h = 0, lty = 2)
Those residuals form a clear arch: the line is too low in the middle and too high at both ends. The errors are tiny (thousandths), but they aren’t random. The real relationship curves slightly.
lm() can fit any linear combination of terms (p. 354), not just \(b + mx\). Adding an \(x^2\) term gives a quadratic, \(y = b_0 + b_1 x + b_2 x^2\). Inside a formula, I(load_thousands^2) tells R to square the values.
quad_fit <- lm(deflection ~ load_thousands + I(load_thousands^2), data = pontius)
coef(quad_fit) (Intercept) load_thousands I(load_thousands^2)
6.735658e-04 7.320592e-04 -3.160819e-09 plot(pontius$load_thousands, resid(quad_fit), pch = 16, col = "darkgreen",
xlab = "Load (thousands)", ylab = "Residual",
main = "Residuals from the quadratic")
abline(h = 0, lty = 2)
The arch is gone, and the residuals are about ten times smaller. What’s left looks like measurement noise. The \(x^2\) coefficient is tiny, but it matters at this level of precision. Notice that we found it by looking at the residuals, not at the original plot.
Our next data relate the temperature of a carbon lamp filament (in thousands of kelvins) to the energy it radiates per square centimeter per second.
danwood <- read.csv("danwood.csv")
danwood temperature energy
1 1.309 2.138
2 1.471 3.421
3 1.490 3.597
4 1.565 4.340
5 1.611 4.882
6 1.680 5.660plot(danwood$temperature, danwood$energy, pch = 16, cex = 1.3,
col = "darkorange", xlab = "Temperature (1000 K)",
ylab = "Energy radiated", main = "Filament data (NIST DanWood)")
lines(range(danwood$temperature), range(danwood$energy), lty = 2, col = "gray50")
The dashed line through the first and last points shows that the data curve upward (concave up). The book’s idea (pp. 357–360): find a power \(p\) so that plotting \(y\) against \(x^p\) looks straight. Raising values greater than 1 to a power stretches the larger ones more, which straightens a concave-up pattern. Then fit a line to \((x^p, y)\).
Rather than trying powers by hand, let a loop do it. For each power, fit the line and record the sum of squared residuals (smaller is better):
powers <- seq(1, 5, by = 0.5)
ssr <- numeric(length(powers))
for (k in seq_along(powers)) {
z <- danwood$temperature^powers[k]
fit <- lm(danwood$energy ~ z)
ssr[k] <- sum(resid(fit)^2)
}
data.frame(power = powers, sum_sq_residuals = round(ssr, 5)) power sum_sq_residuals
1 1.0 0.07461
2 1.5 0.04687
3 2.0 0.02567
4 2.5 0.01106
5 3.0 0.00300
6 3.5 0.00142
7 4.0 0.00618
8 4.5 0.01709
9 5.0 0.03390The power 3.5 gives the smallest sum, about 50 times smaller than a straight line (\(p = 1\)). It’s the same power the textbook picks by looking at plots. (Try a finer grid between 3 and 4. Does the best power move much?) Let’s use it:
danwood$temp_power <- danwood$temperature^3.5
power_fit <- lm(energy ~ temp_power, data = danwood)
coef(power_fit)(Intercept) temp_power
-0.3931313 0.9881859 These match the book: \(y = -0.393 + 0.988\,z\) with \(z = x^{3.5}\). Substituting back gives the empirical model
\[ f(x) = -0.393 + 0.988\,x^{3.5}. \]
Always check the model against the original, untransformed data. Here’s the curve with the measurements:
temp_grid <- seq(1.25, 1.75, by = 0.005)
plot(danwood$temperature, danwood$energy, pch = 16, cex = 1.3,
col = "darkorange", xlab = "Temperature (1000 K)",
ylab = "Energy radiated", main = "Empirical model f(x) = -0.393 + 0.988 x^3.5")
lines(temp_grid, -0.393 + 0.988 * temp_grid^3.5, col = "firebrick", lwd = 2)
A physics aside: a law from physics (the Stefan–Boltzmann law) says an ideal glowing body radiates energy in proportion to \(T^4\). Our data alone suggested a power near 3.5, close to but not equal to 4. A real filament isn’t an ideal body, and six points can’t pin the exponent down precisely. This is the empirical approach’s strength and weakness at once: it follows these data, without knowing why.
Biochemical oxygen demand (BOD) measures how much oxygen bacteria use while breaking down waste in a water sample. These six measurements track BOD over 10 days of incubation.
boxbod <- read.csv("boxbod.csv")
boxbod days bod
1 1 109
2 2 149
3 3 149
4 5 191
5 7 213
6 10 224With six points, a polynomial of degree 5 has six coefficients, so it can pass through every point exactly. Its residuals are all zero (up to rounding). Compare it with a quadratic:
quad_bod <- lm(bod ~ days + I(days^2), data = boxbod)
poly5_bod <- lm(bod ~ days + I(days^2) + I(days^3) + I(days^4) + I(days^5),
data = boxbod)
round(resid(poly5_bod), 6)1 2 3 4 5 6
0 0 0 0 0 0 day_grid <- seq(0.5, 11, by = 0.05)
plot(boxbod$days, boxbod$bod, pch = 16, cex = 1.3, col = "purple",
xlim = c(0.5, 11), ylim = c(0, 400),
xlab = "Incubation time (days)", ylab = "BOD (mg/L)",
main = "Two empirical models for BOD")
lines(day_grid, predict(quad_bod, data.frame(days = day_grid)),
col = "steelblue", lwd = 2)
lines(day_grid, predict(poly5_bod, data.frame(days = day_grid)),
col = "firebrick", lwd = 2)
legend("topleft", c("Quadratic", "Degree 5"), col = c("steelblue", "firebrick"),
lwd = 2, bty = "n")
new_days <- data.frame(days = c(8, 12, 20))
data.frame(days = new_days$days,
quadratic = round(predict(quad_bod, new_days)),
degree_5 = round(predict(poly5_bod, new_days))) days quadratic degree_5
1 8 219 177
2 12 218 1577
3 20 83 121395The degree-5 polynomial is “perfect” on the data and absurd everywhere else. It predicts BOD dropping to 177 on day 8 (between measurements of 213 and 224) and over 100,000 on day 20. It has memorized six numbers rather than learning a trend. The quadratic misses the points a bit, but behaves sensibly inside the data. Even it turns downward by day 20, though, which the chemistry says shouldn’t happen: BOD rises and levels off as the food for the bacteria runs out.
Which model would you report? What shape of function would match “rises, then levels off”? (Hint: think back to constrained growth.)