Code
library(rms)
library(ggplot2)
library(splines)
library(lspline)
library(car)
library(kableExtra)
library(plotrix)
tryCatch(source('pander_registry.R'), error = function(e) invisible(e))ANOVA and Model Choice
Lecture 12
April 6, 2026
1 Overview
There are many different ways we can model predictors
Consider what alternative models will make scientific sense
What is the impact of letting the data drive the selection of a model
I am going to discuss in terms of a clinical trial where we have replicates at dose levels
- We can find dose-specific means and compare modeling approaches
2 ANOVA versus Linear Continuous Models
Compare power of linear continuous models versus ANOVA as a function
of trend in means AND
standard errors within groups
ANOVA (dummy variables)
Uses indicator variables for every dose (group) level
Again, I am thinking about “dose” in a general sense that could include covariates like age, cholesterol, blood pressure, etc.
Traditionally, dose would just be dose of some treatment
Fits group means exactly (saturated model)
One way ANOVA: One categorical predictor
Two way ANOVA: Two categorical predictors
- Fit with the interactions to get group means exactly
Saturated models do not mix random error with systematic error
Systematic error: Error due to differences from sample means from predicted means
Random error: Error that cannot be explained after controlling for dose
ANOVA ignores the ordering of the groups, so it gains no power from trends
e.g. does not assume that the difference between dose=15 and dose=30 group is similar to the difference between the dose=30 and dose=45 groups
In fact, the same level of significance is gained no matter what permutation of dose groups is used
Linear continuous models
Borrows information across groups
- Accurate and efficient if the model is correct
If model is incorrect, mixes random and systematic error
- Will have some systematic error because the means are not predicted exactly
Can gain power from ordering of groups in order to detect a trend
- But, no matter how low the standard error is, if there is no trend in the mean, there is no statistical significance
Code
n <- 10
dose.levels <- c(0,20,40,60,80)
dose <- rep(dose.levels, n)
set.seed(37)
y1 <- 300 + 5*dose + rnorm(5*n, 0, 200)
y2 <- 300 + 200*(dose==20) + 100*(dose==40) - 100*(dose==60) + 0*(dose==80) + rnorm(5*n, 0, 200)
m1 <- lm(y1 ~ dose)
m2 <- lm(y1 ~ factor(dose))
m4 <- lm(y2 ~ factor(dose))
p1 <- predict(m1, newdata=data.frame(dose=dose.levels), se.fit=TRUE)
p2 <- predict(m2, newdata=data.frame(dose=dose.levels), se.fit=TRUE)
p4 <- predict(m4, newdata=data.frame(dose=dose.levels), se.fit=TRUE)Code
par(mfrow=c(2,2))
plotCI(x=dose.levels, y=p1$fit, p1$se.fit, ylab="Response", xlab="Dose", ylim=c(0,1000), main="Linear: High power; ANOVA: High Power")
plotCI(x=dose.levels, y=p1$fit, 3*p1$se.fit, ylab="Response", xlab="Dose", ylim=c(0,1000), main="Linear: Mod power; ANOVA: Low Power")
plotCI(x=dose.levels[c(1,5,4,2,3)], y=p1$fit, p1$se.fit, ylab="Response", xlab="Dose", ylim=c(0,1000), main="Linear: No power; ANOVA: High Power")
plotCI(x=dose.levels[c(1,5,4,2,3)], y=p1$fit, 3*p1$se.fit, ylab="Response", xlab="Dose", ylim=c(0,1000), main="Linear: No power; ANOVA: Low Power")
Other options for modeling continuous predictors
Combinations of linear trends and indicator variables
Splines
Polynomials
etc.
3 Choice of Transformation
The exact form used to model predictors should be based on scientific (first) and statistical (second) criteria
Scientific issues
The form used to model predictors must address the specific scientific question
Should be the next logical step in the process of investigating the overall goal
First, establish some sort of an association
Second, detect a first order trend
Third, detecting specific forms of non-linearities
Threshold effects?
U- or S-shaped trends?
Finally, more complex models
When the scientific question relates to prediction, it is imperative that the regression model accurately reflects the true relationship between predictors and the summary measure of response
- Failure to have the correct model will guarantee that some groups may not have the correct predicted response
When the scientific question relates to detection of associations, the importance of having the true model depends on the statistical role of the predictor
With the predictor of interest, the most important issues is to protect the validity of the statistical inference
- Data driven decision will inflate the type I error rate
With precision variables, it is not as crucial that the true relationship be modeled
- An approximate model will provide most of the precision gains
With confounders, failure to accurately model the relationship between the confounder and the response may lead to residual confounding
- Sometimes we will use very flexible models for continuous confounders as there is little cost to doing so and the potential for important gain
As the goal of any analysis is to communicate findings to the greater scientific community, it is also important that modeling of predictors is easy to understand
This is an issue that matters most for your predictor of interest
We are generally not worried about making inference about precision variables or confounders
Statistical issues
The greatest statistical precision will be gained when the model reflects the true relationship between the predictor and the response
Accurate modeling of the relationship will avoid introducing systematic error in the estimates of the standard errors
Parsimony: Using the fewest parameters to model the relationship will allow greater precision
Precision is a trade-off between parsimony and increased accuracy from including more parameters
We should select the form of modeling the predictor before looking at the data
Data driven selection of transformations will tend to lead to inaccurate (anti-conservative) statistical inference
Overfitting of the data leads to spuriously low estimates of the within group variability
Thus standard errors estimates are too low
Type-I errors are inflated
Confidence interval are too narrow (inaccurate coverage probabilities)
Data-driven model selection will also lead to coefficient estimates that are biased away from the null (leading you to overstate your scientific effects)
4 Example: Beta Carotene Supplements
4.1 Overview
Before doing large scale clinical trials, it is important to understand the pharmacokinetics of a drug
Phase II prevention trials often administer a drug in various doses to volunteers, and pertinent plasma levels are then measured at regular intervals
Of particular interest is how dose level affects the build up of drug in the plasma over time, as well as how the dose level might affect other blood chemistry
Forty-six (46) volunteers were randomly assigned to receive one of five doses of beta-carotene (0, 15, 30, 45, or 60 mg/day) for 9 months in a double blind fashion
The specific aim was to determine how different dose levels affected the serum beta-carotene levels after 9 months
Other measured variables available in this data set include subject age, sex, weight, body mass index, percent body fat, and serum cholesterol level at baseline
Code
carot <- read.csv("data/carot-nomissing.csv")
res1 <- aggregate(carot3 ~ dose, data = carot,
FUN = function(x) c(n = length(x), mean = mean(x), sd = sd(x),
min = min(x), q25 = quantile(x, 0.25),
median = median(x), q75 = quantile(x, 0.75),
max = max(x)))
res1_df <- data.frame(Dose = res1$dose, res1$carot3, check.names = FALSE)
kable(res1_df, digits=1,
col.names = c("Dose", "N", "Mean", "SD", "Min", "Q25", "Median", "Q75", "Max"),
caption = "Summary of 9-month plasma beta-carotene by dose group") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Dose | N | Mean | SD | Min | Q25 | Median | Q75 | Max |
|---|---|---|---|---|---|---|---|---|
| 0 | 8 | 177.6 | 84.9 | 84.5 | 123.6 | 144.5 | 218.3 | 323.0 |
| 15 | 10 | 1167.1 | 537.1 | 576.8 | 729.6 | 1056.7 | 1521.8 | 2018.8 |
| 30 | 10 | 1536.7 | 462.9 | 849.3 | 1161.2 | 1528.9 | 1836.5 | 2248.5 |
| 45 | 8 | 1732.0 | 538.3 | 950.2 | 1457.3 | 1760.6 | 2227.7 | 2310.4 |
| 60 | 10 | 1695.8 | 703.4 | 59.7 | 1642.9 | 1836.3 | 1916.2 | 2855.0 |
In this randomized trial, we can consider several potential response variables
Plasma level at the end of treatment
Change in plasma level over the treatment period
Either of the above adjusted for baseline plasma (ANCOVA model)
Accounting for baseline
Dose group \(i\), subject \(j\), time \(t\)
\(Y_{ijt} \sim (\mu_{it}, \sigma^2)\); \(\textrm{corr}(Y_{ij0}, Y_{ij9}) = \rho\) \[\begin{aligned} \overline{Y}_{i\cdot9} & \sim & \left(\mu_{i9}, \sigma^2/n \right) \\ \overline{Y}_{i\cdot9} - \overline{Y}_{i\cdot0} & \sim & \left(\mu_{i9} - \mu_{i0}, 2\sigma^2(1-\rho)/n \right) \\ \overline{Y}_{i\cdot9} - \rho \overline{Y}_{i\cdot0} & \sim & \left(\mu_{i9} - \rho \mu_{i0}, \sigma^2(1-\rho^2)/n \right)\end{aligned}\]
Compared variances of the above three equations
When are the variances equal, smaller, larger
Which is always smallest
By randomization, there will be equal means at baseline
- \(\mu_{T,0} = \mu_{P,0}\) where \(T\) is any of the treatment doses and \(P\) is placebo
Contrast across dose groups \[\begin{aligned} \overline{Y}_{T,\cdot9} - \overline{Y}_{P,\cdot9} & \sim & \left(\mu_{T,9} - \mu_{P,9}, 2\sigma^2/n \right) \\ \left(\overline{Y}_{T,\cdot9} - \overline{Y}_{T, \cdot0}\right) - \left(\overline{Y}_{P,\cdot9} - \overline{Y}_{P, \cdot0}\right) & \sim & \left(\mu_{T,9} - \mu_{P,9}, 4\sigma^2(1-\rho)/n \right) \\ \left(\overline{Y}_{T,\cdot9} - \rho \overline{Y}_{T,\cdot0}\right) - \left(\overline{Y}_{P,\cdot9} - \rho \overline{Y}_{P,\cdot0}\right) & \sim & \left(\mu_{T,9} - \mu_{P,9}, 2\sigma^2(1-\rho^2)/n \right)\end{aligned}\]
Simple linear regression
Regress \(Y\) on \(X\)
\(Y_i \sim \left(\mu_Y, \sigma^2_Y\right)\); \(X_i \sim \left(\mu_X, \sigma^2_X\right)\)
\(\textrm{corr}\left(Y_i, X_i \right) = \rho\)
Regression model: \(E[Y_i | X_i] = \beta_0 + \beta_1 X_i\)
\(\beta_0 = \mu_Y - \beta_1 \mu_x\)
\(\beta_1 = \rho \frac{\sigma_Y}{\sigma_X}\)
Analysis of Covariance
Dose group \(i\), subject \(j\), time \(t\)
\(Y_{ijt} \sim (\mu_{it}, \sigma^2)\); \(\textrm{corr}(Y_{ij0}, Y_{ij9}) = \rho\)
Regression model: \(E[Y_{ij9} | Y_{i_j0}] = \beta_0 + \beta_1 Y_{ij0}\)
- \(\beta_1 = \rho\)
4.2 Methods for modeling dose response
In a randomized clinical trial, we will tend to have the greatest precision if we adjust for baseline as a predictor in a linear regression model
A wide variety of models may be considered for examining the relationship between dose and plasma levels
Dummy variables where we model each dose level independently, without borrowing information across groups (ANOVA)
Linear continuous predictors (transformed or untransformed)
Dichotomization (at any of several thresholds)
Polynomials, splines, other flexible methods
Combinations of the above
Even more complex models
I will compare possible models
Graphically: Show data and fitted values without adjustment for baseline
Numerically: Show regression estimates and tests after adjustment for baseline
Note that this is an academic exercise and not something you would do in practice to come up with the “best” model
Predicted values
After fitting a regression model, the
predictfunction will provide predicted values for each case- Mathematically, this is just the intercept plus the regression parameters multiplied by the covariates for each case
4.3 ANOVA analysis
Fits each group independently
Does not use the ordering of the dose groups when looking for an effect
- Completely ignores the magnitude and ordering of the \(x\)-axis
A priori, we might expect this is not the most efficient method if the alternative hypothesis is true
We expect larger plasma levels with increasing dose
We will thus have less power to detect a first-order trend
Code
carot$dose.factor <- factor(carot$dose)
m.factor <- lm(carot3 ~ dose.factor + carot0, data=carot)
summary(m.factor)Fitting linear model: carot3 ~ dose.factor + carot0
| Min | 1Q | Median | 3Q | Max |
|---|---|---|---|---|
| -1302.6 | -243.47 | -11.031 | 227.9 | 982.51 |
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | -389.0119 | 235.38396 | -1.6527 | 1.0623e-01 |
| dose.factor15 | 1094.7307 | 222.46410 | 4.9209 | 1.5228e-05 |
| dose.factor30 | 1465.8224 | 222.52693 | 6.5872 | 7.0800e-08 |
| dose.factor45 | 1645.0543 | 233.72692 | 7.0384 | 1.6600e-08 |
| dose.factor60 | 1628.1655 | 222.66634 | 7.3121 | 6.9000e-09 |
| carot0 | 2.0967 | 0.62423 | 3.3589 | 1.7282e-03 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 46 | 464.32 | 0.65401 | 0.61076 |
F-statistic: 15.1219 on 5 and 40 DF, pvalue: 0
Code
# Helper objects used across all dose-response plots
m1 <- lm(carot3 ~ factor(dose), data=carot)
dose.seq <- seq(0, 60, by=15)
group.means <- predict(m1, newdata=data.frame(dose=dose.seq))
plot(jitter(carot$dose, factor=.2), carot$carot3, ylab="Plasma Beta Carotene", xlab="Dose of Supplementation", main="Model: Dummy Variables (ANOVA)")
lines(dose.seq, group.means)
points(dose.seq, group.means, pch="O", cex=2)
points(dose.seq, group.means, pch="X", cex=2)
legend("topleft", c("Fitted", "Group Means"), pch=c("X","O"), bty="n")
Testing for the dose effect
- We must use
linearHypothesisto jointly test the dose coefficients because the model includes the baseline measurement
- We must use
Code
linearHypothesis(m.factor, c("dose.factor15",
"dose.factor30",
"dose.factor45",
"dose.factor60"))| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 44 | 23995659 | ||||
| 40 | 8623858 | 4 | 15371801 | 17.825 | 1.79e-08 |
We would have had the same fitted values (and thus inference) if we had decided to drop a different dose group
- Example: Making my own dummy variables for dose, with dose at 60 being the reference group
Code
carot$dose.new <- relevel(carot$dose.factor, ref="60")
m2.factor <- lm(carot3 ~ dose.new + carot0, data=carot)
summary(m2.factor)Fitting linear model: carot3 ~ dose.new + carot0
| Min | 1Q | Median | 3Q | Max |
|---|---|---|---|---|
| -1302.6 | -243.47 | -11.031 | 227.9 | 982.51 |
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 1239.1536 | 200.11296 | 6.192271 | 2.5360e-07 |
| dose.new0 | -1628.1655 | 222.66634 | -7.312131 | 6.9000e-09 |
| dose.new15 | -533.4348 | 207.65661 | -2.568831 | 1.4041e-02 |
| dose.new30 | -162.3431 | 207.65408 | -0.781796 | 4.3894e-01 |
| dose.new45 | 16.8887 | 220.32243 | 0.076655 | 9.3928e-01 |
| carot0 | 2.0967 | 0.62423 | 3.358942 | 1.7282e-03 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 46 | 464.32 | 0.65401 | 0.61076 |
F-statistic: 15.1219 on 5 and 40 DF, pvalue: 0
Code
linearHypothesis(m2.factor, c("dose.new0",
"dose.new15",
"dose.new30",
"dose.new45"))| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 44 | 23995659 | ||||
| 40 | 8623858 | 4 | 15371801 | 17.825 | 1.79e-08 |
Note that the parameter estimates all will lead to the same fitted values
- e.g. Intercept in the re-leveled model equals the intercept + dose60 coefficient in the original model
Overall F statistics, R-squared, Root MSE all the same
Partial t-tests tend to differ as we are making comparisons to different reference groups
Could also fit the same model with no intercept
Would then have to include all five dose groups
In R, fit a model without an intercept by adding a \(-1\) in the model equation (e.g. \(y \sim -1 + x\))
Not including the intercept changes the overall F statistic and the R-squared measures
Code
m3.factor <- lm(carot3 ~ -1 + dose.factor + carot0, data=carot)
summary(m3.factor)Fitting linear model: carot3 ~ -1 + dose.factor + carot0
| Min | 1Q | Median | 3Q | Max |
|---|---|---|---|---|
| -1302.6 | -243.47 | -11.031 | 227.9 | 982.51 |
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| dose.factor0 | -389.0119 | 235.38396 | -1.6527 | 1.0623e-01 |
| dose.factor15 | 705.7188 | 201.07057 | 3.5098 | 1.1256e-03 |
| dose.factor30 | 1076.8105 | 200.76804 | 5.3635 | 3.6956e-06 |
| dose.factor45 | 1256.0424 | 216.85144 | 5.7922 | 9.2610e-07 |
| dose.factor60 | 1239.1536 | 200.11296 | 6.1923 | 2.5360e-07 |
| carot0 | 2.0967 | 0.62423 | 3.3589 | 1.7282e-03 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 46 | 464.32 | 0.91487 | 0.9021 |
F-statistic: 71.6448 on 6 and 40 DF, pvalue: < 2.22e-16
Correspondence of the no-intercept model compared to previous models
Some textbooks refer to this as a “cell means” coding system
If we didn’t have baseline beta carotene in the model, the dose parameters would correspond directly to the means in each dose group
With baseline beta carotene in the model, the dose parameters are the means when carot0 is 0
In terms of model fit, the model is the same as before
- No intercept means each dose group is compared to a mean of 0
Fitted values will be the same
Test of dose effect will need to test equality of all five dose covariates
This is not a test that these 5 parameters are 0
\(H_0: dose0 = dose15 = dose30 = dose45 = dose60\)
\(H_1:\) at least one of the above is not equal
Code
linearHypothesis(m3.factor, c("dose.factor0=dose.factor15",
"dose.factor0=dose.factor30",
"dose.factor0=dose.factor45",
"dose.factor0=dose.factor60"))| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 44 | 23995659 | ||||
| 40 | 8623858 | 4 | 15371801 | 17.825 | 1.79e-08 |
4.4 Binary dose: Placebo versus Active
Dichotomize into dose 0 versus dose \(>\) 0
Will be an accurate model if all (or virtually all) of the effect is attained at the lowest dose level
Often used when little is know about a treatment, or when dose is difficult to quantify
e.g. Smoking
We are relatively certain of a smoking effect, so our major scientific interest is likely related to the dose-response relationship above the lowest dose
Code
carot$trt <- (carot$dose>0)+0
m.trt <- lm(carot3 ~ trt + carot0, data=carot)
summary(m.trt)Fitting linear model: carot3 ~ trt + carot0
| Min | 1Q | Median | 3Q | Max |
|---|---|---|---|---|
| -1121.1 | -251.97 | 56.167 | 312.41 | 938.71 |
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | -392.2859 | 250.57735 | -1.5655 | 1.2479e-01 |
| trt | 1449.2213 | 195.10662 | 7.4278 | 3.1000e-09 |
| carot0 | 2.1089 | 0.66439 | 3.1741 | 2.7764e-03 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 46 | 494.39 | 0.57833 | 0.55872 |
F-statistic: 29.4876 on 2 and 43 DF, pvalue: 0
Code
m2 <- lm(carot3 ~ dose>0, data=carot)
plot(jitter(carot$dose, factor=.2), carot$carot3, ylab="Plasma Beta Carotene", xlab="Dose of Supplementation", main="Model: Dichotomous Dose")
lines(seq(0,60, by=15), predict(m1, newdata=data.frame(dose=seq(0,60, by=15))), pch="O", cex=2)
points(seq(0,60, by=15), predict(m1, newdata=data.frame(dose=seq(0,60, by=15))), pch="O", cex=2)
points(seq(0,60, by=15), predict(m2, newdata=data.frame(dose=seq(0,60, by=15))), pch="X", cex=2)
legend("topleft", c("Fitted", "Group Means"), pch=c("X","O"), lty=c(NA,1), bty="n")
4.5 Linear, continuous dose
Estimates the best fitting straight line to response
- Accurate if the response is linear
Often used when little is know about the treatment and a general trend is expected
- In this particular application, we are relatively certain of an effect, so our major interest is in modeling the dose response relationship above 0.
Code
m.cont <- lm(carot3 ~ dose + carot0, data=carot)
summary(m.cont)Fitting linear model: carot3 ~ dose + carot0
| Min | 1Q | Median | 3Q | Max |
|---|---|---|---|---|
| -1623.4 | -348.21 | -57.623 | 380.53 | 1220.7 |
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 148.7970 | 228.49026 | 0.65122 | 5.1837e-01 |
| dose | 23.8235 | 3.83976 | 6.20443 | 1.8430e-07 |
| carot0 | 1.7856 | 0.72344 | 2.46817 | 1.7635e-02 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 46 | 542.63 | 0.49203 | 0.46841 |
F-statistic: 20.8256 on 2 and 43 DF, pvalue: 0
- To test the treatment effect, we can use the partial t-test from the output directly as we are only testing one parameter
Code
m3 <- lm(carot3 ~ dose, data=carot)
plot(jitter(carot$dose, factor=.2), carot$carot3, ylab="Plasma Beta Carotene", xlab="Dose of Supplementation", main="Model: Linear Dose")
lines(seq(0,60, by=15), predict(m1, newdata=data.frame(dose=seq(0,60, by=15))), pch="O", cex=2)
points(seq(0,60, by=15), predict(m1, newdata=data.frame(dose=seq(0,60, by=15))), pch="O", cex=2)
points(seq(0,60, by=15), predict(m3, newdata=data.frame(dose=seq(0,60, by=15))), pch="X", cex=2)
legend("topleft", c("Fitted", "Group Means"), pch=c("X","O"), lty=c(NA,1), bty="n")
4.6 Polynomial models of dose
Fit terms involving dose, dose squared
Often used to fit U-shaped trends
In general, a quadratic is a pretty strong assumption in that it assumes constant curvature over dose
Code
m.poly <- lm(carot3 ~ poly(dose,2) + carot0, data=carot)
summary(m.poly)Fitting linear model: carot3 ~ poly(dose, 2) + carot0
| Min | 1Q | Median | 3Q | Max |
|---|---|---|---|---|
| -1272.6 | -259.09 | -8.8187 | 173.8 | 1105.8 |
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 821.367 | 157.05012 | 5.2300 | 5.0306e-06 |
| poly(dose, 2)1 | 3410.018 | 463.88612 | 7.3510 | 4.6000e-09 |
| poly(dose, 2)2 | -1943.138 | 463.05720 | -4.1963 | 1.3744e-04 |
| carot0 | 2.037 | 0.61736 | 3.2996 | 1.9795e-03 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 46 | 460.87 | 0.64209 | 0.61653 |
F-statistic: 25.1162 on 3 and 42 DF, pvalue: 0
Code
m4 <- lm(carot3 ~ poly(dose, 2), data=carot)
plot(jitter(carot$dose, factor=.2), carot$carot3, ylab="Plasma Beta Carotene", xlab="Dose of Supplementation", main="Model: Quadratic Dose")
lines(seq(0,60, by=15), predict(m1, newdata=data.frame(dose=seq(0,60, by=15))), pch="O", cex=2)
points(seq(0,60, by=15), predict(m1, newdata=data.frame(dose=seq(0,60, by=15))), pch="O", cex=2)
points(seq(0,60, by=15), predict(m4, newdata=data.frame(dose=seq(0,60, by=15))), pch="X", cex=2)
legend("topleft", c("Fitted", "Group Means"), pch=c("X","O"), lty=c(NA,1), bty="n")
The partial t-test for the quadratic term can be interpreted as a test for departure from a linear dose response
- It is highly significant, suggesting the relationship is not purely linear
To test the treatment effect, we need to test the two dose covariates
Code
linearHypothesis(m.poly, c("poly(dose, 2)1",
"poly(dose, 2)2"))| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 44 | 23995659 | ||||
| 42 | 8920871 | 2 | 15074787 | 35.487 | 9e-10 |
4.7 Highest order polynomial models versus ANOVA
With 5 discrete dose levels, a 4th degree polynomial will fit the means exactly
Thus, the model will have the same fit as the ANOVA model using dummy variables for each levels of dose
Higher order polynomials are borrowing less information across dose groups
Highest order polynomial borrows no information across dose groups
Code
m.poly4 <- lm(carot3 ~ poly(dose, 4, raw=TRUE) + carot0, data=carot)
summary(m.poly4)Fitting linear model: carot3 ~ poly(dose, 4, raw = TRUE) + carot0
| Min | 1Q | Median | 3Q | Max |
|---|---|---|---|---|
| -1302.6 | -243.47 | -11.031 | 227.9 | 982.51 |
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | -3.8901e+02 | 2.3538e+02 | -1.65267 | 0.1062273 |
| poly(dose, 4, raw = TRUE)1 | 1.1785e+02 | 5.6117e+01 | 2.10015 | 0.0420684 |
| poly(dose, 4, raw = TRUE)2 | -3.8818e+00 | 4.6633e+00 | -0.83240 | 0.4101273 |
| poly(dose, 4, raw = TRUE)3 | 6.5967e-02 | 1.2483e-01 | 0.52845 | 0.6001084 |
| poly(dose, 4, raw = TRUE)4 | -4.4119e-04 | 1.0414e-03 | -0.42366 | 0.6740853 |
| carot0 | 2.0967e+00 | 6.2423e-01 | 3.35894 | 0.0017282 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 46 | 464.32 | 0.65401 | 0.61076 |
F-statistic: 15.1219 on 5 and 40 DF, pvalue: 0
- Compare to the ANOVA model — same fit statistics (R-squared, Root MSE, F-statistic)
Code
summary(m.factor)Fitting linear model: carot3 ~ dose.factor + carot0
| Min | 1Q | Median | 3Q | Max |
|---|---|---|---|---|
| -1302.6 | -243.47 | -11.031 | 227.9 | 982.51 |
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | -389.0119 | 235.38396 | -1.6527 | 1.0623e-01 |
| dose.factor15 | 1094.7307 | 222.46410 | 4.9209 | 1.5228e-05 |
| dose.factor30 | 1465.8224 | 222.52693 | 6.5872 | 7.0800e-08 |
| dose.factor45 | 1645.0543 | 233.72692 | 7.0384 | 1.6600e-08 |
| dose.factor60 | 1628.1655 | 222.66634 | 7.3121 | 6.9000e-09 |
| carot0 | 2.0967 | 0.62423 | 3.3589 | 1.7282e-03 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 46 | 464.32 | 0.65401 | 0.61076 |
F-statistic: 15.1219 on 5 and 40 DF, pvalue: 0
4.8 Threshold at 0 and Linear Term
Threshold at 0 and linear dose
To fit, use a dummy variable for dose0 plus dose (continuous)
Fits dose 0 by its group mean
Fits dose \(>\) 0 by a line (an intercept and slope)
Allows us to address two scientific questions
Is there any effect of dose? (test both slopes)
Is there any additional benefit beyond the lowest dose? (test linear term’s slope)
Code
m.thres <- lm(carot3 ~ trt + dose + carot0, data=carot)
summary(m.thres)Fitting linear model: carot3 ~ trt + dose + carot0
| Min | 1Q | Median | 3Q | Max |
|---|---|---|---|---|
| -1389.4 | -282.23 | -8.3425 | 298.72 | 889.03 |
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | -392.4465 | 234.35829 | -1.6746 | 0.10145283 |
| trt | 1014.1871 | 244.41557 | 4.1494 | 0.00015893 |
| dose | 11.7251 | 4.38259 | 2.6754 | 0.01059449 |
| carot0 | 2.1094 | 0.62139 | 3.3947 | 0.00151081 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 46 | 462.39 | 0.63973 | 0.61399 |
F-statistic: 24.8594 on 3 and 42 DF, pvalue: 0
Code
m5 <- lm(carot3 ~ (dose >0) + dose, data=carot)
plot(jitter(carot$dose, factor=.2), carot$carot3, ylab="Plasma Beta Carotene", xlab="Dose of Supplementation", main="Model: Threshold and Linear Dose")
lines(seq(0,60, by=15), predict(m1, newdata=data.frame(dose=seq(0,60, by=15))), pch="O", cex=2)
points(seq(0,60, by=15), predict(m1, newdata=data.frame(dose=seq(0,60, by=15))), pch="O", cex=2)
points(seq(0,60, by=15), predict(m5, newdata=data.frame(dose=seq(0,60, by=15))), pch="X", cex=2)
legend("topleft", c("Fitted", "Group Means"), pch=c("X","O"), lty=c(NA,1), bty="n")
Testing the effect of treatment
Two variables model dose, so we need to test both
If response increases from dose 0 to lowest dose OR
… response increases as dose increase, THEN
… we will declare an effect of treatment
The partial t-test for the
trtterm tests whether there is a jump from placebo to the active dose groups (the threshold effect)- Here, it is highly significantly different from 0, indicating a large initial effect of treatment
The partial t-test for the
doseterm tests whether there is an additional linear trend among the active dose groupsIt is significantly different from 0 (\(p = 0.01\))
There is a multiple comparison issue here, but many people are comfortable doing this ‘step down’ test after they have already tested from any treatment effect
Code
linearHypothesis(m.thres, c("trt", "dose"))| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 44 | 23995659 | ||||
| 42 | 8979818 | 2 | 15015841 | 35.116 | 1.1e-09 |
4.9 Summary comparison of models
- Comparing all models side by side helps illustrate the trade-off between parsimony and flexibility
Code
dose.seq <- seq(0, 60, by=15)
group.means <- predict(m1, newdata=data.frame(dose=dose.seq))
plot(jitter(carot$dose, factor=.2), carot$carot3,
ylab="Plasma Beta Carotene", xlab="Dose of Supplementation",
main="All Models Compared", col="gray70", pch=16)
points(dose.seq, group.means, pch=16, cex=2)
# Fitted values from each unadjusted model
m2.plot <- lm(carot3 ~ dose, data=carot)
m3.plot <- lm(carot3 ~ poly(dose, 2), data=carot)
m4.plot <- lm(carot3 ~ (dose > 0), data=carot)
lines(dose.seq, predict(m2.plot, newdata=data.frame(dose=dose.seq)), col="blue", lwd=2)
lines(dose.seq, predict(m3.plot, newdata=data.frame(dose=dose.seq)), col="red", lwd=2)
lines(dose.seq, predict(m4.plot, newdata=data.frame(dose=dose.seq)), col="darkgreen", lwd=2, lty=2)
lines(dose.seq, predict(m5, newdata=data.frame(dose=dose.seq)), col="purple", lwd=2, lty=3)
lines(dose.seq, group.means, col="black", lwd=2, lty=1)
legend("topleft",
c("Group Means (ANOVA)", "Linear", "Quadratic", "Dichotomized", "Threshold + Linear"),
col=c("black", "blue", "red", "darkgreen", "purple"),
lwd=2, lty=c(1,1,1,2,3), bty="n", cex=0.85)
Code
# Model comparison table (all adjusted for baseline)
model_names <- c("ANOVA (dummy variables)", "Dichotomized", "Linear", "Quadratic", "Threshold + Linear")
model_list <- list(m.factor, m.trt, m.cont, m.poly, m.thres)
model_df <- min(sapply(model_list, function(m) c("df" = m$df.residual)))
comp_table <- data.frame(
Model = model_names,
`Dose df` = c(4, 1, 1, 2, 2),
`Residual df` = sapply(model_list, function(m) m$df.residual),
`R-squared` = sapply(model_list, function(m) summary(m)$r.squared),
`Root MSE` = sapply(model_list, function(m) summary(m)$sigma),
AIC = sapply(model_list, AIC),
check.names = FALSE
)
kable(comp_table, digits=c(0, 0, 0, 3, 1, 1),
caption = "Comparison of dose-response models (all adjusted for baseline)") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Model | Dose df | Residual df | R-squared | Root MSE | AIC |
|---|---|---|---|---|---|
| ANOVA (dummy variables) | 4 | 40 | 0.654 | 464.3 | 703.0 |
| Dichotomized | 1 | 43 | 0.578 | 494.4 | 706.1 |
| Linear | 1 | 43 | 0.492 | 542.6 | 714.7 |
| Quadratic | 2 | 42 | 0.642 | 460.9 | 700.6 |
| Threshold + Linear | 2 | 42 | 0.640 | 462.4 | 700.9 |
5 Data driven model selection
Suppose we look at a scatterplot before deciding which model we fit and choose a model that can fit the data well
If the data looks like a straight line, choose the model linear in dose
If the data looks like a U, choose a quadratic
If the data is a complicated pattern of differences among groups, we might choose dummy variables or splines
etc.
This approach would tend to mimic the behavior of fitting several different models and choosing the model with the lowest \(p\)-value
- When our eye sees some trend in the data, we would be most likely to pick the model giving the lowest \(p\)-value
5.1 Simulation
Using the 46 subjects in this dataset, I can randomly permute the dose they received
Effectively, randomize subjects to a different dose
But, keep their 9-month and baseline beta carotene levels the same (not permuted)
- Should remove any association between dose and beta carotene
Next, fit each of the five models (linear, quadratic, ANOVA, dichotomized, and dichotomized plus linear)
Repeat the process 5000 times (representing 5000 studies)
- Calculate how often each model rejects the null hypothesis of a dose effect
Code
set.seed(80)
reps <- 5000
n <- length(carot$carot3)
p.vals <- matrix(NA, nrow=reps, ncol=5)
colnames(p.vals) <- c("ANOVA", "Linear", "Quadratic", "Dichotomized", "Dichot + Linear")
y <- carot$carot3
y0 <- carot$carot0
m.restrict <- lm(y ~ y0)
for(i in 1:reps) {
perm <- sample(1:n)
# Permute the dose, but not the outcome and baseline
d <- carot$dose[perm]
m.anova <- lm(y ~ y0 + factor(d))
m.linear <- lm(y ~ y0 + d)
m.quad <- lm(y ~ y0 + poly(d,2))
m.dichot <- lm(y ~ y0 + (d>0))
m.dilin <- lm(y ~ y0 + d + (d>0))
p.vals[i, ] <- c(anova(m.anova)[2,"Pr(>F)"],
anova(m.linear)[2,"Pr(>F)"],
anova(m.quad)[2,"Pr(>F)"],
anova(m.dichot)[2,"Pr(>F)"],
anova(m.restrict, m.dilin)[2,"Pr(>F)"])
}Individual Model Results
- Empirical type I error for each method of analysis individually
Code
sim_results <- data.frame(
Model = colnames(p.vals),
`Emp. Type-I Error` = apply(p.vals < 0.05, 2, mean),
check.names = FALSE)
kable(sim_results, digits=3,
caption = "Empirical type I error rate by model (nominal = 0.05)") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Model | Emp. Type-I Error | |
|---|---|---|
| ANOVA | ANOVA | 0.058 |
| Linear | Linear | 0.047 |
| Quadratic | Quadratic | 0.048 |
| Dichotomized | Dichotomized | 0.047 |
| Dichot + Linear | Dichot + Linear | 0.048 |
Multiple comparison issues
With 5 hypothesis tests at a nominal 0.05 level, experiment-wise error rate is at most \(0.25\) (\(0.05 \times 5\))
Worst-case assumes that all tests are mutually exclusive
e.g. If the linear dose-response model is significant, no other model is more likely to be significant
In fact, the tests will be correlated
How many of the 5000 simulated trials had at least one model with a \(p\)-value \(< 0.05\)?
Code
p.hat <- mean(apply(p.vals < 0.05, 1, any))
ci <- p.hat + c(-1.96, 1.96) * sqrt(p.hat * (1 - p.hat) / reps) - From the simulation, this was 641 or 12.8%
- Note that there is error in this estimate (due to the simulation randomness)
- 95% CI: [11.9%, 13.7%]General statistical issues
The true type 1 error rate for such data driven analyses will depend on several factors
The number of tests performed
The models considered
- Similar models will tend to reject the null hypotheses on the same dataset
The distribution of the data
- In particular, heavy tailed distributions decreases the concordance between the tests
When you have multiple models you are considering, the conclusions are less strong
The p-values (or other metrics) can still be useful in ordering the associations
Among all of the models considered, it appears as if SNP X is the most strongly associated with CVD
- Would be useful to put a CI around this ranking as well
5.2 Post hoc adjustments for multiple comparisons
In frequentist reasoning, we try to ensure that our error rate is held at some level \(\alpha\)
When only considering one decision, this is relatively easy
When making multiple decisions, we must consider the experiment-wise error rate
In the worst case scenario, an error rate of \(\alpha\) on each decision could lead to an experiment-wise error rate that is as high as \(k \times \alpha\)
- Such would be the case if all of our errors were mutually exclusive
If all error were independent of each other, then the experiment-wise error rate is
- \(1 - (1 - \alpha)^k\)
Experiment-wise error rates (\(\alpha = 0.05\) at each decision)
Code
alpha <- 0.05
k <- c(1, 2, 3, 5, 10, 20, 50)
ewer_table <- data.frame(
`Number of Comparisons` = k,
`Worst Case (Bonferroni bound)` = pmin(k * alpha, 1),
`Independent Errors` = round(1 - (1 - alpha)^k, 4),
check.names = FALSE
)
kable(ewer_table, digits=4,
caption = "Experiment-wise error rates at nominal alpha = 0.05") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Number of Comparisons | Worst Case (Bonferroni bound) | Independent Errors |
|---|---|---|
| 1 | 0.05 | 0.0500 |
| 2 | 0.10 | 0.0975 |
| 3 | 0.15 | 0.1426 |
| 5 | 0.25 | 0.2262 |
| 10 | 0.50 | 0.4013 |
| 20 | 1.00 | 0.6415 |
| 50 | 1.00 | 0.9231 |
When making multiple comparison which all tend address the same scientific question, we may adjust our level of significance to protect the experiment-wise error rate
- The problem with this approach is does not adjust for any bias in parameter estimates
Bonferroni Correction
Assumes the worst case scenario
When making \(k\) comparisons, either
Tests individual \(p\)-values against \(\frac{\alpha}{k}\)
Multiply \(p\)-values by \(k\) and compare to \(\alpha\) (keeping the \(p\)-values \(< 1\))
Bonferroni is easy and it can be applied in all settings
- Extremely conservative when the statistics from various tests are positively correlated
Many other varieties of adjusting after performing multiple comparisons
Tukey, Scheffe, etc.
None are great
Did they really adjust for all of the comparisons they made? Probably not.
My strong preference is to avoid multiple comparisons in the first place
If there was some model fitting involved to get to the final model, acknowledge that fact in the paper
Understand the science
Avoid data-driven approaches when you care about correct statistical inference (CIs and p-values)
6 Practical Guidelines
Pre-specify your model before looking at the data
Use scientific knowledge and prior studies to guide model choice
Document the analysis plan (ideally in a protocol or statistical analysis plan)
When in doubt, use flexible but pre-specified approaches
Restricted cubic splines with a pre-specified number of knots provide a good balance between flexibility and parsimony
The number of knots can be chosen based on the sample size, not the data pattern
If you must explore multiple models, be transparent
Report all models considered, not just the “winner”
Use information criteria (AIC, BIC) rather than p-values for model selection when the goal is prediction
Acknowledge the exploratory nature of the analysis
Remember the hierarchy of goals
- Address the scientific question
- Protect the validity of inference
- Maximize precision (power)
In Lecture 13, we will see how model choice decisions directly impact predictive accuracy; both for continuous outcomes (prediction intervals) and binary outcomes (classification, ROC curves)