Mixture models with multiple assignments per mixture

[David Westergaard]

Hi,

I would like to ask for some model guidance. I apologise in advance if my problem seems trivial, but my experience with Bayesian models is limited.

I am trying to identify mixtures in a population, in which each data point has its own profile. The tricky part is that each data point may be a member of multiple mixtures. From what I have found so far, this can be modelled using the so-called "Indian Buffet Process" (IBP) (Griffiths et al., 2011). However, as I have been able to read in the reference and on the board, Stan does not (yet) handle these kind of infinite models. From other threads on this forum, I have seen people have been mimicking the Chinese Restaurant Process using a very large number of mixtures. Can, or have, something similar been done for the IBP?

I am not even sure this is the best model to use, but I do not know of any better. To go a bit more into detail, the data is a NxM matrix, in which there are N data points and M measurements. The measurements are continous and (approximately) normally distributed (on log scale). Each data point has a lot of measurements with a zero value, simply because nothing was detected.

I have tried to come up with a complete description, but if anything is missing, please do ask!

Thanks in advance.

Cheers,

David

[Bob Carpenter]

Trivial's a matter of perspective --- almost everything's
hard the first time and trivial after the tenth time.

It doesn't sound like you have an infinite mixture if you have
a fixed number of columns. From Griffith's and Gharamani's
abstract ( https://cocosci.berkeley.edu/tom/papers/indianbuffet.pdf ):

The Indian buffet process is a stochastic process defining
a probability distribution over equiva- lence classes of sparse
binary matrices with a finite number of rows and an unbounded number
of columns. This distribution is suitable for use as a prior
in probabilistic models that represent objects using a potentially
infinite array of features, or that involve bipartite graphs in
which the size of at least one class of nodes is unknown.

The thing to do when building a Bayesian model is to consider how
the data is generated (here it sounds like some kind of measurement
process). Could you say a bit more about what the
colummns are and how the zero entries come about if the measurements
are approximately lognormal (0 isn't going to come out of a lognormal
observation)? For example, is it a censoring process that truncates
low measurements to zero? If so, you probably just want a model to
account for the censoring. Is there any relation among the columns
or are they all different kinds of measurements?

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

[David Westergaard]

It could be me that read it wrong. I was thinking of it as mixture model with an infinite number of mixtures, in which each data point could be member of multiple mixtures. 

Yes, it is indeed a measurement process. Each "feature" is the same measurement, measured under slightly different conditions. What we expect is that feature 1 is correlated to feature 2, and feature 2 to 3, ..., approximately until we measure 0 (absence). The standard in my field is to work with log-values, so I use the log(value+1) transformation to avoid the log(0) errors, and this is (approximately) normally distributed, excluding the 0 measurements. 

To give an example:

A = [4.5, 5, 6, 5.5, 8, 0, 0, 0, 0, 3, 2.3, 4.5]

B = [4.5, 5, 6, 5.5, 8, 0, 0, 0, 0, 0, 0, 0]

C = [0, 0, , 0, 0, 0, 0, 0, 0, 2.9, 2.2, 4.4]

So A would have the profile [4.5, 5, 6] and [3, 2.3, 4.5]. What I would like to identify is that point A has two mixtures (1 and 2), B has one mixture (1) and C has one mixture (2). Thus mixture 1 consists of A+B, and mixture 2 consists of A+C.

I hope the explanation was clear. 

Cheers,

David

[Bob Carpenter]

> On Mar 4, 2015, at 12:01 AM, David Westergaard <da...@harsk.dk> wrote:
>
> It could be me that read it wrong. I was thinking of it as mixture model with an infinite number of mixtures, in which each data point could be member of multiple mixtures.

It's for infinitely fat mixtures. If you have finitely
many columns and model each one, then there's only finitely
many. It's the "unbounded number of columns" part that
your data's missing.

> Yes, it is indeed a measurement process. Each "feature" is the same measurement, measured under slightly different conditions. What we expect is that feature 1 is correlated to feature 2, and feature 2 to 3, ..., approximately until we measure 0 (absence).

It would probably help if you could be more specific (less abstract)
about what these measurements are.

A mixture would indeed give you that correlation.

> The standard in my field is to work with log-values, so I use the log(value+1) transformation to avoid the log(0) errors, and this is (approximately) normally distributed, excluding the 0 measurements.
>
> To give an example:
>
> A = [4.5, 5, 6, 5.5, 8, 0, 0, 0, 0, 3, 2.3, 4.5]
> B = [4.5, 5, 6, 5.5, 8, 0, 0, 0, 0, 0, 0, 0]
> C = [0, 0, , 0, 0, 0, 0, 0, 0, 2.9, 2.2, 4.4]
>
> So A would have the profile [4.5, 5, 6] and [3, 2.3, 4.5].

Is that the result of log(measurement + 1) or the raw measurements
that you'll add 1 to and take the log?

Is there some physical correspondence to a profile here?

> What I would like to identify is that point A has two mixtures (1 and 2), B has one mixture (1) and C has one mixture (2). Thus mixture 1 consists of A+B, and mixture 2 consists of A+C.

Does the order of the columns matter? For instance, are
they measurements over time? Clearly the behavior isn't
continuous, because you go from 8 down to 0 in one step. But
the values close to each other look related in A. That would
matter in terms of how you'd want to model it.

Then the question is what do you expect the relationship to
be between A's and B's measurement. It looks like B has
one of the components of A and C has one of the components
of A.

If so, you might want something like a change point
time-series model. The data for A is probably three mixtures,
if you include the zeroes. You want to model the entire
data generating process, including generating the zeroes.

So a simple model might assume some number of mixture components K,
each with means mu[k] and std devs sigma[k], and then each
row would start in one of the comonents, continue to draw from
that mixture for some number of observations, then drop to 0 or
another mixture.

- Bob

[David Westergaard]

> On 03 Mar 2015, at 14:28, Bob Carpenter <ca...@alias-i.com> wrote:
>
>
>> On Mar 4, 2015, at 12:01 AM, David Westergaard <da...@harsk.dk> wrote:
>>
>> It could be me that read it wrong. I was thinking of it as mixture model with an infinite number of mixtures, in which each data point could be member of multiple mixtures.
>
> It's for infinitely fat mixtures. If you have finitely
> many columns and model each one, then there's only finitely
> many. It's the "unbounded number of columns" part that
> your data's missing.
>

I see. I might have to re-read the paper. It seems I misunderstood the concept.

>
>> Yes, it is indeed a measurement process. Each "feature" is the same measurement, measured under slightly different conditions. What we expect is that feature 1 is correlated to feature 2, and feature 2 to 3, ..., approximately until we measure 0 (absence).
>
> It would probably help if you could be more specific (less abstract)
> about what these measurements are.
>

Ok, I will try to exemplify the process in more detail, still keeping it a bit abstract.

Imagine a lane. On the lane, there are different obstacles, varying in size (getting larger as one moves down the lane). Now along comes two components A and B, which are physically connected. They move down the lane until they cannot proceed any further due to obstacles. Next comes A and C, which are also physically connected. C is, however, slightly smaller than B, and the A-C component can move further down the lane. The number A-B and A-C are typically in the millions, and what we count is not the actual number of A-B components, but A and B individually (typically these are in millions++, so we use the log(count+1) transformation). Not all of the components stops the exact same place (hence they have a profile).

I hope that made it a bit more clear!

> A mixture would indeed give you that correlation.
>
>> The standard in my field is to work with log-values, so I use the log(value+1) transformation to avoid the log(0) errors, and this is (approximately) normally distributed, excluding the 0 measurements.
>>
>> To give an example:
>>
>> A = [4.5, 5, 6, 5.5, 8, 0, 0, 0, 0, 3, 2.3, 4.5]
>> B = [4.5, 5, 6, 5.5, 8, 0, 0, 0, 0, 0, 0, 0]
>> C = [0, 0, , 0, 0, 0, 0, 0, 0, 2.9, 2.2, 4.4]
>>
>> So A would have the profile [4.5, 5, 6] and [3, 2.3, 4.5].
>
> Is that the result of log(measurement + 1) or the raw measurements
> that you'll add 1 to and take the log?
> Is there some physical correspondence to a profile here?
>

Depends what you mean by physical. Each entry is an actual count of the component. I.e. A[0], B[0], and C[0] (using array index notation) are from the exam same place on the lane.

>> What I would like to identify is that point A has two mixtures (1 and 2), B has one mixture (1) and C has one mixture (2). Thus mixture 1 consists of A+B, and mixture 2 consists of A+C.
>
> Does the order of the columns matter? For instance, are
> they measurements over time? Clearly the behavior isn't
> continuous, because you go from 8 down to 0 in one step. But
> the values close to each other look related in A. That would
> matter in terms of how you'd want to model it.
>

Yes, these values should indeed be related.

> Then the question is what do you expect the relationship to
> be between A's and B's measurement. It looks like B has
> one of the components of A and C has one of the components
> of A.

Yes, that was exactly the point I was trying to make with the example. And this is why I was looking into infinite mixture models. Essentially we don’t know how many mixtures there are, or how many mixtures each component can be part of. We want to be informed by the data about this (with some prior assumptions).

>
> If so, you might want something like a change point
> time-series model. The data for A is probably three mixtures,
> if you include the zeroes. You want to model the entire
> data generating process, including generating the zeroes.
>
> So a simple model might assume some number of mixture components K,
> each with means mu[k] and std devs sigma[k], and then each
> row would start in one of the comonents, continue to draw from
> that mixture for some number of observations, then drop to 0 or
> another mixture.

I am not quite sure I understand the last part. Wouldn’t I then run into the problem that each row would only be associated with one mixture?

Cheers,
David


> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/dNyt1uXvaYI/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

OK, I think I see this a bit better. What you probably want
to be looking at in the way of non-parametrics is a Dirichlet
process prior on the number of mixture components. That would
let the K in my earlier comment be given a prior. If you
have a reasonable upper bound, this can be well approximated with
a finite upper bound on K.

And no, you wouldn't associate each row with a single mixture
component. You'd have something like this for row n:

x[n,1], x[n,2], ..., x[n,M]

and maybe some kind of process for generating a sequence of
components where lambda[n,m] in 1:K is the mixture from which
x[n,m] is drawn,

lambda[n,1], ..., lambda[n,M]

Then you'd want to put a model on how the lambda values transition.
An HMM would be a general way to do it, but they're not very good
at modeling length --- you only get exponential decay beyond a certain
finite memory limit. That's what I was calling the change-point structure.

Thanks for the extra detail on the process --- I'll have to
think a bit more about that. The point's just that you want
to model that process as closely as you can in the statistical
model.

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

[Michael Betancourt]

The IBP is just a certain implementation of a Dirichlet process
prior. Because we don’t handle discrete parameters you’d have
to think about the collapsed IBP, but then you still have to worry
about the unbounded mixture components. As Bob alluded
to in his first email, the only way this will work in Stan is if
you threshold the number of mixture components a priori
and hope that your posterior mass concentrates away from
that threshold.

[David Westergaard]

On 03 Mar 2015, at 16:55, Michael Betancourt <betan...@gmail.com> wrote:

The IBP is just a certain implementation of a Dirichlet process
prior.  Because we don’t handle discrete parameters you’d have
to think about the collapsed IBP, but then you still have to worry
about the unbounded mixture components.  As Bob alluded
to in his first email, the only way this will work in Stan is if
you threshold the number of mixture components a priori
and hope that your posterior mass concentrates away from
that threshold.

Yes, that is also what I deduced from previous threads on the board. Is there a process that would shrink the number of mixtures the model use, so to speak? Or how would I go about checking the posterior mass concentrates away?

On Mar 3, 2015, at 3:35 PM, Bob Carpenter <ca...@alias-i.com> wrote:

OK, I think I see this a bit better.  What you probably want
to be looking at in the way of non-parametrics is a Dirichlet
process prior on the number of mixture components.  That would
let the K in my earlier comment be given a prior.  If you
have a reasonable upper bound, this can be well approximated with
a finite upper bound on K.

And no, you wouldn't associate each row with a single mixture
component.  You'd have something like this for row n:

x[n,1], x[n,2], ..., x[n,M]

and maybe some kind of process for generating a sequence of
components where lambda[n,m] in 1:K is the mixture from which
x[n,m] is drawn,

lambda[n,1], ..., lambda[n,M]

Then you'd want to put a model on how the lambda values transition.

That is actually a very good and interesting idea. Can you recommend any literature on this (apart from the Stan reference)?

An HMM would be a general way to do it, but they're not very good

at modeling length --- you only get exponential decay beyond a certain
finite memory limit.  That's what I was calling the change-point structure.

Not sure I follow here. My knowledge of HMMs are limited to the very basics. 

Cheers,

David

[Michael Betancourt]

Yes, that is also what I deduced from previous threads on the board. Is there a process that would shrink the number of mixtures the model use, so to speak?

You would tune alpha to be small.

Or how would I go about checking the posterior mass concentrates away?

Just histogram the number of active components.

[David Westergaard]

This might be a very basic question, but how do I "know" if a component is active? Simply if the probability > 0, or is there a lower-bound that is considered to be ~0 for practical purposes?

I tried to make a simple mixture model in stan (guided mostly by the Stan reference and wikipedia), but right now it is running incredibly slow. After a day it had barely done 25 iterations, so I am guessing it is mixing very poorly for a reason I do not understand. I have put a copy of the model at https://gist.github.com/demodw/f5e39f08ebb5f49e2e85#file-gistfile1-txt. There is a lot of (potentially) unnecessary comments, but these are mostly for myself. 

Cheers,

David

 

[Michael Betancourt]

> This might be a very basic question, but how do I "know" if a component is active? Simply if the probability > 0, or is there a lower-bound that is considered to be ~0 for practical purposes?

You don’t.

When you marginalize out over the mixture components you’re assigning at least some prior
probability to each component which can never be forced exactly to zero. A proper Bayesian
analysis never makes a hard decision to exclude components so there is no canonical way
of trying to impose a “significance” threshold. You have to define a criteria for defining one
yourself.

For example, one criteria might be that for a mixture component to be excluded the posterior
must shrink towards zero relative to the prior. Of course how you define “shrink” and “shrink
enough to be excluded” are open questions.

> I tried to make a simple mixture model in stan (guided mostly by the Stan reference and wikipedia), but right now it is running incredibly slow. After a day it had barely done 25 iterations, so I am guessing it is mixing very poorly for a reason I do not understand. I have put a copy of the model at https://gist.github.com/demodw/f5e39f08ebb5f49e2e85#file-gistfile1-txt. There is a lot of (potentially) unnecessary comments, but these are mostly for myself.

Your model is not implemented correctly — you need to compute the density of each
component and then sum them up for each datum,

real ps[K];

for (n in 1:N) {
for (k in 1:K) {
ps[k] <- log(theta[k]) + normal_log(y[n], mu[k], sigma[k]);
}
increment_log_prob(log_sum_exp(ps));
}

[David Westergaard]

Den mandag den 9. marts 2015 kl. 16.40.45 UTC+1 skrev Michael Betancourt:

> This might be a very basic question, but how do I "know" if a component is active? Simply if the probability > 0, or is there a lower-bound that is considered to be ~0 for practical purposes?

You don’t.

When you marginalize out over the mixture components you’re assigning at least some prior
probability to each component which can never be forced exactly to zero.  A proper Bayesian
analysis never makes a hard decision to exclude components so there is no canonical way
of trying to impose a “significance” threshold.  You have to define a criteria for defining one
yourself.

For example, one criteria might be that for a mixture component to be excluded the posterior
must shrink towards zero relative to the prior.  Of course how you define “shrink” and “shrink
enough to be excluded” are open questions.

I see. On the other hand, I am not particularly interested in which of the components are active or not, but modelling which of the samples come from the same mixture(s) (I think, but I guess it is all correlated to some degree). Bayesian inference is still a new thing to me, so I am still struggling a bit with defining the generating process for the problem. :-)

