# Looking for a stan program critique

**URL:** <https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865>\
**Category:** Modeling\
**Tags:** techniques\
**Created:** [April 13, 2018, 6:36pm UTC](https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865 "2018-04-13T18:36:15Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![statwonk](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/statwonk/32/937_2.png) [@statwonk](https://discourse.mc-stan.org/u/statwonk)\
**Post date:** [April 13, 2018, 6:36pm UTC](https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865/1 "2018-04-13T18:36:15Z")

</div>

[model.stan](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/3/39c4f4399cc0ffc2cf369f19516ef07c643cec34.stan) (2.2 KB)

Hi there! :wave:

I’m new to stan and I could use a bit of feedback on this model. Is there anything big I’m overlooking? Does something stand out to you?

```
data {
  int<lower=1> J; // number of apps
  int<lower=3> L; // number of years
  int<lower=1> N; // number of observations, the rows
  int K; // number of features, the columns
  int<lower=0> y[N]; // dependent variable index by the number of rows
  matrix[N, K] X; // trend covariates
  matrix[N, K] X2; // action covariates
  matrix[N, L] X3; // year covariates
  int IDs[N]; // apps
}
parameters {
  real mu_alpha;
  real<lower=0> tau_alpha; // hyper parameter for alpha variance
  // I split out the intercept from other paraters in the design matrix;
  vector[J] alpha; // the intercepts indexed by J apps
  vector[K] beta; // a feature parameter K in length
  vector[L] year; // a feature parameter L years in length
  real mu_year;
  real<lower=0> tau_year;
  vector<lower=0>[K] action_slope; // a feature parameter K in length
  vector<lower=0>[J] theta;
  real<lower=0> mu_theta;
  real<lower=0> tau_theta;
  real mu_action_slope;
  real<lower=0> tau_action_slope;
  real mu_beta;
  real<lower=0> tau_beta;
}
model {
  mu_year ~ normal(0.3, 0.5);
  tau_year ~ normal(0.8, 0.2);
  year ~ normal(mu_year, tau_year);

  mu_theta ~ normal(1, 5);
  tau_theta ~ normal(3.6, 0.25);
  theta ~ normal(mu_theta, tau_theta);

  mu_alpha ~ normal(1, 2);
  tau_alpha ~ normal(2.3, 0.5);
  alpha ~ normal(mu_alpha, tau_alpha);

  mu_beta ~ normal(0.07, 0.03);
  tau_beta ~ normal(0.07, 0.04);
  for(k in 1:K) {
    beta[k] ~ normal(mu_beta, tau_beta);
  }

  mu_action_slope ~ normal(0, 0.1);
  tau_action_slope ~ normal(0.12, 0.8);
  for(k in 1:K) {
    action_slope[k] ~ normal(mu_action_slope, tau_action_slope);
  }

  y ~ neg_binomial_2_log(alpha[IDs] + X * beta + X2 * action_slope + X3 * year, theta[IDs]);
}
generated quantities {
  real predicted_tasks[N];
  real marginal_tasks[N];
  matrix[N, K] X2_extra_actions;
  for(n in 1:N) {
    for(k in 1:K) {
      if(X2[n, k] > 0) {
        X2_extra_actions[n, k] = X2[n, k] + 1;
      } else {
        X2_extra_actions[n, k] = 0;
      }
    }
    predicted_tasks[n] = exp(alpha[IDs[n]] + X[n] * beta + (X2[n]) * action_slope + X3[n] * year);
    marginal_tasks[n] = exp(alpha[IDs[n]] + X[n] * beta + (X2_extra_actions[n]) * action_slope + X3[n] * year);
  }
}

```

A few things I’ve learned so far:

- Shorten iter to 1200 and set chain = 1 when prototyping
- Divergences can be caused by non-sensical priors
- Divergencies can occur when you need a hyperparameter but instead have hard-coded a parameter
- hyperparameters can be defined over a for loop of vectors
- use\_cache = FALSE with pars to avoid really slow wait times for print via summary.
- If there are divergences look for the params with low n\_eff and R\_hat \> 1 then focus on their prior, also consider if a hyperparameter is needed.
- bayesplot::ppc\_ecdf\_overlay is also a good way to find gross model errors that likely lead to divergences.
- One of the most important things to do is “center” or “scale” your data. Basically do what you can center and scaling wise such that the posterior can be as close to normal(0, 1) as possible. You want the posterior to be easy to explore and don’t want to be estimating any really large parameters like beta(1, 1e6) for example.

If I’ve got any of the above wrong corrections would be appreciated. Thank you.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [April 13, 2018, 8:07pm UTC](https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865/2 "2018-04-13T20:07:04Z")

</div>

> [@statwonk](#):
>
> - Divergences can be caused by non-sensical priors
> - Divergencies can occur when you need a hyperparameter but instead have hard-coded a parameter

Those are both instances of model misspecification.

> [@statwonk](#):
>
> One of the most important things to do is “center” or “scale” your data.

It’s the _parameters_ that should be as standard normal as possible in the posterior. You can sometimes do that by standardizing predictors and outcomes in a regression. You are not doing that in your model, which contains non-unit scale statements like this

```no-highlight
mu_beta ~ normal(0.07, 0.03);

```

To center and scale that, you’d write

```no-highlight
parameters {
  real mu_beta_raw;
  ...
transformed parameters {
  real mu_beta = 0.07 + 0.03 * mu_beta_raw;
  ...
model {
  mu_beta_raw ~ normal(0, 1);
  ...

```

You also want to vectorize those sampling statements, e.g.,

```no-highlight
beta ~ normal(mu_beta, tau_beta);

```

without the loop.

I take it `alpha[id]` is an exposure term. It can be more useful to use a time-series prior if you have data evolving over years than a regression on the year. The regression assumes the time effect is monotonic.

---

<div class="post-metadata">

**Author:** ![statwonk](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/statwonk/32/937_2.png) [@statwonk](https://discourse.mc-stan.org/u/statwonk)\
**Post date:** [April 13, 2018, 10:09pm UTC](https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865/3 "2018-04-13T22:09:43Z")

</div>

Thank you Bob. These are all very helpful clues and give me new avenues to learn / improve the model.

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [April 13, 2018, 11:25pm UTC](https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865/4 "2018-04-13T23:25:57Z")

</div>

> [@statwonk](#):
>
> I’m new to stan

Welcome on board! Great start with your model. I agree with Bob’s tips.

On the question of scaling or standardizing, I prefer the former since scaling by a constant doesn’t assume the mean and sd are known a priori.

In `rstanarm` the default is to do some rescaling of the prior standard deviations for regression coefficients based on observed standard deviations, but that’s just so we can make the default priors not terrible.

It’s basically an approximation to making the model hierarchical and learning the values from the data instead of assuming them. That can work well sometimes, but just be aware of the assumptions

Hope that helps!

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [April 14, 2018, 12:18am UTC](https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865/5 "2018-04-14T00:18:25Z")

</div>

> [@jonah](#):
>
> It’s basically an approximation to making the model hierarchical and learning the values from the data instead of assuming them.

Neat—hadn’t thought of it that way before.

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [April 14, 2018, 1:15am UTC](https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865/6 "2018-04-14T01:15:25Z")

</div>

In case it’s of interest, here’s a thread that has some other views on related topics:

> [@Standardizing predictors and outputs in a hierarchical model](http://discourse.mc-stan.org/t/standardizing-predictors-and-outputs-in-a-hierarchical-model/2974):
>
> Hello everyone, First of all thanks for the StanCon—it was a great experience. Here is a question that I have been discussing there and many suggested that I post it on the Forum. My question is about Manual Section 28.11 Standardizing Predictors and Outputs (v. 2.17.0), where it is suggested that we standardize our inputs as part of optimization. The model written is a very simple one predictor population variable model of y\_{n} = alpha + beta \* x\_{n} + epsilon\_{n}. I was wondering how to ba…

---

<div class="post-metadata">

**Author:** ![statwonk](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/statwonk/32/937_2.png) [@statwonk](https://discourse.mc-stan.org/u/statwonk)\
**Post date:** [April 14, 2018, 2:43pm UTC](https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865/7 "2018-04-14T14:43:47Z")

</div>

Thank you both! I’d say so far the parameter centering and choosing priors was the hardest part for me to figure out.

This final comment is really interesting: [Standardizing predictors and outputs in a hierarchical model](http://discourse.mc-stan.org/t/standardizing-predictors-and-outputs-in-a-hierarchical-model/2974/12)

Okay so my understanding is that my original model with the non-centered parameters given no divergences would yield appropriate inference as well? The main difference being computational time?

I ask because when I switched to centered I got a cryptic communication error with the sockets. Plan to debug next week. It smells like I blew out memory or something (need to check). If could trust the unscaled version I’d have a nice fallback.

---

<div class="post-metadata">

**Author:** ![statwonk](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/statwonk/32/937_2.png) [@statwonk](https://discourse.mc-stan.org/u/statwonk)\
**Post date:** [April 14, 2018, 4:31pm UTC](https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865/8 "2018-04-14T16:31:56Z")

</div>

So strange I’m hitting this error when I up the dataset using the `transform parameters` centering trick (it’s still fairly smallish N ~ 15k),

On Ubuntu 14.04, I get `Error in unserialize(socklist[[n]]) : error reading from connection` — sounds related to: [Unserialize error](http://discourse.mc-stan.org/t/unserialize-error/3670/5)

I wasn’t running into the same issue before centering the parameters.

---

<div class="post-metadata">

**Author:** ![statwonk](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/statwonk/32/937_2.png) [@statwonk](https://discourse.mc-stan.org/u/statwonk)\
**Post date:** [April 15, 2018, 2:11pm UTC](https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865/9 "2018-04-15T14:11:39Z")

</div>

I was able to verify that that the serialization error above comes is due to being out of memory. I used the `free` command to log memory usage and sure enough I ran out,

`nohup watch -n 1 'free -mh' > mem.log &`

I noticed in another thread that either Bob or Ben said `rstan` allocates memory up front while `CmdStan` does not (it streams?). I think I’ll give that a shot. I also might give the `parallel::mcapply` and `sampling()` route a shot: [https://github.com/stan-dev/rstan/issues/243#issuecomment-164560369](https://github.com/stan-dev/rstan/issues/243#issuecomment-164560369)

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [April 16, 2018, 2:53pm UTC](https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865/10 "2018-04-16T14:53:51Z")

</div>

You can also stream out in R. What I’ve never understood is whether you can turn off saving in memory.

You can also do fewer draws (rarely does anyone need more than a couple hundred effective sample size) or save fewer parameters.

---

<div class="post-metadata">

**Author:** ![statwonk](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/statwonk/32/937_2.png) [@statwonk](https://discourse.mc-stan.org/u/statwonk)\
**Post date:** [April 17, 2018, 5:12pm UTC](https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865/11 "2018-04-17T17:12:38Z")

</div>

I don’t think that’s true. I tried the `sample_file` arg on `stan` but memory still ballooned.

Since standardizing in `transform parameters` adds a significant amount of overhead and most of my priors are centered between 0 - 5 with standard deviations of +/- 2 does standardizing buy me much?

Could I move these back to the model block and specify non-centered priors with practically no loss given my scales aren’t too out of whack here? I know that if I do this I lose the nice interpretability of the standardized parameters.

I think if I can safely do that it could let me fit the model on the full larger dataset.

As an aside, I tried out an AR(1) process. I specified a hierarchical rho parameter bounded between -1 and 1. I was surprised that when the model returned the `n_eff` had dropped significantly for the rest of the covariates. I’m not yet exactly sure how to diagnose this type of issue but it’s one that I’m going to read up on.

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [April 17, 2018, 6:12pm UTC](https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865/12 "2018-04-17T18:12:48Z")

</div>

> [@statwonk](#):
>
> I don’t think that’s true. I tried the sample\_file arg on stan but memory still ballooned.

What if you use `sample_file` in combination with `pars=""`?

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [April 17, 2018, 6:48pm UTC](https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865/13 "2018-04-17T18:48:42Z")

</div>

> [@statwonk](#):
>
> Since standardizing in transform parameters adds a significant amount of overhead

Do it in the model block to have less overhead. You can do the same thing in transformed parameters or in the model block and only difference is whether the values are stored or not.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [April 17, 2018, 7:52pm UTC](https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865/14 "2018-04-17T19:52:23Z")

</div>

> [@statwonk](#):
>
> I tried out an AR(1) process. I specified a hierarchical rho parameter bounded between -1 and 1. I was surprised that when the model returned the n\_eff had dropped significantly for the rest of the covariates.

When you don’t have the constraint, is there any probability mass past the boundaries? If so, the hard boundary will likely introduce difficult geometry that will drive down the overall `n_eff` rate per iteration. With a hard boundary, it pushes the unconstrained lower or upper bound to plus or minus infinity, which is more difficult computationally to adapt for step size. It will also bias the fit compared to the unconstrained parameter, as all that mass that would be beyond the boundary is now piled up at the boundary, so any expectation will be further away from the boundary.

---

<div class="post-metadata">

**Author:** ![statwonk](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/statwonk/32/937_2.png) [@statwonk](https://discourse.mc-stan.org/u/statwonk)\
**Post date:** [April 17, 2018, 10:04pm UTC](https://discourse.mc-stan.org/t/looking-for-a-stan-program-critique/3865/15 "2018-04-17T22:04:04Z")

</div>

@jonah Using `sample_file` and `pars = ""`, I get the error message,

`no parameter ; sampling not done`.

@avehtari’s advice was a big help! I moved the centering process to the `model` block it cut down the memory footprint a good bit as expected.

> You can do the same thing in transformed parameters or in the model block and only difference is whether the values are stored or not.