Week 2 Linear Models

The second week covers Chapter 4 (Geocentric Models).

2.1 Lectures

Lecture 3:

Lecture 4:

2.2 Exercises

2.2.1 Chapter 4

4E1. In the model definition below, which line is the likelihood? \[\begin{align} y_i &\sim \text{Normal}(\mu,\sigma) \\ \mu &\sim \text{Normal}(0,10) \\ \sigma &\sim \text{Exponential}(1) \end{align}\]

The likelihood is given by the first line, \(y_i \sim \text{Normal}(\mu,\sigma)\). The other lines represent the prior distributions for \(\mu\) and \(\sigma\).

4E2. In the model definition just above, how many parameters are in the posterior distribution?

There are two parameters, \(\mu\) and \(\sigma\).

4E3. Using the model definition above, write down the appropriate form of Bayes’ theorem that includes the proper likelihood and priors.

\[ \text{Pr}(\mu,\sigma|y) = \frac{\prod_i\text{Normal}(y_i|\mu,\sigma)\text{Normal}(\mu|0,10)\text{Exponential}(\sigma|1)}{\int\int\prod_i\text{Normal}(y_i|\mu,\sigma)\text{Normal}(\mu|0,10)\text{Exponential}(\sigma|1)d \mu d \sigma} \]

4E4. In the model definition below, which line is the linear model? \[\begin{align} y_i &\sim \text{Normal}(\mu,\sigma) \\ \mu_i &= \alpha + \beta x_i \\ \alpha &\sim \text{Normal}(0,10) \\ \beta &\sim \text{Normal}(0,1) \\ \sigma &\sim \text{Exponential}(2) \end{align}\]

The linear model is the second line, \(\mu_i = \alpha + \beta x_i\).

4E5. In the model definition just above, how many parameters are in the posterior distribution?

There are now three model parameters: \(\alpha\), \(\beta\), and \(\sigma\). The mean, \(\mu\) is no longer a parameter, as it is defined deterministically, as a function of other parameters in the model.

4M1. For the model definition below, simulate observed y values from the prior (not the posterior). \[\begin{align} y_i &\sim \text{Normal}(\mu,\sigma) \\ \mu &\sim \text{Normal}(0,10) \\ \sigma &\sim \text{Exponential}(1) \end{align}\]

library(tidyverse)

sim <- tibble(mu = rnorm(n = 10000, mean = 0, sd = 10),
       sigma = rexp(n = 10000, rate = 1)) %>%
  mutate(y = rnorm(n = 10000, mean = mu, sd = sigma))

ggplot(sim, aes(x = y)) +
  geom_density() +
  labs(x = "y", y = "Density")

4M2. Trasnlate the model just above into a quap formula.

flist <- alist(
  y ~ dnorm(mu, sigma),
  mu ~ dnorm(0, 10),
  sigma ~ dexp(1)
)

4M3. Translate the quap model formula below into a mathematical model definition.

y ~ dnorm( mu , sigma ),
mu <- a + b*x,
a ~ dnorm( 0 , 10 ),
b ~ dunif( 0 , 1 ),
sigma ~ dexp( 1 )

\[\begin{align} y_i &\sim \text{Normal}(\mu_i,\sigma) \\ \mu_i &= a + bx_i \\ a &\sim \text{Normal}(0, 10) \\ b &\sim \text{Uniform}(0, 1) \\ \sigma &\sim \text{Exponential}(1) \end{align}\]

4M4. A sample of students is measured for height each year for 3 years. After the third year, you want to fit a linear regression predicting height using year as a predictor. Write down the mathematical model definition for this regression, using any variable names and priors you choose. Be prepared to defend your choice of priors.

\[\begin{align} h_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta(x_i - \bar{x}) \\ \alpha &\sim \text{Normal}(169, 20) \\ \beta &\sim \text{Normal}(0, 8) \\ \sigma &\sim \text{Exponential}(1) \end{align}\]

Because height is centered, \(\alpha\) represents the average height in the average year (i.e., year 2). The prior of \(\text{Normal}(165,20)\) was chosen assuming that height is measured in centimeters and that that sample is of adults (e.g., college students). The mean of 169cm represents the average height people aged 20 and older in the United States (Fryar et al., 2016). Further, Fryar et al. (2016) report the 5th and 95th percentiles for females and males as [150,174] and [163,188], respectively. A standard deviation of 20 presumes 95% of the probability between 169 \(\pm\) 40, or [129,209], easily encapsulating the ranges seen in previous observed data.

At this point, we have no reason to believe that year would have any effect on height, as we’re assuming that all students are adults, and therefore likely done growing. Accordingly, we choose a prior centered on zero. The standard deviation of the prior of 5 represents plausible growth (or shrinkage). During growth spurts, height growth averages 6–13 cm/year. The standard deviation of 8 encompasses the range we might expect to see if growth were occurring at a high rate.

Prior predictive simulations for \(\alpha\) and \(\beta\) also appear to give reasonably plausible regression lines, given our current assumptions.

n <- 1000
tibble(group = seq_len(n),
       alpha = rnorm(n, 169, 20),
       beta = rnorm(n, 0, 8)) %>%
  expand(nesting(group, alpha, beta), year = c(1, 3)) %>%
  mutate(height = alpha + beta * (year - 2)) %>%
  ggplot(aes(x = year, y = height, group = group)) +
  geom_line(alpha = 1/10) +
  geom_hline(yintercept = c(0, 272), linetype = 2:1, color = "red") +
  annotate(geom = "text", x = 1, y = 0, hjust = 0, vjust = 1,
           label = "Embryo") +
  annotate(geom = "text", x = 1, y = 272, hjust = 0, vjust = 0,
           label = "World's tallest person (272cm)") +
  coord_cartesian(ylim = c(-25, 300)) +
  labs(x = "Year", y = "Height")

Finally, the exponential prior on \(\sigma\) assumes an average deviation of 1.

4M5. Now suppose I remind you that every student got taller each year. Does this information lead you to change your choice of priors? How?

Yes. Because we know that an increase in year will always lead to increased height, we know that \(\beta\) will be positive. Therefore, our prior should reflect this by using, for example, a log-normal distribution.

\[ \beta \sim \text{Log-Normal}(2,0.5) \]

This log-normal prior has an 89% highest density interval between about 2.3 and 14.1cm growth per year.

library(tidybayes)

samples <- rlnorm(1e8, 2, 0.5)
bounds <- median_hdi(samples, .width = 0.89)