> I tried to make a simple mixture model in stan (guided mostly by the Stan reference and wikipedia), but right now it is running incredibly slow. After a day it had barely done 25 iterations, so I am guessing it is mixing very poorly for a reason I do not understand. I have put a copy of the model at https://gist.github.com/demodw/f5e39f08ebb5f49e2e85#file-gistfile1-txt. There is a lot of (potentially) unnecessary comments, but these are mostly for myself.

Your model is not implemented correctly — you need to compute the density of each
component and then sum them up for each datum,

real ps[K];

for (n in 1:N) {
  for (k in 1:K) {
    ps[k] <- log(theta[k]) + normal_log(y[n], mu[k], sigma[k]);
  }
  increment_log_prob(log_sum_exp(ps));
}

I had my suspicion about this line. I thought it was a "smart" idea to vectorize it, instead of looping over the components, but I guess I have misunderstood what increment_log_prob does.

Cheers,

David

 

[Bob Carpenter]

increment_log_prob(x) just adds x to the result of the log
probability function.

It's probably log_sum_exp that is the problem. Just expand
the definitions out and look at your density --- that's not
a Bayesian issue --- it works the same way if you're doing
optimization for penalized maximum likelihood.

- Bob

[David Westergaard]

Bob,

Following on your previous suggestion a few mails back, how would I get a value for each sample, feature (as you write: x[n,1], x[n,2], ..., x[n,M]). Are you referring to simply moving the x variable in my model from the data to the parameter block? Or using the estimated parameters in the generated quantities block?

