Loo with brms custom family

General
1

Hi all
There’s something wrong with my R code, I’m trying to fit a custom family in brms to
a right censored dataset and do loo(). When I try to compute the log likelihood, I run out of memory. It says

Error: vector memory exhausted (limit reached?)

Below the reproducible example:

set.seed(343)
dt.sim <- rnbinom(10125, size = 1, prob = rbeta(n = 10125, shape1 = 3, 25))
cens <- sample(1:10125, size = 1700, replace = F)
dt.sim[cens] <- dt.sim[cens] - 5
dt.sim <- ifelse(dt.sim < 0, 0, dt.sim)
censoring <- rep(0, 10125)
censoring[cens] <- 1
test.data <- data.frame(d = dt.sim, censored = censoring)

discrete_weibull <- custom_family(
  "discrete_weibull", dpars = c("mu", "delta"),
  lb = c(0, 0), ub = c(NA, NA), links = c('log', 'log'), 
  type = "int")

stan_funs <- "
  real discrete_weibull_lpmf(int y, real mu, real delta) {
    return log(((1/exp(mu))^y^delta) - (1/exp(mu))^(y+1)^delta);
  }
  real discrete_weibull_lccdf(int y, real mu, real delta){
    return log((1/exp(mu))^(y+1)^delta);
  }
    real discrete_weibull_lcdf(int y, real mu, real delta){
    return log(1 - ((1/exp(mu))^(y+1)^delta));
  }
"

stanvars <- stanvar(scode = stan_funs, block = "functions")

fit.dwei.brms <- brm(form=bf(d|cens(censored)~1, mu ~ 1, delta ~ 1,
                             family=discrete_weibull),
                     stanvars = stanvars,
                     data= test.data, 
                     chains=2,
                     iter = 2000,
                     prior = c(prior('normal(0.1, 0.5)', class = Intercept),
                               prior('normal(1, 0.5)', class = Intercept, dpar = delta)),
                     cores=2,
                     control = list(adapt_delta =0.9))


log_lik_discrete_weibull <- function(i, prep) {
  mu <- prep$dpars$mu[,i]
  delta <- prep$dpars$delta
  y <- prep$data$Y[i]
  log(((1/exp(mu))^y^delta) - (1/exp(mu))^(y+1)^delta)
}

loo(fit.dwei.brms)

My session info:
R version 4.0.0 (2020-04-24)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.4

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages:
[1] brms_2.13.0 Rcpp_1.0.4.6

2

As always, things became clearer when posting the issue. So because there is censored data,
I might need some lccdf function in the log_lik… However, also without using censoring, the error persists.

3

The line

delta <- prep$dpars$delta

is the problem as you predict delta as well. Please replace by

delta <- prep$dpars$delta[, i]
1 Like
4

Thanks, I saw it as well just now.