ggplot() +
  stat_function(data = tibble(x = c(0, 30)), mapping = aes(x = x),
                geom = "line", fun = dlnorm,
                args = list(meanlog = 2, sdlog = 0.5)) +
  geom_ribbon(data = tibble(x = seq(bounds$ymin, bounds$ymax, 0.01)),
              aes(x = x, ymin = 0, ymax = dlnorm(x, 2, 0.5)),
              alpha = 0.8) +
  labs(x = expression(beta), y = "Density")

Prior predictive simulations of plausible lines using this new log-normal prior indicate that these priors still represent plausible values, while also constraining all lines to have a positive slope.

n <- 1000
tibble(group = seq_len(n),
       alpha = rnorm(n, 169, 20),
       beta = rlnorm(n, 2, 0.75)) %>%
  expand(nesting(group, alpha, beta), year = c(1, 3)) %>%
  mutate(height = alpha + beta * (year - 2)) %>%
  ggplot(aes(x = year, y = height, group = group)) +
  geom_line(alpha = 1/10) +
  geom_hline(yintercept = c(0, 272), linetype = 2:1, color = "red") +
  annotate(geom = "text", x = 1, y = 0, hjust = 0, vjust = 1,
           label = "Embryo") +
  annotate(geom = "text", x = 1, y = 272, hjust = 0, vjust = 0,
           label = "World's tallest person (272cm)") +
  coord_cartesian(ylim = c(-25, 300)) +
  labs(x = "Year", y = "Height")

4M6. Now suppose I tell you that the variance among heights for students of the same age is never more than 64cm. How does this lead you to revise your priors?

A variance of 64cm corresponds to a standard deviation of 8cm. Our current prior of \(\sigma \sim \text{Exponential}(1)\) gives very little probability mass to values greater than 8. However, we can also see that the vast majority of the probability mass for this prior is less than 5cm. Therefore, if we think the standard deviation could get close to 8cm, we might consider reducing the rate of the exponential prior to allow for slightly larger values.

ggplot(tibble(x = c(0, 20)), aes(x = x)) +
  stat_function(fun = dexp, args = list(rate = 1)) +
  labs(x = expression(sigma), y = "Density")

4M7. Refit model m4.3 from the chapter, but omit the mean weight xbar this time. Compare the new model’s posterior to that of the original model. In particular, look at the covariance among the parameters. What is different? Then compare the posterior predictions of both models.

First, we’ll reproduce m4.3 using quap(), and then again using brm(). The covariance matrix from both is the same, as we would expect!

library(rethinking)
library(brms)

data(Howell1)
how_dat <- Howell1 %>%
  filter(age >= 18) %>%
  mutate(weight_c = weight - mean(weight))

# first, duplicate model with `quap`
m4.3 <- quap(alist(height ~ dnorm(mu, sigma),
                   mu <- a + b * (weight_c),
                   a ~ dnorm(178, 20),
                   b ~ dlnorm(0, 1),
                   sigma ~ dunif(0, 50)),
             data = how_dat)

round(vcov(m4.3), 3)
#>           a     b sigma
#> a     0.073 0.000 0.000
#> b     0.000 0.002 0.000
#> sigma 0.000 0.000 0.037

# and then with brms
b4.3 <- brm(height ~ 1 + weight_c, data = how_dat, family = gaussian,
            prior = c(prior(normal(178, 20), class = Intercept),
                      prior(lognormal(0, 1), class = b, lb = 0),
                      prior(uniform(0, 50), class = sigma)),
            iter = 28000, warmup = 27000, chains = 4, cores = 4, seed = 1234,
            file = here("fits", "chp4", "b4.3-0"))

posterior_samples(b4.3) %>%
  select(-lp__) %>%
  cov() %>%
  round(digits = 3)
#>             b_Intercept b_weight_c sigma
#> b_Intercept       0.074      0.000 0.000
#> b_weight_c        0.000      0.002 0.000
#> sigma             0.000      0.000 0.034

Now, let’s using the non-centered parameterization. This time, we’ll only use brm().

b4.3_nc <- brm(height ~ 1 + weight, data = how_dat, family = gaussian,
               prior = c(prior(normal(178, 20), class = Intercept),
                         prior(lognormal(0, 1), class = b, lb = 0),
                         prior(uniform(0, 50), class = sigma)),
               iter = 28000, warmup = 27000, chains = 4, cores = 4, seed = 1234,
               file = here("fits", "chp4", "b4.3_nc"))

posterior_samples(b4.3_nc) %>%
  select(-lp__) %>%
  cov() %>%
  round(digits = 3)
#>             b_Intercept b_weight  sigma
#> b_Intercept       3.634   -0.079 -0.016
#> b_weight         -0.079    0.002  0.000
#> sigma            -0.016    0.000  0.037

We now see non-zero covariances between the parameters. Lets compare the posterior predictions. We’ll generate hypothetical outcome plots which are animated to show the uncertainty in estimates (see Hullman et al., 2015; Kale et al., 2019). Here we’ll just animate the estimated regression line using {gganimate} (Pedersen & Robinson, 2020). We can see that the predictions from the two models are nearly identical.

library(gganimate)

weight_seq <- tibble(weight = seq(25, 70, length.out = 100)) %>%
  mutate(weight_c = weight - mean(how_dat$weight))

predictions <- bind_rows(
  predict(b4.3, newdata = weight_seq) %>%
    as_tibble() %>%
    bind_cols(weight_seq) %>%
    mutate(type = "Centered"),
  predict(b4.3_nc, newdata = weight_seq) %>%
    as_tibble() %>%
    bind_cols(weight_seq) %>%
    mutate(type = "Non-centered")
)

fits <- bind_rows(
  weight_seq %>%
    add_fitted_draws(b4.3) %>%
    mutate(type = "Centered"),
  weight_seq %>%
    add_fitted_draws(b4.3_nc) %>%
    mutate(type = "Non-centered")
) %>%
  ungroup()

bands <- fits %>%
  group_by(type, weight) %>%
  median_qi(.value, .width = c(.67, .89, .97))

lines <- fits %>%
  filter(.draw %in% sample(unique(.data$.draw), size = 50))

ggplot(lines, aes(x = weight)) +
  facet_wrap(~type, nrow = 1) +
  geom_ribbon(data = predictions, aes(ymin = Q2.5, ymax = Q97.5), alpha = 0.3) +
  geom_lineribbon(data = bands, aes(y = .value, ymin = .lower, ymax = .upper),
                  color = NA) +
  scale_fill_brewer(palette = "Blues", breaks = c(.67, .89, .97)) +
  geom_line(aes(y = .value, group = .draw)) +
  geom_point(data = how_dat, aes(y = height), shape = 1, alpha = 0.7) +
  labs(x = "Weight", y = "Height", fill = "Interval") +
  theme(legend.position = "bottom") +
  transition_states(.draw, 0, 1)

