Need some advice on how to reduce or eliminate the divergent transitions with this model. I’m not sure if they even matter as the results are fairly sensible. This is a simple logistic regression model with the outcome whether a health insurance claim was deemed pre-existing or not(approved)PEC_notrauma.txt (34.5 KB) ,and the predictors are age,gender,type of policy and diagnosis.
I can minimise the divergency by cranking up the adapt_delta to 0.999 but this doesn’t seem reasonable.
Tried various priors and not sure how to “re-parameterize” the model.
stanDat <- list(a = as.integer(factor(PEC_notrauma$gender)),
b = as.integer(factor(PEC_notrauma$policy)),
g = as.integer(factor(PEC_notrauma$Diag)),
y=PEC_notrauma$outcome,
N = nrow(PEC_notrauma),
a_N=length(unique(PEC_notrauma$gender)),
b_N=length(unique(PEC_notrauma$policy)),
g_N=length(unique(PEC_notrauma$Diag)),
c=as.matrix(PEC_notrauma$age),#have to store this as a matrix to use this code
#could be handy if more than one continuous predictor
c_N=1)
model_string<-
"data {
int N;
int a_N; // number of categories in gender
int b_N; // number of categories in policy
int g_N; // number of categories in diagnosis
int c_N; //number of continuous covariates(only 1 here)
int<lower=1, upper=a_N> a[N]; // number applied to gender
int<lower=1, upper=b_N> b[N]; //number applied to policy
int<lower=1,upper=g_N> g[N]; //number applied to diagnosis
matrix[N,c_N] c; //observed ages
int<lower=0, upper=1> y[N]; //observed outcome(deemed PEC or approved)
}
parameters {
real intercept;
vector[a_N] a_coeff;
vector[b_N] b_coeff;
vector[g_N] g_coeff;
vector[c_N] c_coeff;
real<lower=0> sigma_a;
real<lower=0> sigma_b;
real<lower=0> sigma_g;
real<lower=0> sigma_c;
real mu_a;
real mu_b;
real mu_g;
real mu_c;
}
transformed parameters {
vector[a_N] normalized_a_coeff;
vector[b_N] normalized_b_coeff;
vector[g_N] normalized_g_coeff;
normalized_a_coeff = a_coeff - a_coeff[1];
normalized_b_coeff = b_coeff - b_coeff[2];
normalized_g_coeff = g_coeff - g_coeff[1];
}
model {
vector[N] p; // probabilities
// priors
intercept ~ normal(0,10);//was normal(0,100)
sigma_a ~ normal(0,1); //all the sigmas were cauchy(0,5)
sigma_b ~ normal(0,1);
sigma_g ~ normal(0,1);
sigma_c ~ normal(0,1);
mu_a ~ normal(0,1); // can leave these out and just have 0 as mean of coefficients
mu_b ~ normal(0,1);
mu_g ~ normal(0,1);
mu_c ~ normal(0,1);
// level 1
a_coeff ~ normal(mu_a, sigma_a);
b_coeff ~ normal(mu_b, sigma_b);
g_coeff ~ normal(mu_g, sigma_g);
c_coeff ~ normal(mu_c, sigma_c);
// level 2
for (i in 1:N) {
p[i] = a_coeff[a[i]] + b_coeff[b[i]] + g_coeff[g[i]];
}
y ~ bernoulli_logit(intercept + p + c * c_coeff);//can use poisson_log here;get different parameter estimates
}
generated quantities{
vector[N] y_rep;
for ( n in 1:N ) {
y_rep[n]=bernoulli_logit_rng(intercept + a_coeff[a[n]] + b_coeff[b[n]] + g_coeff[g[n]]+ c[n] * c_coeff);
}
}
"
stan_samples <- stan(model_code = model_string, data = stanDat,warmup=1000,iter=4000,chains=4,
control = list(adapt_delta = 0.95,stepsize=0.01,max_treedepth=15))
I removed one of the diagnostic categories which was a perfect predictor(trauma,all approved).
Regards
Chris
