Will STAN address these modelling and sampling limitations in WinBugs

[Josh Korman]

Hello,

 

I am a WinBugs user that is very excited about using STAN based on what I have read in the manual. Before I jump in, I'm wondering if STAN can deal with two problems I have encountered in WinBugs due to limitations of its modelling environment:

 

1) 'multiple definitions of node'. This error occurs when I calculate a parameter as a function of data and other parameters as a logical expression, but also assume the parameter is part of a stochastic distribution (say the order of a binomial). I think this limitation is due to the DAG structure of WinBugs, and I'm hoping STAN eliminates this issue.

 

 

2) 'unable to choose update method'. According to WB documentation, this occurs if more than one binomial response is assumed to share the same order parameter. From the documentation, this has something to do with the sampler, and I'm hoping the substantially revised sampler in STAN will no longer have this limitation.

 

Will STAN get around these problems?

 

Any thoughts much appreciated,

 

Josh Korman

 

[Ben Goodrich]

On Thursday, October 18, 2012 12:22:11 PM UTC-4, Josh Korman wrote:

I am a WinBugs user that is very excited about using STAN based on what I have read in the manual. Before I jump in, I'm wondering if STAN can deal with two problems I have encountered in WinBugs due to limitations of its modelling environment:

I haven't used WinBUGS much in a long time, so I am guessing. But I am guessing that these are not going to be a problem in Stan.

 

 1) 'multiple definitions of node'. This error occurs when I calculate a parameter as a function of data and other parameters as a logical expression, but also assume the parameter is part of a stochastic distribution (say the order of a binomial). I think this limitation is due to the DAG structure of WinBugs, and I'm hoping STAN eliminates this issue.

Pretty sure this is not going to be a problem. See appendix C of the Stan manual. The code in the .stan file is parsed to a C++ file and works as would be expected. So, you can define a variable and then redefine it later, but you can't use it before it is defined and allocated. Stan doesn't do DAG structures; everything hinges on the textual representation of the model in the .stan file.

 

2) 'unable to choose update method'. According to WB documentation, this occurs if more than one binomial response is assumed to share the same order parameter. From the documentation, this has something to do with the sampler, and I'm hoping the substantially revised sampler in STAN will no longer have this limitation.

We only have 3 variations on the HMC sampler, but the default one is usually good enough because it is adaptive. You won't get an error message like this, but whether you will get a large number of effective samples from that posterior distribution remains to be seen. 

It is hard to say anything for sure given the particularities of WinBUGS and your data. Just read the Stan manual and the FAQ thread, give Stan a try, and ask here if you have more questions.

Ben

[Bob Carpenter]

On 10/18/12 12:22 PM, Josh Korman wrote:
> Hello,
> I am a WinBugs user that is very excited about using STAN based on what I have read in the manual. Before I jump in, I'm
> wondering if STAN can deal with two problems I have encountered in WinBugs due to limitations of its modelling environment:

> 1) 'multiple definitions of node'. This error occurs when I calculate a parameter as a function of data and other
> parameters as a logical expression, but also assume the parameter is part of a stochastic distribution (say the order of
> a binomial). I think this limitation is due to the DAG structure of WinBugs, and I'm hoping STAN eliminates this issue.

It would help if you gave us a concrete example of what
you want to do.

You can usually work around this problem indirectly in WinBUGS itself
by defining auxiliary nodes.

Stan does not permit assignments to parameters.

Transformed parameters may be assigned and reassigned.
Local variables can also be assigned and reassigned.
So it would be OK in Stan to do this:

model {
real sum;

sum <- 0;
for (n in 1:N) sum <- sum + x[n];
}

Stan also lets you put anything on the left-hand side of
a sampling statement, but then it's up to you to add the
log Jacobian to lp__.

model {
sqrt(1/tau) ~ uniform(0,100);
lp__ <- lp__ + log(

}


> 2) 'unable to choose update method'. According to WB documentation, this occurs if more than one binomial response is
> assumed to share the same order parameter.