4M8. In the chapter, we used 15 knots with the cherry blossom spline. Increase the number of know and observe what happens to the resulting spline. Then adjust also the width of the prior on the weights—change the standard deviation of the prior and watch what happens. What do you think the combination of know number and the prior on the weights controls?

Because {brms} estimates splines a little differently than {rethinking} (see here), I’ll use {rethinking} for estimating the spline models. First lets look again at the 15-knot spline model, and then double the number of knots and play with the prior. As expected, increasing the number of knots increases the “wiggly-ness” of the spline. We can also see that tightening the prior on the weights takes away some of the “wiggly-ness”.

library(splines)
library(tidybayes.rethinking)
library(colorblindr)

data(cherry_blossoms)
cb_dat <- cherry_blossoms %>%
  drop_na(doy)

# original m4.7 model
knots_15 <- quantile(cb_dat$year, probs = seq(0, 1, length.out = 15))
B_15 <- bs(cb_dat$year, knots = knots_15[-c(1, 15)],
        degree = 3, intercept = TRUE)

m4.7 <- quap(alist(D ~ dnorm(mu, sigma),
                   mu <- a + B_15 %*% w,
                   a ~ dnorm(100, 10),
                   w ~ dnorm(0, 10),
                   sigma ~ dexp(1)),
             data = list(D = cb_dat$doy, B = B_15),
             start = list(w = rep(0, ncol(B_15))))

# double the number of knots
knots_30 <- quantile(cb_dat$year, probs = seq(0, 1, length.out = 30))
B_30 <- bs(cb_dat$year, knots = knots_30[-c(1, 30)],
        degree = 3, intercept = TRUE)

m4.7_30 <- quap(alist(D ~ dnorm(mu, sigma),
                      mu <- a + B_30 %*% w,
                      a ~ dnorm(100, 10),
                      w ~ dnorm(0, 10),
                      sigma ~ dexp(1)),
                data = list(D = cb_dat$doy, B_30 = B_30),
                start = list(w = rep(0, ncol(B_30))))

# and modify the prior
m4.7_30_tight_prior <- quap(alist(D ~ dnorm(mu, sigma),
                                  mu <- a + B_30 %*% w,
                                  a ~ dnorm(100, 10),
                                  w ~ dnorm(0, 2),
                                  sigma ~ dexp(1)),
                            data = list(D = cb_dat$doy, B_30 = B_30),
                            start = list(w = rep(0, ncol(B_30))))

# create plot data
spline_15 <- fitted_draws(m4.7, newdata = cb_dat) %>%
  group_by(.row) %>%
  slice_sample(n = 1000) %>%
  group_by(year) %>%
  median_hdci(.value, .width = 0.89) %>%
  mutate(knots = "15 knots")

spline_30 <- fitted_draws(m4.7_30, newdata = cb_dat) %>%
  group_by(.row) %>%
  slice_sample(n = 1000) %>%
  group_by(year) %>%
  median_hdci(.value, .width = 0.89) %>%
  mutate(knots = "30 knots")

spline_30_tight_prior <- fitted_draws(m4.7_30_tight_prior, newdata = cb_dat) %>%
  group_by(.row) %>%
  slice_sample(n = 1000) %>%
  group_by(year) %>%
  median_hdci(.value, .width = 0.89) %>%
  mutate(knots = "30 knots; Tight prior")

all_splines <- bind_rows(spline_15, spline_30, spline_30_tight_prior)

ggplot(cb_dat, aes(x = year)) +
  geom_point(aes(y = doy), alpha = 0.2) +
  geom_ribbon(data = all_splines, aes(ymin = .lower, ymax = .upper),
              alpha = 0.7, fill = palette_OkabeIto[1]) +
  facet_wrap(~knots, ncol = 1) +
  labs(x = "Year", y = "Day in Year")

:::question > 4H1. The weights listed below were recored in the !Kung census, but heights were not recorded for these individuals. Provide predicted heights and 89% intervals for each of these individuals. That is, fill in the table, below, using model-based predictions.

Individual weight expected height 89% interval
1 47.0
2 43.7
3 64.8
4 32.6
5 54.6

:::

tibble(individual = 1:5,
       weight = c(46.95, 43.72, 64.78, 32.59, 54.63)) %>%
  mutate(weight_c = weight - mean(how_dat$weight)) %>%
  add_predicted_draws(b4.3) %>%
  group_by(individual, weight) %>%
  mean_qi(.prediction, .width = 0.89) %>%
  mutate(range = glue("[{sprintf('%0.1f', .lower)}--",
                      "{sprintf('%0.1f', .upper)}]"),
         .prediction = sprintf("%0.1f", .prediction)) %>%
  select(individual, weight, exp = .prediction, range) %>%
  kbl(align = "c", booktabs = TRUE,
      col.names = c("Individual", "weight", "expected height", "89% interval"))
Individual weight expected height 89% interval
1 47.0 156.3 [148.0–164.4]
2 43.7 153.4 [144.9–161.7]
3 64.8 172.5 [164.4–180.9]
4 32.6 143.3 [135.2–151.6]
5 54.6 163.3 [155.4–171.3]

4H2. Select out all the rows in the Howell1 data with ages below 18 years of age. If you do it right, you should end up with a new data frame with 192 rows in it.

young_how <- Howell1 %>%
  filter(age < 18) %>%
  mutate(weight_c = weight - mean(weight))
nrow(young_how)
#> [1] 192
  1. Fit a linear regression to these data, using quap. Present and interpret the estimates. For every 10 units of increase in weight, how much taller does the model predict a child gets?

We’ll use brms::brm() for fitting the model. We’ll use the same priors as for the model of adults, except for the weight for \(\beta\). Because children are shorter than adults, we’ll use a prior of \(\text{Normal}(138,20)\), based on the data reported in Fryar et al. (2016). Based on these estimates, an increase 10 units of weights corresponds to an average increase in height of 27.2 centimeters.

b4h2 <- brm(height ~ 1 + weight_c, data = young_how, family = gaussian,
            prior = c(prior(normal(138, 20), class = Intercept),
                      prior(lognormal(0, 1), class = b, lb = 0),
                      prior(exponential(1), class = sigma)),
            iter = 4000, warmup = 2000, chains = 4, cores = 4, seed = 1234,
            file = here("fits", "chp4", "b4h2.rds"))

