This lab will introduce two ways of inspecting the "fit" of a model: posterior mean plots and posterior predictive plots.
Start by loading the data on 2016 annual income, and taking sample of 150 individuals to make our posterior distributions less precise.
d <- read.csv("https://soci620.netlify.com/data/USincome_subset_2016.csv")
# take a random sample of 150 rows
d <- d[sample(nrow(d),150),]
d
id | family_size | children_in_hh | female | age | marital_status | times_married | race | hispanic | hispanic_origin | health_ins | health_ins_priv | health_ins_employer | in_school | ed_highschool | ed_college | income | poverty | occ_income | log_income | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
<int> | <int> | <int> | <int> | <int> | <chr> | <int> | <chr> | <int> | <chr> | <int> | <int> | <int> | <int> | <int> | <int> | <int> | <lgl> | <int> | <dbl> | |
606 | 571365 | 3 | 0 | 0 | 36 | Never married/single | 0 | White | 0 | Not Hispanic | 0 | 0 | 0 | 0 | 1 | 0 | 25800 | FALSE | 24 | 10.158130 |
2477 | 2387903 | 2 | 0 | 1 | 49 | Married, spouse present | 1 | Black/African American/Negro | 0 | Not Hispanic | 1 | 0 | 0 | 0 | 1 | 0 | 13700 | FALSE | NA | 9.525151 |
1105 | 1099431 | 1 | 0 | 0 | 75 | Divorced | 2 | White | 0 | Not Hispanic | 1 | 0 | 0 | 0 | 1 | 0 | 56600 | FALSE | NA | 10.943764 |
88 | 78464 | 4 | 2 | 0 | 30 | Never married/single | 0 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 30600 | FALSE | 41 | 10.328755 |
1800 | 1746745 | 1 | 0 | 1 | 60 | Never married/single | 0 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 18100 | FALSE | 23 | 9.803667 |
3134 | 2990863 | 4 | 2 | 0 | 43 | Married, spouse present | 1 | Other Asian or Pacific Islander | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 85000 | FALSE | 33 | 11.350407 |
2986 | 2869773 | 4 | 2 | 0 | 43 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 0 | 0 | 1 | 0 | 110000 | FALSE | 26 | 11.608236 |
2563 | 2461937 | 2 | 0 | 0 | 64 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 21000 | FALSE | 25 | 9.952278 |
3227 | 3077613 | 4 | 2 | 1 | 30 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 31000 | FALSE | 38 | 10.341742 |
980 | 961705 | 2 | 0 | 1 | 46 | Married, spouse present | 1 | Black/African American/Negro | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 1 | 125000 | FALSE | 33 | 11.736069 |
2466 | 2380693 | 1 | 0 | 0 | 85 | Widowed | 1 | White | 0 | Not Hispanic | 1 | 0 | 0 | 0 | 1 | 0 | 23300 | NA | NA | 10.056209 |
2504 | 2409271 | 2 | 0 | 1 | 62 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 24000 | FALSE | 25 | 10.085809 |
2289 | 2210649 | 1 | 0 | 0 | 60 | Divorced | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 1 | 38000 | FALSE | 27 | 10.545341 |
3119 | 2974225 | 2 | 0 | 1 | 53 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 36000 | FALSE | 25 | 10.491274 |
1911 | 1856605 | 2 | 0 | 0 | 36 | Married, spouse present | 1 | White | 0 | Mexican | 1 | 1 | 1 | 0 | 1 | 0 | 30000 | FALSE | 34 | 10.308953 |
1600 | 1541948 | 2 | 0 | 0 | 58 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 46800 | FALSE | 19 | 10.753638 |
70 | 63251 | 2 | 0 | 0 | 65 | Divorced | 1 | White | 0 | Not Hispanic | 1 | 0 | 0 | 0 | 1 | 1 | 412000 | FALSE | 62 | 12.928779 |
2344 | 2271390 | 3 | 0 | 1 | 25 | Never married/single | 0 | American Indian or Alaska Native | 0 | Not Hispanic | 1 | 0 | 0 | 0 | 1 | 0 | 18800 | FALSE | 13 | 9.841612 |
1289 | 1250200 | 2 | 0 | 0 | 48 | Never married/single | 0 | White | 0 | Not Hispanic | 1 | 0 | 0 | 0 | 1 | 0 | 13800 | FALSE | NA | 9.532424 |
2965 | 2848590 | 2 | 0 | 1 | 67 | Married, spouse present | 2 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 24400 | FALSE | 22 | 10.102338 |
747 | 736019 | 6 | 3 | 1 | 26 | Never married/single | 0 | Black/African American/Negro | 0 | Not Hispanic | 1 | 0 | 0 | 0 | 1 | 0 | 5500 | FALSE | 22 | 8.612503 |
846 | 828071 | 4 | 2 | 0 | 43 | Never married/single | 0 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 60000 | FALSE | 33 | 11.002100 |
1151 | 1137032 | 5 | 2 | 1 | 48 | Married, spouse present | 3 | White | 0 | Not Hispanic | 1 | 0 | 0 | 0 | 1 | 0 | 15800 | FALSE | 23 | 9.667765 |
1013 | 997966 | 4 | 2 | 1 | 31 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 95000 | FALSE | 21 | 11.461632 |
435 | 409170 | 1 | 0 | 0 | 56 | Divorced | 1 | White | 0 | Not Hispanic | 1 | 0 | 0 | 0 | 1 | 0 | 5000 | TRUE | 20 | 8.517193 |
2870 | 2775303 | 4 | 0 | 0 | 49 | Never married/single | 0 | White | 0 | Mexican | 1 | 0 | 0 | 0 | 1 | 0 | 8800 | FALSE | NA | 9.082507 |
612 | 577818 | 1 | 0 | 0 | 62 | Married, spouse absent | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 1 | 114500 | FALSE | 42 | 11.648330 |
583 | 548296 | 1 | 0 | 0 | 93 | Widowed | 1 | White | 0 | Not Hispanic | 1 | 1 | 0 | 0 | 1 | 0 | 14100 | FALSE | NA | 9.553930 |
457 | 428496 | 1 | 0 | 1 | 28 | Never married/single | 0 | Two major races | 0 | Not Hispanic | 0 | 0 | 0 | 0 | 1 | 0 | 28800 | FALSE | 42 | 10.268131 |
197 | 183666 | 2 | 0 | 1 | 46 | Divorced | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 85000 | FALSE | 33 | 11.350407 |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
2293 | 2216204 | 2 | 0 | 1 | 76 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 0 | 0 | 7000 | FALSE | NA | 8.853665 |
3295 | 3134474 | 1 | 0 | 1 | 18 | Never married/single | 0 | White | 0 | Not Hispanic | 1 | 1 | 1 | 1 | 1 | 0 | 2200 | NA | 13 | 7.696213 |
2923 | 2813298 | 1 | 0 | 1 | 40 | Divorced | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 1 | 70000 | FALSE | 27 | 11.156251 |
1284 | 1247595 | 5 | 3 | 1 | 33 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 1 | 1 | 1 | 46000 | FALSE | 27 | 10.736397 |
3012 | 2891885 | 2 | 0 | 1 | 69 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 12000 | FALSE | NA | 9.392662 |
2396 | 2318689 | 2 | 0 | 1 | 64 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 31800 | FALSE | 33 | 10.367222 |
320 | 305240 | 1 | 0 | 0 | 27 | Never married/single | 0 | Other race, nec | 0 | Mexican | 1 | 1 | 1 | 0 | 1 | 0 | 90000 | FALSE | 42 | 11.407565 |
1231 | 1208900 | 3 | 1 | 1 | 30 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 0 | 0 | 15000 | FALSE | 22 | 9.615805 |
999 | 982895 | 6 | 3 | 1 | 46 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 33800 | FALSE | 43 | 10.428216 |
1774 | 1719562 | 4 | 2 | 1 | 42 | Married, spouse present | 1 | Japanese | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 103000 | FALSE | 42 | 11.542484 |
2667 | 2577590 | 1 | 0 | 0 | 49 | Never married/single | 0 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 39000 | FALSE | 38 | 10.571317 |
69 | 62895 | 2 | 0 | 0 | 51 | Married, spouse present | 2 | White | 0 | Not Hispanic | 0 | 0 | 0 | 0 | 1 | 0 | 30000 | FALSE | 25 | 10.308953 |
1000 | 984363 | 4 | 2 | 0 | 48 | Married, spouse present | 1 | Other Asian or Pacific Islander | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 132000 | FALSE | 33 | 11.790557 |
446 | 421992 | 1 | 0 | 1 | 62 | Divorced | 1 | Black/African American/Negro | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 40000 | FALSE | 16 | 10.596635 |
2742 | 2649201 | 4 | 2 | 0 | 44 | Married, spouse present | 1 | Black/African American/Negro | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 13000 | FALSE | 42 | 9.472705 |
2897 | 2791523 | 6 | 1 | 1 | 55 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 0 | 0 | 0 | 1 | 0 | 7000 | FALSE | NA | 8.853665 |
1167 | 1155428 | 1 | 0 | 0 | 89 | Married, spouse absent | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 1 | 34000 | NA | NA | 10.434116 |
2353 | 2279738 | 1 | 0 | 0 | 38 | Divorced | 1 | White | 0 | Not Hispanic | 0 | 0 | 0 | 0 | 1 | 0 | 26000 | FALSE | 23 | 10.165852 |
3009 | 2886864 | 1 | 0 | 0 | 27 | Never married/single | 0 | Other Asian or Pacific Islander | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 1 | 22000 | FALSE | 42 | 9.998798 |
1940 | 1882346 | 1 | 0 | 0 | 33 | Never married/single | 0 | Black/African American/Negro | 0 | Not Hispanic | 0 | 0 | 0 | 0 | 1 | 0 | 28000 | FALSE | 22 | 10.239960 |
1586 | 1529038 | 5 | 3 | 1 | 34 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 0 | 0 | 1 | 0 | 24000 | FALSE | 36 | 10.085809 |
2197 | 2115506 | 1 | 0 | 1 | 21 | Never married/single | 0 | White | 0 | Not Hispanic | 1 | 1 | 1 | 1 | 1 | 0 | 9000 | TRUE | 24 | 9.104980 |
2311 | 2231587 | 4 | 2 | 1 | 54 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 22000 | FALSE | 25 | 9.998798 |
2144 | 2071582 | 2 | 1 | 0 | 55 | Divorced | 2 | White | 0 | Not Hispanic | 1 | 0 | 0 | 0 | 1 | 0 | 25000 | FALSE | NA | 10.126631 |
19 | 19751 | 3 | 1 | 0 | 33 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 1 | 75000 | FALSE | 42 | 11.225243 |
395 | 379665 | 2 | 1 | 1 | 88 | Widowed | 1 | White | 0 | Not Hispanic | 1 | 1 | 0 | 0 | 1 | 0 | 18800 | FALSE | NA | 9.841612 |
3276 | 3118806 | 3 | 1 | 0 | 56 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 22900 | FALSE | 19 | 10.038892 |
736 | 722825 | 4 | 2 | 1 | 46 | Married, spouse present | 1 | Other Asian or Pacific Islander | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 33000 | FALSE | 24 | 10.404263 |
2714 | 2621194 | 3 | 1 | 0 | 32 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 0 | 22000 | FALSE | 16 | 9.998798 |
1227 | 1207810 | 2 | 0 | 0 | 70 | Married, spouse present | 1 | White | 0 | Not Hispanic | 1 | 1 | 1 | 0 | 1 | 1 | 43800 | FALSE | 80 | 10.687389 |
We can start by building a simple model that predicts a person's log-income using just their age.
library(rethinking)
m1 <- alist(
log_income ~ dnorm(mu,sigma),
mu <- a + b1*age,
a ~ dnorm(8,3),
b1 ~ dnorm(0,2),
sigma ~ dunif(0,5)
)
fit1 <- quap(m1,data=d)
# look at the estimates
precis(fit1)
Loading required package: rstan Loading required package: StanHeaders Loading required package: ggplot2 rstan (Version 2.21.2, GitRev: 2e1f913d3ca3) For execution on a local, multicore CPU with excess RAM we recommend calling options(mc.cores = parallel::detectCores()). To avoid recompilation of unchanged Stan programs, we recommend calling rstan_options(auto_write = TRUE) Loading required package: parallel rethinking (Version 2.13) Attaching package: ‘rethinking’ The following object is masked from ‘package:stats’: rstudent
mean | sd | 5.5% | 94.5% | |
---|---|---|---|---|
<dbl> | <dbl> | <dbl> | <dbl> | |
a | 9.895820675 | 0.258221894 | 9.483132215 | 10.30850914 |
b1 | 0.006979848 | 0.005069483 | -0.001122166 | 0.01508186 |
sigma | 1.132539930 | 0.065388349 | 1.028036719 | 1.23704314 |
If the estimate for age
seems small, that is because we are measuring age in years and the range of ages in the data is large.
Let's re-estimate using standardized age.
# just to look at the age range in the data
range(d$age)
# construct standardized age
d$age.s <- scale(d$age)
range(d$age.s)
# recreate fit1, using standardized age
fit1 <- quap(alist(
log_income ~ dnorm(mu,sigma),
mu <- a + b1*age.s,
a ~ dnorm(8,3),
b1 ~ dnorm(0,2),
sigma ~ dunif(0,5)
),data=d)
# look at the estimates
precis(fit1)
exp(fit1@coef[2])
mean | sd | 5.5% | 94.5% | |
---|---|---|---|---|
<dbl> | <dbl> | <dbl> | <dbl> | |
a | 10.2279588 | 0.09247971 | 10.08015832 | 10.3757592 |
b1 | 0.1222032 | 0.09273354 | -0.02600287 | 0.2704093 |
sigma | 1.1331757 | 0.06547991 | 1.02852616 | 1.2378252 |
We can start to look at the predictions of our model by plotting our "point estimate" of $\mu$ as a function of standardized age.
# set up plot window for projector visibility
options(repr.plot.width=12, repr.plot.height=8)
# first plot the data points
plot(d$age.s, d$log_income)
# note: here we used the plot(x,y) form of the plot function,
# which plots x on the horizontal axis against y on the vertical.
# There is another way to make the same scatter plot using the
# 'formula' form of plot:
# plot(log_income~age.s,data=d)
# This takes the form of plot(y~x, data=d). The two forms create
# identical output.
# get the coefficients
cf1 <- coef(fit1)
cf1
a_fit <- cf1['a']
b_fit <- cf1['b1']
# create a grid of hypothetical ages to predict
n_grid <- 400
grid <- seq(min(d$age.s),max(d$age.s),length.out=n_grid)
# plot the "deterministic" portion of the model:
# mu <- a + b1*age.s
lines(grid , a_fit + b_fit*grid)
This shows us the overall trend predicted by the data, but ignores the uncertainty in our estimates of $\mu$ that comes from uncertainty in a
and b1
.
One simple way to try to express this uncertainty is to draw the same line as above multiple times, using draws from our posterior rather than the maximum a-posteriori estimates.
# get a sample of coefficients
n_samp <- 300
s <- extract.samples(fit1,n_samp)
head(s)
# plot the data
plot(d$age.s,d$log_income)
# create a grid of hypothetical ages to predict
n_grid <- 400
grid <- seq(min(d$age.s),max(d$age.s),length.out=n_grid)
# create a loop to plot each of the lines
for(i in 1:n_samp){
# each time we go through this loop, i will take
# on one of the values between 1 and n_samp (30).
# get the coefficients from the i'th row of the sample s
a_fit <- s$a[i]
b_fit <- s$b1[i]
# plot the line, but make it slightly "translucent" for clarity
lines(grid , a_fit + b_fit*grid , col="#00000011")
}
a | b1 | sigma | |
---|---|---|---|
<dbl> | <dbl> | <dbl> | |
1 | 10.32448 | 0.17737059 | 1.1066351 |
2 | 10.29114 | 0.03338418 | 1.0486046 |
3 | 10.26172 | 0.11108619 | 1.1285676 |
4 | 10.19481 | 0.02294440 | 1.0948543 |
5 | 10.27887 | 0.15726755 | 0.9869483 |
6 | 10.13689 | 0.34857763 | 1.1873645 |
This gives us an idea of the spread of uncertainty in our model, but can be hard to read or interpret. A common approach here is to replace the pile of spaghetti in this plot with a shaded region that will show, say, the 90% credible interval for our posterior.
The procedure for this is straightforward:
grid
, draw a large number of posterior samples (say 1000) for a
and b1
.grid
.There are lots of ways to do this in R, and the rethinking
package has some tools to make it easier. But we will first do it "manually" to make it clear what is going on.
# create a grid of hypothetical ages to predict
n_grid <- 400
grid <- seq(min(d$age.s),max(d$age.s),length.out=n_grid)
# how many samples at each grid point?
n_samp <- 1000
s1 <- extract.samples(fit1,n_samp)
head(s1)
# use the `sapply()` function to calculate the 90% credible
# interval at each point x in the grid
mu_ci <- sapply(grid,function(x){
# this is like a loop, where x will take on every value
# in the grid
quantile(s1$a + s1$b1 * x , probs=c(.05,.95))
})
# see what the output of this looks like
mu_ci
# plot the data
plot(d$age.s,d$log_income,pch=16)
# use the `polygon` function to plot a shaded region
polygon(
x = c(grid,rev(grid)), # rev() returns its argument in reverse order
y = c(mu_ci[1,],rev(mu_ci[2,])),
border = NA , col = '#00000055')
a | b1 | sigma | |
---|---|---|---|
<dbl> | <dbl> | <dbl> | |
1 | 10.31780 | 0.04336469 | 1.223364 |
2 | 10.28770 | 0.05917803 | 1.002291 |
3 | 10.34150 | 0.15617497 | 1.207649 |
4 | 10.22377 | 0.30798533 | 1.060152 |
5 | 10.24130 | -0.07191926 | 1.036488 |
6 | 10.27759 | 0.15282044 | 1.147524 |
5% | 9.741618 | 9.744198 | 9.745957 | 9.748336 | 9.75127 | 9.754089 | 9.756498 | 9.758661 | 9.760387 | 9.761957 | ⋯ | 10.13678 | 10.13619 | 10.13560 | 10.13475 | 10.13386 | 10.13297 | 10.13207 | 10.13118 | 10.13029 | 10.12940 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
95% | 10.296917 | 10.297433 | 10.297591 | 10.297695 | 10.29797 | 10.298574 | 10.298795 | 10.298737 | 10.298675 | 10.298585 | ⋯ | 10.89463 | 10.89696 | 10.89928 | 10.90159 | 10.90390 | 10.90621 | 10.90855 | 10.91125 | 10.91395 | 10.91664 |
Looks good! But a bit of a pain. Let's try the same thing using some functions from the rethinking
package:
link()
will draw samples from posterior distribution of $\mu$ and calculated predicted values from the data data.PI()
calculates the percentile interval on a vector.shade()
is just a convenient wrapper around polygon()
that we used above.# create a grid of hypothetical ages to predict
n_grid <- 400
# create a data frame listing our grid as possible values of age.s
pred_dat <- data.frame(
age.s = seq(min(d$age.s),max(d$age.s),length.out=n_grid)
)
head(pred_dat)
# how many samples at each grid point?
n_samp <- 1000
# use the 'link' function to take n_samp posterior samples at each point in pred_dat
mu_post <- link(fit1, data=pred_dat, n=n_samp)
# mu_post has one row for each of our posterior samples, and one column
# for each point in our grid.
dim(mu_post)
age.s | |
---|---|
<dbl> | |
1 | -1.677729 |
2 | -1.667288 |
3 | -1.656847 |
4 | -1.646406 |
5 | -1.635964 |
6 | -1.625523 |
# calculate the credible interval for every column of mu_post
# `apply()` here will call PI(mu[,i],prob=0.9) for every column i.
mu_ci <- apply(mu_post, 2, PI, prob=0.9)
mu_ci
5% | 9.725457 | 9.728678 | 9.731687 | 9.734335 | 9.736982 | 9.739629 | 9.742277 | 9.744924 | 9.747572 | 9.750231 | ⋯ | 10.14056 | 10.14083 | 10.1411 | 10.14136 | 10.14163 | 10.14189 | 10.14216 | 10.14242 | 10.14269 | 10.14246 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
95% | 10.307084 | 10.307032 | 10.307035 | 10.307017 | 10.307045 | 10.307072 | 10.307054 | 10.307005 | 10.306957 | 10.306909 | ⋯ | 10.89348 | 10.89654 | 10.8996 | 10.90217 | 10.90457 | 10.90695 | 10.90915 | 10.91148 | 10.91415 | 10.91657 |
# plot the data
plot(d$age.s,d$log_income,pch=16)
# use the `shade()` function from rethinking to draw the shaded region
shade(mu_ci, pred_dat$age.s)