You won't get that error in Stan because we can
sample from anything that you can express in the
language. WinBUGS has this problem in some cases
if one of its sampling strategies (see the manual
for the four or five that are used in different cases)
cannot be applied to the structure of the problem.

I'm again not sure what exactly you are trying to
do from this description. An example would help.

Stan does not have discrete parameters. So you
can't have n as a parameter in:

n ~ binomial(N,theta);

- Bob

[Sebastian Weber]

Hello,

so stan could for example model modeled parameters? For example, suppose I fit a fixed intercept, varying slope model to J groups (J large). I would then be interested in regressing the slope as a logistic regression of some other covariate. In my case, I fit longitudinal data and summarize the kinetics of subjects by a slope. These subjects are under a treatment of different concentrations of whatever and I want to find out how the treatment concentration varies with the treatment effect which I define as slope larger than a threshold. So in short, I want to regress estimated parameters from a model.

Up to now I simply estimate the kinetics parameters, draw posterior samples, and fit again. But this is kind of awkward to me, there should be a direct way of doing so. In BUGS I always ran into the DAG problem mentioned of redefining nodes again...

Best,

Sebastian

> --
>
>

--
Sebastian Weber
Group of Cell Communication and Control
Freiburg Institute for Advanced Studies - FRIAS
School of Life Sciences - LIFENET
Albert-Ludwigs-Universität Freiburg
Albertstr. 19
79104 Freiburg
T.: +49-761-203-97237
Fax:+49-761-203-97334

[Josh Korman]

Hi Bob,

Thanks for your reply.

I have included a simplified example that gets at my multiple definition
problem. Here I want to estimate abundance from two population (N[]) and
assume that each of these estimates comes from a hyper distribution with a
fixed value for tau. The objective is to estimate the overall abundance
(SuperN) as an information weighted average of the two estimates conditional
on the fixed value of tau. This works fine using the uncommented lines
(expu<- and u~dpois) below. Problem is that in my more complex model, N is a
function of other estimated parameters so need to do something like the two
commented-out lines, which won't work in WinBugs because N is defined twice
(as a logical expression and as a stochastic node).

As for as you second point, you say that STAN cannot estimate discrete
parameters. Would N from u~dbin(theta, N) be estimable if you allow N to be
continuous rather than discrete? Or is the order of a binomial (e.g. N in
this example) just not estimable in STAN?

Thanks,

Josh

Example CODE

r, M, and u are data, p,N, and superN are parameters to be estimated.

model
{
mu~dnorm(0,1.0E-6) #hyper prior on log of single
escapement estimate

for (i in 1:2)
{
p[i]~dbeta(1,1) #uninformative prior n p
r[i]~dbin(p[i],M[i]) #binomial likelihood to estimate
capture probability p given recaps (r) and marks release (M)

logN[i]~dnorm(mu,tau) #log N for each population comes
from a common (super) distribution)
N[i]<-exp(logN[i]) #transform

#The following two lines work as expected
expu[i]<-N[i]*p[i] #expected number of unmarked fish
captured on 2nd pass given N and p
u[i]~dpois(expu[i]) #This makes N[i] a stochastic node
which allows HBM to work

#Replacing the previous 2 lines with the following 2 does
not work due to multiple definitions of node N[]
#u[i]~dbin(p[i],N[i])
#N[i]<-u[i]/p[i]

}
SuperN<-exp(mu) #This is the single estimate based
on information across populations

--

[Bob Carpenter]

On 10/18/12 1:21 PM, Sebastian Weber wrote:
> Hello,
>
> so stan could for example model modeled parameters? For example, suppose I fit a fixed intercept, varying slope model to J groups (J large). I would then be interested in regressing the slope as a logistic regression of some other covariate. In my case, I fit longitudinal data and summarize the kinetics of subjects by a slope. These subjects are under a treatment of different concentrations of whatever and I want to find out how the treatment concentration varies with the treatment effect which I define as slope larger than a threshold. So in short, I want to regress estimated parameters from a model.

> Up to now I simply estimate the kinetics parameters, draw posterior samples, and fit again. But this is kind of awkward to me, there should be a direct way of doing so. In BUGS I always ran into the DAG problem mentioned of redefining nodes again...

