What happens when the iteration progress stops/freezes?


#1

I was running a Stan model (on AWS). It got up to 400/2000 iterations and then stopped there for 30+ hours. What could cause this?

Thanks!


#2

Most likely, because this is during adaptation, the stepsize was adjusted to be << 1e-6 and it just looks like the model froze. The best way to avoid this is to re-parameterize the model but there’s not generic advice for re-parameterizing that always leads to success.

You can check this by writing adaptation output to a file (the save_warmup=1 parameter in CmdStan, something similar in rstan), and loading the .csv file that gets saved. It should have a stepsize__ column.


#3

Thanks! I should mention, when I ran the model the first time I used the defaults and got this warning:
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
so I increased to max_treedepth = 15.

The model I’m running is a slight modification of the one in this vignette:
http://mc-stan.org/users/documentation/case-studies/icar_stan.html


#4

could you share the modified code? there’s a few models in that vignette…


#5

Sure!

We want to estimate rates of y = 1 at the level of HSA (hospital service area). HSAs are nested into HRRs (hospital referral regions).

Thanks so much!

p.s. what’s the best way to format code in these message boards so it isn’t awful to look at? Mitzi answered!

data {
  // control flags
  int<lower=0, upper=1> t_prior; // 0=normal, 1=t
  int<lower=0, upper=1> has_spatial; // 0=false, 1=true
  int<lower=0, upper=1> has_covar; // 0=false, 1=true
  
  // array sizes
  int N_prac;
  int N_hsa;
  int N_hrr;
  int N_covar;
  int N_edges_hrr;
  int N_edges_hsa;
  
  // index information for group membership
  int i_hsa[N_prac];
  int i_hrr[N_prac];
  
  // data
  matrix[N_prac * has_covar, N_covar] X;
  int y[N_prac];
  
  // spatial/adjacency
  int node_hsa1[N_edges_hsa * has_spatial];
  int node_hsa2[N_edges_hsa * has_spatial];
  int node_hrr1[N_edges_hrr * has_spatial];
  int node_hrr2[N_edges_hrr * has_spatial];
}
transformed data {
  matrix[N_prac * has_covar, N_covar] Q_ast;
  matrix[N_covar, N_covar] R_ast;
  matrix[N_covar, N_covar] R_ast_inverse;
  if (has_covar) {
    // thin and scale the QR decomposition
    Q_ast = qr_Q(X)[, 1:N_covar] * sqrt(N_prac - 1);
    R_ast = qr_R(X)[1:N_covar, ] / sqrt(N_prac - 1);
    R_ast_inverse = inverse(R_ast);
  }
}
parameters {
  real alpha; // intercept
  vector[N_covar * has_covar] theta; // slopes of Q_ast (QR decomposition)
  vector[N_hsa] a_j_raw; // random hsa effects
  vector[N_hrr] b_k_raw; // random hrr effects
  vector[N_hsa * has_spatial] phi_j_std; // spatial hsa effects
  // vector[N_hrr * has_spatial] psi_k_std; // spatial hrr effects
  real<lower=0> tau_a; // precision of random hsa effects
  real<lower=0> tau_b; // precision of random hrr effects
  real<lower=0> tau_phi; // precision of spatial hsa effects
  //real<lower=0> tau_psi; // precision of spatial hrr effects
}
transformed parameters {
  vector[N_prac] eta_i; // linear predictor
  vector[N_hsa] a_j; 
  vector[N_hrr] b_k;  
  vector[N_hsa * has_spatial] phi_j;
  // vector[N_hrr * has_spatial] psi_k;
  real<lower=0> sigma_a;
  real<lower=0> sigma_b;
  real<lower=0> sigma_phi;
  // real<lower=0> sigma_psi;
  
  // convert precisions to sigmas
  sigma_a = inv(sqrt(tau_a));
  sigma_b = inv(sqrt(tau_b));
  sigma_phi = inv(sqrt(tau_phi));
  // sigma_psi = inv(sqrt(tau_psi));
  
  // matt trick random effects
  a_j = a_j_raw * sigma_a;
  b_k = b_k_raw * sigma_b;
  // matt trick spatial effects
  phi_j = phi_j_std * sigma_phi;
  // psi_k = psi_k_std * sigma_psi;
    
  // linear predictor
  eta_i = alpha + a_j[i_hsa] + b_k[i_hrr];
  if (has_covar)
    eta_i += Q_ast * theta;
  if (has_spatial)  
    eta_i += phi_j[i_hsa]; // += psi_k[i_hrr] + phi_j[i_hsa];
}
model {
  for (i in 1:N_prac) y[i] ~ bernoulli(inv_logit(eta_i[i])); // could vectorize?
  
  // weakly-informative priors
  alpha ~ normal(0,1);
  theta ~ normal(0,1);
  
  // raw for random effects
  if (t_prior == 0) {
    a_j_raw ~ normal(0, 1);
    b_k_raw ~ normal(0, 1);
  } else if (t_prior == 1) {
    a_j_raw ~ student_t(4, 0, 1);
    b_k_raw ~ student_t(4, 0, 1);
  }
  
  // Carlin WinBUGS priors on spatial precision hyperparameters
  tau_phi ~ gamma(1, 1);
  // tau_psi ~ gamma(1, 1);
  
  // spatial effects
  if (has_spatial) {
    // ICAR prior
    target += -0.5 * dot_self(phi_j_std[node_hsa1] - phi_j_std[node_hsa2]);
    // target += -0.5 * dot_self(psi_k_std[node_hrr1] - psi_k_std[node_hrr2]);
    // soft sum-to-zero constraint
    sum(phi_j_std) ~ normal(0, 0.001 * N_hsa);
    // sum(psi_k_std) ~ normal(0, 0.001 * N_hrr);
    
    // Carlin WinBUGS priors on random effects precision hyperparameters
    tau_a ~ gamma(3.2761, 1.81);
    tau_b ~ gamma(3.2761, 1.81);
  } else {
    tau_a ~ normal(0,1);
    tau_b ~ normal(0,1);
  }
}
generated quantities {
  vector[N_covar] beta;
  if (has_covar) {
    beta = R_ast_inverse * theta; // coefficients on X
  }
}

