Linear Regression Example
The stan_lm function, which has its own vignette, fits regularized linear models using a
novel means of specifying priors for the regression coefficients. Here
we focus using the stan_glm function, which can be used to
estimate linear models with independent priors on the regression
coefficients.
To illustrate the usage of stan_glm and some of the
post-processing functions in the rstanarm package we’ll
use a simple example from Chapter 3 of Gelman and Hill
(2007):
We shall fit a series of regressions predicting cognitive test scores of three- and four-year-old children given characteristics of their mothers, using data from a survey of adult American women and their children (a subsample from the National Longitudinal Survey of Youth).
Using two predictors – a binary indicator for whether the mother has
a high-school degree (mom_hs) and the mother’s score on an
IQ test (mom_iq) – we will fit four contending models. The
first two models will each use just one of the predictors, the third
will use both, and the fourth will also include a term for the
interaction of the two predictors.
For these models we’ll use the default weakly informative priors for
stan_glm, which are currently set to
normal(0,10) for the intercept and normal(0,5)
for the other regression coefficients. For an overview of the many other
available prior distributions see
help("prior", package = "rstanarm").
library(rstanarm)
data(kidiq)
post1 <- stan_glm(kid_score ~ mom_hs, data = kidiq,
family = gaussian(link = "identity"),
seed = 12345)
post2 <- update(post1, formula = . ~ mom_iq)
post3 <- update(post1, formula = . ~ mom_hs + mom_iq)
(post4 <- update(post1, formula = . ~ mom_hs * mom_iq))stan_glm
family: gaussian [identity]
formula: kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq
observations: 434
predictors: 4
------
Median MAD_SD
(Intercept) -10.1 13.7
mom_hs 49.5 15.3
mom_iq 1.0 0.1
mom_hs:mom_iq -0.5 0.2
Auxiliary parameter(s):
Median MAD_SD
sigma 18.0 0.6
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregFollowing Gelman and Hill’s example, we make some plots overlaying the estimated regression lines on the data.
base <- ggplot(kidiq, aes(x = mom_hs, y = kid_score)) +
geom_point(size = 1, position = position_jitter(height = 0.05, width = 0.1)) +
scale_x_continuous(breaks = c(0,1), labels = c("No HS", "HS"))
base + geom_abline(intercept = coef(post1)[1], slope = coef(post1)[2],
color = "skyblue4", size = 1)
There several ways we could add the uncertainty in our estimates to
the plot. One way is to also plot the estimated regression line at each
draw from the posterior distribution. To do this we can extract the
posterior draws from the fitted model object using the
as.matrix or as.data.frame methods:
draws <- as.data.frame(post1)
colnames(draws)[1:2] <- c("a", "b")
base +
geom_abline(data = draws, aes(intercept = a, slope = b),
color = "skyblue", size = 0.2, alpha = 0.25) +
geom_abline(intercept = coef(post1)[1], slope = coef(post1)[2],
color = "skyblue4", size = 1)
For the second model we can make the same plot but the x-axis will
show the continuous predictor mom_iq:
draws <- as.data.frame(as.matrix(post2))
colnames(draws)[1:2] <- c("a", "b")
ggplot(kidiq, aes(x = mom_iq, y = kid_score)) +
geom_point(size = 1) +
geom_abline(data = draws, aes(intercept = a, slope = b),
color = "skyblue", size = 0.2, alpha = 0.25) +
geom_abline(intercept = coef(post2)[1], slope = coef(post2)[2],
color = "skyblue4", size = 1)
For the third and fourth models, each of which uses both predictors,
we can plot the continuous mom_iq on the x-axis and use
color to indicate which points correspond to the different
subpopulations defined by mom_hs. We also now plot two
regression lines, one for each subpopulation:
reg0 <- function(x, ests) cbind(1, 0, x) %*% ests
reg1 <- function(x, ests) cbind(1, 1, x) %*% ests
args <- list(ests = coef(post3))
kidiq$clr <- factor(kidiq$mom_hs, labels = c("No HS", "HS"))
lgnd <- guide_legend(title = NULL)
base2 <- ggplot(kidiq, aes(x = mom_iq, fill = relevel(clr, ref = "HS"))) +
geom_point(aes(y = kid_score), shape = 21, stroke = .2, size = 1) +
guides(color = lgnd, fill = lgnd) +
theme(legend.position = "right")
base2 +
stat_function(fun = reg0, args = args, aes(color = "No HS"), size = 1.5) +
stat_function(fun = reg1, args = args, aes(color = "HS"), size = 1.5)
reg0 <- function(x, ests) cbind(1, 0, x, 0 * x) %*% ests
reg1 <- function(x, ests) cbind(1, 1, x, 1 * x) %*% ests
args <- list(ests = coef(post4))
base2 +
stat_function(fun = reg0, args = args, aes(color = "No HS"), size = 1.5) +
stat_function(fun = reg1, args = args, aes(color = "HS"), size = 1.5)
Model comparison
One way we can compare the four contending models is to use an
approximation to Leave-One-Out (LOO) cross-validation, which is
implemented by the loo function in the loo
package:
# Compare them with loo
loo1 <- loo(post1, cores = 1)
loo2 <- loo(post2, cores = 1)
loo3 <- loo(post3, cores = 1)
loo4 <- loo(post4, cores = 1)
(comp <- loo_compare(loo1, loo2, loo3, loo4)) elpd_diff se_diff
post4 0.0 0.0
post3 -3.5 2.7
post2 -6.2 4.1
post1 -42.4 8.7 In this case the fourth model is preferred as it has the highest
expected log predicted density (elpd_loo) or, equivalently,
the lowest value of the LOO Information Criterion (looic).
The fourth model is preferred by a lot over the first model
loo_compare(loo1, loo4) elpd_diff se_diff
post4 0.0 0.0
post1 -42.4 8.7 because the difference in elpd is so much larger than
the standard error. However, the preference of the fourth model over the
others isn’t as strong:
loo_compare(loo3, loo4) elpd_diff se_diff
post4 0.0 0.0
post3 -3.5 2.7 loo_compare(loo2, loo4) elpd_diff se_diff
post4 0.0 0.0
post2 -6.2 4.1 The posterior predictive distribution
The posterior predictive distribution is the distribution of the outcome implied by the model after using the observed data to update our beliefs about the unknown parameters. When simulating observations from the posterior predictive distribution we use the notation \(y^{\rm rep}\) (for replicate) when we use the same observations of \(X\) that were used to estimate the model parameters. When \(X\) contains new observations we use the notation \(\tilde{y}\) to refer to the posterior predictive simulations.
Simulating data from the posterior predictive distribution using the observed predictors is useful for checking the fit of the model. Drawing from the posterior predictive distribution at interesting values of the predictors also lets us visualize how a manipulation of a predictor affects (a function of) the outcome(s).
Graphical posterior predictive checks
The pp_check function generates a variety of plots
comparing the observed outcome \(y\) to
simulated datasets \(y^{\rm rep}\) from
the posterior predictive distribution using the same observations of the
predictors \(X\) as we used to fit the
model. He we show a few of the possible displays. The documentation at
help("pp_check.stanreg", package = "rstanarm") has details
on all of the pp_check options.
First we’ll look at a plot directly comparing the distributions of
\(y\) and \(y^{\rm rep}\). The following call to
pp_check will create a plot juxtaposing the histogram of
\(y\) and histograms of five \(y^{\rm rep}\) datasets:
pp_check(post4, plotfun = "hist", nreps = 5)
The idea is that if the model is a good fit to the data we should be able to generate data \(y^{\rm rep}\) from the posterior predictive distribution that looks a lot like the observed data \(y\). That is, given \(y\), the \(y^{\rm rep}\) we generate should be plausible.
Another useful plot we can make using pp_check shows the
distribution of a test quantity \(T(y^{\rm
rep})\) compared to \(T(y)\),
the value of the quantity in the observed data. When the argument
plotfun = "stat" is specified, pp_check will
simulate \(S\) datasets \(y_1^{\rm rep}, \dots, y_S^{\rm rep}\), each
containing \(N\) observations. Here
\(S\) is the size of the posterior
sample (the number of MCMC draws from the posterior distribution of the
model parameters) and \(N\) is the
length of \(y\). We can then check if
\(T(y)\) is consistent with the
distribution of \(\left(T(y_1^{\rm yep}),
\dots, T(y_S^{\rm yep})\right)\). In the plot below we see that
the mean of the observations is plausible when compared to the
distribution of the means of the \(S\)
\(y^{\rm rep}\) datasets:
pp_check(post4, plotfun = "stat", stat = "mean")
Using plotfun="stat_2d" we can also specify two test
quantities and look at a scatterplot:
pp_check(post4, plotfun = "stat_2d", stat = c("mean", "sd"))
Generating predictions
The posterior_predict function is used to generate
replicated data \(y^{\rm rep}\) or
predictions for future observations \(\tilde{y}\). Here we show how to use
posterior_predict to generate predictions of the outcome
kid_score for a range of different values of
mom_iq and for both subpopulations defined by
mom_hs.
IQ_SEQ <- seq(from = 75, to = 135, by = 5)
y_nohs <- posterior_predict(post4, newdata = data.frame(mom_hs = 0, mom_iq = IQ_SEQ))
y_hs <- posterior_predict(post4, newdata = data.frame(mom_hs = 1, mom_iq = IQ_SEQ))
dim(y_hs)[1] 4000 13We now have two matrices, y_nohs and y_hs.
Each matrix has as many columns as there are values of
IQ_SEQ and as many rows as the size of the posterior
sample. One way to show the predictors is to plot the predictions for
the two groups of kids side by side:
par(mfrow = c(1:2), mar = c(5,4,2,1))
boxplot(y_hs, axes = FALSE, outline = FALSE, ylim = c(10,170),
xlab = "Mom IQ", ylab = "Predicted Kid IQ", main = "Mom HS")
axis(1, at = 1:ncol(y_hs), labels = IQ_SEQ, las = 3)
axis(2, las = 1)
boxplot(y_nohs, outline = FALSE, col = "red", axes = FALSE, ylim = c(10,170),
xlab = "Mom IQ", ylab = NULL, main = "Mom No HS")
axis(1, at = 1:ncol(y_hs), labels = IQ_SEQ, las = 3)