Hey @StaffanBetner, I’ve been practicing with your code set for cumulative ordinal models, both with setting phi to either mean(Intercept) or to 0. Either way, the fit seems substantially worse than the default method. Where am I going wrong?
Here’s the model set up:
# load packages
library(tidyverse)
library(brms)
library(patchwork)
# define dirichlet_prior for the probit approach
dirichlet_prior <- "
real induced_dirichlet_lpdf(vector c, vector alpha, real phi) {
int K = num_elements(c) + 1;
vector[K - 1] anchoredcutoffs = phi - c;
// sigma is the CDF of the latent distribution
// logit:
// vector[K - 1] sigma = inv_logit(anchoredcutoffs);
// probit:
vector[K - 1] sigma = Phi(anchoredcutoffs);
// or:
// vector[K - 1] sigma = Phi_approx(anchoredcutoffs);
vector[K] p;
matrix[K, K] J = rep_matrix(0, K, K);
// Induced ordinal probabilities
p[1] = 1 - sigma[1];
for (k in 2:(K - 1))
p[k] = sigma[k - 1] - sigma[k];
p[K] = sigma[K - 1];
// Baseline column of Jacobian
for (k in 1:K) J[k, 1] = 1;
// Diagonal entries of Jacobian
for (k in 2:K) {
// rho is the PDF of the latent distribution
// logit:
// real rho = sigma[k - 1] * (1 - sigma[k - 1]);
// probit:
real rho = exp(std_normal_lpdf(anchoredcutoffs));
J[k, k] = - rho;
J[k - 1, k] = rho;
}
return dirichlet_lpdf(p | alpha)
+ log_determinant(J);
}
"
# define the stanvar()
dirichlet_prior_stanvar <- stanvar(scode = dirichlet_prior, block = "functions")
# fit the conventional cumulative probit model
test_standard <- brm(
formula = cyl ~ 1 + mpg,
data = mtcars,
family = cumulative(probit),
cores = 4, seed = 4
)
# fit the model with induced_dirichlet() and phi = mean(Intercept)
test_dirichlet_mean_intercept <- brm(
formula = cyl ~ 1 + mpg,
data = mtcars,
family = cumulative(probit),
stanvars = dirichlet_prior_stanvar,
prior = set_prior("induced_dirichlet(rep_vector(1, nthres+1), mean(Intercept))", class = "Intercept"),
cores = 4, seed = 4
)
# fit the model with induced_dirichlet() and phi = 0
test_dirichlet_zero <- brm(
formula = cyl ~ 1 + mpg,
data = mtcars,
family = cumulative(probit),
stanvars = dirichlet_prior_stanvar,
prior = set_prior("induced_dirichlet(rep_vector(1, nthres+1), 0)", class = "Intercept"),
cores = 4, seed = 4
)
A few quick plots suggest the default model is much better.
ce <- conditional_effects(test_standard)
p1 <- plot(ce, points = TRUE, plot = F)[[1]] +
labs(subtitle = "default t prior")
ce <- conditional_effects(test_dirichlet_mean_intercept)
p2 <- plot(ce, points = TRUE, plot = F)[[1]] +
labs(subtitle = "dirichlet w/ phi = mean(Intercept)")
ce <- conditional_effects(test_dirichlet_zero)
p3 <- plot(ce, points = TRUE, plot = F)[[1]] +
labs(subtitle = "dirichlet w/ phi = 0")
# combine
p1 | p2 | p3
The LOO confirms the plots.
test_standard <- add_criterion(test_standard, criterion = c("loo", "waic"))
test_dirichlet_mean_intercept <- add_criterion(test_dirichlet_mean_intercept, criterion = c("loo", "waic"))
test_dirichlet_zero <- add_criterion(test_dirichlet_zero, criterion = c("loo", "waic"))
loo_compare(test_standard, test_dirichlet_mean_intercept, test_dirichlet_zero) %>% print(simplify = F)
elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic se_looic
test_standard 0.0 0.0 -14.8 2.6 2.4 0.5 29.6 5.2
test_dirichlet_mean_intercept -19.4 3.6 -34.2 4.2 2.2 0.3 68.4 8.4
test_dirichlet_zero -25.2 3.1 -40.0 3.3 2.3 0.4 79.9 6.6
Where have I gone wrong?
