/
Assignment08.Rmd
199 lines (166 loc) · 10.5 KB
/
Assignment08.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
---
title: "Assignment08"
author: "Shane"
date: "14/03/2024"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Libraries:
```{r, warning=FALSE,, message=FALSE}
#install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
library(lme4)
library(brms)
options(brms.backend = "cmdstanr")
library(broom.mixed)
library(purrr)
library(dplyr)
library(tidybayes)
library(bayesplot)
library(bayestestR)
library(ggplot2); theme_set(theme_bw())
library(see)
options(ggplot2.discrete.colour= scale_color_okabeito)
library(ggrastr)
library(cowplot)
```
Note, I included all packages used in the example in class: I am sure some are not completely necessary here. I have also suppressed errors and warning messages as I was getting info about compiling and libraries loading additional packages.
## Load Data
```{r}
dd <- readRDS("PupBirthdaysClean.rds")
summary(dd)
```
## Formula to Use for Priors
```{r}
form01 <- Julian.Date ~ 1 + Group
get_prior(form01, dd)
```
Here, we have the prior of Group (*GroupWildMCaught*), the Intercept with parameters - as set by Student t distribution with a df = 3, a mean equal to the median of Julian Date (168.5) and a standard deviation equal to the mean absolute deviation. Note, this intercept is somewhat irrational within Bayes Theorm: The function *get_prior* is estimating the intercept by examining the data and estimates the mean and SD based on observations within the data; however, this method is suitable to provide estimates. The df = 3 allows the model relaxes the restrictions of the Bayesian model while maintaining reasonable estimates. **BMB**: "relaxes the restrictions of the Bayesian model" sounds a bit loose/vague, it sounds like there's something non-Bayesian about it (what it does is to allow for *heavy tails* in the prior distribution ...
```{r}
b_prior <- c(set_prior("normal(160, 2.5)", "Intercept"),
set_prior("normal(0, 10)", "b"))
b_prior_02 <- c(set_prior("normal(20, 2.5)", "Intercept"),
set_prior("normal(0, 10)", "b"))
```
Here, I have constrained the intercept and slope to reasonable values (I believe). I have set up a second prior (*b_prior_02*) with, what I think, is an unreasonable value - Setting the mean to 20 would imply that most bats give birth in the winter, which is an incredibly rare occurrence. For learning purposes, I am curious how this will affect the model.
Note, the only thing that differs between these two models are the means set for each prior where b_prior is an expected prior and b_prior_02 is an unreasonable prior given what is known about gestation and birth in bats. Again, I included it purely out of curiosity and to learn (i..e test limits) of the Bayesian Model. **Reasonable priors give rise to reasonable outputs!**
These priors were set one at a time, assuming independence (this may not be good). However, this is the method that proves reliable. **BMB**: yes, awkward to do anything else.
```{r}
test_prior <- function(p) {
## https://discourse.mc-stan.org/t/suppress-all-output-from-brms-in-markdown-files/30117/3
capture.output(
b <- brm(form01, dd, prior = p,
seed = 101, ## reproducibility
sample_prior = 'only', ## for prior predictive sim
chains = 1, iter = 500, ## very short sample for convenience
silent = 2, refresh = 0 ## be vewy vewy quiet ...
)
)
p_df <- dd |> add_predicted_draws(b)
## 'spaghetti plot' of prior preds
gg0 <- ggplot(p_df,aes(Group, .prediction)) + #, group=interaction(Subject,.draw))) +
geom_jitter(height=0) + geom_violin(fill = "red", alpha = 0.4)
print(gg0)
invisible(b) ## return without printing
}
test_prior(b_prior)
# Note, this is code used in the example given in class.
```
This figure does not look like what was present in class. I believe it is because my predictor variable is a two-level factor (Group). These model diagnostics serve to see if the model is reasonable (i.e. would it be able to reproduce what is seen in my observed data).
**BMB**: you confused predictor and response!
```{r, message=FALSE}
model.lm <- lm(form01, dd)
model.brms <- brm(form01, dd, prior = b_prior, # This is the model of interest
seed = 101,
control = list(adapt_delta = 0.95))
model.brms.02 <- brm(form01, dd, prior = b_prior_02, # Note, I do not expect this model to be good.
seed = 101,
control = list(adapt_delta = 0.95))
model.brms.def <- brm(form01, dd, # This model is the default, with no priors set
seed = 101)
```
## Diagnose
```{r}
print(bayestestR::diagnostic_posterior(model.brms), digits = 4)
```
Here, we see our parameters of Group (i.e. b_GroupWildMCaught) and Intercept (i.e. b_Intercept). To my understanding, the measures of Rhat, Effective Sampling Size (ESS) and Monte Carlo Standard Error (MCSE) are are diagnostics of whether the set priors are reasonable. The Rhat values should hover around 1, which for both parameters in this model, they do. The ESS should be well beyond 400 (Note, this value is a recommendation), which is true for both parameters in this model (Group = 3825 and Intercept = 3525). and the MCSE should be "small enough" for scientific purposes. In this case, the MCSE values are incredibly small for both parameters.
**BMB**: interesting that Rhat<1. I'm not sure if this is theoretically impossible (i.e. a computational issue) or just unlikely. (Not something to worry about)
## Trace Plots
```{r}
mcmc_trace(model.brms, regex_pars = "b_|sd_")
```
Shown in this plot are the diagnostics of four (4) chains displayed in slightly different shades of blue. The Figure displays history of samples bouncing around the mean for Intercept (Left) and Group (Right). These plots appear steady (i.e. no obscurities that stand out) and the chains appear to be following similar patterns (steady) as one another.
**BMB**: the 'new' version of trace plots is recommended (but this is fine)
## Looking at the Results of the Model
```{r}
summary(model.brms)
```
```{r}
brms_modlist <- list(brms_default = model.brms.def, brms_reg = model.brms, brms_irreg = model.brms.02)
res_bayes <- (brms_modlist
|> purrr::map_dfr(~ tidy(., conf.int = TRUE), .id = "model")
)
## need to do separately - different conf.method choices
res_lm <- suppressMessages(model.lm
|> tidy(conf.int = TRUE, conf.method = "profile")
|> mutate(model = "lm", .before = 1)
)
res <- (bind_rows(res_lm, res_bayes)
|> select(-c(std.error, statistic, component, group))
|> filter(term != "(Intercept)")
|> mutate(facet = ifelse(grepl("^cor", term), "cor",
ifelse(grepl("Days", term), "Days",
"int")))
|> mutate(across(model, ~ factor(., levels = c("lm", names(brms_modlist)))))
)
ggplot(res, aes(estimate, term, colour = model, shape = model)) +
geom_pointrange(aes(xmin = conf.low, xmax = conf.high),
position = position_dodge(width = 0.5)) +
facet_wrap(~ facet, scales = "free", ncol = 1) +
guides(colour = guide_legend(reverse=TRUE),
shape = guide_legend(reverse=TRUE))
```
This figure is comparing the estimates determined from the different models: (1) brms_reg - which uses the manually set priors from model.brms, which is a normal distribution with a mean of 150 and SD of 2.5 for the intercept and a normal distribution with a mean of 0 and SD of 10 for b (*See b_prior above*), (2) brms_irreg - which uses the manually set priors from model.brms.02, which is a normal distribution with a mean of 20 and SD of 2.5 for the intercept and a normal distribution with a mean of 0 and SD of 10 for b (*See b_prior_02 above*), and (3) brms_default - which uses the default set parameters.
As expected, the irregular model does not fit well to the others, as it is obscure relative to what is expected (i.e. the priors are incorrect / inappropriate). I will remove this model so we can see with clarity how the regular model and the default model compare:
```{r}
brms_modlist_clean <- list(brms_default = model.brms.def, brms_reg = model.brms)
res_bayes <- (brms_modlist_clean
|> purrr::map_dfr(~ tidy(., conf.int = TRUE), .id = "model")
)
## need to do separately - different conf.method choices
res_lm <- suppressMessages(model.lm
|> tidy(conf.int = TRUE, conf.method = "profile")
|> mutate(model = "lm", .before = 1)
)
res <- (bind_rows(res_lm, res_bayes)
|> select(-c(std.error, statistic, component, group))
|> filter(term != "(Intercept)")
|> mutate(facet = ifelse(grepl("^cor", term), "cor",
ifelse(grepl("Days", term), "Days",
"int")))
|> mutate(across(model, ~ factor(., levels = c("lm", names(brms_modlist)))))
)
ggplot(res, aes(estimate, term, colour = model, shape = model)) +
geom_pointrange(aes(xmin = conf.low, xmax = conf.high),
position = position_dodge(width = 0.5)) +
facet_wrap(~ facet, scales = "free", ncol = 1) +
guides(colour = guide_legend(reverse=TRUE),
shape = guide_legend(reverse=TRUE))
```
This looks much better.
Note, the estimates are very similar between the two models: The two models, despite having different priors set, are modelling the data well. I am not sure why the lm model is on a separate line for the Group estimate and not present for the sd_Observation estimate. However, the Group estimate for the lm model seems to align closely with the brms_reg and brms_default models.
**BMB**: this is just a difference in naming (not quite sure where that happened)
## Posterior Predictive Simulations , Compare with Data
```{r}
post_df1 <- dd |> add_predicted_draws(model.brms)
gg1 <- ggplot(post_df1,
aes(Group, .prediction)) +
geom_violin(fill = "gray") +
stat_sum(data=dd, aes(y=Julian.Date), col = "red") +
labs(y = "Day") +
scale_size(breaks = c(1, 5, 10))
print(gg1 + labs(title = "Priors (weird)"))
```
Note, this plot differs from the example showed in class in that it does not have the black lines showcasing the sampling from the Bayesian model. I *think* this is because my model uses a two-level factor as the predictor variable as opposed to a continuous variable, as used in the example. I believe this figure is showing random selections of possible estimates (i.e. **red dots**) for each group as determined from the model. I am not certain to what this figure is meant to display or what information can be gleaned from it.
**BMB**: see adjusted version. Mark: 2.2