Are these jags model features possible in stan

[Stephen]

Dear Stan users,

Apologies in advance for my question, which is really more of a clarification.

I have a model written for jags that runs well but slowly. I want to use the model on a large (150,000+ cases) dataset and it will probably take too long with jags.

My question is, is there anything obvious in the model below that would not translate to stan?

---
model{
 
for(i in 1:N){

  w[i] ~ dbern(p[i])
  logit(p[i]) <- a0 + (b0 * x1[i])
  test[i] <- step(w[i] - 1) + 1
 
}

for(i in 1:N){

  y[i] ~ dnorm(y.hat[i], tau[test[i]])
  y.hat[i] <- a1[test[i]] + (b1[test[i]] * x2[i])
   
}
   
a0 ~ dnorm(0, 1e-3)
b0 ~ dnorm(0, 1e-3)
a1[2] ~ dnorm(0, 1e-3)
a1[1] ~ dnorm(0, 1e-3)
b1[2] ~ dnorm(0, 1e-3)
b1[1] ~ dnorm(0, 1e-3)
tau[1] <- pow(sigma[1], -2)
sigma[1] ~ dunif(0, 100)
tau[2] <- pow(sigma[2], -2)
sigma[2] ~ dunif(0, 100)

}
---

Specifically, the jags model is designed to allow for:

• missing values in the response variable w; all other variables are complete
  
• a step function or ifelse statement to switch between different coefficients in the linear regression
• interactions among categorical variables, with the possibility to set a contrast as 0

Are all these features available in stan? My reading suggests that they are but I need to be sure before attempting to rewrite my code.

Thanks in advance for any insight,

Stephen

[Bob Carpenter]

It helps to know what's data and what's parameters to figure out
how a JAGS model translates to Stan.

Stan doesn't have integer parameters, so you can't implement missing
w directly.

There is an if_else() in Stan, but also full conditional statements,
which are usually easier to read. And local variables. If you see
step() being used in a Stan model, it's usually just a straight translation
from BUGS/JAGS that can be improved.

Stan gives you a logit-scaled Bernoulli. So rather than writing

for (i in 1:N)
y[i] ~ bernoulli(inv_logit(theta[i]));

you can write

for (i in 1:N)
y[i] ~ bernoulli_logit(theta[i]);

and then once you've done that you can vectorize to

y ~ bernoulli_logit(theta);

Unless you expect parameters on the scale of 100, use tighter
priors than precision=1e-3 (sd=100). For comparison, Stan's
normal uses the standard deviation/scale parameterization, not
precision or variance (both precision and covariance matrices
available for multinormals).

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Stephen]

Hi Bob,

Thanks for your quick response. I attach the R code to run the model in jags.

The data in the model are: w, x1 and x2.

• w is the Bernoulli response variable and is missing for approximately 10% of cases
• x1 and x2 are covariates, which may be categorical or continuous and are complete

RE missing w: how might I go about accounting for missing w?

Thanks for all the other tips.

All the best,
Stephen

[Bob Carpenter]

For missing w it depends what you want to do with it.

If you just need to impute values, you can generate those in
the generated_quantities block, though it's more efficient
statistically to work in terms of expectations. It works like
this on the logit scale:

w_missing[i] <- bernoulli_logit_rng(...);

You'll also want to define test in the generated quantities block
as nothing depends on it in the model. It's more efficient and
clearer in terms of intent that way. The generated quantities block is
where predictive quantities go that don't impact estimation.

If you want the missing w to impact estimation, then you need to
marginalize out the w explicitly in the model. You basically get
a mixture model over the two outcomes, w[i] = 0 and w[i] = 1. There's
a chapter in the manual, and a whole new chapter with more examples coming
in Stan 2.5 (out real soon, as in hopefully the next month).

That is, where you currently have

PROD_i bernoulli(w[i] | inv_logit(a0 + b0 * x1[i]))

you need to marginalize out the w. But I don't see any reason you'd
want to do that as it won't impact the estimates in the model you have
because the two blocks in your model are completely separated and w doesn't
appear anywhere else.

- Bob

> <bern_norm.R>

[Stephen]

Hi Bob,

Thanks for your explanation. Are you referring to chapter 9 on Mixture Modeling in the Stan 2.4.0 manual?

Actually, I want estimation of w to impact the estimation of y: I want uncertainty in w estimation (particularly missing w) to propagate to uncertainty in y estimation. I'm not convinced I achieve this in the bugs model I shared, which returns very similar estimates to the "combined block" model below.

I'm want to apply this model to animal length data where w is a variable indicating whether the animal is a young of the year (w == 0) or not (w == 1) and y is animal length as a function of several possible environmental drivers.

I have two questions:
    - age is usually determined by length and so length is a predictor of age (w): is this OK? Perhaps there is circularity here that complicates estimation?
    - if the above is OK, then is mixture modeling using Stan (or MCMC generally) a suitable approach to this problem?
   
Thanks again for your help,
Stephen

## "combined block" model that returns very similar estimates to a model separating the dbern and dnorm estimations:

model{
  for(i in 1:N){
   

    Age[i] ~ dbern(p[i])
    logit(p[i]) <- a0 + b0[River[i]] + (c0 * Lgth[i])
   
    test[i] <- step(Age[i] - 1) + 1
   
#   }
#
#   for(i in 1:N){
   
    Length[i] ~ dnorm(y.hat[i], tau[test[i]])
    y.hat[i] <- a1[test[i]] + b1[test[i], River[i]]
   
  }

  a0 ~ dnorm(0, 1e-3)
  b0[1] <- 0
  c0 ~ dnorm(0, 1e-3)
  b1[1, 1] <- 0
  b1[2, 1] <- 0
  for(r in 2:N.rivers){
    b0[r] ~ dnorm(0, 1e-3)
    b1[1, r] ~ dnorm(0, 1e-3)
    b1[2, r] ~ dnorm(0, 1e-3)

  }
  a1[1] ~ dnorm(0, 1e-3)
  a1[2] ~ dnorm(0, 1e-3)

[Bob Carpenter]

Yes, that's a mixture model where w is the "indicator" variable
for the mixture. In Stan, w will have to be marginalized out
of the log probability. You will be able to recover the posterior
expectation of w for each animal, but you won't get samples.

I'm attaching the new manual chapter, which should provide some
more examples. It'll be out in Stan 2.5 and all of the models
in it should work in 2.4 as is.

Rather than

step(Age[i] - 1)

you can write

(Age[i] > 1)

in Stan.

As to the first question, you will indeed run into problems if you try
to do this:

u ~ regression based on v
v ~ regression based on u

To get that effect, you need to formulate a consistent joint model of (u,v).
The usual approach here, if I understand what you're trying to do, is
to just draw the mixture component from a Bernoulli and then based on its
outcome do a regression based on it as you did below. You don't want to also
do a regression for w based on the length. The effect of the length on the
choice of w still comes in indirectly because the nodes are connected in
the graphical model, so they're not independent in the posterior.

- Bob