Sorry to bug you, but have you had any chance to consider if there is a better way/more appropriate way to model this process?

This is the model I have so far: https://gist.github.com/demodw/ed5a6159a183679a1d52 (it’s pretty much a replica of the one from the reference.)

Cheers,

David

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/dNyt1uXvaYI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

I'm afraid I'm not sure what exactly what process you're
talking about. I answer a whole lot of these questions
and need a bit more reminder of context.

When I tried to look up your gist,

https://gist.github.com/demodw/ed5a6159a183679a1d52

I get a redirected 404 message:

Whoops. We seem to have missed the gist of that gist you were looking for.

- Bob

[David Westergaard]

Bob,

Sorry - it is quite far back in the thread. Just to copy my (pseudo) explanation of the generative process from earlier: 

"Imagine a lane. On the lane, there are different obstacles, varying in size (getting larger as one moves down the lane). Now along comes two components A and B, which are physically connected. They move down the lane until they cannot proceed any further due to obstacles. Next comes A and C, which are also physically connected. C is, however, slightly smaller than B, and the A-C component can move further down the lane. The number A-B and A-C are typically in the millions, and what we count is not the actual number of A-B components, but A and B individually (typically these are in millions++, so we use the log(count+1) transformation). Not all of the components stops the exact same place (hence they have a profile).”

