Hello everyone, I am trying to fit a nester mixed effects model in R. I have tried to fit the model on my personal PC and on the PC at my working place. When I tried to run the MCMC with the “$sample” on the PC of my working place i got the error: “this user is not authorized to execute BIN/STANC.EXE”. Do you know why is this happening? How can someone in academia fit this model with the Cmdstan package if extra admin rights are required? When I run the code on my personal PC I do not get this error. Instead I get the warnings:
Init values were only set for a subset of parameters.
Missing init values for the following parameters:
- chain 1: y_gamma, y_z_beta, b_sd, z_b_vec, z_u_mat, u_sd, u_cholesky
- chain 2: y_gamma, y_z_beta, b_sd, z_b_vec, z_u_mat, u_sd, u_cholesky
Warning: 1 of 2000 (0.0%) transitions ended with a divergence.
See Runtime warnings and convergence problems for details.
and after fitting the model i get this message:
** Warning messages:**
1: The ESS has been capped to avoid unstable estimates.
2: The ESS has been capped to avoid unstable estimates.
3: The ESS has been capped to avoid unstable estimates
The version of the cmdstan is 2.29.2.
My Stan code is:
// Longitudinal sub-model for multicenter
/********************************************************/
functions {
vector evaluate_eta(matrix X, vector[] Z_u, int[] U_id, int[] C_id, real gamma,
vector beta, vector bVec, matrix uMat) {
int N = rows(X); // num rows in design matrix
int K = rows(beta); // num predictors
vector[N] eta;
if (K > 0) eta = X * beta + gamma;
else eta = rep_vector(0.0, N)+gamma;
for (n in 1:N)
eta[n] = eta[n] + bVec[C_id[n]] + uMat[U_id[n],1] * Z_u[1,n] + uMat[U_id[n],2] * Z_u[2,n];
return eta;
}
/**
* Get the indices corresponding to the lower tri of a square matrix
*
* @param dim The number of rows in the square matrix
* @return A vector of indices
*/
int[] lower_tri_indices(int dim) {
int indices[dim + choose(dim, 2)];
int mark = 1;
for (r in 1:dim) {
for (c in r:dim) {
indices[mark] = (r - 1) * dim + c;
mark = mark + 1;
}
}
return indices;
}
}
data {
//----- Longitudinal submodels
// population level dimensions
int<lower=0> y_N; // num observations
int<lower=0> y_K; // num predictors
// population level data
// NB these design matrices are evaluated at the observation times
vector[y_N] y; // response vectors
matrix[y_N,y_K] y_X; // fe design matrix
vector[y_K] y_Xbar; // predictor means
// group level dimensions
int<lower=0> b_N; // num center
int<lower=0> b_K; // total num params
int<lower=0> u_N; // num patients
int<lower=0> u_K; // total num params
int<lower=0> y_C_id[y_N]; // center ids
int<lower=0> y_U_id[y_N]; // patients ids
vector[y_N] y_Z[u_K]; // re design matrix
//----- Hyperparameters for prior distributions
// scale for priors
vector<lower=0>[y_K] y_prior_scale;
real<lower=0> y_prior_scale_for_intercept;
real<lower=0> y_prior_scale_for_aux;
// lkj prior stuff
real<lower=0> b_prior_scale;
vector<lower=0>[u_K] u_prior_scale;
vector<lower=0>[u_K] u_prior_df;
real<lower=0> u_prior_regularization;
}
transformed data {
// indexing used to extract lower tri of RE covariance matrix
int u_cov_idx[u_K + choose(u_K, 2)];
if (u_K > 0)
u_cov_idx = lower_tri_indices(u_K);
}
parameters {
real y_gamma; // intercept in long. submodel
vector[y_K] y_z_beta; // primitive coefs in long. submodels
real<lower=0> y_aux_unscaled; // unscaled residual error SDs
// group level params
real<lower=0> b_sd; // center level sds
vector [b_N] z_b_vec; // unscaled patient level params
matrix[u_K,u_N] z_u_mat; // unscaled group level params
vector<lower=0>[u_K] u_sd; // group level sds
cholesky_factor_corr[u_K > 1 ? u_K : 0]
u_cholesky; // cholesky factor of corr matrix
}
transformed parameters {
vector[y_K] y_beta; // population level params for long. submodels
real<lower=0> y_aux; // residual error SDs for long. submodels
matrix[u_N, u_K] u_mat; // group level params
vector[b_N] b_vec; //center level params
// coefs for long. submodels
y_beta = y_z_beta .* y_prior_scale;
// residual error SDs for long. submodels
y_aux = y_aux_unscaled * y_prior_scale_for_aux;
// group level params
// group level params
b_vec= b_sd * z_b_vec;
if (u_K == 1)
u_mat = (u_sd[1] * z_u_mat)';
else if (u_K > 1)
u_mat = (diag_pre_multiply(u_sd, u_cholesky) * z_u_mat)';
}
model {
//---- Log-lik for longitudinal submodels
vector[y_N] y_eta;
// declare linear predictors
y_eta = evaluate_eta(y_X, y_Z, y_U_id, y_C_id, y_gamma, y_beta, b_vec, u_mat);
// increment the target with the log-lik
target += normal_lpdf(y | y_eta, y_aux);
//----- Log-priors
// intercepts for long. submodels
target += normal_lpdf(y_gamma | 0, y_prior_scale_for_intercept);
// coefficients for long. submodels
target += normal_lpdf(y_z_beta | 0, 1);
// residual error SDs for long. submodels
target += normal_lpdf(y_aux_unscaled | 0, 1);
// group level terms
// sds
target += normal_lpdf(z_b_vec | 0, 1);
target += cauchy_lpdf(b_sd | 0, b_prior_scale);
target += student_t_lpdf(u_sd | u_prior_df, 0, u_prior_scale);
target += normal_lpdf(to_vector(z_u_mat) | 0, 1);
// corr matrix
if (u_K > 1)
target += lkj_corr_cholesky_lpdf(u_cholesky | u_prior_regularization);
}
generated quantities {
real y_alpha; // transformed intercepts for long. submodels
vector[size(u_cov_idx)] u_cov; // var-cov for REs
// Transformed intercepts for long. submodels
y_alpha = y_gamma - dot_product(y_Xbar,y_beta);
// Transform variance-covariance matrix for REs
if (u_K == 1)
u_cov[1] = u_sd[1] * u_sd[1];
else
u_cov = to_vector(quad_form_diag(
multiply_lower_tri_self_transpose(u_cholesky), u_sd))[u_cov_idx];
}
And the R code for fitting this model to my data is:
################################################################################
#Objective: 1) Generate data list; 2a) Fit simulated data by cmdstanr;
################################################################################
#Functions
# Function to return the range or SD of the predictors, used for scaling the priors
get_scale_value <- function(x) {
num.categories <- n_distinct(x)
x.scale <- 1
if (num.categories == 2) {
x.scale <- diff(range(x))
} else if (num.categories > 2) {
x.scale <- sd(x)
}
return(x.scale)
}
# Split the random effects part of a model formula into
# - the formula part (ie. the formula on the LHS of "|"), and
# - the name of the grouping factor (ie. the variable on the RHS of "|")
split_at_bars <- function(x) {
terms <- strsplit(deparse(x, 500), "\\s\\|\\s")[[1L]]
if (!length(terms) == 2L)
stop2("Could not parse the random effects formula.")
re_form <- formula(paste("~", terms[[1L]]))
group_var <- terms[[2L]]
nlist(re_form, group_var)
}
# Create a named list using specified names or, if names are omitted, using the
# names of the objects in the list
nlist <- function(...) {
m <- match.call()
out <- list(...)
no_names <- is.null(names(out))
has_name <- if (no_names) FALSE else nzchar(names(out))
if (all(has_name))
return(out)
nms <- as.character(m)[-1L]
if (no_names) {
names(out) <- nms
} else {
names(out)[!has_name] <- nms[!has_name]
}
return(out)
}
################################################################################
#1525 obs. with 300 patients who come from 6 centers
datLong=readRDS("U://simDatLong.rds")
################################################################################
library(dplyr)
library(lme4)
#Generate datalist----
formulaLong=y~time+I(time^2)+(1|centerID)+(1+time|id)
lmer.fit=lme4::lmer(formula=formulaLong, data=datLong)
model_frame <- stats::model.frame(lme4::subbars(formulaLong),datLong)
x_form <- lme4::nobars(formulaLong)
x <- model.matrix(x_form,model_frame)
predictors.X <- x[,-1] #drop intercept
y_Xbar <- colMeans(predictors.X)
y_X.scl <- sweep(predictors.X, 2, y_Xbar, FUN = "-")
y_prior_scale0=2.5*sd(datLong$y)
min_prior_scale=1e-12
y_prior_scale=pmax(min_prior_scale, y_prior_scale0 / apply(predictors.X, 2L, get_scale_value))
u_prior_scale0=b_prior_scale0=10*sd(datLong$y)
bars <- lme4::findbars(formulaLong)
z_parts <- lapply(bars,split_at_bars)
z_forms <- lapply(z_parts, `[[`, 're_form')
z <- lapply(z_forms,model.matrix,model_frame)
names(z) <- lapply(z_parts,`[[`,'group_var')
u_prior_scale=pmax(min_prior_scale, u_prior_scale0 / apply(z[['id']], 2L, get_scale_value))
b_prior_scale=10 #manually assign by better performance (tested)
standata.long=list(
y_N=nrow(datLong),
y_K=ncol(predictors.X), # exclude intercept
y=as.vector(datLong$y),
y_X=y_X.scl, # exclude intercept+scale
y_Xbar=y_Xbar,
b_N=length(unique(datLong$centerID)), #number of centers
b_K=ncol(z$centerID), #intercept
u_N=length(unique(datLong$id)), # number of patients
u_K=ncol(z$id), #intercept + slope
y_C_id=datLong$centerID,
y_U_id=datLong$id,
y_Z=t(z$id), # b_K times y_N matrix
y_prior_scale=y_prior_scale, # exclude intercept
y_prior_scale_for_intercept=10*sd(datLong$y),
y_prior_scale_for_aux=5*sd(datLong$y),
b_prior_scale=b_prior_scale,
u_prior_scale=u_prior_scale,
u_prior_df=rep(1,ncol(z$id)),
u_prior_regularization=2
)
#Fit by cmdstanr----
##Install cmdstanr
#install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
library(cmdstanr)
check_cmdstan_toolchain()
library(posterior)
library(bayesplot)
install_cmdstan(cores = 2)
file.long <- file.path("U://LONG.stan")
mod.long <- cmdstan_model(file.long)
#Default: num_sample=1000 per chain post-warmup (num_warmup=1000)
fit.long <- mod.long$sample(
data = standata.long,
chains = 2,
save_warmup = FALSE,
parallel_chains = 2,
refresh = 500,
adapt_delta=0.95,
max_treedepth=12,
seed=2022,
init=function() list(y_aux_unscaled=0.2)
)```