I may be missing something here, but why not just build another level
into the model?

The regression on the slope predictors would give you a mean around which to
draw an individual's slope.

// u[j] predictor for group j's slope

slope[j] ~ normal(alpha * u[j] + beta, sigma_slope); // slope for group j


// jj[n] is group for prediction n, ii[n] is index for predictor for predictor n

y[n] ~ normal(slope[jj[n]] * x[ii[n]] + gamma, sigma_y);

- Bob

[Bob Carpenter]

On 10/18/12 2:22 PM, Josh Korman wrote:
...

> I have included a simplified example that gets at my multiple definition
> problem. Here I want to estimate abundance from two population (N[]) and
> assume that each of these estimates comes from a hyper distribution with a
> fixed value for tau. The objective is to estimate the overall abundance
> (SuperN) as an information weighted average of the two estimates conditional
> on the fixed value of tau. This works fine using the uncommented lines
> (expu<- and u~dpois) below. Problem is that in my more complex model, N is a
> function of other estimated parameters so need to do something like the two
> commented-out lines, which won't work in WinBugs because N is defined twice
> (as a logical expression and as a stochastic node).
>
> As for as you second point, you say that STAN cannot estimate discrete
> parameters. Would N from u~dbin(theta, N) be estimable if you allow N to be
> continuous rather than discrete? Or is the order of a binomial (e.g. N in
> this example) just not estimable in STAN?

