Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

prior_summary() incorrect after updating a cumulative() fit using fewer response categories #1380

Closed
fweber144 opened this issue Jul 18, 2022 · 3 comments
Labels
Milestone

Comments

@fweber144
Copy link
Contributor

fweber144 commented Jul 18, 2022

When update()-ing a cumulative() fit with a new dataset that contains fewer response categories than the original one, prior_summary() still displays all the thresholds from the initial model fit:

options(mc.cores = parallel::detectCores(logical = FALSE))

# Data --------------------------------------------------------------------

nobsv <- 150L
ncat <- 4L
yunq <- paste0("ycat", seq_len(ncat))
nthres <- ncat - 1L
thres <- seq(-1, 1, length.out = nthres)
npreds_cont <- 5L
coefs_cont <- seq(-4, 4, length.out = npreds_cont)
set.seed(3738)
dat_sim <- matrix(rnorm(nobsv * npreds_cont, sd = 0.5), nrow = nobsv,
                  dimnames = list(NULL, paste0("Xcont", seq_len(npreds_cont))))
eta <- drop(dat_sim %*% coefs_cont)
dat_sim <- as.data.frame(dat_sim)
thres_eta <- sapply(thres, function(thres_k) {
  thres_k - eta
})
yprobs <- brms:::inv_link_cumulative(thres_eta, link = "logit")
dat_sim$Y <- factor(sapply(seq_len(nobsv), function(i) {
  sample(yunq, size = 1, replace = TRUE, prob = yprobs[i, ])
}), ordered = TRUE)

# Model fit ---------------------------------------------------------------

stopifnot(identical(levels(dat_sim$Y), yunq))
bfit <- brms::brm(
  formula = Y ~ Xcont1 + Xcont2 + Xcont3 + Xcont4 + Xcont5,
  data = dat_sim,
  family = brms::cumulative(),
  prior = brms::prior(normal(0, 1), class = "Intercept"),
  seed = 878232
)
brms::prior_summary(bfit)

## Update model fit -------------------------------------------------------

# Update using a dataset which lacks a response category:
dat_upd <- dat_sim[dat_sim$Y != "ycat3", , drop = FALSE]
dat_upd$Y <- droplevels(dat_upd$Y)
stopifnot(identical(levels(dat_upd$Y), setdiff(yunq, "ycat3")))
bupd <- update(
  bfit,
  formula. = Y ~ Xcont1 + Xcont2 + Xcont3 + Xcont4 + Xcont5,
  newdata = dat_upd,
  ### Actually not needed (only added to check if this makes a difference, but
  ### it doesn't):
  prior = brms::prior(normal(0, 1), class = "Intercept"),
  ###
  seed = 54543
)
brms::prior_summary(bupd)

where the last line gives

       prior     class   coef group resp dpar nlpar lb ub       source
      (flat)         b                                         default
      (flat)         b Xcont1                             (vectorized)
      (flat)         b Xcont2                             (vectorized)
      (flat)         b Xcont3                             (vectorized)
      (flat)         b Xcont4                             (vectorized)
      (flat)         b Xcont5                             (vectorized)
normal(0, 1) Intercept                                            user
normal(0, 1) Intercept      1                             (vectorized)
normal(0, 1) Intercept      2                             (vectorized)
normal(0, 1) Intercept      3                             (vectorized)

just like for the initial model fit.

I tested this with brms v2.17.0. As shown in the update() call above, I tried to specify argument prior explicitly a second time, but that didn't help.

Apart from this prior_summary() thing, all other post-processing functions seem to work correctly, at least as far as I have checked them.

@paul-buerkner
Copy link
Owner

This is because priors are not updated automatically by update in the sense that they are not recomputed from newdata. This of course only concerns the priors that are data-dependent. I will see what options we have in the case you presented.

@paul-buerkner paul-buerkner added this to the brms 2.17.0++ milestone Jul 18, 2022
paul-buerkner added a commit that referenced this issue Aug 12, 2022
@paul-buerkner
Copy link
Owner

Thank you for opening this issue. Should now be fixed.

@fweber144
Copy link
Contributor Author

Yes, thank you. The example from above works now as expected.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

2 participants