more

avehtari · Jan 3, 2024 · b84addd · b84addd
1 parent 2eec955
commit b84addd
Show file tree

Hide file tree

Showing 2 changed files with 97 additions and 5 deletions.
diff --git a/demos_rstan/brms_demo.R b/demos_rstan/brms_demo.R
@@ -535,13 +535,18 @@ data_lin |>
   theme(legend.position="none")+
   scale_x_continuous(breaks=seq(1950,2030,by=10))
 
+#' LOO-PIT check is good for checking whether the normal distribution
+#' is well describing the variation. LOO-PIT ploty looks good.
+pp_check(fit_lin, type='loo_pit_qq', ndraws=4000)
+
 #' # Linear Student's $t$ model
 #' 
 #' The temperatures used in the above analyses are averages over three
 #' months, which makes it more likely that they are normally
 #' distributed, but there can be extreme events in the feather and we
 #' can check whether more robust Student's $t$ observation model would
-#' give different results.
+#' give different results (although LOO-PIT check did already indicate
+#' that the normal would be good).
 #' 
 #+  results='hide'
 fit_lin_t <- brm(temp ~ year, data = data_lin, family = student(),
@@ -1000,13 +1005,25 @@ theta |>
   labs(x='Odds-ratio', y='Treatment', title='Pooled over studies, separate over treatments') +
   geom_vline(xintercept=1, linetype='dashed')
 
-
 #' We see a big variation between treatments and two treatments seem
 #' to be harmful, which is suspicious. Looking at the data we see that
 #' not all studies included all treatments, and thus if some of the
 #' studies had more events, then the above estimates can be wrong.
 #' 
 
+#' Posterior predictive checking with kernel density estimates for the
+#' data and 10 posterior predictive replicates shows clear discrepancy.
+pp_check(fit_pooled, type='dens_overlay')
+
+#' Posterior predictive checking with PIT values and ECDF difference
+#' plot with envelope shows clear discrepancy.
+pp_check(fit_pooled, type='pit_ecdf', ndraws=4000)
+
+#' Posterior predictive checking with LOO-PIT values show clear discrepancy.
+pp_check(fit_pooled, type='loo_pit_qq', ndraws=4000) +
+  geom_abline() +
+  ylim(c(0,1))
+
 #' The second model uses a hiearchical model both for treatment
 #' effects and study effects.
 #+ results='hide'
@@ -1026,6 +1043,23 @@ loo_compare(loo(fit_pooled), loo(fit_hier))
 #' between studies. Clearly there are many highly influential
 #' observations.
 #'
+#' Posterior predictive checking with kernel density estimates for the
+#' data and 10 posterior predictive replicates looks good (although
+#' with this many parameters, this check is likely to be optimistic).
+pp_check(fit_hier, type='dens_overlay')
+
+#' Posterior predictive checking with PIT values and ECDF difference
+#' plot with envelope looks good (although with this many parameters,
+#' this check is likely to be optimistic).
+pp_check(fit_hier, type='pit_ecdf', ndraws=4000)
+
+#' Posterior predictive checking with LOO-PIT values look good
+#' (alhough as there are Pareto-khat warnings, it is possible that
+#' this diagnostic is optimistic).
+pp_check(fit_hier, type='loo_pit_qq', ndraws=4000) +
+  geom_abline() +
+  ylim(c(0,1))
+
 #' Treatment effect posteriors have now much less variation.
 fit_hier |>
   spread_rvars(b_Intercept, r_treatment[treatment,]) |>
@@ -1034,7 +1068,7 @@ fit_hier |>
   stat_halfeye() +
   labs(x='theta', y='Treatment', title='Hierarchical over studies, hierarchical over treatments')  
 
-#' Study effect posteriors show the expected high variation
+#' Study effect posteriors show the expected high variation.
 fit_hier |>
   spread_rvars(b_Intercept, r_study[study,]) |>
   mutate(theta_study = rfun(plogis)(b_Intercept + r_study)) |>
@@ -1061,6 +1095,8 @@ theta |>
 #' similar uncertainty about the overall theta due to high variation
 #' between studies.
 
+
+
 #' The third model includes interaction so that the treatment can depend on study.
 #+ results='hide'
 fit_hier2 <- brm(exac | trials(total) ~ (1 | treatment) + (treatment | study),

diff --git a/demos_rstan/brms_demo.Rmd b/demos_rstan/brms_demo.Rmd
@@ -710,13 +710,21 @@ data_lin |>
   scale_x_continuous(breaks=seq(1950,2030,by=10))
 ```
 
+LOO-PIT check is good for checking whether the normal distribution
+is well describing the variation. LOO-PIT ploty looks good.
+
+```{r}
+pp_check(fit_lin, type='loo_pit_qq', ndraws=4000)
+```
+
 # Linear Student's $t$ model
 
 The temperatures used in the above analyses are averages over three
 months, which makes it more likely that they are normally
 distributed, but there can be extreme events in the feather and we
 can check whether more robust Student's $t$ observation model would
-give different results.
+give different results (although LOO-PIT check did already indicate
+that the normal would be good).
 
 
 ```{r results='hide'}
@@ -1329,6 +1337,28 @@ to be harmful, which is suspicious. Looking at the data we see that
 not all studies included all treatments, and thus if some of the
 studies had more events, then the above estimates can be wrong.
 
+Posterior predictive checking with kernel density estimates for the
+data and 10 posterior predictive replicates shows clear discrepancy.
+
+```{r}
+pp_check(fit_pooled, type='dens_overlay')
+```
+
+Posterior predictive checking with PIT values and ECDF difference
+plot with envelope shows clear discrepancy.
+
+```{r}
+pp_check(fit_pooled, type='pit_ecdf', ndraws=4000)
+```
+
+Posterior predictive checking with LOO-PIT values show clear discrepancy.
+
+```{r}
+pp_check(fit_pooled, type='loo_pit_qq', ndraws=4000) +
+  geom_abline() +
+  ylim(c(0,1))
+```
+
 The second model uses a hiearchical model both for treatment
 effects and study effects.
 
@@ -1356,6 +1386,32 @@ hierarchical model is much better and there is high variation
 between studies. Clearly there are many highly influential
 observations.
 
+Posterior predictive checking with kernel density estimates for the
+data and 10 posterior predictive replicates looks good (although
+with this many parameters, this check is likely to be optimistic).
+
+```{r}
+pp_check(fit_hier, type='dens_overlay')
+```
+
+Posterior predictive checking with PIT values and ECDF difference
+plot with envelope looks good (although with this many parameters,
+this check is likely to be optimistic).
+
+```{r}
+pp_check(fit_hier, type='pit_ecdf', ndraws=4000)
+```
+
+Posterior predictive checking with LOO-PIT values look good
+(alhough as there are Pareto-khat warnings, it is possible that
+this diagnostic is optimistic).
+
+```{r}
+pp_check(fit_hier, type='loo_pit_qq', ndraws=4000) +
+  geom_abline() +
+  ylim(c(0,1))
+```
+
 Treatment effect posteriors have now much less variation.
 
 ```{r}
@@ -1367,7 +1423,7 @@ fit_hier |>
   labs(x='theta', y='Treatment', title='Hierarchical over studies, hierarchical over treatments')  
 ```
 
-Study effect posteriors show the expected high variation
+Study effect posteriors show the expected high variation.
 
 ```{r}
 fit_hier |>