Yes of course. I used the Stan models that @fabio provided earlier in this thread.
library(rstan)
library(projpred) # 2.0.0 devtools::install_github("stan-dev/projpred@develop")
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
set.seed(123)
N <- 500 # number of observations
M <- 30 # number of predictors
K <- 5 # number of ordinal categories
c <- sort(runif(n = K - 1, min = -5, max = 5))
generative <- stan_model("./generate_data.stan")
gen <- sampling(generative,
data = list(N = N, M = M, K = K, c = c),
seed = 123, algorithm = "Fixed_param", iter = 1, chains = 1
)
beta <- as.numeric(extract(gen, pars = "beta")$beta)
X <- matrix(
data = as.numeric(extract(gen, pars = "X")$X),
nrow = N, ncol = M
)
y <- as.integer(extract(gen, pars = "y")$y)
model <- stan_model("./ordinal_regression.stan")
fit <- sampling(model,
data = list(N = N, M = M, K = K, X = X, y = y),
seed = 123, chains = 4
)
mu <- extract(fit, pars = "mu")
Emu <- as.matrix(colMeans(mu$mu))
beta_out <- extract(fit, pars = "beta")$beta
beta <- as.matrix(colMeans(beta_out))
predfun <- function(fit, newdata = NULL) {
if (is.null(newdata)) {
newdata <- data
}
return(as.matrix(newdata[, -1]) %*% t(beta_out))
}
## compute empirical sd; too large
sigma <- apply(mu$mu, 1, sd)
colnames(X) <- paste0("x", seq_len(NCOL(X)))
data <- as.data.frame(cbind(y = Emu, x = X))
formula <- as.formula(paste("y ~", paste(colnames(X), collapse = " + ")))
.extract_model_data <- function(object, newdata = NULL, wrhs = NULL,
orhs = NULL, resp_form = NULL) {
if (is.null(newdata)) {
newdata <- data
}
if (inherits(wrhs, "formula")) {
weights <- eval_rhs(wrhs, newdata)
} else if (is.null(wrhs)) {
weights <- rep(1, NROW(newdata))
}
if (inherits(orhs, "formula")) {
offset <- eval_rhs(orhs, newdata)
} else if (is.null(orhs)) {
offset <- rep(0, NROW(newdata))
}
if (inherits(resp_form, "formula")) {
y <- eval_rhs(resp_form, newdata)
} else {
y <- NULL
}
return(nlist(y, weights, offset))
}
extract_model_data <- function(object, newdata = NULL, wrhs = NULL,
orhs = NULL) {
resp_form <- lhs(formula)
args <- nlist(object, newdata, wrhs, orhs, resp_form)
return(do_call(.extract_model_data, args))
}
ref <- init_refmodel(fit, data, formula, gaussian(),
ref_predfun = predfun, div_minimizer = linear_mle,
proj_predfun = linear_proj_predfun,
extract_model_data = extract_model_data, dis = sigma
)
vs <- varsel(ref, method = "forward")
cvs <- cv_varsel(ref, method = "forward")
## fix sigma
sigma <- rep(1.6, 4000)
ref <- init_refmodel(fit, data, formula, gaussian(),
ref_predfun = predfun, div_minimizer = linear_mle,
proj_predfun = linear_proj_predfun,
extract_model_data = extract_model_data, dis = sigma
)
vs <- varsel(ref, method = "forward")
cvs <- cv_varsel(ref, method = "forward")