Notebook for Assignment 8

Author

Aki Vehtari et al.

1 General information

This assignment relates to Lectures 8-9 and Chapter 7 of BDA3

We recommend using JupyterHub (which has all the needed packages pre-installed).

Reading instructions:

General Instructions for Answering the Assignment Questions
  • Questions below are exact copies of the text found in the MyCourses quiz and should serve as a notebook where you can store notes and code.
  • We recommend opening these notebooks in the Aalto JupyterHub, see how to use R and RStudio remotely.
  • For inspiration for code, have a look at the BDA R Demos and the specific Assignment code notebooks
  • Recommended additional self study exercises for each chapter in BDA3 are listed in the course web page. These will help to gain deeper understanding of the topic.
  • Common questions and answers regarding installation and technical problems can be found in Frequently Asked Questions (FAQ).
  • Deadlines for all assignments can be found on the course web page and in MyCourses.
  • You are allowed to discuss assignments with your friends, but it is not allowed to copy solutions directly from other students or from internet.
  • Do not share your answers publicly.
  • Do not copy answers from the internet or from previous years. We compare the answers to the answers from previous years and to the answers from other students this year.
  • Use of AI is allowed on the course, but the most of the work needs to by the student, and you need to report whether you used AI and in which way you used them (See points 5 and 6 in Aalto guidelines for use of AI in teaching).
  • All suspected plagiarism will be reported and investigated. See more about the Aalto University Code of Academic Integrity and Handling Violations Thereof.
  • If you have any suggestions or improvements to the course material, please post in the course chat feedback channel, create an issue, or submit a pull request to the public repository!
if (!require(brms)) {
    install.packages("brms")
    library(brms)
}
Loading required package: brms
Loading required package: Rcpp
Loading 'brms' package (version 2.22.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').

Attaching package: 'brms'
The following object is masked from 'package:stats':

    ar
if(!require(cmdstanr)){
    install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
    library(cmdstanr)
}
Loading required package: cmdstanr
This is cmdstanr version 0.8.1
- CmdStanR documentation and vignettes: mc-stan.org/cmdstanr
- CmdStan path: /home/runner/.cmdstan/cmdstan-2.35.0
- CmdStan version: 2.35.0
cmdstan_installed <- function(){
  res <- try(out <- cmdstanr::cmdstan_path(), silent = TRUE)
  !inherits(res, "try-error")
}
if(!cmdstan_installed()){
    install_cmdstan()
}

2 Sleep deprivation

The dataset sleepstudy is available by using the command data(sleepstudy, package = "lme4")

Below is some code for fitting a brms model. This model is a simple pooled model. You will need to fit a hierarchical model as explained in the assignment, but this code should help getting started.

Load the dataset

data(sleepstudy, package = "lme4")
Error in find.package(package, lib.loc, verbose = verbose): there is no package called 'lme4'

Specify the formula and observation family:

sleepstudy_pooled_formula <- bf(
  Reaction ~ 1 + Days,
  family = "gaussian",
  center = FALSE
)

We can then specify the priors:

(sleepstudy_pooled_priors <- c(
  prior(
    normal(250, 100),
    class = "b",
    coef = "Intercept"
  ),
  prior(
    normal(0, 20),
    class = "b",
    coef = "Days"
  ),
  prior(
    normal(0, 100),
    class = "sigma"
  )
))
            prior class      coef group resp dpar nlpar   lb   ub source
 normal(250, 100)     b Intercept                       <NA> <NA>   user
    normal(0, 20)     b      Days                       <NA> <NA>   user
   normal(0, 100) sigma                                 <NA> <NA>   user

And then fit the model:

sleepstudy_pooled_fit <- brm(
  formula = sleepstudy_pooled_formula,
  prior = sleepstudy_pooled_priors,
  data = sleepstudy
)
Error: object 'sleepstudy' not found

We can add the leave-one-out CV criterion to the fit object, and then view the output

sleepstudy_pooled_fit <- add_criterion(
  sleepstudy_pooled_fit,
  criterion = "loo"
)
Error: object 'sleepstudy_pooled_fit' not found
loo(sleepstudy_pooled_fit)
Error: object 'sleepstudy_pooled_fit' not found

3 Sleep deprivation extension

Here is an example of a pooled model with spline

sleepstudy_pooled_spline_formula <- bf(
  Reaction ~ 1 + Days + s(Days),
  family = "gaussian",
  center = FALSE
)


sleepstudy_pooled_spline_priors <- c(
  prior(
    normal(250, 100),
    class = "b",
    coef = "Intercept"
  ),
  prior(
    normal(0, 20),
    class = "b",
    coef = "Days"
  ),
  prior(
    normal(0, 100),
    class = "sigma"
  ),
  prior(
    normal(0, 20),
    class = "b",
    coef = "sDays_1"
  )
)

sleepstudy_pooled_spline_fit <- brm(
  formula = sleepstudy_pooled_spline_formula,
  prior = sleepstudy_pooled_spline_priors,
  data = sleepstudy,
  family = "gaussian"
)
Error: object 'sleepstudy' not found

Adding the loo criterion:

sleepstudy_pooled_spline_fit <- add_criterion(
  sleepstudy_pooled_spline_fit,
  criterion = "loo"
)
Error: object 'sleepstudy_pooled_spline_fit' not found

Compare to the other model

loo_compare(
  loo(sleepstudy_pooled_fit),
  loo(sleepstudy_pooled_spline_fit)
)
Error: object 'sleepstudy_pooled_fit' not found

If there are high Pareto-k values

First try moment matching:

sleepstudy_pooled_spline_fit <- add_criterion(
  sleepstudy_pooled_spline_fit,
  criterion = "loo",
  moment_match = TRUE,
  overwrite = TRUE
)
Error: object 'sleepstudy_pooled_spline_fit' not found

Then if there are a few remaining (and refitting is not going to take forever), try moment matching and reloo

sleepstudy_pooled_spline_fit <- add_criterion(
  sleepstudy_pooled_spline_fit,
  criterion = "loo",
  moment_match = TRUE,
  reloo = TRUE,
  overwrite = TRUE
)
Error: object 'sleepstudy_pooled_spline_fit' not found