Using Loo with brms

Interfaces brms
1
  • Operating System: Mac OS X El Capitan Version 10.11.6
  • brms Version: 2.3.1

Rstan version 2.17.3

Hello,

I’m trying to do some model checking and comparison with loo through brms, and I’m running into difficulties. I’ve been trying to do two things.

The first thing is I’ve been trying to create the Loo-PIT Overlay Plot in the Marginal Posterior Predictive Checks section of the “Using the Loo Package” Vignette (https://cran.r-project.org/web/packages/loo/vignettes/loo2-example.html) for my data and models.

loo_model7_3l <- loo(model7_3l, reloo = TRUE)     # Run loo on the model. 

loo_model7_3l    # Show loo output 

Computed from 4000 by 521 log-likelihood matrix

         Estimate   SE
elpd_loo    463.3 17.5
p_loo        94.7  5.8
looic      -926.6 35.0
------
Monte Carlo SE of elpd_loo is NA.

Pareto k diagnostic values:
                         Count Pct.    Min. n_eff
(-Inf, 0.5]   (good)     487   93.5%   574       
 (0.5, 0.7]   (ok)        33    6.3%   110       
   (0.7, 1]   (bad)        1    0.2%   578       
   (1, Inf)   (very bad)   0    0.0%   <NA>      
See help('pareto-k-diagnostic') for details.

plot(loo_model7_3l)    # Plotting Pareto k diagnostics 

![Loo_Plot|654x360](upload://94H5YGRWVvfJpfkn79F2dh6bOTx.png)

So everything is running fine until I get to the code for the Loo-PIT Overlay Plot

yrep <- posterior_predict(model7_3l)  

ppc_loo_pit_overlay(
  y = rwbl_filter$log_MeHg,
  yrep = yrep,
  lw = weights(loo_model7_3l$psis_object)
)

# I get the following error message: 

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

![PIT_check|690x363](upload://uvfchmcoEYSMQGvF7tYrH4Yt9Dq.png)

So I don’t understand why I’m getting a blank plot since the only changes I’ve made to the code is to input my own models and objects. In the vignette, the examples are based on using rstanarm. Do I have to tweak the code to make it run with brms?

The second thing is that I’m having issues running loo with my other models. In this example, I run loo on another model, model8_3l.

loo_model8_3l <- loo(model8_3l)
loo_model8_3l

Computed from 4000 by 521 log-likelihood matrix

         Estimate   SE
elpd_loo    463.4 17.7
p_loo        94.4  6.0
looic      -926.7 35.3
------
Monte Carlo SE of elpd_loo is NA.

Pareto k diagnostic values:
                         Count Pct.    Min. n_eff
(-Inf, 0.5]   (good)     497   95.4%   416       
 (0.5, 0.7]   (ok)        20    3.8%   192       
   (0.7, 1]   (bad)        4    0.8%   37        
   (1, Inf)   (very bad)   0    0.0%   <NA>      
See help('pareto-k-diagnostic') for details.

From the above I can see that I have 4 observations with a pareto_k > 0.7, so I try to run it again with reloo = TRUE. This is where I’m running into issues.

loo_model8_3l <- loo(model8_3l, reloo = TRUE)

Error: New factor levels are not allowed.
Levels allowed: '17-RNA-03_K', '17-RNA-09_B', '17-RNA-09_G', '17-RNA-09_K', '17-RNA-09_R', '17-RNA-5_B', '17-RNA-5_G', '17-RNA-5_K', '17-RNA-5_R', '17-TJT-1_K', '17-TJT-1_R', '17-TJT-100_B', '17-TJT-100_K', '17-TJT-100_R', '17-TJT-102_K', '17-TJT-105_B', '17-TJT-105_K', '17-TJT-105_R', '17-TJT-107_B', '17-TJT-107_G', '17-TJT-107_K', '17-TJT-107_R', '17-TJT-108_B', '17-TJT-108_K', '17-TJT-110_B', '17-TJT-110_K', '17-TJT-111_B', '17-TJT-111_K', '17-TJT-111_R', '17-TJT-115_B', '17-TJT-115_K', '17-TJT-115_R', '17-TJT-117_B', '17-TJT-117_R', '17-TJT-121_B', '17-TJT-121_K', '17-TJT-121_R', '17-TJT-127_B', '17-TJT-127_G', '17-TJT-127_K', '17-TJT-127_R', '17-TJT-128_B', '17-TJT-128_K', '17-TJT-128_R', '17-TJT-13_B', '17-TJT-13_K', '17-TJT-131_B', '17-TJT-131_K', '17-TJT-131_R', '17-TJT-134_B', '17-TJT-134_K', '17-TJT-134_R', '17-TJT-135_B', '17-TJT-135_K', '17-TJT-136_B', '17-TJT-136_K', '17-TJT-139_B', '17-TJT-139_G', '17-TJT-139_K', '17-TJT-139_R',

I’m not sure why this model is having issues when the other wasn’t. Does anyone have an idea of what’s going on here?

Thanks.

2

Sorry, don’t know why the photos didn’t upload. Here’s the plot for the Pareto k diagnostics.

Loo_Plot

And here is the blank Loo-PIT Overlay plot I’m trying to figure out.

PIT_check
PIT_check.png1100×580 3.12 KB

Thanks again.

3

can you provide a minimal reproducible example for me to check out whats happening?

4

Hi Paul,

I would be happy to.

Forgive me, but I’m still a little new to the Github community and sharing code and data. What would you need to see in a minimal reproducible example?

Thanks again,

Tom

5

Provide some code that I can copy and paste into my editor to run it and see the error you see.

6

Code for Loo-PIT Overlay Plot

loo_model7_3l <- loo(model7_3l, reloo = TRUE)

yrep <- posterior_predict(model7_3l)  

ppc_loo_pit_overlay(
  y = rwbl_filter$log_MeHg,
  yrep = yrep,
  lw = weights(loo_model7_3l$psis_object)
)

Code for running loo

loo_model8_3l <- loo(model8_3l)

loo_model8_3l <- loo(model8_3l, reloo = TRUE)

Will this do, or do you need more code?

Thanks again.

7

I can’t run it because I don’t have your model objects

8

I see. Sorry about that!

Code for Loo-PIT Overlay Plot

model7_fm_3l  <- bf(log_MeHg ~ 0 + intercept + Pond_s +
                      (1|NestID/ID))

get_prior(model7_fm_3l, data = rwbl_filter)

priors_7 <- c(prior("student_t(3,1.3,1)", class = "b", coef = "intercept"),
              prior("student_t(3,0,1)", class = "b", coef = "Pond_s"),
              prior("student_t(3,0,1)", class = "sd"),
              prior("student_t(3,0,1)", class = "sigma"))

model7_3l <- brm(model7_fm_3l, data = rwbl_filter, prior = priors_7,
                 iter = 2000, chains = 4, cores = 4, 
                 control = list(adapt_delta = 0.99))

loo_model7_3l <- loo(model7_3l, reloo = TRUE)

yrep <- posterior_predict(model7_3l)  

ppc_loo_pit_overlay(
  y = rwbl_filter$log_MeHg,
  yrep = yrep,
  lw = weights(loo_model7_3l$psis_object)
)

And here’s the code for when I try to use loo on another model

model8_fm_3l  <- bf(log_MeHg ~ 0 + intercept + Grass_s +
                  (1|NestID/ID))

get_prior(model8_fm_3l, data = rwbl_filter)

priors_8 <- c(prior("student_t(3,1.3,1)", class = "b", coef = "intercept"),
              prior("student_t(3,0,1)", class = "b", coef = "Grass_s"),
              prior("student_t(3,0,1)", class = "sd"),
              prior("student_t(3,0,1)", class = "sigma")) #Priors will be the same for both models 


model8_3l <- brm(model8_fm_3l, data = rwbl_filter, prior = priors_8,
                 iter = 2000, chains = 4, cores = 4, 
                 control = list(adapt_delta = 0.99))

loo_model8_3l <- loo(model8_3l)

loo_model8_3l <- loo(model8_3l, reloo = TRUE)

Is there anything else I should include?

Thanks again!

9

To make it fully reproducible, I need some fake data or otherwise I won’t be able to run the above code.

10

Hi Paul,

Here’s code for fake data for my first question

set.seed(25)

a <- 1.23
sigma_nest <- 0.24
sigma_ind <- 0.01
sigma_y <- 0.09
beta <- 0.02      # For Pond_s 

pond_area <- runif(105, min = -2, max = 3)   # Pond area around each nest 

hg_nest <- rstudent_t(105, df = 3, mu = a + beta*pond_area, sigma = sigma_nest) # Hg level of each nest

df_1 <- data.frame(Pond_s = pond_area, NestID = 1:length(hg_nest), hg_nest = hg_nest)  # Joining Pond area, NestID, and Hg level of nests

df_2 <- df_1 %>%
  group_by(NestID) %>%
  mutate(Number = round(rnorm(1, mean = 3, sd = 1))) # Generating numbers for the number of nestlings in each nest

nestlings <- rep(df_2$NestID, times = df_2$Number) # Repeating NestID by the number of nestlings in each nest

df_3 <- data.frame(NestID = nestlings)  # Make nestlings into a dataframe so I can join it with df_1

df_4 <- df_1 %>%
  left_join(df_3, by = "NestID") %>%
  mutate(ID = 1:nrow(df_3)) %>%
  group_by(ID) %>%
  mutate(hg_ind = rstudent_t(n = 1, df = 3, mu = hg_nest, sigma = sigma_ind)) %>%
  ungroup() %>%
  group_by(NestID) %>%
  mutate(number_of_samples = round(runif(1, min = 1, max = 3)))    # Dataframe with Hg level in individuals and the number of samples taken from each nestling

df_5 <- data.frame(ID = rep(df_4$ID, times = df_4$number_of_samples)) # Dataframe for the number of times each nestling was sampled. To be joined with df_4

df_6 <- df_4 %>%
  left_join(df_5, by = "ID") %>%
  ungroup() %>%
  mutate(y = 1:nrow(df_5)) %>%
  group_by(y) %>%
  mutate(log_MeHg = ifelse(number_of_samples == 1, hg_ind,
                           rstudent_t(n = 1, df = 3, mu = hg_ind, sigma = sigma_y))) %>%
  ungroup() %>%
  select(NestID, ID, Pond_s, log_MeHg) # Final data frame with NestID, ID, Pond_s, and log_MeHg

df_6$NestID <- as.character(df_6$NestID)
df_6$ID <- as.character(df_6$ID)           # In my data, NestID and ID are character vectors. Not sure if that's important, but I tried to stay true to form. 

model7_fd_fm  <- bf(log_MeHg ~ 0 + intercept + Pond_s +
                      (1|NestID/ID))

get_prior(model7_fd_fm, data = df_6)

priors_7 <- c(prior("student_t(3,1.3,1)", class = "b", coef = "intercept"),
              prior("student_t(3,0,1)", class = "b", coef = "Pond_s"),
              prior("student_t(3,0,1)", class = "sd"),
              prior("student_t(3,0,1)", class = "sigma"))

model7_fd <- brm(model7_fd_fm, data = df_6, prior = priors_7,
                 iter = 2000, chains = 4, cores = 4, 
                 control = list(adapt_delta = 0.99))

loo_model7_fd <- loo(model7_fd, reloo = TRUE, cores = 4, save_psis = TRUE)  # Loo isn't recognizing the save_psis = TRUE agrument, which I need for the Loo-Pit Overlay Graph

yrep <- posterior_predict(model7_fd)  

ppc_loo_pit_overlay(
  y = df_6$log_MeHg,
  yrep = yrep,
  lw = weights(loo_model7_fd$psis_object)
)

The code for loo runs fine, but I am having difficulty saving the psis object, which is part of the code for the Loo-PIT Overlay graph.

Here’s a fake data simulation for my second model

set.seed(25)

a <- 1.23
sigma_nest <- 0.24
sigma_ind <- 0.01
sigma_y <- 0.09
beta <- -0.07      # For Pond_s 

Grass_area <- runif(105, min = -2, max = 2)   # Grass area around each nest 

hg_nest <- rstudent_t(105, df = 3, mu = a + beta*Grass_area, sigma = sigma_nest) # Hg level of each nest

df_1 <- data.frame(Grass_s = Grass_area, NestID = 1:length(hg_nest), hg_nest = hg_nest)  # Joining Grass area, NestID, and Hg level of nests

df_2 <- df_1 %>%
  group_by(NestID) %>%
  mutate(Number = round(rnorm(1, mean = 3, sd = 1))) # Generating numbers for the number of nestlings in each nest

nestlings <- rep(df_2$NestID, times = df_2$Number) # Repeating NestID by the number of nestlings in each nest

df_3 <- data.frame(NestID = nestlings)  # Make nestlings into a dataframe so I can join it with df_1

df_4 <- df_1 %>%
  left_join(df_3, by = "NestID") %>%
  mutate(ID = 1:nrow(df_3)) %>%
  group_by(ID) %>%
  mutate(hg_ind = rstudent_t(n = 1, df = 3, mu = hg_nest, sigma = sigma_ind)) %>%
  ungroup() %>%
  group_by(NestID) %>%
  mutate(number_of_samples = round(runif(1, min = 1, max = 3)))    # Dataframe with Hg level in individuals and the number of samples taken from each nestling

df_5 <- data.frame(ID = rep(df_4$ID, times = df_4$number_of_samples)) # Dataframe for the number of times each nestling was sampled. To be joined with df_4

df_6 <- df_4 %>%
  left_join(df_5, by = "ID") %>%
  ungroup() %>%
  mutate(y = 1:nrow(df_5)) %>%
  group_by(y) %>%
  mutate(log_MeHg = ifelse(number_of_samples == 1, hg_ind,
                           rstudent_t(n = 1, df = 3, mu = hg_ind, sigma = sigma_y))) %>%
  ungroup() %>%
  select(NestID, ID, Grass_s, log_MeHg) # Final data frame with NestID, ID, Grass_s, and log_MeHg

df_6$NestID <- as.character(df_6$NestID)
df_6$ID <- as.character(df_6$ID)           # In my data, NestID and ID are character vectors. Not sure if that's important, but I tried to stay true to form. 

model8_fd_fm  <- bf(log_MeHg ~ 0 + intercept + Grass_s +
                      (1|NestID/ID))

get_prior(model8_fd_fm, data = df_6)

priors_8 <- c(prior("student_t(3,1.3,1)", class = "b", coef = "intercept"),
              prior("student_t(3,0,1)", class = "b", coef = "Grass_s"),
              prior("student_t(3,0,1)", class = "sd"),
              prior("student_t(3,0,1)", class = "sigma"))

model8_fd <- brm(model8_fd_fm, data = df_6, prior = priors_8,
                 iter = 2000, chains = 4, cores = 4, 
                 control = list(adapt_delta = 0.99))

loo_model8_fd <- loo(model8_fd, reloo = TRUE, cores = 4)

I can get loo to run for this model, which is strange, because I was hoping for the same warning messages when I run the real model 8. As a reminder, that warning message looks like:

Error: New factor levels are not allowed.
Levels allowed: '17-RNA-03_K', '17-RNA-09_B', '17-RNA-09_G', '17-RNA-09_K', '17-RNA-09_R', '17-RNA-5_B', '17-RNA-5_G', '17-RNA-5_K', '17-RNA-5_R', '17-TJT-1_K', '17-TJT-1_R', '17-TJT-100_B', '17-TJT-100_K', '17-TJT-100_R', '17-TJT-102_K', '17-TJT-105_B', '17-TJT-105_K', '17-TJT-105_R', '17-TJT-107_B', '17-TJT-107_G', '17-TJT-107_K', '17-TJT-107_R', '17-TJT-108_B', '17-TJT-108_K', '17-TJT-110_B', '17-TJT-110_K', '17-TJT-111_B', '17-TJT-111_K', '17-TJT-111_R', '17-TJT-115_B', '17-TJT-115_K', '17-TJT-115_R', '17-TJT-117_B', '17-TJT-117_R', '17-TJT-121_B', '17-TJT-121_K', '17-TJT-121_R', '17-TJT-127_B', '17-TJT-127_G', '17-TJT-127_K', '17-TJT-127_R', '17-TJT-128_B', '17-TJT-128_K', '17-TJT-128_R', '17-TJT-13_B', '17-TJT-13_K', '17-TJT-131_B', '17-TJT-131_K', '17-TJT-131_R', '17-TJT-134_B', '17-TJT-134_K', '17-TJT-134_R', '17-TJT-135_B', '17-TJT-135_K', '17-TJT-136_B', '17-TJT-136_K', '17-TJT-139_B', '17-TJT-139_G', '17-TJT-139_K', '17-TJT-139_R',

Thanks again for your help, and please let me know if you need anything else.

11

Thanks. Now I see that your brms version is outdated. You use 2.3.1, but the release is 2.6.0. You should definitely update it. Perhaps, the 2. problem goes away after updating your brms version.

A suggestion for the 1. problem with loo_pit. Try out instead

pp_check(model7_fd, type = "loo_pit_overlay")

With this, brms is able to figure out that it needs lw and will add it internally.

For your benefit: Minimal reproducible examples should be minimal in the sense that they only contain the code necessary without all the other stuff around. And they need to produce the error you see to help developers figure out what’s wrong. Else you give the full burden to the developer to figure out what was producing your error which may take a considerable amount of time on top of the time it takes to fix the problem.

12

I am also running into the same issue where pp_check(fit1, type=“loo_pit_overlay”) gives a message:
“Coordinate system already present. Adding new coordinate system, which will replace the existing one.”
and then the resulting plot looks empty, just like @tthalhuber 's plot. I am using the first example in the brm documentation, i.e.:

bprior1 ← prior(student_t(5,0,10), class = b) +
prior(cauchy(0,2), class = sd)
fit1 ← brm(count ~ log_Age_c + log_Base4_c * Trt + (1|patient),
data = epilepsy, family = poisson(), prior = bprior1)

sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.1

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] bindrcpp_0.2.2 brms_2.6.0 Rcpp_1.0.0 patchwork_0.0.1 glue_1.3.0 lme4_1.1-18-1 Matrix_1.2-15
[8] reshape2_1.4.3 readxl_1.1.0 forcats_0.3.0 stringr_1.3.1 dplyr_0.7.7 purrr_0.2.5 readr_1.1.1
[15] tidyr_0.8.2 tibble_1.4.2 ggplot2_3.1.0 tidyverse_1.2.1

loaded via a namespace (and not attached):
[1] nlme_3.1-137 matrixStats_0.54.0 xts_0.11-2 lubridate_1.7.4 threejs_0.3.1
[6] httr_1.3.1 rprojroot_1.3-2 rstan_2.18.2 tools_3.5.1 backports_1.1.2
[11] R6_2.3.0 DT_0.5 lazyeval_0.2.1 colorspace_1.3-2 withr_2.1.2
[16] prettyunits_1.0.2 processx_3.2.0 tidyselect_0.2.5 gridExtra_2.3 Brobdingnag_1.2-6
[21] compiler_3.5.1 cli_1.0.1 rvest_0.3.2 xml2_1.2.0 shinyjs_1.0
[26] labeling_0.3 colourpicker_1.0 scales_1.0.0 dygraphs_1.1.1.6 mvtnorm_1.0-8
[31] callr_3.0.0 ggridges_0.5.1 qwraps2_0.3.0 StanHeaders_2.18.0 digest_0.6.18
[36] minqa_1.2.4 rmarkdown_1.10 base64enc_0.1-3 pkgconfig_2.0.2 htmltools_0.3.6
[41] htmlwidgets_1.3 rlang_0.3.0.1 rstudioapi_0.8 shiny_1.2.0 bindr_0.1.1
[46] zoo_1.8-4 jsonlite_1.5 crosstalk_1.0.0 gtools_3.8.1 inline_0.3.15
[51] magrittr_1.5 loo_2.0.0 bayesplot_1.6.0 munsell_0.5.0 abind_1.4-5
[56] stringi_1.2.4 yaml_2.2.0 MASS_7.3-51.1 pkgbuild_1.0.2 plyr_1.8.4
[61] grid_3.5.1 parallel_3.5.1 promises_1.0.1 crayon_1.3.4 miniUI_0.1.1.1
[66] lattice_0.20-38 haven_1.1.2 splines_3.5.1 hms_0.4.2 ps_1.2.1
[71] knitr_1.20 pillar_1.3.0 igraph_1.2.2 markdown_0.8 shinystan_2.5.0
[76] codetools_0.2-15 stats4_3.5.1 rstantools_1.5.1 evaluate_0.12 modelr_0.1.2
[81] nloptr_1.2.1 httpuv_1.4.5 cellranger_1.1.0 gtable_0.2.0 assertthat_0.2.0
[86] mime_0.6 xtable_1.8-3 broom_0.5.0 coda_0.19-2 later_0.7.5
[91] rsconnect_0.8.8 shinythemes_1.1.2 bridgesampling_0.6-0

13

There appears to be a bug in bayesplot not in brms, which merely passes stuff to the former. Would you mind opening an issue at https://github.com/stan-dev/bayesplot describing this problem?

14

Sure - thanks for your quick reply.

15

This has been solved in the dev version of bayesplot.

See: Issue 167

16

Great!

17

Hi Kyle,

Thanks for your input. Have you gotten your graphs to work now? I updated brms like Paul recommended and tried pp_check(model7_fd, type = "loo_pit_overlay"), but I’m still getting the same warning and the same blank plot. How do I take advantage of the bug fix? Do I need to uninstall and then reinstall bayesplot?

Thanks again,

Tom

18

Run

if (!require("devtools")) {
  install.packages("devtools")
}
devtools::install_github("stan-dev/bayesplot", dependencies = TRUE, build_vignettes = TRUE)
19

Hi Paul,

That worked! Thank you so much for your patience and assistance!

Tom