Re: change point models?

[Andrew Gelman]

Hi all. I generally don’t like change point models because it’s my impression that real life does not usually have sudden changes—and when sudden changes exist, we usually know when they happen.

But, setting those qualms aside, the model seems well-enough specified. I see the potential HMC problem, though, as there is a discontinuous derivative. One way to handle this would be to allow the change to be continuous, which in some way is an application of the folk theorem in that I expect actual changes really are continuous.

Beyond this, we (Stan-dev, that is) might want to think more generally about how to patch models that have bits of discontinuity so they work ok in HMC. Maybe we can map to some latent continuous model?

A

> Has anybody fitted a change point model in Stan, please?
>
> We are interested to fit a change point (broken stick) model in RStan to model longitudinal continuous outcomes measured about three to ten times (at baseline and follow-ups) for each of about hundreds of subjects (i.e. patients in a clinical trial); (So this is different from a time series on a single outcome for a single unit of observation, say a NYSE stock.). The timing of the outcome observations is different for each individual (the patients come in for appointments at different, irregularly), as are the slopes before and after the change point, the intercept and the actual change point itself.
>
> The model:
>
> y[i,k] = beta[pre, i] * (t[i, k] – tau[i]) * I(t[i, k] < tau[i]) + beta[intercept, i]+ beta[post, i] * (t[i,k] – tau[i]) * I(t[i,k] > tau[i]) + epsilon[i,j]
>
> where the error term epsilon [i,k] ~MVN (0, Sigma), iid;
>
> y[i,j,t] is the continuous outcome of the ith subject at the kth observation at time t[i,k]; beta[pre,i] is the random slope for each subject before, beta[post,i] is the random slope for each subject after the subject specific change point tau. We think a reasonable prior for beta may be ~MVN, iid, but the slopes and intercept are certainly not iid.
>
> I( t[i,k] > tau[i]) is the indicator function used to define a linear spline: (t[i, k] – tau[i]) * I(t[i, k] < tau[i]) = t[i,k] - tau[k] iff t[i,k] < tau[i] ; 0 otherwise;
>
> tau [i] is the random change point (to be estimated for each subject i); tau [i] is the time of the ‘kink’ in the spline measured before say a diagnosis of dementia.;
>
> See for example attached: Hall CB, Ying J, Kuo L, Sliwinski M, Buschke H, Katz M, Lipton RB. Estimation of bivariate measurements having different change points, with application to cognitive ageing. Stat Med. 2001 Dec 30;20(24):3695-714. PubMed PMID: 11782027.
>
> We are interested to learn from similar models other users may have fitted or know about or to share related functions, now possible in Rstan.
>
> Thank you.
>
> Michael Andreae

[Bob Carpenter]

Is there a discontinuous derivative if we marginalize out
the change point?

There was some discussion of moving to a continuous change
point, but I think that's where you get the discontinuous derivative
right? I don't think that helps --- it would just be one of those
cases where we devolved to a truly random walk (might not be
too bad since there are only two directions in which to walk here).

- Bob

On Jul 2, 2014, at 3:34 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

> Hi all. I generally don't like change point models because it's my impression that real life does not usually have sudden changes--and when sudden changes exist, we usually know when they happen.

>
> But, setting those qualms aside, the model seems well-enough specified. I see the potential HMC problem, though, as there is a discontinuous derivative. One way to handle this would be to allow the change to be continuous, which in some way is an application of the folk theorem in that I expect actual changes really are continuous.
>
> Beyond this, we (Stan-dev, that is) might want to think more generally about how to patch models that have bits of discontinuity so they work ok in HMC. Maybe we can map to some latent continuous model?
>
> A
>
>> Has anybody fitted a change point model in Stan, please?
>>
>> We are interested to fit a change point (broken stick) model in RStan to model longitudinal continuous outcomes measured about three to ten times (at baseline and follow-ups) for each of about hundreds of subjects (i.e. patients in a clinical trial); (So this is different from a time series on a single outcome for a single unit of observation, say a NYSE stock.). The timing of the outcome observations is different for each individual (the patients come in for appointments at different, irregularly), as are the slopes before and after the change point, the intercept and the actual change point itself.
>>
>> The model:
>>

>> y[i,k] = beta[pre, i] * (t[i, k] - tau[i]) * I(t[i, k] < tau[i]) + beta[intercept, i]+ beta[post, i] * (t[i,k] - tau[i]) * I(t[i,k] > tau[i]) + epsilon[i,j]