The function signatures for the distributions are in the
manual -- binomial(N,theta) requires N to be an integer.
(We're not using the R notation dbin() for pmfs and pdfs).

Unfortunately, the discretizing transforms like floor()
lose derivative information because they're flat step
functions. So while we allow them as transforms, they don't
play into the log probability (we should make that clearer
in the manual).


Stan also uses different parameterizations than BUGS for
some distros (e.g., deviation rather than precision for
normals).

More below.

> Example CODE
>
> r, M, and u are data, p,N, and superN are parameters to be estimated.
>
> model
> {
> mu~dnorm(0,1.0E-6) #hyper prior on log of single
> escapement estimate
>
> for (i in 1:2)
> {
> p[i]~dbeta(1,1) #uninformative prior n p
> r[i]~dbin(p[i],M[i]) #binomial likelihood to estimate
> capture probability p given recaps (r) and marks release (M)
>
> logN[i]~dnorm(mu,tau) #log N for each population comes
> from a common (super) distribution)
> N[i]<-exp(logN[i]) #transform
>
> #The following two lines work as expected
> expu[i]<-N[i]*p[i] #expected number of unmarked fish
> captured on 2nd pass given N and p
> u[i]~dpois(expu[i]) #This makes N[i] a stochastic node
> which allows HBM to work
>
> #Replacing the previous 2 lines with the following 2 does
> not work due to multiple definitions of node N[]
> #u[i]~dbin(p[i],N[i])
> #N[i]<-u[i]/p[i]

I don't understand what you'd like this to do.

Stan is imperative, so what would happen in Stan is
that you'd define N[i] the second time and never use it.

If you just want to sample u[i] twice because its
probability is a product, that's OK in Stan:

u[i] ~ poisson(exp(logN[i]) * p[i]);

u[i] ~ binomial(floor(exp(logN[i])), p[i]);

That'll result in the following effect on the probability
function being defined:

p(u,...) propto poisson(u[i]|...) * binomial(u[i]|...)

What's the log probability function supposed to be?

Then if you want to look at something, just define a new
generated quantity for it.

- Bob

[Sebastian Weber]

Hmm, maybe something like that is possible. However, for now each group J has Nj subjects allocated and I intended to model the mean number of subjects which have a slope larger than a constant, depending on the treatment of group J. So something like

# get the number of subjects in group j which have a treatment to be counted
Ntreat[j] = sum_i step(slope[j,i] > const)

Ntreat[j] ~ binomial(Nj, inv_logit(alpha + beta * conc[j]))

Then Ntreat needs first to be calculated from thresholded slopes. The problem is the thresholding, but maybe Stan can handle this as a simple transformation step.

Sebastian

[Ben Goodrich]

On Thursday, October 18, 2012 2:59:46 PM UTC-4, Sebastian Weber wrote:

Hmm, maybe something like that is possible. However, for now each group J has Nj subjects allocated and I intended to model the mean number of subjects which have a slope larger than a constant, depending on the treatment of group J. So something like

# get the number of subjects in group j which have a treatment to be counted
Ntreat[j] = sum_i step(slope[j,i] > const)

Ntreat[j] ~ binomial(Nj, inv_logit(alpha + beta * conc[j]))

Then Ntreat needs first to be calculated from thresholded slopes. The problem is the thresholding, but maybe Stan can handle this as a simple transformation step.

I think something like this can be done, but you would have to do it inside the model{} block to work around some integer vs. double issues.

model {

  int Ntreat[J];

  for (j in 1:J) {

    Ntreat[j] = 0;

    for (i in 1:I) {

      Ntreat[j] <- Ntreat[j] + step(slope[i,j] - const);

    }

    Ntreat[j] ~ binomial(N[j], inv_logit(alpha + beta * conc[j]);

  }

}

But, I wonder what that does to the geometry? Usually, we say that if you are going to put a function of parameters on the left-hand side of a sampling statement (with the ~), then you have to make a manual Jacobian adjustment. And that doesn't make sense for a discrete function of continuous parameters.

Ben

[Josh Korman]

Hi,

Thanks for hanging in there for this. The example model I provided estimates
N (abundance) for each population and also SuperN, which is an
information-weighted average abundance across the populations.

So, log(N[i)]~dnorm(mu,sd)
SuperN<-exp(mu)

The problem in my more complex model is that N is actually a function of
other estimated parameters. Say N[i]=a/b. So what I would like Stan to be
able to do is:

1) First calculate N as a function of other model parameters (N[i]=a/b)
2) Then estimate mu based on the values of N via log(N[i])~dnorm(mu,sd)

With this structure the parameter estimates of a and b would adjust so that
the N[i]'s would be consistent, depending on what one assumed for sd. A
small value of sd would result in N[1]~=N[2], and a large value would yield
estimates of N[i] that would be equivalent to values if they were
independently estimated.

If I follow your message below, the first definition of N[i] (~dnorm(mu,sd))
would be ignored and replaced by the expression N[i]=a/b, thus not having
the effect I am after. But perhaps there is a way in STAN to have some sort
of penalty function that increases as N[i]'s diverge, thus forcing N[i]'s to
end up near the same value (depending on the information in that penalty
function).

If this still isn't clear than I'll just code this up in STAN and sort this
out!

Thanks,

Josh

-----Original Message-----
From: stan-...@googlegroups.com [mailto:stan-...@googlegroups.com] On
Behalf Of Bob Carpenter
Sent: October-18-12 11:59 AM
To: stan-...@googlegroups.com
Subject: Re: [stan-users] Will STAN address these modelling and sampling
limitations in WinBugs

--

[Bob Carpenter]

The short answer is that discrete random
quantities don't play nicely with Hamiltonian
Monte Carlo. Geometrically there's a discontinuity due
to the step() function, which is bad news.

In practice, Stan loses the connection between slope[j,i]
and the log probability because the step function flattens
it so the derivative is zero where it's defined. You want
that binomial to have some effect on the slopes, I'd think.

I should talk about the dangers of step() in the manual.
I added yet another to-do item to that effect.

- Bob

> --
>
>

[Bob Carpenter]

On 10/18/12 3:57 PM, Josh Korman wrote:
> Hi,
>
> Thanks for hanging in there for this. The example model I provided estimates
> N (abundance) for each population and also SuperN, which is an
> information-weighted average abundance across the populations.
>
> So, log(N[i)]~dnorm(mu,sd)
> SuperN<-exp(mu)
>
> The problem in my more complex model is that N is actually a function of
> other estimated parameters. Say N[i]=a/b. So what I would like Stan to be
> able to do is:
>
> 1) First calculate N as a function of other model parameters (N[i]=a/b)
> 2) Then estimate mu based on the values of N via log(N[i])~dnorm(mu,sd)