summary(b4h2)
#>  Family: gaussian 
#>   Links: mu = identity; sigma = identity 
#> Formula: height ~ 1 + weight_c 
#>    Data: young_how (Number of observations: 192) 
#> Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
#>          total post-warmup samples = 8000
#> 
#> Population-Level Effects: 
#>           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept   108.34      0.61   107.14   109.54 1.00     7838     5731
#> weight_c      2.72      0.07     2.58     2.85 1.00     7926     5234
#> 
#> Family Specific Parameters: 
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma     8.36      0.42     7.59     9.21 1.00     6965     5695
#> 
#> Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
  1. Plot the raw data, with height on the vertical axis and weight on the horizontal axis. Superimpose the MAP regression line and 89% interval for the mean. Also superimpose the 89% interval for predicted heights.
library(modelr)

mod_fits <- young_how %>%
  data_grid(weight = seq_range(weight, 100)) %>%
  mutate(weight_c = weight - mean(young_how$weight)) %>%
  add_fitted_draws(b4h2) %>%
  group_by(weight) %>%
  mean_qi(.value, .width = 0.89)

mod_preds <- young_how %>%
  data_grid(weight = seq_range(weight, 100)) %>%
  mutate(weight_c = weight - mean(young_how$weight)) %>%
  add_predicted_draws(b4h2) %>%
  group_by(weight) %>%
  mean_qi(.prediction, .width = 0.89)

ggplot(young_how, aes(x = weight)) +
  geom_point(aes(y = height), alpha = 0.4) +
  geom_ribbon(data = mod_preds, aes(ymin = .lower, ymax = .upper),
              alpha = 0.2) +
  geom_lineribbon(data = mod_fits,
                  aes(y = .value, ymin = .lower, ymax = .upper),
                  fill = "grey60", size = 1) +
  labs(x = "Weight", y = "Height")

  1. What aspects of the model fit concern you? Describe the kinds of assumptions you would change, if any, to improve the model. You don’t have to write any new code. Just explain what the model appears to be doing a bad job of, and what you hypothesize would be a better model.

The model consistently over-estimates the height for individuals with weight less than ~13 and great than ~35. The model is also consistently underestimating the height of individuals with weight between ~13-35. Thus, our data appears to have a curve that our assumption of a straight line is violating. If we wanted to improve the model, we should relax the assumption of a straight line.

4H3. Suppose a colleauge of yours, who works on allometry, glances at the practice problems just above. Your colleague exclaims, “That’s silly. Everyone knows that it’s only the logarithm of body weight that scales with height!” Let’s take your colleague’s advice and see what happens.

  1. Model the relationship between height (cm) and the natural logarithm of weight (log-kg). Use the entire Howell1 data frame, all 544 rows, adults and non-adults. Can you interpret the resulting estimates?
full_how <- Howell1 %>%
  mutate(log_weight = log(weight),
         log_weight_c = log_weight - mean(log_weight))

b4h3 <- brm(height ~ 1 + log_weight_c, data = full_how, family = gaussian,
            prior = c(prior(normal(158, 20), class = Intercept),
                      prior(lognormal(0, 1), class = b, lb = 0),
                      prior(exponential(1), class = sigma)),
            iter = 4000, warmup = 2000, chains = 4, cores = 4, seed = 1234,
            file = here("fits", "chp4", "b4h3"))
summary(b4h3)
#>  Family: gaussian 
#>   Links: mu = identity; sigma = identity 
#> Formula: height ~ 1 + log_weight_c 
#>    Data: full_how (Number of observations: 544) 
#> Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
#>          total post-warmup samples = 8000
#> 
#> Population-Level Effects: 
#>              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept      138.26      0.22   137.83   138.71 1.00     7945     5700
#> log_weight_c    47.08      0.39    46.32    47.83 1.00     8210     5909
#> 
#> Family Specific Parameters: 
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma     5.13      0.16     4.84     5.44 1.00     8458     5616
#> 
#> Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

Conditional on the data and the model, the intercept estimate of 138.3 represents the predicted average height for an individual with an average log-weight (log-kg). The \(\beta\) estimate of 47.1 represents the average expected increase in height associated with an one-unit increase in weight (log-kg).

  1. Begin with this plot: plot( height ~ weight , data = Howell1 ). Then use samples fro the quadratic approximate posterior of the model in (a) to superimpose on the plot: (1) the predicted mean height as a function of weight, (2) the 97% interval for the mean, and (3) the 97% interval for predicted heights.
how_fits <- full_how %>%
  data_grid(weight = seq_range(weight, 100)) %>%
  mutate(log_weight = log(weight),
         log_weight_c = log_weight - mean(full_how$log_weight)) %>%
  add_fitted_draws(b4h3) %>%
  group_by(weight) %>%
  mean_qi(.value, .width = 0.97)

how_preds <- full_how %>%
  data_grid(weight = seq_range(weight, 100)) %>%
  mutate(log_weight = log(weight),
         log_weight_c = log_weight - mean(full_how$log_weight)) %>%
  add_predicted_draws(b4h3) %>%
  group_by(weight) %>%
  mean_qi(.prediction, .width = 0.97)

ggplot(full_how, aes(x = weight)) +
  geom_point(aes(y = height), alpha = 0.4) +
  geom_ribbon(data = how_preds, aes(ymin = .lower, ymax = .upper),
              alpha = 0.2) +
  geom_lineribbon(data = how_fits,
                  aes(y = .value, ymin = .lower, ymax = .upper),
                  fill = "grey60", size = 1) +
  labs(x = "Weight", y = "Height")

4H4. Plot the prior predictive distribution for the parabolic polynomial regression model in the chapter. You can modify the code that plots the linear regression prior predictive distribution. Can you modify the prior distributions of \(\alpha\), \(\beta_1\), and \(\beta_2\) so that the prior predictions stay within the biologically reasonable outcome space? That is to say: Do not try to fit the data by hand. But do try to keep the curves consistent with what you know about height and weight, before seeing these exact data.

The polynomial model from the chapter is defined as:

\[\begin{align} h_i &\sim \text{Normal}(\mu_i,\sigma) \\ \mu_i &= \alpha + \beta_1x_i + \beta_2x_i^2 \\ \alpha &\sim \text{Normal}(178,20) \\ \beta_1 &\sim \text{Log-Normal}(0,1) \\ \beta_2 &\sim \text{Normal}(0,1) \\ \sigma &\sim \text{Uniform}(0,50) \end{align}\]