>>
>> where the error term epsilon [i,k] ~MVN (0, Sigma), iid;
>>
>> y[i,j,t] is the continuous outcome of the ith subject at the kth observation at time t[i,k]; beta[pre,i] is the random slope for each subject before, beta[post,i] is the random slope for each subject after the subject specific change point tau. We think a reasonable prior for beta may be ~MVN, iid, but the slopes and intercept are certainly not iid.
>>

>> I( t[i,k] > tau[i]) is the indicator function used to define a linear spline: (t[i, k] - tau[i]) * I(t[i, k] < tau[i]) = t[i,k] - tau[k] iff t[i,k] < tau[i] ; 0 otherwise;

>>
>> tau [i] is the random change point (to be estimated for each subject i); tau [i] is the time of the 'kink' in the spline measured before say a diagnosis of dementia.;
>>
>> See for example attached: Hall CB, Ying J, Kuo L, Sliwinski M, Buschke H, Katz M, Lipton RB. Estimation of bivariate measurements having different change points, with application to cognitive ageing. Stat Med. 2001 Dec 30;20(24):3695-714. PubMed PMID: 11782027.
>>
>> We are interested to learn from similar models other users may have fitted or know about or to share related functions, now possible in Rstan.
>>
>> Thank you.
>>
>> Michael Andreae
>
>

> --
> 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.
>
>

[Andrew Gelman]

The change point, if unknown, should be continuous—but, even so, I think discontinuity arises if that change point is a random variable. It makes sense that a mixture formulation would resolve the problem.

[Michael Betancourt]

The closest continuous model is given by integrating the change points out, yielding

a highly multimodal model that will be nearly impossible to fit in its own right.  Unless

someone can figure out a way to isolate a single mode?

[Andrew Gelman]

When the model is highly multimodal, in practice one can deal with this crudely via starting points or more systematically via informative priors.  (or ultimately via thermodynamic sampling!)

[Michael Andreae]

Update: (thanks to a hint by Dr. Goodrich, with most credit to Jonah Sol Gabry):

Is this a way around the discontinuous derivative issue?

Problem: Observational cohort shows an reduction in a disease outcome y after initiation of treatment:

tau, the single change point where this change occurs, is CONTINUOUS

y, the outcomes is measured k =5 times at irregular intervals x, for each subject n=100

data are in long format, we use a inv-wishart as prior for the beta intercept, beta pre slope and beta post change point slope.

Preliminary model below. 