That you can do.

> With this structure the parameter estimates of a and b would adjust so that
> the N[i]'s would be consistent, depending on what one assumed for sd.

That's what should happen with the model coded like that in Stan.

You don't even need to define a variable N[i], you can just
write:

log(a/b) ~ normal(mu,sd);

The only issue is that it's up to you to do the Jacobian adjustment
for the change of variables. I keep meaning to write a longer section
on how to do this in the manual.

> A
> small value of sd would result in N[1]~=N[2], and a large value would yield
> estimates of N[i] that would be equivalent to values if they were
> independently estimated.

Right. Small values of sd will "shrink" log(a/b) toward mu.

> If I follow your message below, the first definition of N[i] (~dnorm(mu,sd))
> would be ignored and replaced by the expression N[i]=a/b, thus not having
> the effect I am after. But perhaps there is a way in STAN to have some sort
> of penalty function that increases as N[i]'s diverge, thus forcing N[i]'s to
> end up near the same value (depending on the information in that penalty
> function).

Hopefully the above should work.

I think the above should work and you don't need to define N[i] even
once, much less twice :-)

> If this still isn't clear than I'll just code this up in STAN and sort this
> out!

Good luck.

- Bob

[Sebastian Weber]

So using a bernoulli distribution on the per subject slopes is numerically better? Is that what you are saying or should I stay with my current approach of basically fitting the posterior?

[Bob Carpenter]

On 10/18/12 4:29 PM, Sebastian Weber wrote:
>
> So using a bernoulli distribution on the per subject slopes is numerically better? Is that what you are saying or should I stay with my current approach of basically fitting the posterior?

I'm saying that it's not going to work to use
a step function to discretize a continuous parameter
that then contributes to the log probability.

The issue is that we follow the gradient of the log
probabilty function and applying the step function
effectively cuts the tie.

- Bob

[Ben Goodrich]

On Thursday, October 18, 2012 3:57:34 PM UTC-4, Josh Korman wrote:

Thanks for hanging in there for this. The example model I provided estimates
N (abundance) for each population and also SuperN, which is an
information-weighted average abundance across the populations.

So, log(N[i)]~dnorm(mu,sd)
SuperN<-exp(mu)

The problem in my more complex model is that N is actually a function of
other estimated parameters. Say N[i]=a/b. So what I would like Stan to be
able to do is:

1) First calculate N as a function of other model parameters (N[i]=a/b)
2) Then estimate mu based on the values of N via log(N[i])~dnorm(mu,sd)

Stan has the lognormal(mu,sd) distribution, which may make this easier. It seems to me like what you are saying would look like this in Stan

data {
  real<lower=0> sd; // fixed by assumption
  // other stuff
}
parameters {
  real<lower=0> a[2];
  real<lower=0> b[2];
  real mu[2];
}
transformed parameters {
  real<lower=0> N[2];
  N[1] <- a[1] / b[1];
  N[2] <- a[2] / b[2];
}
model {
  N ~ lognormal(mu, sd);
  lp__ <- lp__ + JacobianAdjustment;
  // more stuff but not priors on a or b
}

Bob mentioned that this is legal Stan code, but the JacobianAdjustment is one of the easiest things to mess up in Stan. By defining N as a transformed parameter and then putting it on the left-hand side of a sampling statement (with a ~), that implies a some belief about a and b.

f(a,b | mu, sd) \propto abs(det(J{a,b -> N[1], N[2]})) * f(N[1],N[2] | mu, sd)

where J{a,b -> N[1], N[2]} means the Jacobian matrix for the transformation from a,b to N[1],N[2]. That Jacobian matrix would be 2x4 but needs to be square to calculate its determinant. In other words, more than one combination of a[1] and b[1] implies the same N[1] and more than one combination of a[2] and b[2] implies N[2]. So, the posterior is improper I think. If, for example, b were fixed, then the Jacobian matrix would be 2x2.

If this still isn't clear than I'll just code this up in STAN and sort this
out! 

That would be good. Maybe in your actual example the transformation is one-to-one?

Ben