Assignment 9

Decision analysis

1 General information

This is the template for assignment 9. You can download the qmd-file or copy the code from this rendered document after clicking on </> Code in the top right corner.

Please replace the instructions in this template by your own text, explaining what you are doing in each exercise.

Setup

This block will only be visible in your HTML output, but will be hidden when rendering to PDF with quarto for the submission. Make sure that this does not get displayed in the PDF!

The following loads several needed packages:

library(bayesplot)

This is bayesplot version 1.10.0

- Online documentation and vignettes at mc-stan.org/bayesplot

- bayesplot theme set to bayesplot::theme_default()

   * Does _not_ affect other ggplot2 plots

   * See ?bayesplot_theme_set for details on theme setting

library(cmdstanr)

This is cmdstanr version 0.5.3

- CmdStanR documentation and vignettes: mc-stan.org/cmdstanr

- CmdStan path: /root/.cmdstan/cmdstan-2.31.0

- CmdStan version: 2.31.0

A newer version of CmdStan is available. See ?install_cmdstan() to install it.
To disable this check set option or environment variable CMDSTANR_NO_VER_CHECK=TRUE.

library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(ggplot2)
library(ggdist) # for stat_dotsinterval
library(posterior)

This is posterior version 1.4.0

Attaching package: 'posterior'

The following object is masked from 'package:bayesplot':

    rhat

The following objects are masked from 'package:stats':

    mad, sd, var

The following objects are masked from 'package:base':

    %in%, match

library(brms)

Loading required package: Rcpp

Loading 'brms' package (version 2.19.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 objects are masked from 'package:ggdist':

    dstudent_t, pstudent_t, qstudent_t, rstudent_t

The following object is masked from 'package:bayesplot':

    rhat

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

    ar

# Globally specfiy cmdstan backend for brms
options(brms.backend="cmdstanr")
# Tell brms to cache results if possible
options(brms.file_refit="on_change")

# Set more readable themes with bigger font for plotting packages
ggplot2::theme_set(theme_minimal(base_size = 14))
bayesplot::bayesplot_theme_set(theme_minimal(base_size = 14))

2 Escaping from the chicken coop

2.1 (a)

A simple GP model

This block will only be visible in your HTML output, but will be hidden when rendering to PDF with quarto for the submission. Make sure that this does not get displayed in the PDF!

The below fits a GP model to the chicken growth curves. It may take a few minutes to fit, but you can also download the fit .rds-file and work with that fit object.

fit <- brm(
  weight ~ gp(Time) + (0+Time|Diet) + (0+Time|Chick),
  data = ChickWeight,
  family = "lognormal",
  file="additional_files/assignment9/gp_chicken_fit"
)
brms::pp_check(fit, type = "intervals_grouped", group = "Diet")

Using all posterior draws for ppc type 'intervals_grouped' by default.

# Useful r functions: 
#   rep(..., each=...), cbind, colMeans, 
#   posterior_predict(..., newdata=..., allow_new_levels=TRUE, sample_new_levels="gaussian")
#   ggplot, geom_line, aes(..., group=..., color=...)

Chickenwise probability of escape function

This block will only be visible in your HTML output, but will be hidden when rendering to PDF with quarto for the submission. Make sure that this does not get displayed in the PDF!

bump <- function(x, loc=0, scale=1){
  xi = (x - loc) / scale
  ifelse(abs(xi) < 1, exp(-1/(1-xi^2)), 0.)
}
daily_probability_of_escape <- function(day, weight){
  # Expects a day and a weight and computes the daily probability of escape 
  bump(day, 30, 10) * (1e-2 + bump(weight, 200, 150)+bump(weight, 700, 150))
}
chickenwise_probability_of_escape <- function(weights){
  # Expects a vector of daily weights from day 1 to N and computes the probability of 
  # escape at the end of the time series
  prob_of_failure = 1
  for(day in 1:length(weights)){
    prob_of_failure = prob_of_failure * (1-daily_probability_of_escape(day, weights[day]))
  }
  return(1 - prob_of_failure)
}

2.2 (b)

# Useful r functions: chickenwise_probability_of_escape (see above)
# rep(..., each=...), apply, 
# ggplot, stat_dotsinterval

2.3 (c)

# Useful r functions: chickenwise_probability_of_escape (see above)
# apply, aggregate,