First, let’s generate some prior predictive checks from the original priors.

n <- 1000
tibble(group = seq_len(n),
       alpha = rnorm(n, 178, 20),
       beta1 = rlnorm(n, 0, 1),
       beta2 = rnorm(n, 0, 1)) %>%
  expand(nesting(group, alpha, beta1, beta2),
         weight = seq(25, 70, length.out = 100)) %>%
  mutate(height = alpha + (beta1 * weight) + (beta2 * (weight ^ 2))) %>%
  ggplot(aes(x = weight, y = height, group = group)) +
  geom_line(alpha = 1 / 10) +
  geom_hline(yintercept = c(0, 272), linetype = 2:1, color = "red") +
  annotate(geom = "text", x = 25, y = 0, hjust = 0, vjust = 1,
           label = "Embryo") +
  annotate(geom = "text", x = 25, y = 272, hjust = 0, vjust = 0,
           label = "World's tallest person (272cm)") +
  coord_cartesian(ylim = c(-25, 300)) +
  labs(x = "Weight", y = "Height")

Clearly there is room for improvement here. However, because it’s not intuitive how exactly each parameter effects the parabolic curve, finding a good prior distribution is really hard! After much trial and error and playing with parabola calculators online, here is what I ended up with:

\[\begin{align} h_i &\sim \text{Normal}(\mu_i,\sigma) \\ \mu_i &= \alpha + \beta_1x_i + \beta_2x_i^2 \\ \alpha &\sim \text{Normal}(-190,5) \\ \beta_1 &\sim \text{Normal}(13,0.2) \\ \beta_2 &\sim \text{Uniform}(-0.13,-0.10) \\ \sigma &\sim \text{Uniform}(0,50) \end{align}\]

Which has the following prior predictive distribution.

n <- 1000
tibble(group = seq_len(n),
       alpha = rnorm(n, -190, 5),
       beta1 = rnorm(n, 13, 0.2),
       beta2 = runif(n, -0.13, -0.1)) %>%
  expand(nesting(group, alpha, beta1, beta2),
         weight = seq(25, 70, length.out = 100)) %>%
  mutate(height = alpha + (beta1 * weight) + (beta2 * (weight ^ 2))) %>%
  ggplot(aes(x = weight, y = height, group = group)) +
  geom_line(alpha = 1 / 10) +
  geom_hline(yintercept = c(0, 272), linetype = 2:1, color = "red") +
  annotate(geom = "text", x = 25, y = 0, hjust = 0, vjust = 1,
           label = "Embryo") +
  annotate(geom = "text", x = 25, y = 272, hjust = 0, vjust = 0,
           label = "World's tallest person (272cm)") +
  coord_cartesian(ylim = c(-25, 300)) +
  labs(x = "Weight", y = "Height")

4H5. Return to data(cherry_blossoms) and model the association between blossom date (doy) and March temperature (temp). Note that there are many missing values in both variables. You may consider a linear model, a polynomial, or a spline on temperature. How well does temperature trend predict the blossom trend?

We’ll try each type of model: linear, polynomial, and spline. For each, we’ll fit the model, and then visualize the predictions with the observed data. Because we’re no longer trying to replicate a spline model from the text, I’ll use brm() to estimate this spline model. Overall the predictions from each model are remarkably similar. Therefore, I would go with a linear model, as that is the simplest of the models.

cb_temp <- cherry_blossoms %>%
  drop_na(doy, temp) %>%
  mutate(temp_c = temp - mean(temp),
         temp_s = temp_c / sd(temp),
         temp_s2 = temp_s ^ 2,
         temp_s3 = temp_s ^ 3)

# linear model
lin_mod <- brm(doy ~ 1 + temp_c, data = cb_temp, family = gaussian,
               prior = c(prior(normal(100, 10), class = Intercept),
                         prior(normal(0, 10), class = b),
                         prior(exponential(1), class = sigma)),
               iter = 4000, warmup = 2000, chains = 4, cores = 4, seed = 1234,
               file = here("fits", "chp4", "b4h5-linear"))

# quadratic model
qad_mod <- brm(doy ~ 1 + temp_s + temp_s2, data = cb_temp, family = gaussian,
               prior = c(prior(normal(100, 10), class = Intercept),
                         prior(normal(0, 10), class = b, coef = "temp_s"),
                         prior(normal(0, 1), class = b, coef = "temp_s2"),
                         prior(exponential(1), class = sigma)),
               iter = 4000, warmup = 2000, chains = 4, cores = 4, seed = 1234,
               file = here("fits", "chp4", "b4h5-quadratic"))

# cubic model
cub_mod <- brm(doy ~ 1 + temp_s + temp_s2 + temp_s3, data = cb_temp,
               family = gaussian,
               prior = c(prior(normal(100, 10), class = Intercept),
                         prior(normal(0, 10), class = b, coef = "temp_s"),
                         prior(normal(0, 1), class = b, coef = "temp_s2"),
                         prior(normal(0, 1), class = b, coef = "temp_s3"),
                         prior(exponential(1), class = sigma)),
               iter = 4000, warmup = 2000, chains = 4, cores = 4, seed = 1234,
               file = here("fits", "chp4", "b4h5-cubic"))

# spline model
spl_mod <- brm(doy ~ 1 + s(temp, bs = "bs"), data = cb_temp, family = gaussian,
               prior = c(prior(normal(100, 10), class = Intercept),
                         prior(normal(0, 10), class = b),
                         prior(normal(0, 10), class = sds),
                         prior(exponential(1), class = sigma)),
               iter = 10000, warmup = 8000, chains = 4, cores = 4, seed = 1234,
               control = list(adapt_delta = 0.99),
               file = here("fits", "chp4", "b4h5-spline"))

Now let’s visualize the predictions from each model.

grid <- cb_temp %>%
  data_grid(temp = seq_range(temp, 100)) %>%
  mutate(temp_c = temp - mean(cb_temp$temp),
         temp_s = temp_c / sd(cb_temp$temp),
         temp_s2 = temp_s ^ 2,
         temp_s3 = temp_s ^ 3)

