Problem with MCMC execution ($sample) of cmdstanr and admin rights

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)
)```

Some workplaces place limitations on the ability of users to run arbitrary executables, only allowing them to run white-listed programs. I believe this is the case here.

You’ll have to either contact your IT department to request permissions for running executables, or you’ll have to change your code to use rstan instead of cmdstanr, since that doesn’t require calling any programs outside of R

1 Like