model{

  for (i in 1:n) {

          y [i] ~ normal(beta[1]*100,sigma);                               //intercept

     for (j in 2:k) {

       if (x [ i + (n* (j-1)) ] <= tau )                            //pre tau point

          y [i + (n* (j-1))] ~ normal (y[i] 

                                      +beta [2] * ( x [i + (n* (j-1))] - x[i] ), 

                                      sigma);

       

       else                                                         //post tau point

          y [i + (n* (j-1))] ~ normal( y[i] 

                                       +( beta [2] * (tau - x[i]) )  

                                       +( beta [3] * (x [i + (n* (j-1))] - tau)), 

                                       sigma);

     }

  }

beta ~ multi_normal(means, SIGMA); 

means ~ normal(0, 1000);

SIGMA ~ inv_wishart(3.0,W);

[Bob Carpenter]

You only provided the model part of the Stan program. It
really helps to understand the parameter and data definitions,
too.

We have alternative recommendations for priors on covariance
matrices, but that's a second-order detail.

I'm also not clear on why beta[1] is multiplied by 100 on the
first line.

And wouldn't it be easier to use 2D arrays than fiddling indexes?
The latter's very error prone, at least when I'm coding.

- Bob

[Sebastian Weber]

Hi!

Well, not really a change-point model, but something very similar are so-called Emax models in PK/PD. The hill function can be written in terms of the log-concentration like

hill(conc) = inv_logit(h * (lconc - lEC50) )

here the "l" denotes log. I experienced some issues in Stan when fitting this type of model with h larger than 1. I suppose the derivatives can become quite steep in that case. So if such a function could be stabilized by providing a specialized function, that would be great. The behavior would be such that with h >> 1 the function essentially becomes the step function.

Best,

Sebastian

[Luc Coffeng]

As discussed elsewhere on the Stan user list, functions with step-like behaviour cause ridges in the posterior (like you say, areas where the gradient is extreme), which are problematic as they may confine the Hamiltonian particle to an area around the posterior mode (or one of the modes), unless the random momentum is just exactly right (not too big to overshoot and prevent efficient sampling, not too small to get stuck around the (or a particular) mode). Correct me if I'm wrong, but the random momentum is calibrated to an optimal values during adaptation, but this is probably optimized such that the sampler explores the area around the mode efficiently, and not the whole posterior. As a result, the sampler will poorly explore tails of the posterior. I remember the Stan dev team pointing out that Stan's samplers work best with posteriors with a curvature close to a normal distribution. Maybe therefore, Stan is currently not very good at fitting non-linear models. Hopefully, Riehmannian HMC will make it possible to better explore extremely curved posteriors (what's the planning for that?). Until that's available, I'm considering using some brute-force methods such as Bayesian melding and incremental mixed importance sampling (IMIS; there's an R package) to fit non-linear models as an alternative to Stan.

Best,
Luc

[Bob Carpenter]

Adapatation does a global estimation of the mass matrix (diagonal
by default). It only uses the second half of the adaptation iterations
to do this, so if the sampler has converged, it will do this in the
high-mass volume of the posterior, aka the typical set, which is
not necessarily very close (especially in high dimensions) to the mode.

For the "step-like" function Sebastian's talking about, if the
multiplier (like an IRT discrimation value) isn't too high, then
the derivatives will still be stable. You could put a prior on
that multiplier to give it effectively step-like behavior, but
still be stably differentiable, I would think. Just be sure to
use the log_inv_logit variant for arithmetic stability on the log scale.

I'm not sure what you mean by fitting non-linear models.
Often the posteriors, which is the curvature that matters, can
be normal even when the models don't look particularly normal.
But as Luc says, if there's very high curvature, the Hamiltonian
dynamics simulator (aka the leapfrog integrator) will struggle
and need small step sizes (and hence more time).

- Bob

[Luc Coffeng]

Is there any way to have Stan perform the estimate the optimal number of leap-frog steps and the mass matrix, while telling it to take a certain step-size (presumably a very small step-size to deal with extreme curvature)? Or would you have to perform a test run to estimate tuning parameters, and then manually re-enter an adapted version on the next run?

[Bob Carpenter]

You probably always want to use NUTS rather than estimating
an optimal number of leapfrog steps, because it's not necessarily
the same for every position and random momentum.

You can fix the step size.

But the easiest thing to do is just set a higher target acceptance
rate for adaptation (the delta parameter), which forces a smaller
step size to ensure a higher acceptance rate.

Michael's been recommending validation runs that crank down the
acceptance rate target substantially to make sure runs with higher
acceptance rates are doing the right things in the tails or other hard
to read parts of the posterior.

- Bob

On Jul 18, 2014, at 12:28 AM, Luc Coffeng <lucco...@gmail.com> wrote:

> Is there any way to have Stan perform the estimate the optimal number of leap-frog steps and the mass matrix, while telling it to take a certain step-size (presumably a very small step-size to deal with extreme curvature)? Or would you have to perform a test run to estimate tuning parameters, and then manually re-enter an adapted version on the next run?
>

[Sebastian Weber]

Hi Bob!

It certainly makes sense to calculate the log_inv_logit with respect to smooth derivatives - but at some point in time I simply have to back transform it with exp and here I am a bit confused, i.e. at what point should I do that?

I mean the hill function does not live in the log-space, but I need it on the natural scale at some point in time in order to formulate the likelihood.

Best,

Sebastian

[Michael Betancourt]

Right — if you’re worried that the adaptation is not finding a small enough step size
then you can rerun with a larger target acceptance probability, or turn off adaptation
altogether and run with a constant step size. The ultimate issue here is one of
spatially-varying curvature — the integrator needs to have a step size adapted to the
largest curvatures and when there are huge variations this inevitably leads to
inefficient samples in regions of low curvature.

As Bob noted, _always_ run a second chain with a larger target acceptance
probability to be sure that you didn’t miss anything. And make sure you look at
the n_divergent statistic! If that’s non-zero you’re already having problems.

If you’re already seeing problems there’s not much you can do. In fact, such a
pathological model will be hard to fit using almost _any_ algorithm. Importance
sampling is going to fail just as hard because the hard edges mean that the
importance weights are concentrated on a very small filament that will be hard to
hit with unguided samples.

RHMC will definitely help here, but it won’t be a complete panacea. I’ve really
bought into the folk theorem in the last year — these kinds of pathological models
typically indicate that the model isn’t quite what you want it to be, either in
terms of the actual specification or it’s explicit parameterization.

[Sebastian Weber]

Hmm... RHMC sounds like it is going to solve many problems, but this beast is also very costly computational wise. Is it possible to run RHMC on a sub-space of the parameters space and use HMC for the rest?

Sorry for already asking for features for RHMC which isn't even in Stan yet... just a thought...

[Michael Betancourt]

Hmm... RHMC sounds like it is going to solve many problems, but this beast is also very costly computational wise. Is it possible to run RHMC on a sub-space of the parameters space and use HMC for the rest?

Anytime you try to partition parameters you’re going to face a huge drop is performance and

this would be no different.  There are various approximate RHMC implementations that can

be much faster if the pathologies are limited (see the SoftAbs paper), but in the end the kinds 

of pathologies we’re talking about increase exponentially with dimension, N, and the N^{3} 

cost of RHMC isn’t that bad in comparison.

[Luc Coffeng]

Thanks, I'll keep that in mind! Not sure what the n_divergent statistic is though. And yes, importance sampling fails hard in case of highly curved posteriors and unguided samples, but that's actually what incremental mixed importance sampling tries to address. The algorithm iteratively resamples highly curved areas of the posterior at a higher resolution (and reweights those sample appropriately) until some satisfactory number of different samples has been acquired. I haven't tried it yet, but it sounded promising.

[Michael Betancourt]

Thanks, I'll keep that in mind! Not sure what the n_divergent statistic is though.

It's discussed in the manual and http://arxiv.org/abs/1312.0906, and is displayed in the bin/print output.

 

And yes, importance sampling fails hard in case of highly curved posteriors and unguided samples, but that's actually what incremental mixed importance sampling tries to address. The algorithm iteratively resamples highly curved areas of the posterior at a higher resolution (and reweights those sample appropriately) until some satisfactory number of different samples has been acquired. I haven't tried it yet, but it sounded promising.

None of these adaptive algorithms scale well to high dimensions.  If you're lucky enough to find a point with high curvature, sampling around it in any naive way is even more unlikely to find neighboring points

with high curvature.  High dimensions are hard and importance sampling is never the answer, but I'll leave that for you to find empirically.

[Ben Goodrich]

On Friday, July 18, 2014 5:59:12 AM UTC-4, Michael Betancourt wrote:

Thanks, I'll keep that in mind! Not sure what the n_divergent statistic is though.

It's discussed in the manual and http://arxiv.org/abs/1312.0906, and is displayed in the bin/print output.

For RStan users, it is a bit harder but can be displayed with get_sampler_params(), as in:

dogs <- stan_demo("dogs")
sampler_params <- get_sampler_params(dogs)
names(sampler_params) <- colnames(dogs)
lapply(sampler_params, FUN = function(x) table(x[,"n_divergent__"]))
$`chain:1`

   0    1
1991    9

$`chain:2`

   0    1
1991    9

$`chain:3`

   0    1
1991    9

$`chain:4`

   0    1
1993    7

> lapply(sampler_params, FUN = function(x) which(x[,"n_divergent__"] > 0))
$`chain:1`
[1]   2   3   4  10  48 251 252 452 952

$`chain:2`
[1]   2   3   4 102 252 452 518 952 987

$`chain:3`
[1]   2   3   4 102 152 252 452 532 952

$`chain:4`
[1]   2   3 102 152 251 252 952

Not too bad. Appears to be more likely to diverge after a window switch.

Ben

 

[Avraham Adler]

Could that be built as a test/function in rstan itself for ease of use?

Avi

[Luc Coffeng]

Thanks Michael, I'll take your word for it! Regarding the n_divergent statistics, what manual do you mean? I can't seem to find it in the reference manual.

Thanks Ben, as an rstan used, I had indeed never come across this statistic.

-Luc

[Michael Betancourt]

Looks like it got moved to the cmdstan manual, but then it's only a quick definition and not much exposition.  I'd go through the paper I referenced for an attempt at explanation.

This is all related to theoretical guarantees on the performance of HMC there we're actively working on right now.  Once we (myself and my colleagues at Warkwick/UCL) have 

started drafting papers we'll start incorporating the material into the manual.

--

[Bob Carpenter]

But Stan only needs log probability, never linear probabilities.
The more you can do things on the log scale, the more stable it
tends to be.

- Bob