fits <- bind_rows(
  add_fitted_draws(grid, lin_mod) %>%
    mean_qi(.width = c(0.67, 0.89, 0.97)) %>%
    mutate(model = "Linear"),
  add_fitted_draws(grid, qad_mod) %>%
    mean_qi(.width = c(0.67, 0.89, 0.97)) %>%
    mutate(model = "Quadratic"),
  add_fitted_draws(grid, cub_mod) %>%
    mean_qi(.width = c(0.67, 0.89, 0.97)) %>%
    mutate(model = "Cubic"),
  add_fitted_draws(grid, spl_mod) %>%
    mean_qi(.width = c(0.67, 0.89, 0.97)) %>%
    mutate(model = "Spline")
) %>%
  ungroup() %>%
  mutate(model = factor(model, levels = c("Linear", "Quadratic", "Cubic",
                                          "Spline")))

preds <- bind_rows(
  add_predicted_draws(grid, lin_mod) %>%
    mean_qi(.width = 0.89) %>%
    mutate(model = "Linear"),
  add_predicted_draws(grid, lin_mod) %>%
    mean_qi(.width = 0.89) %>%
    mutate(model = "Quadratic"),
  add_predicted_draws(grid, lin_mod) %>%
    mean_qi(.width = 0.89) %>%
    mutate(model = "Cubic"),
  add_predicted_draws(grid, lin_mod) %>%
    mean_qi(.width = 0.89) %>%
    mutate(model = "Spline")
) %>%
  ungroup() %>%
  mutate(model = factor(model, levels = c("Linear", "Quadratic", "Cubic",
                                          "Spline")))

ggplot(cb_temp, aes(x = temp)) +
  facet_wrap(~model, nrow = 2) +
  geom_point(aes(y = doy), alpha = 0.2) +
  geom_ribbon(data = preds, aes(ymin = .lower, ymax = .upper),
              alpha = 0.2) +
  geom_lineribbon(data = fits, aes(y = .value, ymin = .lower, ymax = .upper),
                  size = .6) +
  scale_fill_brewer(palette = "Blues", breaks = c(0.67, 0.89, 0.97)) +
  labs(x = "March Temperature", y = "Day in Year") +
  theme(legend.position = "bottom")

4H6. Simulate the prior predictive distribution for the cherry blossom spline in the chapter. Adjust the prior on the weights and observe what happens. What do you think the prior on the weights is doing?

As a reminder, here is the cherry blossom spline model from the chapter:

\[\begin{align} D_i &\sim \text{Normal}(\mu_i,\sigma) \\ \mu_i &= \alpha + \sum_{k=1}^Kw_kB_{k,i} \\ \alpha &\sim \text{Normal}(100, 10) \\ w_k &\sim \text{Normal}(0,10) \\ \sigma &\sim \text{Exponential}(1) \end{align}\]

We’ll also need to recreate the basis functions that that model uses:

cb_dat <- cherry_blossoms %>%
  drop_na(doy)

num_knots <- 15
knot_list <- quantile(cb_dat$year, probs = seq(0, 1, length.out = num_knots))

B <- bs(cb_dat$year,
        knots = knot_list[-c(1, num_knots)],
        degree = 3, intercept = TRUE)

Finally, we can generate data from the priors, and combine those parameters with the basis functions to get the prior predictive distributions.

n <- 1000
tibble(.draw = seq_len(n),
       alpha = rnorm(n, 100, 10),
       w = purrr::map(seq_len(n),
                      function(x, knots) {
                        w <- rnorm(n = knots + 2, 0, 10)
                        return(w)
                      },
                      knots = num_knots)) %>%
  mutate(mu = map2(alpha, w,
                   function(alpha, w, b) {
                     res <- b %*% w
                     res <- res + alpha
                     res <- res %>%
                       as_tibble(.name_repair = ~".value") %>%
                       mutate(year = cb_dat$year, .before = 1)
                     return(res)
                   },
                   b = B)) %>%
  unnest(cols = mu) %>%
  ggplot(aes(x = year, y = .value)) +
  geom_line(aes(group = .draw), alpha = 0.1) +
  geom_vline(xintercept = knot_list, color = palette_OkabeIto[2]) +
  labs(x = "Year", y = "Day in Year")

Now let’s tighten the prior on w to \(\text{Normal}(0,2)\), as we used for exercise 4M8. Now the lines are much less wiggly, which is consistent with what we found in the previous exercise, which used the observed data.

n <- 1000
tibble(.draw = seq_len(n),
       alpha = rnorm(n, 100, 10),
       w = purrr::map(seq_len(n),
                      function(x, knots) {
                        w <- rnorm(n = knots + 2, 0, 2)
                        return(w)
                      },
                      knots = num_knots)) %>%
  mutate(mu = map2(alpha, w,
                   function(alpha, w, b) {
                     res <- b %*% w
                     res <- res + alpha
                     res <- res %>%
                       as_tibble(.name_repair = ~".value") %>%
                       mutate(year = cb_dat$year, .before = 1)
                     return(res)
                   },
                   b = B)) %>%
  unnest(cols = mu) %>%
  ggplot(aes(x = year, y = .value)) +
  geom_line(aes(group = .draw), alpha = 0.1) +
  geom_vline(xintercept = knot_list, color = palette_OkabeIto[2]) +
  labs(x = "Year", y = "Day in Year")

4H8. (sic; there is no 4H7 in the text, so I’ve kept this labeled 4H8 to be consistent with the book) The cherry blossom spline in the chapter used an intercept \(\alpha\), but technically it doesn’t require one. The first basis functions could substitute for the intercept. Try refitting the cherry blossom spline without the intercept. What else about the model do you need to change to make this work?

We can remove the intercept by removing the a parameter from the model.

no_int <- quap(alist(D ~ dnorm(mu, sigma),
                     mu <- B %*% w,
                     w ~ dnorm(0, 10),
                     sigma ~ dexp(1)),
               data = list(D = cb_dat$doy, B = B),
               start = list(w = rep(0, ncol(B))), control = list(maxit = 5000))

fitted_draws(no_int, cb_dat) %>%
  mean_qi(.width = 0.89) %>%
  ggplot(aes(x = year, y = doy)) +
  geom_point(alpha = 0.2) +
  geom_hline(yintercept = mean(cb_dat$doy), linetype = "dashed") +
  geom_lineribbon(aes(y = .value, ymin = .lower, ymax = .upper),
                  alpha = 0.5, fill = "grey20") +
  labs(x = "Year", y = "Day in Year")

This looks a lot like our original model, except the left hand side of the spline is pulled down. This is likely due to the prior on w. The prior is centered on 0, but that assumes an intercept is present (i.e., the curves of the spline average a deviation of 0 from the mean). However, with the intercept, the prior drags the line down to actual zero when the first basis function in non-zero. By changing the prior of w to the same prior originally used for the intercept, we get almost our original model back.