#6

have your tried using different priors? the Carlin WINBUGS priors help Gibbs samplers (?) but aren’t appropriate for HMC.


#7

set the code by putting a line consisting of three backticks “```” before and after your code.
here’s the formatted code:

data {
  // control flags
  int<lower=0, upper=1> t_prior; // 0=normal, 1=t
  int<lower=0, upper=1> has_spatial; // 0=false, 1=true
  int<lower=0, upper=1> has_covar; // 0=false, 1=true

  // array sizes
  int N_prac;
  int N_hsa;
  int N_hrr;
  int N_covar;
  int N_edges_hrr;
  int N_edges_hsa;

  // index information for group membership
  int i_hsa[N_prac];
  int i_hrr[N_prac];

  // data
  matrix[N_prac * has_covar, N_covar] X;
  int y[N_prac];

  // spatial/adjacency
  int node_hsa1[N_edges_hsa * has_spatial];
  int node_hsa2[N_edges_hsa * has_spatial];
  int node_hrr1[N_edges_hrr * has_spatial];
  int node_hrr2[N_edges_hrr * has_spatial];
}
transformed data {
  matrix[N_prac * has_covar, N_covar] Q_ast;
  matrix[N_covar, N_covar] R_ast;
  matrix[N_covar, N_covar] R_ast_inverse;
  if (has_covar) {
    // thin and scale the QR decomposition
    Q_ast = qr_Q(X)[, 1:N_covar] * sqrt(N_prac - 1);
    R_ast = qr_R(X)[1:N_covar, ] / sqrt(N_prac - 1);
    R_ast_inverse = inverse(R_ast);
  }
}
parameters {
  real alpha; // intercept
  vector[N_covar * has_covar] theta; // slopes of Q_ast (QR decomposition)
  vector[N_hsa] a_j_raw; // random hsa effects
  vector[N_hrr] b_k_raw; // random hrr effects
  vector[N_hsa * has_spatial] phi_j_std; // spatial hsa effects
  // vector[N_hrr * has_spatial] psi_k_std; // spatial hrr effects
  real<lower=0> tau_a; // precision of random hsa effects
  real<lower=0> tau_b; // precision of random hrr effects
  real<lower=0> tau_phi; // precision of spatial hsa effects
  //real<lower=0> tau_psi; // precision of spatial hrr effects
}
transformed parameters {
  vector[N_prac] eta_i; // linear predictor
  vector[N_hsa] a_j;
  vector[N_hrr] b_k;
  vector[N_hsa * has_spatial] phi_j;
  // vector[N_hrr * has_spatial] psi_k;
  real<lower=0> sigma_a;
  real<lower=0> sigma_b;
  real<lower=0> sigma_phi;
  // real<lower=0> sigma_psi;

  // convert precisions to sigmas
  sigma_a = inv(sqrt(tau_a));
  sigma_b = inv(sqrt(tau_b));
  sigma_phi = inv(sqrt(tau_phi));
  // sigma_psi = inv(sqrt(tau_psi));

  // matt trick random effects
  a_j = a_j_raw * sigma_a;
  b_k = b_k_raw * sigma_b;
  // matt trick spatial effects
  phi_j = phi_j_std * sigma_phi;
  // psi_k = psi_k_std * sigma_psi;

  // linear predictor
  eta_i = alpha + a_j[i_hsa] + b_k[i_hrr];
  if (has_covar)
    eta_i += Q_ast * theta;
  if (has_spatial)
    eta_i += phi_j[i_hsa]; // += psi_k[i_hrr] + phi_j[i_hsa];
}
model {
  for (i in 1:N_prac) y[i] ~ bernoulli(inv_logit(eta_i[i])); // could vectorize?

  // weakly-informative priors
  alpha ~ normal(0,1);
  theta ~ normal(0,1);

  // raw for random effects
  if (t_prior == 0) {
    a_j_raw ~ normal(0, 1);
    b_k_raw ~ normal(0, 1);
  } else if (t_prior == 1) {
    a_j_raw ~ student_t(4, 0, 1);
    b_k_raw ~ student_t(4, 0, 1);
  }

  // Carlin WinBUGS priors on spatial precision hyperparameters
  tau_phi ~ gamma(1, 1);
  // tau_psi ~ gamma(1, 1);

  // spatial effects
  if (has_spatial) {
    // ICAR prior
    target += -0.5 * dot_self(phi_j_std[node_hsa1] - phi_j_std[node_hsa2]);
    // target += -0.5 * dot_self(psi_k_std[node_hrr1] - psi_k_std[node_hrr2]);
    // soft sum-to-zero constraint
    sum(phi_j_std) ~ normal(0, 0.001 * N_hsa);
    // sum(psi_k_std) ~ normal(0, 0.001 * N_hrr);

    // Carlin WinBUGS priors on random effects precision hyperparameters
    tau_a ~ gamma(3.2761, 1.81);
    tau_b ~ gamma(3.2761, 1.81);
  } else {
    tau_a ~ normal(0,1);
    tau_b ~ normal(0,1);
  }
}
generated quantities {
  vector[N_covar] beta;
  if (has_covar) {
    beta = R_ast_inverse * theta; // coefficients on X
  }
}

working through your model, hope to have further comments soon.


#8

I haven’t yet, no! What kind of priors do you recommend for HMC?

Thanks so much!


#9

Ah, I realized that we haven’t yet tried the BYM2 parametrization (again, see Mitzi’s awesome write-up: http://mc-stan.org/users/documentation/case-studies/icar_stan.html). Mitzi, do you think that would help a lot here?


#10

That’s a good point. I wonder how hard it would be to change the code to have a minimum stepsize that you can’t go below. That way instead of this “freezing” or slowind down to a crawl behavior happening you just get divergences or lower accept rates. I think this would “fail faster” and be a more informative clue to the user that they should re-parameterize.

I should look into how hard that would be to do. I think @Bob_Carpenter mentioned that one of the upcoming releases of Stan will allow you to specify the stepsize manually.


#11

My preference is to deal with the output s.t. we can make online monitoring easier. I have an abstract due so I didn’t get a chance to work on it this weekend like I wanted to … :( Would be pretty easy to fire a warning if the stepsize drops, but I’ve also had models recover (b/c of mass matrix scaling adaptation, I think). Even currently you can load a partial .csv file and plotting log10(stepsize__) is a solid way of spotting the problem… we just don’t have tooling in place to make it easy for users.