My initial idea was to use a mixture model (with a dirichlet process prior for the mixture - in this case just a reasonable high number of mixtures) to estimate which components came from the same mixture. That is, identify A-B and A-C as physically connected. My question was, what is the most convenient way in stan to get the probability of each sample, feature (x[n,m) being in mixture 1 .. K. Normally I would just compute the density under the estimated gaussian distribution, but I am not quite sure what is the most optimal way to do this in Stan. 

Using those probabilities, my (very naive) idea to estimate which components come from the same mixture would simply be to apply some lower cutoff, and if the probability exceeds that they are from the same mixture. Following this, then simply just identify which components have the same active mixture(s). But again, this is a very naive idea. You suggested a HMM instead. I am not quite sure how I would apply a HMM in this case. I guess another approach could also be to construct a binary matrix. I.e. b[n,m] == 1 if the mixture is active, otherwise zero, and then apply some latent allocation procedure. 

Sorry for the broken link! https://gist.github.com/demodw/dec9724b96fbf8f8e428 is the correct one. 

Cheers,

David

[Bob Carpenter]

Thanks --- I remember now. Though I probably don't understand
the problem well enough to be of much help.

What are the numbers in y[n]? I didn't understand what the A, B,
and C were or what was observed.

You didn't even declare y in the
gist, so it can't possibly compile.

The model you have (modulo not defining y) is just a standard
mixture of K normals.

If you want to approximate a Dirichlet process, you have to
let the Dirichlet basis alpha shrink to alpha/K as you increase
K. At least I'm pretty sure that's how it'd work, though I have
no idea if it'll work in Stan.

Do you really think there will be 1000 mixture components? It's
likely to lead to very slow going and potentially arithmetic instability.

Yes, it's a good idea to put a prior on everything. Choose something
consistent with what you expect the posterior to be given your
knowledge of the problem and the scale of the variables involved.

If your data's all positive, normal can still be a good approximation.
For example, it's a good approximation of a binomial as the number
of observations grows (just the central limit theorem). And in that
case, you wouldn't want to log transform the counts.

- Bob

[David Westergaard]

On 12 Mar 2015, at 15:46, Bob Carpenter <ca...@alias-i.com> wrote:

Thanks --- I remember now.  Though I probably don't understand
the problem well enough to be of much help.  

You have been of much help already. This is one of my first times actually working with Bayesian models, instead of just reading about them, so I appreciate all the advice. 

The problem is explained a bit more in scientific jargon in the introduction of http://bioinformatics.oxfordjournals.org/content/early/2014/10/08/bioinformatics.btu623.full.pdf+html (paywall’ed). Here they use a simple hierarchical clustering scheme to solve the problem. Others have attempted to solve it using deconvolution of Gaussian distributions (i.e. non-linear regression). The reason I thought of a mixture model was figure 1b. It looks very much, to me, like a mixture of distributions. 

What are the numbers in y[n]?  I didn't understand what the A, B,
and C were or what was observed.

You didn't even declare y in the
gist, so it can't possibly compile.

Sorry, I usually don’t make this many mistakes when copy/pasting. y is a vector[M] y[N] (https://gist.github.com/demodw/8a3f0963b7b2cad36508). A, B and C are physical entities (proteins), that are observed to co-migrate.

Technically the numbers in y are counts. But the counts are in the billions (I actually suspect Stan is overflowing when I feed the numbers — it gives mysterious errors that does not occur when I log2 transform the data). The log2 values typically range between 20-50. 

The model you have (modulo not defining y) is just a standard
mixture of K normals.  

If you want to approximate a Dirichlet process, you have to
let the Dirichlet basis alpha shrink to alpha/K as you increase
K.  At least I'm pretty sure that's how it'd work, though I have
no idea if it'll work in Stan.

Do you really think there will be 1000 mixture components?  It's
likely to lead to very slow going and potentially arithmetic instability.

1,000 is a somewhat conservative estimate. It is indeed very slow when working with the full data, and identifiability across chains is a big problem. I’m wondering if using a EM optimisation is better in this situation, given the large number of mixtures. 

Yes, it's a good idea to put a prior on everything.  Choose something
consistent with what you expect the posterior to be given your
knowledge of the problem and the scale of the variables involved.

If your data's all positive, normal can still be a good approximation.
For example, it's a good approximation of a binomial as the number
of observations grows (just the central limit theorem).  And in that
case, you wouldn't want to log transform the counts.

That is one of my problems, also related to identifiability. I know that most likely each of the 1,000 centers will be in the range of ~20-50ish, but I suspect that won’t really make the prior strong enough to solve the problem of identifiability. 

Cheers,

David

[Bob Carpenter]

Bayesian models aren't that different form classical ones, it's
the inference where things tend to be different. You can use
Stan's optimization, which is an alternative to the EM algorithm.

I still don't understand the problem and I'm afraid I don't have time to
read a bio journal. It takes a long time for me to understand a problem if
it's not broken down into standard statistical terminology.

So my question is: what are the entries in the y[n] vectors for n in 1:N?
How big is N that you want 1000 mixtures?

If they're counts, what are they counts of? I can't connect the dots between
counts and migrating proteins. I'm guessing that you have N different
proteins, M places they can get stuck in an experiment, and y[n] represents
the counts of the places they get stuck over billions of replications for each
n in 1:N. So sum(y[n]) is in the billions, but how big is N and how big is M?

And the goal's then to cluster the proteins by the distributions of where
they get stuck? To what end? Are you sure you want to run a Bayesian
model here?

I don't see how regression is related to clustering.

- Bob

> On Mar 13, 2015, at 2:28 AM, David Westergaard <da...@harsk.dk> wrote:
>
>
>
>> On 12 Mar 2015, at 15:46, Bob Carpenter <ca...@alias-i.com> wrote:
>>
>> Thanks --- I remember now. Though I probably don't understand
>> the problem well enough to be of much help.
>
> You have been of much help already. This is one of my first times actually working with Bayesian models,

> instead of just reading about them, so I appreciate all the advice.
>
> The problem is explained a bit more in scientific jargon in the introduction of http://bioinformatics.oxfordjournals.org/content/early/2014/10/08/bioinformatics.btu623.full.pdf+html (paywall’ed). Here they use a simple hierarchical clustering scheme to solve the problem.

> Others have attempted to solve it using deconvolution of Gaussian distributions (i.e. non-linear regression). The reason I thought of a mixture model was figure 1b. It looks very much, to me, like a mixture of distributions.
>
>>
>> What are the numbers in y[n]? I didn't understand what the A, B,
>> and C were or what was observed.
>>
>> You didn't even declare y in the
>> gist, so it can't possibly compile.
>
> Sorry, I usually don’t make this many mistakes when copy/pasting. y is a vector[M] y[N] (https://gist.github.com/demodw/8a3f0963b7b2cad36508). A, B and C are physical entities (proteins), that are observed to co-migrate.
>
> Technically the numbers in y are counts. But the counts are in the billions (I actually suspect Stan is overflowing when I feed the numbers — it gives mysterious errors that does not occur when I log2 transform the data). The log2 values typically range between 20-50

[David Westergaard]

On 12 Mar 2015, at 16:55, Bob Carpenter <ca...@alias-i.com> wrote:

Bayesian models aren't that different form classical ones, it's
the inference where things tend to be different.  You can use
Stan's optimization, which is an alternative to the EM algorithm.

I still don't understand the problem and I'm afraid I don't have time to
read a bio journal.  It takes a long time for me to understand a problem if
it's not broken down into standard statistical terminology.

That is understandable. 

So my question is:  what are the entries in the y[n] vectors for n in 1:N? 
How big is N that you want 1000 mixtures?

If they're counts, what are they counts of?  I can't connect the dots between 
counts and migrating proteins.  I'm guessing that you have N different
proteins, M places they can get stuck in an experiment, and y[n] represents
the counts of the places they get stuck over billions of replications for each
n in 1:N.  So sum(y[n]) is in the billions, but how big is N and how big is M?

Yes, exactly. Except y[n] is is not the count of the places, but the count of the number of proteins in a particular place. In figure 1a from the paper, there are 60 such places. In each of the 60 places, there are either 0 proteins (the protein is not found in the given place or could not be detected), or a number in the billions if the protein is present. Each y[n] thus holds a vector of counts for the given place, so that y[n,m] is the count of protein n in place m. 

And the goal's then to cluster the proteins by the distributions of where
they get stuck?  To what end?  Are you sure you want to run a Bayesian
model here?  

What we want to identify are the proteins that has the same profile / is drawn from the same mixture. If two proteins are co-migrating, they are most likely physically interacting in a complex. At least that is what we hypothesise. The idea running the Bayesian model would be to have priors on the number of mixtures, as well as infer which proteins are co-migrating. It could be there are better ways of doing it, but our initial thought was the Bayesian framework as it is quite flexible. 

I don't see how regression is related to clustering.

Definitely, it is not. It was just another approach used in a prior paper. 

Cheers,

David

[Bob Carpenter]

> On Mar 13, 2015, at 3:19 AM, David Westergaard <da...@harsk.dk> wrote:
>
> Yes, exactly. Except y[n] is is not the count of the places, but the count of the number of proteins in a particular place. In figure 1a from the paper, there are 60 such places. In each of the 60 places, there are either 0 proteins (the protein is not found in the given place or could not be detected), or a number in the billions if the protein is present. Each y[n] thus holds a vector of counts for the given place, so that y[n,m] is the count of protein n in place m.

Thanks --- that's what I was trying to say. I can't see
the paper (and please don't send a pirated copy to our list).

> What we want to identify are the proteins that has the same profile / is drawn from the same mixture. If two proteins are co-migrating, they are most likely physically interacting in a complex. At least that is what we hypothesise. The idea running the Bayesian model would be to have priors on the number of mixtures, as well as infer which proteins are co-migrating. It could be there are better ways of doing it, but our initial thought was the Bayesian framework as it is quite flexible.

Thanks --- now I get it.

So what I'd do from a Bayesian perspective is calculate the posterior
probability that two proteins, indexed n and n', are from the same
cluster given the observed data, i.e.,

Pr[protein[n] drawn from same cluster as protein[n'] | y]

You can do this with Stan.

The basic idea is that you have the mixture model you have, but
use a lognormal if you are going to log transform, then calculate
the above event probabilities in the generated quantities block:

generated quantities {
matrix[N,N] Pr_same_cluster;

for (m in 1:N)
for (n in (n+1):N)
Pr_same_cluster[m,n] <- pr_same(m,n,y,theta,mu,sigma);
}

And use a function pr_same() that calculates the expectation above.
For a mixture model, you can calculate using the same ps calculation
you use for the mixture,

log Pr[n in cluster k | y, theta, mu, sigma]

= ps[k] - log_sum_exp(ps)

For example, when you do log_sum_exp(ps), that's the unnormalized log
probabilities, so you need to divide each of the ps[iby the sum, which is another
so you just need to do



Then calculate

Pr[n and m in same cluster | y]

= SUM_k Pr[n in cluseter k | y, theta, mu, sigma]
* Pr[m in cluster k | y, theta, mu, sigma]

or

log Pr[...]

= log_sum_exp_{k in 1:K} (log Pr[n in k | ...] + log Pr[m in k | ...] );

Then I'd probably try to add some multivariate structure to the clustering,
but first things first.

And I don't think you'll need 1000 clusters --- how many clusters do you think
are actually operating?

The nice part about this is that the Pr[n and m in same cluster | ...] calcs
are not sensitive to label switching.

- Bob

[David Westergaard]

I think it somehow truncated this part, but I am assuming you are referring to the softmax function, yes?

Then calculate

 Pr[n and m in same cluster | y]

   = SUM_k Pr[n in cluseter k | y, theta, mu, sigma]
           * Pr[m in cluster k | y, theta, mu, sigma]

or

 log Pr[...]

    = log_sum_exp_{k in 1:K} (log Pr[n in k | ...] + log Pr[m in k | ...] );

Sometimes I forget the basics, especially when working in these higher-level frameworks. This is very close to what we are trying to infer. However, If we calculated the probability of two proteins coming from the same complex, how would we combine this information to assign multiple memberships per protein? Re-reading the formulation I gave in the prior email, I think I oversimplified it a bit. The co-migration is not restricted to two proteins, it can be anywhere between 2-50 proteins. So I guess the probability of interest is something in the line of,

  Pr[n in same cluster as n’, n’’, … | y]

can this be calculated downstream from the Pr[n in same cluster as n’| y]? Calculating it for all given combinations seems unfeasible. 

Digging deeper in the google scholar pages, I actually stumped on this article http://www.cs.cmu.edu/~sguennem/publications/ICDM2013a.pdf. It describes a mixed membership mixture model (I am assuming it is a legal copy, since it comes from the main authors web site.) Of course, this does not capture the generative process exactly, as each of the places are physically connected sequentially (i.e. 1-2-3-4-5, …). But practically, that might be OK. In theory, I should be able to implement the model in Stan using the same summarisation as for a normal mixture model, right?

Then I'd probably try to add some multivariate structure to the clustering,
but first things first.

And I don't think you'll need 1000 clusters --- how many clusters do you think
are actually operating?

If each of the clusters represent a physical complex (a set of proteins), and assuming each protein can be a member of multiple clusters, 1,000 is not an unreasonable high number. 

The nice part about this is that the Pr[n and m in same cluster | ...] calcs
are not sensitive to label switching.

David

[Bob Carpenter]

I probably just forgot to finish it. You probably don't need
all of softmax, which computes this value for each k.

>>
>> Then calculate
>>
>> Pr[n and m in same cluster | y]
>>
>> = SUM_k Pr[n in cluseter k | y, theta, mu, sigma]
>> * Pr[m in cluster k | y, theta, mu, sigma]
>>
>> or
>>
>> log Pr[...]
>>
>> = log_sum_exp_{k in 1:K} (log Pr[n in k | ...] + log Pr[m in k | ...] );
>>
>
>
> Sometimes I forget the basics, especially when working in these higher-level frameworks. This is very close to what we are trying to infer. However, If we calculated the probability of two proteins coming from the same complex, how would we combine this information to assign multiple memberships per protein?

The model implicitly assumes every protein comes from a single
mixture component. If you want to relax this, you need to go to
some kind of mixed-membership model like LDA.

> Re-reading the formulation I gave in the prior email, I think I oversimplified it a bit. The co-migration is not restricted to two proteins, it can be anywhere between 2-50 proteins. So I guess the probability of interest is something in the line of,
>
> Pr[n in same cluster as n’, n’’, … | y]
>
> can this be calculated downstream from the Pr[n in same cluster as n’| y]? Calculating it for all given combinations seems unfeasible.

Yes, anything combinatoric like this is intractable.
And there's no way of getting around that no matter what you
do. People often try to fudge this with something like k-means
or LDA, but those are so highly multimodal in their posteriors,
you can't actually do the proper inference. Every initial point
leads to a different local maximum and different inferences.

> Digging deeper in the google scholar pages, I actually stumped on this article http://www.cs.cmu.edu/~sguennem/publications/ICDM2013a.pdf. It describes a mixed membership mixture model (I am assuming it is a legal copy, since it comes from the main authors web site.)

That's usually not a good assumption legally. Computer
scientists are notorious IP pirates as a rule.

> Of course, this does not capture the generative process exactly, as each of the places are physically connected sequentially (i.e. 1-2-3-4-5, …). But practically, that might be OK. In theory, I should be able to implement the model in Stan using the same summarisation as for a normal mixture model, right?

I'm afraid I'm really not going to have time to do the background
reading here.

>> Then I'd probably try to add some multivariate structure to the clustering,
>> but first things first.
>>
>> And I don't think you'll need 1000 clusters --- how many clusters do you think
>> are actually operating?
>
> If each of the clusters represent a physical complex (a set of proteins), and assuming each protein can be a member of multiple clusters, 1,000 is not an unreasonable high number.

Then you're probably going to need some customized approximate
software to stand a chance at this problem. The combinatorics
will be well beyond what you can apply full Bayesian inference to.

- Bob

[David Westergaard]

Hi,

Sorry for bumping an old thread. I was not sure if I should make a new thread or not, so I just used this.

Anyway, in continuation of the discussion on mixture models, I was considering if it was possible to discard the principle of exchangeability. 

In the standard model, the features are exchangeable and independant:

p(y|theta, mu, sigma) = \sigma_i^K theta_i * N(y|mu_k, sigma_k)

But would it be possible to make the current feature dependant on the N last? Written in very poor math,

p(y|theta, mu, sigma, y_{t-1}) = \sigma_i^K theta_i * N(y|y_{t-1}, mu_k, sigma_k)

I found one article expanding the LDA model with a first order Markov property for the sequence of words (https://labs.yahoo.com/_c/uploads/ECMLPKDD2013MLJ.pdf). This model is, however, slightly different in that it assumes discrete features coming from a multinomial distribution with a parameter depending on the previous word. I am not sure how this translates into Gaussian distributions. 

I tried creating a very simple model, where I just add the previous observation (https://gist.github.com/demodw/b07a2eb962898e8703da). In Stan code,

log(theta[k]) + normal_log(y[n, m] + y[n, m-1], mu[k], sigma[k])

But I am not sure if that even makes sense. I guess my biggest problem is I find it difficult to formulate the likelihood function mathematically.

Thanks in advance.

David

[Bob Carpenter]

Yes, you add dependency on previous (latent) state --- you
get a time series or state-space model. See below for more info.

> On Sep 8, 2015, at 3:33 PM, David Westergaard <da...@harsk.dk> wrote:
>
> Hi,
>
> Sorry for bumping an old thread. I was not sure if I should make a new thread or not, so I just used this.
>
> Anyway, in continuation of the discussion on mixture models, I was considering if it was possible to discard the principle of exchangeability.
>
> In the standard model, the features are exchangeable and independant:
>
> p(y|theta, mu, sigma) = \sigma_i^K theta_i * N(y|mu_k, sigma_k)
>
> But would it be possible to make the current feature dependant on the N last? Written in very poor math,
>
> p(y|theta, mu, sigma, y_{t-1}) = \sigma_i^K theta_i * N(y|y_{t-1}, mu_k, sigma_k)
>

Indeed --- let's stick to formulas that make sense, and N(y | mu, sigma) only has
three arguments, so I don't know what to make of the last one. But see below.

> I found one article expanding the LDA model with a first order Markov property for the sequence of words (https://labs.yahoo.com/_c/uploads/ECMLPKDD2013MLJ.pdf). This model is, however, slightly different in that it assumes discrete features coming from a multinomial distribution with a parameter depending on the previous word. I am not sure how this translates into Gaussian distributions.
>
> I tried creating a very simple model, where I just add the previous observation (https://gist.github.com/demodw/b07a2eb962898e8703da). In Stan code,
>
> log(theta[k]) + normal_log(y[n, m] + y[n, m-1], mu[k], sigma[k])

You almost certainly would want to modify mu instead of y to stick to
a coherent generative Bayesian model.

To move to an AR(1)-like structure (see the manual sections on time series),
that'd look like this:

normal_log(y[n,m] | mu[k] + alpha[k] * y[n, m-1], sigma[k])

> But I am not sure if that even makes sense. I guess my biggest problem is I find it difficult to formulate the likelihood function mathematically.

Yes, that's where you need to start with Stan.

- Bob

[David Westergaard]

That’s very helpful. Thanks a lot. 

Did you intentionally use alpha as the coefficient for the AR(1) component? In my code I use alpha as a prior on the mixture components. I figure it was just an arbitrary choice, but just to be sure… What kind of priors would be “normal” to use? The uniform seems improper, the goto for linear models is the half-cauchy but I’m not sure it applies in the case of mixture models (or does it?).

Anyway, I plugged it into the Stan model (https://gist.github.com/demodw/3641605a0962dc953fb5), but now I get an error when I try to run it using CmdStan with the Variational inference engine:

Gradient evaluation took 4.34796 seconds

1000 iterations under these settings should take 4347.96 seconds.

Adjust your expectations accordingly!

  iter       ELBO   delta_ELBO_mean   delta_ELBO_med   notes

     0     -4e+05             1.000            1.000

Exception: stan::variational::normal_meanfield::set_mu: Input vector[1] is nan, but must not be nan!

Diagnostic information:

Dynamic exception type: std::domain_error

std::exception::what: stan::variational::normal_meanfield::set_mu: Input vector[1] is nan, but must not be nan!

Any debugging help would be appreciated.

Also, what would actually be the proper way to express this in maths? Reverse engineering your expression I come up with,

p(y | theta, mu, sigma, beta)_t = \sum_i^K theta_i * N(y_t | mu+beta[k]*y_{t-1}, sigma) for t = 2 .. T

Cheers,

David

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/dNyt1uXvaYI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

To post to this group, send email to stan-...@googlegroups.com.

[Bob Carpenter]

> On Sep 8, 2015, at 6:02 PM, David Westergaard <da...@harsk.dk> wrote:
>
> That’s very helpful. Thanks a lot.
>
> Did you intentionally use alpha as the coefficient for the AR(1) component? In my code I use alpha as a prior on the mixture components. I figure it was just an arbitrary choice, but just to be sure…

Sorry --- just an arbitrary choice.

> What kind of priors would be “normal” to use? The uniform seems improper, the goto for linear models is the half-cauchy but I’m not sure it applies in the case of mixture models (or does it?).

This is a question for Andrew or Ben --- if you call out a
more specific subject line, you might get their attention.

> Anyway, I plugged it into the Stan model (https://gist.github.com/demodw/3641605a0962dc953fb5), but now I get an error when I try to run it using CmdStan with the Variational inference engine:

This is a question for Alp or Dustin. I changed the subject here to
try to get their attention too.

>
> Gradient evaluation took 4.34796 seconds
> 1000 iterations under these settings should take 4347.96 seconds.
> Adjust your expectations accordingly!
>
> iter ELBO delta_ELBO_mean delta_ELBO_med notes
> 0 -4e+05 1.000 1.000
>
> Exception: stan::variational::normal_meanfield::set_mu: Input vector[1] is nan, but must not be nan!
> Diagnostic information:
> Dynamic exception type: std::domain_error
> std::exception::what: stan::variational::normal_meanfield::set_mu: Input vector[1] is nan, but must not be nan!
>
> Any debugging help would be appreciated.
>
> Also, what would actually be the proper way to express this in maths? Reverse engineering your expression I come up with,
>
> p(y | theta, mu, sigma, beta)_t = \sum_i^K theta_i * N(y_t | mu+beta[k]*y_{t-1}, sigma) for t = 2 .. T

That looks OK, but I'd put the _t on y, not on p(). And then there's a product
on the outside

p(y[2:T] | theta, mu, sigma, beta)

= PROD_{t in 2:T} SUM_{k in 1:K} (theta[k] * N(y[t] | mu + beta[k] * y[t-1], sigma)

And then you can add in some likelihood for y[1] if you want but it can't depend
on y[0] because there won't be one.

- Bob

[David Westergaard]

Shameless bump.

I tried to fit it with a reduced data set, and I am getting mixed results. Including only 2 rows of the 1457 x 38 matrix it runs without problems. When I add in the third row, for a 3 x 38 matrix, I am getting the same error as below (what I did practically was just duplicate the first row of matrix, but that should not matter). The data I used can be found at https://gist.github.com/demodw/1df3f5141819daea69e2 in R format. I suspect the initialisation values might be too much off, but I am not sure how to come up with a better guesstimate than Stan. I should maybe also add that this data has a lot of zero observations (either the feature was not present, or it was below detection limit), but I don’t suspect that should cause any trouble. 

Any input would be welcome.

Cheers,

David

[Dustin Tran]

Sorry for the delayed response! That informational message in ADVI comes up when it has a numerically unstable sample. The default should allow it to drop the point and continue resampling, so don’t worry too much about the informational message. It generally means that either the current initialization is bad (i.e. the informational message happens all the time in the beginning of the procedure), or the model is simply hard to fit.

If the initialization is bad, and there’s no way to fix it on the user end, then I suggest you up the number of grad_samples to make the gradient updates more reliable. Hope this helps.

Dustin

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.

[Bob Carpenter]

Is there a way to do something like trying smaller step sizes?
That's usually more stable numerically.

Dustin and Alp are working on more stable adaptation,
which should try smaller step sizes automatically.

- Bob

[David Westergaard]

I tried now with steadily increased grad_samples, going from 1 -> 5 -> 10 -> 15 -> 30 -> 50, and it is still failing.

I would not be too worried if it was just a warning, but the program actually stops running after spitting out the error message I pasted below.

Best,

David

[Bob Carpenter]

Have you tried smaller step sizes?

- Bob

[Dustin Tran]

Interesting. The algorithm shouldn't throw an error if it reaches something like that. I'll look into it. Will post possible solution/fix here as well (https://github.com/stan-dev/stan/issues/1610).

Dustin

On Tuesday, March 3, 2015 at 6:03:52 AM UTC-5, David Westergaard wrote:

Hi,

I would like to ask for some model guidance. I apologise in advance if my problem seems trivial, but my experience with Bayesian models is limited.

I am trying to identify mixtures in a population, in which each data point has its own profile. The tricky part is that each data point may be a member of multiple mixtures. From what I have found so far, this can be modelled using the so-called "Indian Buffet Process" (IBP) (Griffiths et al., 2011). However, as I have been able to read in the reference and on the board, Stan does not (yet) handle these kind of infinite models. From other threads on this forum, I have seen people have been mimicking the Chinese Restaurant Process using a very large number of mixtures. Can, or have, something similar been done for the IBP?

I am not even sure this is the best model to use, but I do not know of any better. To go a bit more into detail, the data is a NxM matrix, in which there are N data points and M measurements. The measurements are continous and (approximately) normally distributed (on log scale). Each data point has a lot of measurements with a zero value, simply because nothing was detected.

I have tried to come up with a complete description, but if anything is missing, please do ask!

Thanks in advance.

Cheers,

David