no_int2 <- quap(alist(D ~ dnorm(mu, sigma),
                      mu <- B %*% w,
                      w ~ dnorm(100, 10),
                      sigma ~ dexp(1)),
                data = list(D = cb_dat$doy, B = B),
                start = list(w = rep(0, ncol(B))), control = list(maxit = 5000))

fitted_draws(no_int2, cb_dat) %>%
  mean_qi(.width = 0.89) %>%
  ggplot(aes(x = year, y = doy)) +
  geom_point(alpha = 0.2) +
  geom_hline(yintercept = mean(cb_dat$doy), linetype = "dashed") +
  geom_lineribbon(aes(y = .value, ymin = .lower, ymax = .upper),
                  alpha = 0.5, fill = "grey20") +
  labs(x = "Year", y = "Day in Year")

Session Info

View the session information used to render this week.

devtools::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 4.0.2 (2020-06-22)
#>  os       macOS Catalina 10.15.6      
#>  system   x86_64, darwin17.0          
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  ctype    en_US.UTF-8                 
#>  tz       UTC                         
#>  date     2020-09-19                  
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  ! package              * version     date       lib
#>  P abind                  1.4-5       2016-07-21 [?]
#>  P arrayhelpers           1.1-0       2020-02-04 [?]
#>  P assertthat             0.2.1       2019-03-21 [?]
#>  P backports              1.1.9       2020-08-24 [?]
#>  P base64enc              0.1-3       2015-07-28 [?]
#>  P bayesplot              1.7.2       2020-05-28 [?]
#>  P blob                   1.2.1       2020-01-20 [?]
#>  P bookdown               0.20.3      2020-09-19 [?]
#>  P bridgesampling         1.0-0       2020-02-26 [?]
#>  P brms                 * 2.13.9      2020-09-19 [?]
#>  P Brobdingnag            1.2-6       2018-08-13 [?]
#>  P broom                  0.7.0       2020-07-09 [?]
#>  P callr                  3.4.3       2020-03-28 [?]
#>  P cellranger             1.1.0       2016-07-27 [?]
#>  P checkmate              2.0.0       2020-02-06 [?]
#>  P class                  7.3-17      2020-04-26 [?]
#>  P classInt               0.4-3       2020-04-07 [?]
#>  P cli                    2.0.2       2020-02-28 [?]
#>  P coda                   0.19-3      2019-07-05 [?]
#>  P codetools              0.2-16      2018-12-24 [?]
#>  P colorblindr          * 0.1.0       2020-09-19 [?]
#>  P colorspace           * 1.4-1       2019-03-18 [?]
#>  P colourpicker           1.0         2017-09-27 [?]
#>  P cpp11                  0.2.1       2020-08-11 [?]
#>  P crayon                 1.3.4       2017-09-16 [?]
#>  P crosstalk              1.1.0.1     2020-03-13 [?]
#>  P curl                   4.3         2019-12-02 [?]
#>  P DBI                    1.1.0       2019-12-15 [?]
#>  P dbplyr                 1.4.4       2020-05-27 [?]
#>  P desc                   1.2.0       2018-05-01 [?]
#>  P devtools               2.3.1.9000  2020-09-19 [?]
#>  P digest                 0.6.25      2020-02-23 [?]
#>  P distributional         0.2.0       2020-08-03 [?]
#>  P dplyr                * 1.0.2       2020-08-18 [?]
#>  P DT                     0.15        2020-08-05 [?]
#>  P dygraphs               1.1.1.6     2018-07-11 [?]
#>  P e1071                  1.7-3       2019-11-26 [?]
#>  P ellipsis               0.3.1       2020-05-15 [?]
#>  P evaluate               0.14        2019-05-28 [?]
#>  P extrafont              0.17        2014-12-08 [?]
#>  P extrafontdb            1.0         2012-06-11 [?]
#>  P fansi                  0.4.1       2020-01-08 [?]
#>  P farver                 2.0.3       2020-01-16 [?]
#>  P fastmap                1.0.1       2019-10-08 [?]
#>  P forcats              * 0.5.0       2020-03-01 [?]
#>  P fs                     1.5.0       2020-07-31 [?]
#>  P gdtools                0.2.2       2020-04-03 [?]
#>  P generics               0.0.2       2018-11-29 [?]
#>  P gganimate            * 1.0.6       2020-07-08 [?]
#>  P ggdist                 2.2.0       2020-07-12 [?]
#>  P ggplot2              * 3.3.2       2020-06-19 [?]
#>  P ggridges               0.5.2       2020-01-12 [?]
#>  P gifski                 0.8.6       2018-09-28 [?]
#>  P glue                 * 1.4.2       2020-08-27 [?]
#>  P gridExtra              2.3         2017-09-09 [?]
#>  P gtable                 0.3.0       2019-03-25 [?]
#>  P gtools                 3.8.2       2020-03-31 [?]
#>  P haven                  2.3.1       2020-06-01 [?]
#>  P HDInterval             0.2.2       2020-05-23 [?]
#>  P here                 * 0.1         2017-05-28 [?]
#>  P highr                  0.8         2019-03-20 [?]
#>  P hms                    0.5.3       2020-01-08 [?]
#>  P hrbrthemes           * 0.8.0       2020-03-06 [?]
#>  P htmltools              0.5.0       2020-06-16 [?]
#>  P htmlwidgets            1.5.1       2019-10-08 [?]
#>  P httpuv                 1.5.4       2020-06-06 [?]
#>  P httr                   1.4.2       2020-07-20 [?]
#>  P igraph                 1.2.5       2020-03-19 [?]
#>  P inline                 0.3.15      2018-05-18 [?]
#>  P jsonlite               1.7.0       2020-06-25 [?]
#>  P kableExtra           * 1.2.1       2020-08-27 [?]
#>  P KernSmooth             2.23-17     2020-04-26 [?]
#>  P knitr                  1.29        2020-06-23 [?]
#>  P labeling               0.3         2014-08-23 [?]
#>  P later                  1.1.0.1     2020-06-05 [?]
#>  P lattice                0.20-41     2020-04-02 [?]
#>  P lifecycle              0.2.0       2020-03-06 [?]
#>  P loo                  * 2.3.1       2020-07-14 [?]
#>  P lpSolve                5.6.15      2020-01-24 [?]
#>  P lubridate              1.7.9       2020-06-08 [?]
#>  P magrittr               1.5         2014-11-22 [?]
#>  P markdown               1.1         2019-08-07 [?]
#>  P MASS                   7.3-51.6    2020-04-26 [?]
#>  P Matrix                 1.2-18      2019-11-27 [?]
#>  P matrixStats            0.56.0      2020-03-13 [?]
#>  P memoise                1.1.0       2017-04-21 [?]
#>  P mgcv                   1.8-31      2019-11-09 [?]
#>  P mime                   0.9         2020-02-04 [?]
#>  P miniUI                 0.1.1.1     2018-05-18 [?]
#>  P modelr               * 0.1.8       2020-05-19 [?]
#>  P munsell                0.5.0       2018-06-12 [?]
#>  P mvtnorm                1.1-1       2020-06-09 [?]
#>  P nlme                   3.1-148     2020-05-24 [?]
#>  P pillar                 1.4.6       2020-07-10 [?]
#>  P pkgbuild               1.1.0       2020-07-13 [?]
#>  P pkgconfig              2.0.3       2019-09-22 [?]
#>  P pkgload                1.1.0       2020-05-29 [?]
#>  P plyr                   1.8.6       2020-03-03 [?]
#>  P prettyunits            1.1.1       2020-01-24 [?]
#>  P processx               3.4.3       2020-07-05 [?]
#>  P progress               1.2.2       2019-05-16 [?]
#>  P promises               1.1.1       2020-06-09 [?]
#>  P ps                     1.3.4       2020-08-11 [?]
#>  P purrr                * 0.3.4       2020-04-17 [?]
#>  P R6                     2.4.1       2019-11-12 [?]
#>  P ragg                   0.3.1       2020-07-03 [?]
#>  P ratlas               * 0.0.0.9000  2020-09-19 [?]
#>  P RColorBrewer           1.1-2       2014-12-07 [?]
#>  P Rcpp                 * 1.0.5       2020-07-06 [?]
#>  P RcppParallel           5.0.2       2020-06-24 [?]
#>  P readr                * 1.3.1       2018-12-21 [?]
#>  P readxl                 1.3.1       2019-03-13 [?]
#>  P remotes                2.2.0       2020-07-21 [?]
#>    renv                   0.11.0      2020-06-26 [1]
#>  P reprex                 0.3.0       2019-05-16 [?]
#>  P reshape2               1.4.4       2020-04-09 [?]
#>  P rethinking           * 2.13        2020-09-19 [?]
#>  P rlang                  0.4.7       2020-07-09 [?]
#>  P rmarkdown              2.3.3       2020-09-19 [?]
#>  P rprojroot              1.3-2       2018-01-03 [?]
#>  P rsconnect              0.8.16      2019-12-13 [?]
#>  P rstan                * 2.21.2      2020-07-27 [?]
#>  P rstantools             2.1.1       2020-07-06 [?]
#>  P rstudioapi             0.11        2020-02-07 [?]
#>  P Rttf2pt1               1.3.8       2020-01-10 [?]
#>  P rvest                  0.3.6       2020-07-25 [?]
#>  P scales                 1.1.1       2020-05-11 [?]
#>  P sessioninfo            1.1.1       2018-11-05 [?]
#>  P sf                     0.9-5       2020-07-14 [?]
#>  P shape                  1.4.4       2018-02-07 [?]
#>  P shiny                  1.5.0       2020-06-23 [?]
#>  P shinyjs                1.1         2020-01-13 [?]
#>  P shinystan              2.5.0       2018-05-01 [?]
#>  P shinythemes            1.1.2       2018-11-06 [?]
#>  P StanHeaders          * 2.21.0-6    2020-08-16 [?]
#>  P stringi                1.5.3       2020-09-09 [?]
#>  P stringr              * 1.4.0       2019-02-10 [?]
#>  P svUnit                 1.0.3       2020-04-20 [?]
#>  P systemfonts            0.3.0       2020-09-01 [?]
#>  P testthat               2.99.0.9000 2020-09-19 [?]
#>  P threejs                0.3.3       2020-01-21 [?]
#>  P tibble               * 3.0.3       2020-07-10 [?]
#>  P tidybayes            * 2.1.1.9000  2020-09-19 [?]
#>  P tidybayes.rethinking * 2.0.3.9000  2020-09-19 [?]
#>  P tidyr                * 1.1.2       2020-08-27 [?]
#>  P tidyselect             1.1.0       2020-05-11 [?]
#>  P tidyverse            * 1.3.0       2019-11-21 [?]
#>  P transformr             0.1.3       2020-07-05 [?]
#>  P tweenr                 1.0.1       2018-12-14 [?]
#>  P units                  0.6-7       2020-06-13 [?]
#>  P usethis                1.6.1       2020-04-29 [?]
#>  P V8                     3.2.0       2020-06-19 [?]
#>  P vctrs                  0.3.4       2020-08-29 [?]
#>  P viridisLite            0.3.0       2018-02-01 [?]
#>  P webshot                0.5.2       2019-11-22 [?]
#>  P withr                  2.2.0       2020-04-20 [?]
#>  P xfun                   0.16        2020-07-24 [?]
#>  P xml2                   1.3.2       2020-04-23 [?]
#>  P xtable                 1.8-4       2019-04-21 [?]
#>  P xts                    0.12-0      2020-01-19 [?]
#>  P yaml                   2.2.1       2020-02-01 [?]
#>  P zoo                    1.8-8       2020-05-02 [?]
#>  source                                      
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  Github (rstudio/bookdown@9eb20a8)           
#>  CRAN (R 4.0.2)                              
#>  Github (paul-buerkner/brms@c520a8c)         
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  Github (clauswilke/colorblindr@e6730be)     
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.1)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  Github (r-lib/devtools@df619ce)             
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.1)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  Github (atlas-aai/ratlas@746bf9c)           
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  Github (rmcelreath/rethinking@3b48ec8)      
#>  CRAN (R 4.0.2)                              
#>  Github (rstudio/rmarkdown@204aa41)          
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  Github (r-lib/testthat@b3f2203)             
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  Github (mjskay/tidybayes@17f4dc6)           
#>  Github (mjskay/tidybayes.rethinking@df903c8)
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.1)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#>  CRAN (R 4.0.2)                              
#> 
#> [1] /Users/runner/work/sr2-solutions/sr2-solutions/renv/library/R-4.0/x86_64-apple-darwin17.0
#> [2] /private/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T/RtmpojYdR7/renv-system-library
#> 
#>  P ── Loaded and on-disk path mismatch.