The grading will be done in peergrade. All grading questions and evaluations for this assignment are contained within this document in the collapsible Rubric blocks.
Installing and using CmdStanR:
See the Stan demos on how to use Stan in R (or Python). Aalto JupyterHub has working R and CmdStanR/RStan environment and is probably the easiest way to use Stan. * To use CmdStanR in Aalto JupyterHub: library(cmdstanr) set_cmdstan_path('/coursedata/cmdstan')
The Aalto Ubuntu desktops also have the necessary libraries installed.]{.aalto}
To install Stan on your laptop, run ‘install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))’ in R. If you encounter problems, see additional answers in FAQ. If you don’t succeed in short amount of time, it is probably easier to use Aalto JupyterHub.
If you use Aalto JupyterHub, all necessary packages have been pre-installed. In your laptop, install package cmdstanr. Installation instructions on Linux, Mac and Windows can be found at https://mc-stan.org/cmdstanr/. Additional useful packages are loo, bayesplot and posterior (but you don’t need these in this assignment). For Python users, PyStan, CmdStanPy, and ArviZ packages are useful.
If you edit files ending .stan in RStudio, you can click “Check” in the editor toolbar to make syntax check. This can significantly speed-up writing a working Stan model.
Reporting accuracy
For posterior statistics of interest, only report digits that are not completely random based on the Monte Carlo standard error (MCSE).
Example: If you estimate \(E(\mu) \approx 1.234\) with MCSE(\(E(\mu)\)) = 0.01, then the true expectation is likely to be between \(1.204\) and \(1.264\), it makes sense to report \(E(\mu) \approx 1.2\).
If you want to install R and RStudio locally, download R and RStudio.
There are tons of tutorials, videos and introductions to R and RStudio online. You can find some initial hints from RStudio Education pages.
When working with R, we recommend writing the report using quarto and the provided template. The template includes the formatting instructions and how to include code and figures.
Instead of quarto, you can use other software to make the PDF report, but the the same instructions for formatting should be used.
Report all results in a single, anonymous *.pdf -file and submit it in peergrade.io.
The course has its own R package aaltobda with data and functionality to simplify coding. The package is pre-installed in JupyterHub. To install the package on your own system, run the following code (upgrade="never" skips question about updating other packages):
Many of the exercises can be checked automatically using the R package markmyassignment (pre-installed in JupyterHub). Information on how to install and use the package can be found in the markmyassignment documentation. There is no need to include markmyassignment results in the report.
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.
Deadlines for all assignments can be found on the course web page and in Peergrade. You can set email alerts for the deadlines in Peergrade settings.
You are allowed to discuss assignments with your friends, but it is not allowed to copy solutions directly from other students or from internet.
You can copy, e.g., plotting code from the course demos, but really try to solve the actual assignment problems with your own code and explanations.
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.
Do not submit empty PDFs, almost empty PDFs, copy of the questions, nonsense generated by yourself or AI, as these are just harming the other students as they can’t do peergrading for the empty or nonsense submissions. Violations of this rule will be reported and investigated in the same way was plagiarism.
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!
Rubric S1: General information - 7.5/100 points
Q1:
Can you open the PDF and it’s not blank nor nonsense? If the pdf is blank, nonsense, or something like only a copy of the questions, 1) report it as problematic in Peergrade-interface to get another report to review, and 2) send a message to TAs.
Q2:
Is the report anonymous?
Setup
The following loads several needed packages:
library(aaltobda)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_dotsintervallibrary(posterior)
This is posterior version 1.4.0
Attaching package: 'posterior'
The following object is masked from 'package:bayesplot':
rhat
The following object is masked from 'package:aaltobda':
mcse_quantile
The following objects are masked from 'package:stats':
mad, sd, var
The following objects are masked from 'package:base':
%in%, match
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
# 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))# This registers CmdStan as the backend for compiling cmdstan-chunks.check_cmdstan_toolchain(fix =TRUE, quiet =TRUE)register_knitr_engine(override =FALSE)
2 Hierarchical Model: Chicken Data with Stan (6p)
Do not look at the secret dataset until you have defined your priors in section 2.1.
A secret dataset, included in base r, contains weight measurements from 50 chicks from their birth until the age of 21 days. For this task, we are interested in modeling the weight of a chick at the age of 12 days.
2.1 Choosing a weakly informative prior by intuition
We will first guide you through two processes for choosing weakly informative priors which you can use in your projects as well (though, if expert knowledge is present, we encourage to use that though instead). Many definitions of weak exist, but for our purposes, we intend to set priors in this assignment which don’t overwhelm the likelihood contributions to the posterior while still exerting some amount of regularisation.
In the absence of additional information, we will assume that chick weights at 12 days of age (\(w\)) follow a normal distribution: \[
w\sim\mathcal N (\mu,\sigma).
\] Our task will be to define weakly informative priors for the mean \(\mu\) given observation model standard deviation, \(\sigma\) (a prior distribution for \(\sigma\) will be set in after section 2).
Here, \(\mu\) represents the population mean chick weight at 12 days and \(\sigma\) represents the standard-deviation of the chick weights \(w\) around that population mean. For this exercise, we will be specifying a normal prior for \(\mu\):
\[
\mu\sim\mathcal N (\mu_0,\sigma_0).
\]
\(\mu_0\) represents our prior knowledge of what we believe the population weight to be and \(\sigma_0\) represents our level of certainty in this belief. To specify a weakly-informative prior for \(\mu\) we need to select values of \(\mu_0\) and \(\sigma_0\) that imply the range of plausible values that \(\mu\) could take.
We can do this by our own intuitive prior knowledge (if such intuition exists), or by searching for external references.
Despite the name, weakly informed priors can be quite subjective (see for more theoretical discussion here), such that some justification is always needed. For the subtasks below, you will not be graded on the accuracy or precision of your numbers, but on your justification of them. The numerical choices you make should make sense and be understandable to an external reviewer of your work (even if they may not agree with your choices).
A word of caution on eliciting the priors below
Please note that in the below, we intend to set a prior on \(\mu\) (the mean chick weight), but the intuition we ilicit is based on the weight of individual chicks. We do so to help create intuition about what the mean could be, however, it would be theoretically more accurate to ilicit priors about mean chick weights.
Important
We have made changes to the assignment text and some of the rubrics to make it clearer.
Subtask 2.a)
Based on your own past experience and estimation skills, what would you guess is a typical weight range of a fully grown chicken in grams? Justify your choice with 1-3 sentences.
Tip
Would it make sense if a chicken weighed 1 million grams? 0.005 grams? What about a chicken that weighed more than you? More than a car? Note that is not important to have a very precise or accurate guess here, the key goal is simply to identify the range of plausible (or at least possible) values.
Subrubric 2.a)
Q3:
Does the chosen range meet the following common sense criteria?
The range is always above … and below ….
Q4:
Is the justification based on some sort of logic, even if you may disagree?
Subtask 2.b)
Adjust this range for a 12-day old chick and choose a mean \(\mu_0\) for the weakly informed prior of the parameter \(\mu.\) Justify your adjustment and choice with 1-3 sentences.
Subrubric 2.b)
Q5:
The range is always above … and below ….
Q6:
Is the justification based on some sort of logic, even if you may disagree?
Choosing the prior standard deviation for \(\mu\) requires a little more caution; overconfident (i.e. narrow) priors can have a strong effect on your results, whereas less confident priors are more easily overcome with observed data. Given that overconfidence is a common human bias, a good intuition-based standard deviation should focus on eliminating the impossible values, rather than including the most likely values.
Subtask 2.c)
Choose a conservative lower and upper bound for the weight of any 12-day old chick, with the goal to exclude impossible values. Justify your choice with 1-3 sentences.
Tip
Depending on your choice of range above for a “typical” chick, this range excluding the impossible values should be wider.
Subrubric 2.c)
Q7:
Does the chosen range meet the following common sense criteria?
The range is always above … and below ….
Q8:
Is the justification based on some sort of logic, even if you may disagree?
A common technique to find a weakly informative prior is to have a standard deviation which is an order of magnitude (a factor of 10) larger than a plausible standard deviation of the data.
Subtask 2.d)
What do you think is a plausible standard deviation \(\sigma_\text{plausible}\) of the weight of 12 days-old chicks, based on your ranges stated above? Under the above recommendation, what standard deviation \(\sigma_0\) should you use for your prior for the mean weight \(\mu\)?
Subrubric 2.d)
Q9:
Do the two standard deviations meet the following common sense criteria?
Both values are …
Q10:
Given the choice of mean from above, the interval \((\mu_0 - 3\sigma_0, \mu_0 + 3\sigma_0)\) includes ….
Subtask 2.e)
Write down in mathematical form the final prior for the mean weight \(\mu\) you found using this prior definition technique.
Subrubric 2.e)
Q11:
Does the final prior exist in mathematical notation?
Q12:
Does the prior reflect the choices made above (i.e. \(\mu_0, \sigma_0\))?
2.2 Choosing a weakly informative prior using external references
Next, we’ll use external references to pick the weakly informed prior. This technique is more general and doesn’t assume you would have prior knowledge.
Subtask 2.f)
Consult a trustworthy source on the weight range of farm chickens, e.g. books, articles, a farmer friend. If the recommended values are for a fully grown chicken, make a reasonable adjustment for a 12-day old chick. What is the weight range of a 12-day old chicken? Cite your source and justify any adjustments you make to the reference range.
Subrubric 2.f)
Q13:
Does the reference range meet the following common sense criteria?
The range is always above … and below ….
Q14:
Is there a citation? If an adjustment was made, was it justified by logic, even if you disagree?
Subtask 2.g)
Based on this reference range, what will you choose for the mean of our weakly informed prior? Justify your choice with 1-3 sentences.
Subrubric 2.g)
Q15:
Does the mean value meet the following common sense criteria?
The range is always above … and below ….
Next we choose the standard deviation of the prior. We could use the same technique as before, but we’ll walk you through another common approach. Assume that \(99.7\%\) of all 12-day old chicks fall into the reference range you found. Under our assumption of a normal distribution, this range will encompass values between \(\mu_0 \pm 3\sigma_0\).
Subtask 2.h)
Assuming symmetry, use the mean you chose and either the upper or lower bound \(b\) of your reference range to solve the correct version of the following equations to find your associated choice of \(\sigma_0\): (show your work)
For upper bound \(b_u\): \(Pr(\mu_0 + 3\sigma_0 < b_u) \approx 0.997\), solving for \(\sigma_0\).
For lower bound \(b_l\): \(Pr(\mu_0 - 3\sigma_0 > b_l) \approx 0.997\), solving for \(\sigma_0\).
Subrubric 2.h)
Q16:
Does the calculated \(\sigma_0\) value meet the following common sense criteria?
The value is above … and below ….
Q17:
Did they show their work?
Subtask 2.i)
Write down in mathematical form the final prior for the mean weight \(\mu\) you found using this prior definition technique.
Subrubric 2.i)
Q18:
Does the final prior exist in mathematical notation?
Q19:
Does the prior reflect the choices made above (i.e. \(\mu_0, \sigma_0\))?
2.3 Non-normal priors
The previous steps all assumed we could use a prior which is normally distributed, but this may not always be the correct assumption to make.
Subtask 2.j)
Under what mathematical/statistical circumstances would a normal distributed prior not make sense? List at least one circumstance. Are there values that the normal distribution can take on which would not make sense for some types of variables?
Tip
Consider the nature of any variable you are trying to define a prior over.
Subrubric 2.j)
Q20:
Example cases include variables that … or variables that …
2.4 Modeling diet effects on chicken weight
In addition to chick weights, the data also contains a categorical variable indicating which diet the chick received. In the data file, each column contains the measurements for a single chick at a given point of time.
In addition to the existing diets, we are interested in the quality of another box of feed (the fifth diet), which a farmer happened to find yesterday at a dark corner of his warehouse. To read in the data and select chicks with age of 12 days, and to have a peek at the first 6 rows, just use:
Data inside, don’t peek before you have set your priors!
In the following analysis, we will model the weight of a chick at the age of 12 days. We will use the following three Gaussian models:
a separate model, in which each diet is modeled individually
a pooled model, in which all measurements are combined and there is no distinction between diets
a hierarchical model, which has a hierarchical structure as described in BDA3 Section 11.6
As in the model described in the book, use the same weight standard deviation \(\sigma\) for all the groups in the hierarchical model. In the separate model, however, use separate weight standard deviation \(\sigma_d\) for each diet \(d\). You should use weakly informative priors for all your models.
The provided Stan code below is given as an example of a separate model. Note that the author has left a comment expressing uncertainty about their prior choices. This separate model can be summarized mathematically as \[
\begin{aligned}
\mu_{d} &\sim \pi(\mu_d)&&\text{(diet-wise mean weight),}\\
\sigma_d &\sim \pi(\sigma_d)&&\text{(diet-wise standard deviation of diet-wise chicken weights)},\\
w_{i,d} &\sim N(\mu_d,\sigma_d)&&\text{(diet-wise chicken weights)}\\
\end{aligned}
\tag{1}\]
For the separate and the pooled models, use one of the weakly informative priors you have derived in 2.1) or 2.2) for the diet-wise mean weights \[\mu_d\sim \pi(\mu_d)=N(\mu_0,\sigma_0).\] For the hierarchical model, remember that the parameters in the priors for the diet-wise mean weights itself have to be parameters with their own prior distributions, in our case \[
\begin{aligned}
\mu_d &\sim N(\mu,\tau)&&\text{(diet-wise mean weights)},\\
\mu &\sim \pi(\mu)&&\text{(mean of prior for diet-wise mean weights) and}\\
\tau &\sim \pi(\tau)&&\text{(standard deviation of prior for diet-wise mean weights).}
\end{aligned}
\]
Use the prior you have used for the diet-wise mean weights \(\mu_d\) in the separate and pooled models for the prior on the mean of the prior for the diet-wise mean weights \[\mu \sim \pi(\mu) = N(\mu_0,\sigma_0)\] in the hierarchical model and use \[\tau \sim \pi(\tau) = \mathrm{exponential}(.02)\] for the prior on the standard deviation of the prior for the diet-wise mean weights.
data {int<lower=0> N_observations;int<lower=0> N_diets;array[N_observations] int diet_idx; // Pair observations to their diets.vector[N_observations] weight;}parameters {// Average weight of chicks with a given diet.vector[N_diets] mean_diet;// Standard deviation of weights observed among chicks sharing a diet.vector<lower=0>[N_diets] sd_diet;}model {// Priors// These look bad. I need to think about these again.for (diet in1:N_diets) { mean_diet[diet] ~ normal(0, 10); sd_diet[diet] ~ exponential(.02); }// Likelihoodfor (obs in1:N_observations) { weight[obs] ~ normal(mean_diet[diet_idx[obs]], sd_diet[diet_idx[obs]]); }// Best practice would be to write the likelihood without the for loop as:// weight ~ normal(mean_diet[diet_idx], sd_diet[diet_idx]);}generated quantities {real weight_pred;real mean_five;// The below is just there to make the plotting in the template work with the "wrong model". real sd_diets = sd_diet[4];// Sample from the (posterior) predictive distribution of the fourth diet. weight_pred = normal_rng(mean_diet[4], sd_diet[4]);// Construct samples of the mean of the fifth diet.// We only have the prior... mean_five = normal_rng(0, 10);}
Subtask 2.k)
Describe the models with mathematical notation (as is done for the separate model above). Also describe in words the difference between the model and the other models.
Subrubric 2.k)
Q21:
Are the models described using mathematical notation and the difference to other models described in words?
No equations and no description
Description but no equations
Equations but no description
Equations and description
Subtask 2.l)
Implement the models in Stan and include the code in the report. Use weakly informative priors for all of your models.
Tip
When sampling from the posterior of the hierarchical model, you will very likely get a warning about divergent transitions. This tells you that the sampler is having difficulties in sampling the joint posterior and this may lead to biased inference. We will return to this in the next task.
Subrubric 2.l)
Q22:
Is there a related Stan implementation?
No Stan model implemented
Stan model implemented, but it seems clearly wrong or broken
Seemingly valid Stan model implemented
Subtask 2.m)
Use the provided code in the template to plot the posterior distribution of the mean of the weight measurements of the fourth diet and comment on the possible differences you observe between the models.
Subrubric 2.m)
Q23:
Is there a comparison plotted for the posteriors of the mean of diet 4? Does it look something like the model solution plot?
No comparison plotted
Comparison plotted but it clearly differs from the example
Comparison plotted and it approximately matches the example
Q24:
Separate model: Is the result for the separate model discussed?
The result is not discussed.
There is some discussion, but it is not mentioned that ….
The result is discussed and the separate model is recognised as ….
Q25:
Pooled model: Is the result for the pooled model discussed?
The result is not discussed.
There is some discussion, but it is not mentioned that ….
The result is discussed and the pooled model is recognised as ….
Q26:
Hierarchical model: Is the result for the hierarchical model discussed?
The result is not discussed.
There is some discussion, but it is not mentioned that ….
The result is discussed and the hierarchical model is recognised as ….
Subtask 2.n)
Use the provided code in the template to plot the predictive distribution for another weight measurement from a chick having the fourth diet and comment on the possible differences you observe between the models.
Subrubric 2.n)
Q27:
Is there a comparison plotted for the predictive distributions of the weight of a chick with diet 4? Does it look something like the model solution plot?
No comparison plotted
Comparison plotted but it clearly differs from the example
Comparison plotted and it approximately matches the example
Q28:
Separate model: Is the result for the separate model discussed?
The result is not discussed.
There is some discussion, but it is not mentioned that ….
The result is discussed and the separate model is recognised as ….
Q29:
Pooled model: Is the result for the pooled model discussed?
The result is not discussed.
There is some discussion, but it is not mentioned that ….
The result is discussed and the pooled model is recognised as ….
Q30:
Hierarchical model: Is the result for the hierarchical model discussed?
The result is not discussed.
There is some discussion, but it is not mentioned that ….
The result is discussed and the hierarchical model is recognised as ….
Subtask 2.o)
Use the provided code in the template to plot the posterior distribution of the mean of the weight measurements of a new fifth diet and comment on the possible differences you observe between the models.
Subrubric 2.o)
Q31:
Is there a comparison plotted for the posterior distributions of the mean weight with a new diet? Does it look something like the model solution plot?
No comparison plotted
Comparison plotted but it clearly differs from the example
Comparison plotted and it approximately matches the example
Q32:
Separate model: Is the result for the separate model discussed?
The result is not discussed.
There is some discussion, but it is not mentioned that ….
The result is discussed and the separate model is recognised as ….
In addition to the previous option, it is recognised that ….
Q33:
Pooled model: Is the result for the pooled model discussed?
The result is not discussed.
There is some discussion, but it is not mentioned that ….
The result is discussed and the pooled model is recognised as ….
In addition to the previous option, it is mentioned that ….
Q34:
Hierarchical model: Is the result for the hierarchical model discussed?
The result is not discussed.
There is some discussion, but it is not mentioned that ….
The result is discussed and the hierarchical model is recognised as ….
Subtask 2.p)
For each model, report the posterior expectation for the mean weight of diet 4 with a 90% credible interval.
Q35:
For the separate model, is the posterior 90% credible interval for the mean of the fourth diet close to: … (small/medium deviation is fine).
No or incorrect answer
Answer is only partially correct
Answers look correct
Q36:
For the pooled model, is the posterior 90% credible interval for the mean of the fourth diet close to: … (small/medium deviation is fine).
No answer
Answer is only partially correct
Answer looks correct
Q37:
For the hierarchical model, is the posterior 90% credible interval for the mean of the fourth diet close to: … (small/medium deviation is fine).
No answer
Answer is only partially correct.
Answers look correct.
3 Hierarchical model with BRMS (3p)
Important
We have made changes to the assignment text and some of the rubrics to make it clearer.
The goal of this task is to discuss an alternative parameterisation of hierarchical models and introduce the brms-package.
brms is a high-level interface for Stan providing tools to create a wide range of Bayesian models, including hierarchical ones. In the previous section, you might have noted that with the hierarchical implementation of the model there may be divergent transitions during the sampling from the posterior. This is related to Neal’s funnel, where the sampling algorithm would need to change the step size in order to effectively sample the posterior. One trick to improve sampling is to use an alternative parameterisation for the model. This parameterisation is called the non-centered parameterisation (Equation 3) and often it can help with the divergent transitions.
In short, we sample parameters from a distribution with a potentially easier geometry for the sampler to explore efficiently and then transform those values to the joint distribution of interest.
The noncentered parameterisation for the hierarchical model using the diet-wise helper parameter \(z_d\) becomes
\[
\begin{aligned}
z_d &\sim \operatorname{normal}(0, 1)&&\text{(diet-wise helper parameter)},\\
\mu_d &= \mu + z_d \tau&&\text{(diet-wise mean weight)},\\
\mu &\sim \pi(\mu)&&\text{(mean of prior for diet-wise mean weights),}\\
\tau &\sim \pi(\tau)&&\text{(standard deviation of prior for diet-wise mean weights),}\\
\sigma &\sim \pi(\sigma)&&\text{(standard deviation of diet-wise chicken weights) and},\\
w_{i, d} &\sim \operatorname{normal}(\mu_d, \sigma)&&\text{(chicken weights)}.\\
\end{aligned}
\tag{3}\]
The hierarchical models made with brms use this non-centered parameterisation. Your tasks are the following:
Subtask 3.a)
Use the function in the template to make a scatter plot of the posterior draws given the hierarchical model of the population level standard deviation parameter \(\tau\) and the mean parameter of the fourth diet, \(\mu_4\). The possible divergent transitions are highlighted in red. In the report show the plot and comment if there are any divergent transitions.
Subrubric 3.a)
Q38:
Is there a scatter plot of \(\sigma_0\) and \(\mu_4\) and some comments provided in the report?
No scatter plot or comments
The scatter plot is included, and divergent transitions are visible, but not commented on.
The scatter plot and comments on the divergences are included.
To sample from the posterior given the same model using brms, in the template set the weakly informative priors for the hierarchical model you used in the previous task. brms has an internal convention to name model parameters and the parameter names from Equation 3 map into brms names with the following logic:
\(\mu\) corresponds to class="Intercept"
\(\tau\) corresponds to class="sd",
\(\sigma\) corresponds to class="sigma".
Subtask 3.b)
Replace normal(0,10) in the brm call in the template with your prior for the population mean \[\mu \sim \pi(\mu) = N(\mu_0,\sigma_0)\] you have derived in 2). Then, run the code chunk to sample from the posterior of the hierarchical model using brms.
Subrubric 3.b)
Q39:
Is the brms code shown and do the priors agree with the priors specified above?
Subtask 3.c)
Report the posterior expectation for the mean weight of diet 4 \(\mu_4\) with a 90% credible interval, using the code block in the template. How does it compare with your results from the first task?
Subrubric 3.c)
Q40:
Given the brms model, is the posterior 90% credible interval for the mean of the fourth diet close to: …(small/medium deviation is fine).
No answer
Answer is only partially correct.
Answers look correct.
Subtask 3.d)
Use the code in the template to make a scatter plot of the group standard deviation parameter \(\tau\) and group specific mean parameter.
Subrubric 3.d)
Q41:
Is there a scatter plot of \(\tau\) and \(\mu_4\) and some comments provided in the report?
No scatter plot or comments
The scatter plot is included, and divergent transitions are visible, but not commented on.
The scatter plot and comments on the divergences are included.
Subtask 3.e)
In your report, address the following questions based on the plots you made for both parameterisations:
Which of the parameterisations resulted in fewer divergent transitions?
Comment for both parameterisations: for which kind of values of \(\tau\), do the divergent transitions occur (if there are any)? It also might be that no clear pattern can be seen.
Does centered parameterisation have problems sampling from some specific region of the parameter space compared to the non-centered parameterisation.
Subrubric 3.e)
Q42:
Has it been discussed which parameterisation resulted in fewer divergent transitions?
Q43:
Has it been discussed for which values of \(\tau\) the divergences occurred?
Q44:
Has it been discussed whether/where the centered parameterisation has problems sampling?
4 Overall quality of the report
Rubric S4: Overall quality of the report - 7.5/100 points
Q45:
Does the report include comment on whether AI was used, and if AI was used, explanation on how it was used?
No
Yes
Q46:
Does the report follow the formatting instructions?
Not at all
Little
Mostly
Yes
Q47:
In case the report doesn’t fully follow the general and formatting instructions, specify the instructions that have not been followed. If applicable, specify the page of the report, where this difference is visible. This will help the other student to improve their reports so that they are easier to read and review. If applicable, specify the page of the report, where this difference in formatting is visible.
Q48:
Please also provide feedback on the presentation (e.g. text, layout, flow of the responses, figures, figure captions). Part of the course is practicing making data analysis reports. By providing feedback on the report presentation, other students can learn what they can improve or what they already did well. You should be able to provide constructive or positive feedback for all non-empty and non-nonsense reports. If you think the report is perfect, and you can’t come up with any suggestions how to improve, you can provide feedback on what you liked and why you think some part of the report is better than yours.
Source Code
---title: "Assignment 7"subtitle: "Hierarchical model in Stan"author: "Aki Vehtari et al."format: html: toc: true code-tools: true code-line-numbers: true number-sections: true mainfont: Georgia, serif page-layout: articleeditor: sourcefilters: - includes/assignments.lua - includes/include-code-files.lua ---# General information**The maximum amount of points from this assignment is 9.**We have prepared a **quarto template specific to this assignment ([html](template7.html), [qmd](https://avehtari.github.io/BDA_course_Aalto/assignments/template7.qmd), [pdf](template7.pdf))** to help you get started. :::{.aalto}We recommend you use [jupyter.cs.aalto.fi](https://jupyter.cs.aalto.fi) or the [docker container](docker.html).::::::{.hint} **Reading instructions**:- [**The reading instructions for BDA3 Chapter 5**](../BDA3_notes.html#ch5).{{< include includes/_grading_instructions.md >}}{{< include includes/_cmdstanr.md >}}:::{{< include includes/_reporting_accuracy.md >}}{{< include includes/_general_info.md >}}:::: {.content-hidden when-format="pdf"}::: {.callout-tip collapse=true}## Setup The following loads several needed packages:```{r}#| label: importslibrary(aaltobda)library(bayesplot)library(cmdstanr)library(dplyr)library(ggplot2)library(ggdist) # for stat_dotsintervallibrary(posterior)if(!require(brms)){install.packages("brms")library(brms)}# 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))# This registers CmdStan as the backend for compiling cmdstan-chunks.check_cmdstan_toolchain(fix =TRUE, quiet =TRUE)register_knitr_engine(override =FALSE)```:::::::# Hierarchical Model: Chicken Data with Stan (6p)**Do not look at the secret dataset until you have defined your priors in section 2.1.**A secret dataset, included in base `r`, contains weight measurements from 50 chicks from their birth until the age of 21 days. For this task, we are interested in modeling the weight of a chick at the age of 12 days.## Choosing a weakly informative prior by intuitionWe will first guide you through two processes for choosing *weakly* informative priors which you can use in your projects as well (though, if expert knowledge is present, we encourage to use that though instead). Many definitions of *weak* exist, but for our purposes, we intend to set priors in this assignment which don't overwhelm the likelihood contributions to the posterior while still exerting some amount of regularisation.In the absence of additional information, we will assume that chick weights at 12 days of age ($w$) follow a normal distribution:$$w\sim\mathcal N (\mu,\sigma).$$Our task will be to define weakly informative priors for the mean $\mu$ given observation model standard deviation, $\sigma$ (a prior distribution for $\sigma$ will be set in after section 2).Here, $\mu$ represents the population mean chick weight at 12 days and $\sigma$ represents the standard-deviation of the chick weights $w$ around that population mean. For this exercise, we will be specifying a normal prior for $\mu$:$$\mu\sim\mathcal N (\mu_0,\sigma_0).$$$\mu_0$ represents our prior knowledge of what we believe the population weight to be and $\sigma_0$ represents our level of certainty in this belief. To specify a weakly-informative prior for $\mu$ we need to select values of $\mu_0$ and $\sigma_0$ that imply the range of plausible values that $\mu$ could take.We can do this by our own intuitive prior knowledge (if such intuition exists), or by searching for external references.Despite the name, weakly informed priors can be quite subjective ([see for more theoretical discussion here](http://www.stat.columbia.edu/~gelman/research/published/entropy-19-00555-v2.pdf)), such that some justification is always needed. For the subtasks below, you will not be graded on the accuracy or precision of your numbers, but on your justification of them. The numerical choices you make should make sense and be understandable to an external reviewer of your work (even if they may not agree with your choices).::: {.callout-caution collapse="false"}## A word of caution on eliciting the priors belowPlease note that in the below, we intend to set a prior on $\mu$ (the *mean* chick weight), but the intuition we ilicit is based on the weight of *individual* chicks. We do so to help create intuition about what the mean could be, however, it would be theoretically more accurate to ilicit priors about *mean* chick weights. ::::::{.callout-important}We have made changes to the assignment text and some of the rubrics to make it clearer.::::::{.subtask letter=a}Based on your own past experience and estimation skills, what would you guess is a typical weight range of a fully grown chicken in grams? Justify your choice with 1-3 sentences.::::::{.hint}Would it make sense if a chicken weighed 1 million grams? 0.005 grams? What about a chicken that weighed more than you? More than a car? Note that is not important to have a very precise or accurate guess here, the key goal is simply to identify the range of plausible (or at least possible) values.:::::::{.subrubric}### Fully grown chicken weight range by intuition* Does the chosen range meet the following common sense criteria? * The range is always above **...** and below **...**.* Is the justification based on some sort of logic, even if you may disagree?:::::::{.subtask}Adjust this range for a 12-day old chick and choose a mean $\mu_0$ for the weakly informed prior of the parameter $\mu.$ Justify your adjustment and choice with 1-3 sentences.:::::::{.subrubric}### 12-day old chick weight range by intuition* The range is always above **...** and below **...**.* Is the justification based on some sort of logic, even if you may disagree?:::: Choosing the prior standard deviation for $\mu$ requires a little more caution; overconfident (i.e. narrow) priors can have a strong effect on your results,whereas less confident priors are more easily overcome with observed data.Given that overconfidence is a common human bias, a good intuition-based standard deviation should focus on eliminating the impossible values, rather than including the most likely values. :::{.subtask}Choose a conservative lower and upper bound for the weight of any 12-day old chick, with the goal to exclude impossible values. Justify your choice with 1-3 sentences.::::::{.hint}Depending on your choice of range above for a "typical" chick, this range excluding the impossible values should be wider.:::::::{.subrubric}### Prior standard deviation by intuition* Does the chosen range meet the following common sense criteria? * The range is always above **...** and below **...**.* Is the justification based on some sort of logic, even if you may disagree?::::A common technique to find a weakly informative prior is to have a standard deviation which is an order of magnitude (a factor of 10) larger than a plausible standard deviation of the data.:::{.subtask}What do you think is a plausible standard deviation $\sigma_\text{plausible}$ of the weight of 12 days-old chicks, based on your ranges stated above? Under the above recommendation, what standard deviation $\sigma_0$ should you use for your prior for the mean weight $\mu$?:::::::{.subrubric}#### Prior standard deviation by intuition* Do the two standard deviations meet the following common sense criteria? - Both values are **...*** Given the choice of mean from above, the interval $(\mu_0 - 3\sigma_0, \mu_0 + 3\sigma_0)$ includes **...**.:::::::{.subtask}Write down in mathematical form the final prior for the mean weight $\mu$ you found using this prior definition technique.:::::::{.subrubric}### Prior by intuition* Does the final prior exist in mathematical notation?* Does the prior reflect the choices made above (i.e. $\mu_0, \sigma_0$)?::::## Choosing a weakly informative prior using external referencesNext, we'll use external references to pick the weakly informed prior. This technique is more general and doesn't assume you would have prior knowledge.:::{.subtask}Consult a trustworthy source on the weight range of farm chickens, e.g. books, articles, a farmer friend. If the recommended values are for a fully grown chicken, make a reasonable adjustment for a 12-day old chick. What is the weight range of a 12-day old chicken? Cite your source and justify any adjustments you make to the reference range.:::::::{.subrubric}### Weight range by reference* Does the reference range meet the following common sense criteria? * The range is always above **...** and below **...**.* Is there a citation? If an adjustment was made, was it justified by logic, even if you disagree?:::::::{.subtask}Based on this reference range, what will you choose for the mean of our weakly informed prior? Justify your choice with 1-3 sentences.:::::::{.subrubric}### Weight range by reference* Does the mean value meet the following common sense criteria? * The range is always above **...** and below **...**.::::Next we choose the standard deviation of the prior. We could use the same technique as before, but we'll walk you through another common approach. Assume that $99.7\%$ of all 12-day old chicks fall into the reference range you found. Under our assumption of a normal distribution, this range will encompass values between $\mu_0 \pm 3\sigma_0$. :::{.subtask}Assuming symmetry, use the mean you chose and either the upper or lower bound $b$ of your reference range to solve the correct version of the following equations to find your associated choice of $\sigma_0$: (show your work)- For upper bound $b_u$: $Pr(\mu_0 + 3\sigma_0 < b_u) \approx 0.997$, solving for $\sigma_0$.- For lower bound $b_l$: $Pr(\mu_0 - 3\sigma_0 > b_l) \approx 0.997$, solving for $\sigma_0$.:::::::{.subrubric}### Prior standard deviation by reference* Does the calculated $\sigma_0$ value meet the following common sense criteria? * The value is above **...** and below **...**.* Did they show their work?:::::::{.subtask}Write down in mathematical form the final prior for the mean weight $\mu$ you found using this prior definition technique.:::::::{.subrubric}### Prior by reference* Does the final prior exist in mathematical notation?* Does the prior reflect the choices made above (i.e. $\mu_0, \sigma_0$)?::::## Non-normal priorsThe previous steps all assumed we could use a prior which is normally distributed, but this may not always be the correct assumption to make. :::{.subtask}Under what mathematical/statistical circumstances would a normal distributed prior not make sense? List at least one circumstance. Are there values that the normal distribution can take on which would not make sense for some types of variables?::::::{.hint}Consider the nature of any variable you are trying to define a prior over.:::::::{.subrubric}### Non-normal priors* Example cases include variables that **...** or variables that **...**::::## Modeling diet effects on chicken weightIn addition to chick weights, the data also contains a categorical variable indicating which diet the chick received. In the data file, each column contains the measurements for a single chick at a given point of time.In addition to the existing diets, we are interested in the quality of another box of feed (the fifth diet), which a farmer happened to find yesterday at a dark corner of his warehouse. To read in the data and select chicks with age of 12 days, and to have a peek at the first 6 rows, just use:::::{.callout-important collapse=true}# Data inside, don't peek before you have set your priors!:::{.callout-important collapse=true}# Have you set your priors?```{r}#| message: falsedata("ChickWeight")Chick12 <- ChickWeight |>filter(Time ==12)head(Chick12)```:::::::In the following analysis, we will model the **weight of a chick at the age of 12 days**. We will use the following three Gaussian models:- a separate model, in which each diet is modeled individually- a pooled model, in which all measurements are combined and there is no distinction between diets- a hierarchical model, which has a hierarchical structure as described in [BDA3](https://users.aalto.fi/~ave/BDA3.pdf) Section 11.6As in the model described in the book, use the same weight standard deviation $\sigma$ for all the groups in the hierarchical model. In the separate model, however, use separate weight standard deviation $\sigma_d$ for each diet $d$. You should use weakly informative priors for all your models.The provided Stan code below is given as an example of a separate model. Note that the author has left a comment expressing uncertainty about their prior choices. This separate model can be summarized mathematically as$$\begin{aligned} \mu_{d} &\sim \pi(\mu_d)&&\text{(diet-wise mean weight),}\\ \sigma_d &\sim \pi(\sigma_d)&&\text{(diet-wise standard deviation of diet-wise chicken weights)},\\ w_{i,d} &\sim N(\mu_d,\sigma_d)&&\text{(diet-wise chicken weights)}\\\end{aligned}$$ {#eq-separate}with priors $$\begin{aligned} \mu_{d} &\sim N(0,10)&&\text{(adjust this) and}\\ \sigma_d &\sim\mathrm{exponential}(.02)&&\text{(you can keep this).}\end{aligned}$$ {#eq-separate-priors}For the separate and the pooled models, use one of the weakly informative priors you have derived in 2.1) or 2.2) for the diet-wise mean weights $$\mu_d\sim \pi(\mu_d)=N(\mu_0,\sigma_0).$$For the hierarchical model, remember that the parameters in the priors for the diet-wise mean weights itself have to be parameters with their own prior distributions, in our case $$\begin{aligned}\mu_d &\sim N(\mu,\tau)&&\text{(diet-wise mean weights)},\\\mu &\sim \pi(\mu)&&\text{(mean of prior for diet-wise mean weights) and}\\\tau &\sim \pi(\tau)&&\text{(standard deviation of prior for diet-wise mean weights).}\end{aligned}$$Use the prior you have used for the diet-wise mean weights $\mu_d$ in the separate and pooled models for the prior on the mean of the prior for the diet-wise mean weights $$\mu \sim \pi(\mu) = N(\mu_0,\sigma_0)$$in the hierarchical model and use $$\tau \sim \pi(\tau) = \mathrm{exponential}(.02)$$ for the prior on the standard deviation of the prior for the diet-wise mean weights.```{.stan include="additional_files/assignment7/chickens_separate.stan"}```:::{.subtask}Describe the models with mathematical notation (as is done for the separate model above). Also describe in words the difference between the model and the other models.:::::::{.subrubric}* Are the models described using mathematical notation and the difference to other models described in words? - No equations and no description - Description but no equations - Equations but no description - Equations and description:::::::{.subtask}Implement the models in Stan and include the code in the report. Use weakly informative priors for all of your models.::::::{.hint}When sampling from the posterior of the hierarchical model, you will very likely get a warning about divergent transitions. This tells you that the sampler is having difficulties in sampling the joint posterior and this may lead to biased inference. We will return to this in the next task.:::::::{.subrubric}### All models* Is there a related Stan implementation? - No Stan model implemented - Stan model implemented, but it seems clearly wrong or broken - Seemingly valid Stan model implemented:::::::{.subtask}Use the provided code in the template toplot the **posterior distribution of the mean of the weight measurements of the fourth diet** and comment on the possible differences you observebetween the models.:::::::{.subrubric}* Is there a comparison plotted for the posteriors of the mean of diet 4? Does it look something like the model solution plot? - No comparison plotted - Comparison plotted but it clearly differs from the example - Comparison plotted and it approximately matches the example* Separate model: Is the result for the separate model discussed? - The result is not discussed. - There is some discussion, but it is not mentioned that **...**. - The result is discussed and the separate model is recognised as **...**.* Pooled model: Is the result for the pooled model discussed? - The result is not discussed. - There is some discussion, but it is not mentioned that **...**. - The result is discussed and the pooled model is recognised as **...**.* Hierarchical model: Is the result for the hierarchical model discussed? - The result is not discussed. - There is some discussion, but it is not mentioned that **...**. - The result is discussed and the hierarchical model is recognised as **...**.:::::::{.subtask}Use the provided code in the template toplot the **predictive distribution for another weight measurement from a chick having the fourth diet** and comment on the possible differences you observebetween the models.:::::::{.subrubric}* Is there a comparison plotted for the predictive distributions of the weight of a chick with diet 4? Does it look something like the model solution plot? - No comparison plotted - Comparison plotted but it clearly differs from the example - Comparison plotted and it approximately matches the example* Separate model: Is the result for the separate model discussed? - The result is not discussed. - There is some discussion, but it is not mentioned that **...**. - The result is discussed and the separate model is recognised as **...**.* Pooled model: Is the result for the pooled model discussed? - The result is not discussed. - There is some discussion, but it is not mentioned that **...**. - The result is discussed and the pooled model is recognised as **...**.* Hierarchical model: Is the result for the hierarchical model discussed? - The result is not discussed. - There is some discussion, but it is not mentioned that **...**. - The result is discussed and the hierarchical model is recognised as **...**.::::::: {.subtask}Use the provided code in the template toplot the **posterior distribution of the mean of the weight measurements of a new fifth diet** and comment on the possible differences you observebetween the models.:::::::{.subrubric}* Is there a comparison plotted for the posterior distributions of the mean weight with a new diet? Does it look something like the model solution plot? - No comparison plotted - Comparison plotted but it clearly differs from the example - Comparison plotted and it approximately matches the example* Separate model: Is the result for the separate model discussed? - The result is not discussed. - There is some discussion, but it is not mentioned that **...**. - The result is discussed and the separate model is recognised as **...**. - In addition to the previous option, it is recognised that **...**.* Pooled model: Is the result for the pooled model discussed? - The result is not discussed. - There is some discussion, but it is not mentioned that **...**. - The result is discussed and the pooled model is recognised as **...**. - In addition to the previous option, it is mentioned that **...**.* Hierarchical model: Is the result for the hierarchical model discussed? - The result is not discussed. - There is some discussion, but it is not mentioned that **...**. - The result is discussed and the hierarchical model is recognised as **...**.:::::::{.subtask}For each model, report the **posterior expectation for the mean weight of diet 4 with a 90% credible interval**.::::::{.hint}See the example Stan codes in the demo [Bayesian data analysis - CmdStanR demos: Comparison of $k$ groups with hierarchical models](http://avehtari.github.io/BDA_R_demos/demos_rstan/cmdstanr_demo.html#8_Comparison_of_(k)_groups_with_hierarchical_models) for the comparison of $k$ groups with and without the hierarchical structure.:::::::{.subrubric}* For the separate model, is the posterior 90% credible interval for the mean of the fourth diet close to: **...** (small/medium deviation is fine). - No or incorrect answer - Answer is only partially correct - Answers look correct* For the pooled model, is the posterior 90% credible interval for the mean of the fourth diet close to: **...** (small/medium deviation is fine). - No answer - Answer is only partially correct - Answer looks correct* For the hierarchical model, is the posterior 90% credible interval for the mean of the fourth diet close to: **...** (small/medium deviation is fine). - No answer - Answer is only partially correct. - Answers look correct.::::# Hierarchical model with BRMS (3p):::{.callout-important}We have made changes to the assignment text and some of the rubrics to make it clearer.:::The goal of this task is to discuss an alternative parameterisation of hierarchical models and introduce the `brms`-package.[`brms`](https://paul-buerkner.github.io/brms/) is a high-level interface for Stan providing tools to create a wide range of Bayesian models, including hierarchical ones. In the previous section, you might have noted that with the hierarchical implementation of the model there may be divergent transitions during the sampling from the posterior. This is related to [Neal's funnel](https://mc-stan.org/docs/stan-users-guide/reparameterisation.html), where the sampling algorithm would need to change the step size in order to effectively sample the posterior. One trick to improve sampling is to use an alternative parameterisation for the model. This parameterisation is called the non-centered parameterisation (@eq-non-centered) and often it can help with the divergent transitions. In short, we sample parameters from a distribution with a potentially easier geometry for the sampler to explore efficiently and then transform those values to the joint distribution of interest.The noncentered parameterisation for the hierarchical model using the diet-wise helper parameter $z_d$ becomes$$\begin{aligned}z_d &\sim \operatorname{normal}(0, 1)&&\text{(diet-wise helper parameter)},\\\mu_d &= \mu + z_d \tau&&\text{(diet-wise mean weight)},\\\mu &\sim \pi(\mu)&&\text{(mean of prior for diet-wise mean weights),}\\\tau &\sim \pi(\tau)&&\text{(standard deviation of prior for diet-wise mean weights),}\\\sigma &\sim \pi(\sigma)&&\text{(standard deviation of diet-wise chicken weights) and},\\w_{i, d} &\sim \operatorname{normal}(\mu_d, \sigma)&&\text{(chicken weights)}.\\\end{aligned}$$ {#eq-non-centered}The hierarchical models made with `brms` use this non-centered parameterisation. Your tasks are the following::::{.subtask letter=a}Use the function in the template to make a scatter plot of the posterior draws given the hierarchical model of the population level standard deviation parameter $\tau$ and the mean parameter of the fourth diet, $\mu_4$. The possible divergent transitions are highlighted in red. In the report show the plot and comment if there are any divergent transitions.::::::{.subrubric}* Is there a scatter plot of $\sigma_0$ and $\mu_4$ and some comments provided in the report? - No scatter plot or comments - The scatter plot is included, and divergent transitions are visible, but not commented on. - The scatter plot and comments on the divergences are included.:::To sample from the posterior given the same model using `brms`, in the template set the weakly informative priors for the hierarchical model you used in the previous task. `brms` has an internal convention to name model parameters and the parameter names from @eq-non-centered map into `brms` names with the following logic: - $\mu$ corresponds to `class="Intercept"` - $\tau$ corresponds to `class="sd"`, - $\sigma$ corresponds to `class="sigma"`.:::{.subtask}Replace `normal(0,10)` in the `brm` call in the template with your prior for the population mean $$\mu \sim \pi(\mu) = N(\mu_0,\sigma_0)$$ you have derived in 2). Then, run the code chunk to sample from the posterior of the hierarchical model using `brms`.:::::: {.subrubric}* Is the `brms` code shown and do the priors agree with the priors specified above?::::::{.subtask}Report the **posterior expectation for the mean weight of diet 4 $\mu_4$ with a 90% credible interval**, using the code block in the template. How does it compare with your results from the first task?:::::: {.subrubric}* Given the `brms` model, is the posterior 90% credible interval for the mean of the fourth diet close to: **...**(small/medium deviation is fine). - No answer - Answer is only partially correct. - Answers look correct.::::::{.subtask}Use the code in the template to make a **scatter plot of the group standard deviation parameter $\tau$ and group specific mean parameter**. :::::: {.subrubric}* Is there a scatter plot of $\tau$ and $\mu_4$ and some comments provided in the report? - No scatter plot or comments - The scatter plot is included, and divergent transitions are visible, but not commented on. - The scatter plot and comments on the divergences are included.::::::{.subtask}In your report, address the following questions based on the plots you made for both parameterisations: - Which of the parameterisations resulted in fewer divergent transitions? - Comment for both parameterisations: for which kind of values of $\tau$, do the divergent transitions occur (if there are any)? It also might be that no clear pattern can be seen. - Does centered parameterisation have problems sampling from some specific region of the parameter space compared to the non-centered parameterisation.:::::: {.subrubric}* Has it been discussed which parameterisation resulted in fewer divergent transitions?* Has it been discussed for which values of $\tau$ the divergences occurred?* Has it been discussed whether/where the centered parameterisation has problems sampling?:::{{< include includes/_overall_quality.md >}}