Hi all
I’d like to ask if one needs a further adjustment of the code provided in the custom family vignette
(https://cran.r-project.org/web/packages/brms/vignettes/brms_customfamilies.html)
if one is to evaluate model fit with loo on a dataset with censored observations.
Specifically, do I need to change the log_like_custom function as below to get a proper elpd with brms:loo()?
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[,i]
y <- prep$data$Y[i]
cens <- prep$data$cens[,i]
if(cens == 0){
log(((1/exp(mu))^y^delta) - (1/exp(mu))^(y+1)^delta)
} else{
log((1/exp(mu))^(y+1)^delta)
}
}
loo(fit.dwei.brms)
The log_lik_discrete_weibull function does not work because there is no prep$data$cens, but any comment if this is in principle the way to go is very much appreciated