Help parameterizing a constrained simplex?

[David Manheim]

Hi,

Still working on my dissertation model, and I have a question about how best to parameterize the prior for the epidemiological "contact matrix." (It runs now, but it's not at all elegant, and I'm hoping there are nice solutions to do this.)

This type of matrix matrix has entries c(a,b) that are the number of times an individual in group a is in contact with an individual in group b. It is constrained by the fact that contact is symmetric - though groups sizes are not equal. This means c(b,a) = c(a,b)*size(a)/size(b).

As an additional concern, I had earlier changed the contact matrix to be a simplex that was then multiplied by a fixed value, because allowing the total number of contacts to vary was interacting with the estimated spread variable. (More viral flu with fewer contacts looks exactly like a less viral flu with more contacts.) I found (unsurprisingly) that Stan was arriving at implausible contact matrix simplex posteriors. Using a static prior contact matrix was much faster than having Stan estimate the values, but I'd prefer to allow it to estimate the values based on the provided values as priors.

So, I am attempting to figure out how I can parameterize the model to include a simplex that obeys a constraint. To do this, I use a simplex for the upper triangular part of the matrix, then enforce the constraint on the other half, and then normalize the entire thing all together.

The relevant parts of my attempt (excluding the actual infectious disease model parts,) which seems to work but is very kludge-y:

Data:

  real<lower=0> contact_matrix[groups,groups]; // These are unknown, but depend on the group types. This (approximate) matrix here will be used in the prior.

Transformed Data:

  real Upper_Contact_Matrix[groups*(groups+1)/2];

{

int upper_counter;

upper_counter <- 1;

Contact_Matrix_Normalizer <- 0.0;

for (i in 1:groups) {

for (j in 1:groups) {

                if (i<=j){

Upper_Contact_Matrix[upper_counter]<-contact_matrix[i,j];

upper_counter <- upper_counter + 1;

}

Contact_Matrix_Normalizer <- Contact_Matrix_Normalizer + contact_matrix[i,j]; //This is based on the total

}

}

}

Parameters:

  simplex[groups*(groups+1)/2] contact_matrix_vector; //Upper triangular portion only.

Transformed Parameters:

  real upper_simplex_pct; //The percentage of the total in the upper simplex - used for normalizing later.

  real<lower=0> Lower_Triangular_Simplex[groups*(groups-1)/2];

  {

  real contact_matrix_temp[groups,groups];

  int upper_counter;

  int lower_counter;

  upper_counter <- 1;

  lower_counter <- 1;

  for (i in 1:groups) {

    for (j in 1:groups) {

      real inverse_value;

 inverse_value <- (-1); //If my logic is bad, this will violate constraints later.

      if (i>j) {

   contact_matrix_temp[i,j] <- contact_matrix_temp[j,i]*pop_sizes[j]/pop_sizes[i]; //This creates the contact matrix identity.

   Lower_Triangular_Simplex[lower_counter] <- contact_matrix_temp[i,j];

        lower_counter <- lower_counter + 1;

 }

 if(i<=j) {

 contact_matrix_temp[i,j] <- contact_matrix_vector[upper_counter];

 upper_counter <- upper_counter + 1;

 }

    }

  }

  upper_simplex_pct <- sum(contact_matrix_vector)/(sum(Lower_Triangular_Simplex)+sum(contact_matrix_vector));

  } //End contact matrix creation.

model:

  {

  real contact_matrix_temp[groups,groups];

  int theta_counter;

  int matrix_counter;

  theta_counter <- 1;

  matrix_counter <- 1;

  for (i in 1:groups) {

for (j in 1:groups) {

        contact_matrix_temp[i,j] <- if_else(i<=j,contact_matrix_vector[matrix_counter],contact_matrix_temp[j,i]*pop_sizes[j]/pop_sizes[i]);

        if (i<=j){

          matrix_counter <- matrix_counter + 1;

        }

theta[y, theta_counter] <- contact_matrix_temp[i,j]*Contact_Matrix_Normalizer;

        theta_counter<-theta_counter+1;

}

}

  }

contact_matrix_vector ~ dirichlet(to_vector(Upper_Contact_Matrix)/(Contact_Matrix_Normalizer*upper_simplex_pct));

[Bob Carpenter]

> On Jul 28, 2016, at 2:30 PM, David Manheim <davidm...@gmail.com> wrote:
>
>
> Hi,
>
> Still working on my dissertation model, and I have a question about how best to parameterize the prior for the epidemiological "contact matrix." (It runs now, but it's not at all elegant, and I'm hoping there are nice solutions to do this.)
>
> This type of matrix matrix has entries c(a,b) that are the number of times an individual in group a is in contact with an individual in group b. It is constrained by the fact that contact is symmetric - though groups sizes are not equal. This means c(b,a) = c(a,b)*size(a)/size(b).

I couldn't follow this. Is c(a, b) a count in {0, 1, 2, ...} or a
probability in (0, 1)?

> As an additional concern, I had earlier changed the contact matrix to be a simplex that was then multiplied by a fixed value, because allowing the total number of contacts to vary was interacting with the estimated spread variable. (More viral flu with fewer contacts looks exactly like a less viral flu with more contacts.) I found (unsurprisingly) that Stan was arriving at implausible contact matrix simplex posteriors. Using a static prior contact matrix was much faster than having Stan estimate the values, but I'd prefer to allow it to estimate the values based on the provided values as priors.
>
> So, I am attempting to figure out how I can parameterize the model to include a simplex that obeys a constraint. To do this, I use a simplex for the upper triangular part of the matrix, then enforce the constraint on the other half, and then normalize the entire thing all together.

Just because the empirical probabilities will obey
those constraints doesn't mean that the underlying parameters
have to. What is controlling size(a) and size(b)?

> The relevant parts of my attempt (excluding the actual infectious disease model parts,) which seems to work but is very kludge-y:
>

I'm afraid I couldn't follow the code, but then I'm not
exactly sure what you're trying to do.

- 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[David Manheim]

Hi Bob,

Some answers below - sorry this isn't clear.

On Sunday, July 31, 2016 at 6:31:42 PM UTC-4, Bob Carpenter wrote:

> On Jul 28, 2016, at 2:30 PM, David Manheim <davidm...@gmail.com> wrote:
>
>
> Hi,
>
> Still working on my dissertation model, and I have a question about how best to parameterize the prior for the epidemiological "contact matrix." (It runs now, but it's not at all elegant, and I'm hoping there are nice solutions to do this.)
>
> This type of matrix matrix has entries c(a,b) that are the number of times an individual in group a is in contact with an individual in group b. It is constrained by the fact that contact is symmetric - though groups sizes are not equal. This means c(b,a) = c(a,b)*size(a)/size(b).

I couldn't follow this.  Is c(a, b) a count in {0, 1, 2, ...} or a
probability in (0, 1)?

c(a,b) is the matrix entry in position a,b. It is a count, not a probability - but I end up transforming it to a simplex times a fixed value because I want the total size constrained. (Perhaps I should do that more directly, but I'm unsure how without using a simplex.)

 

> As an additional concern, I had earlier changed the contact matrix to be a simplex that was then multiplied by a fixed value, because allowing the total number of contacts to vary was interacting with the estimated spread variable. (More viral flu with fewer contacts looks exactly like a less viral flu with more contacts.) I found (unsurprisingly) that Stan was arriving at implausible contact matrix simplex posteriors. Using a static prior contact matrix was much faster than having Stan estimate the values, but I'd prefer to allow it to estimate the values based on the provided values as priors.
>
> So, I am attempting to figure out how I can parameterize the model to include a simplex that obeys a constraint. To do this, I use a simplex for the upper triangular part of the matrix, then enforce the constraint on the other half, and then normalize the entire thing all together.

Just because the empirical probabilities will obey
those constraints doesn't mean that the underlying parameters
have to.  What is controlling size(a) and size(b)?

size(a) and size(b) are data elements. (In the below code, I neglected to include "pop_sizes" as a data element.)

I need the parameters to follow some constraints because without doing so, they don't converge well (the n_eff is small, and the values are unreasonable - like group 1 contacts group 3 17 times daily, group 3 contacts group 1 0.0001 times daily. This is partly because using a simplex without constraints is making the total fixed, so, I would think, mass gets assigned to whichever item helps most locally, away from other ones.)

 

> The relevant parts of my attempt (excluding the actual infectious disease model parts,) which seems to work but is very kludge-y:
>

I'm afraid I couldn't follow the code, but then I'm not
exactly sure what you're trying to do.

I'm trying to figure out how to parameterize a prior over the contact matrix that obeys these constraints. I did so, below, by making the upper triangular elements be parameters, then trying to enforce the constraint on the matrix in the transformed data block. To get the total number of contacts size constraint to work, I need to track what percentage of the total is in the upper triangular portion of the matrix, so that I can multiply by the simplex by the correct amount in the model block.

[Bob Carpenter]

I still don't get it, but that's probably just me. Maybe
someone else has a better idea of what you're up to and can
help.

But if it runs the way you want, I don't see a better way to
code it up. I'd check with simulated data to make sure you
can recover parameters within posterior quantiles well.

- Bob

[andre.p...@googlemail.com]

Hi David,

I'm with Bob. Not sure I understand you right.

What I get is that you want to constraint a matrix containing counts.

This is what you mentioned with simplex. The simplex construct in 

STAN represents a special vector of probabilities that sum to 1. 

And now as I hope understand right - you have a constraint not

necessarily 1 and its entries are real values. 

For clarity, the following is not highly optimized STAN code

This could be done by STANs softmax function plus

vectorization. 

Lets create a parameter theta which creates a multinomial

'simplex' pr.theta, each 0 <= pr.theta[i] <= 1

and sum(pr.theta) == 1.

sigma.theta ~ cauchy(0, 2.5);

theta[1] ~ normal(0, sigma.theta)
theta[2] ~ normal(0, sigma.theta)
theta[3] ~ normal(0, sigma.theta)

pr.theta[1] <- 1/(1+exp(theta[1])+exp(theta[2])+exp(theta[3]))
pr.theta[2] <- exp(theta[1]) * pr.theta[1]
pr.theta[3] <- exp(theta[2]) * pr.theta[1]
pr.theta[4] <- exp(theta[3]) * pr.theta[1]

If you now want to change the upper limit to something else, 

just change

 pr.theta[1] <- 1/(1+exp(theta[1])+exp(theta[2])+exp(theta[3]))

by

 pr.theta[1] <- my_fancy_upper_limit /(1+exp(theta[1])+exp(theta[2])+exp(theta[3]))

Take a look here to:
http://www.mathworks.com/help/stats/multinomial-models-for-nominal-responses.html

Another option you have is to specify a Dirichlet distribution upon a simplex in STAN and 

scale it. 

http://stackoverflow.com/questions/24349692/dirichlet-multinomial-winbugs-code

Andre

[Michael Betancourt]

Think of the multinomial distribution — you have a collection of counts

in bins, n_{i}, with the constraint that \sum_{i = 1}^{I} n_{i} = N.  In the

limit of each n_{i} being large, then this can be reasonably approximated

by taking a simplex and scaling each value by N (which was originally

suggested and is equivalent to the softmax construction in the previous

email).  If _any_ of the n_{i} are small, however, then this approximation

falls apart and bad things happen.

Ultimately, however, this is a common mistake of not thinking generatively.

Counts are _data_, and we don’t fit data.  We fit a generative model of

how the data were produced.  So the relevant “contact” matrix isn’t over

counts, it’s over probabilities that could induce said counts.  And those

probabilities can often be easily modeled with simplices.

For example, each column could be a simplex with the corresponding

likelihood a mutlinomial distribution over the counts corresponding to

that column.  Same could be true for the rows.  These are known as

left and right stochastic matrices.

Things become more difficult if you have a double stochastic matrix,

where the rows _and_ columns have to sum to one.  There the matrix

doesn’t decouple into a collection of simplices and you’d have to be

more careful in the implementation (not sure how to do it off the top

of my head).

On Aug 2, 2016, at 2:20 AM, andre.pfeuffer via Stan users mailing list <stan-...@googlegroups.com> wrote:

Hi David,

I'm with Bob. Not sure what  I understand you right.

[David Manheim]

OK, so I am trying something like the suggested approaches, and it seems to keep the constrain, but seemingly samples much faster. (I'm running a larger sample now to see if that is really true, but so far, it's much better.)

I think it's essentially a left stochastic matrix, from what Michael described. Instead of having a matrix, I have a simplex of the percentage of total contacts per group, and a simplex of contacts in each group to each other group (i.e. the percentage of that group's total). I multiply these together, and multiply the product by the sum of the number of contacts.

This enforces only the total constraint. For the contact matrix constraint, I don't use some of those entries - I only use the upper triangular matrix entries, and replace the ones corresponding with the remainder with the values implied by the constraint. Now, I have a bunch of free variables as parameters, which seems non-ideal, but it works. (I checked the resulting matrix, and it follows my constraints.)

(Of course, I'd love to do something closer to what Michael suggested, thinking generatively, about how the contact matrix is a function of some underlying process. This kind-of does that, but I'm unsure how I would decompose my prior about relative propensity to interact, instead of simply providing the prior.)

Thanks,

David

[Bob Carpenter]

> On Aug 2, 2016, at 4:31 AM, Michael Betancourt <betan...@gmail.com> wrote:
>
> Think of the multinomial distribution — you have a collection of counts
> in bins, n_{i}, with the constraint that \sum_{i = 1}^{I} n_{i} = N. In the
> limit of each n_{i} being large, then this can be reasonably approximated
> by taking a simplex and scaling each value by N (which was originally
> suggested and is equivalent to the softmax construction in the previous
> email). If _any_ of the n_{i} are small, however, then this approximation
> falls apart and bad things happen.
>
> Ultimately, however, this is a common mistake of not thinking generatively.
> Counts are _data_, and we don’t fit data. We fit a generative model of
> how the data were produced. So the relevant “contact” matrix isn’t over
> counts, it’s over probabilities that could induce said counts. And those
> probabilities can often be easily modeled with simplices.

That's exactly what I was thinking and trying to suggest, but
I couldn't quite follow what was data.

> For example, each column could be a simplex with the corresponding
> likelihood a mutlinomial distribution over the counts corresponding to
> that column. Same could be true for the rows. These are known as
> left and right stochastic matrices.

I think of these as transition matrices for Markov processes.

- Bob

[Bob Carpenter]

I'm still not sure your decomposition makes sense in that
if you have these symmetry constraints, it's not underlyingly
independent simplexes per item.

Could you take a step back and explain what this contact
matrix is and what it represents? Is it data? If so, what
is its shape and what do the entries represent?

- Bob

[David Manheim]

Hi Bob,
Sure.

Disease transmission happens when an infected person and a noninfected person come in contact with one another, with some probability p that transmission occurs during that contact. (That probability is a key parameter in infectious disease modeling, which I need to estimate.)

In simple compartmental Epi models, we have a set of groups, or compartments, and we assume that these groups are homogenous. We also assume that the rate at which people contact one another is independent of disease status, so we find the underlying contact rate for different groups. The number of infected individuals in group A at t+1 is the sum over all groups of: ( Percentage of individuals in group i who are sick at time t, multiplied by percentage of individuals in group A that are susceptible at time t, multiplied by the number of contacts between group i and group A per time period, multiplied by the probability p that transmission occurs per contact.)

The contact rate in basic Epi models is frequently assumed or given, but is uncertain. I am attempting to estimate it, based on a prior about what these rates are. The contact rate between each pair of groups is summarized as a matrix - that's what I'm trying to find right now.

2 Constraints exist; 

The matrix is governed by a simple constraint that contacts are symmetric; every contact in group 1 with someone in group 2 is also an contact between the person in group 2 with the person in group 1 - but because the groups are potentially different sizes, the two contact rates are not equal. Instead, you have the reciprical constraint; contact(i,j) = contact(j,i) * size(j)/size(i). 

The matrix also interacts with the probability of transmission, p, since, for example, doubling one and halving the other leaves the outcome the same. In order to control for that, I keep the sum of contacts constant in my model.

I need to model the constraints in order to find a posterior estimate of both the contact matrix and the probability of transmission.

Hope this is clearer,

David Manheim

[Bob Carpenter]

Thanks. That helps a lot. But it doesn't tell me what
data it is that you are trying to model, and that's what we
need to know to think of this whole process generatively.

Is it the number of infected individuals per group at each time t
or some measure on subgroups of each group?

- Bob

[Bob Carpenter]

I thought about this some more and can't make things add up.

Shouldn't we be be modeling the chance that an uninfected
group member gets infected and the chance that an infected group
member gets over the infection? Then you don't need to worry
about symmetry because only infected individuals contacting uninfected
individuals matters for creating new infections.

Shouldn't we adjust for possible double counting where an
individual has contact with more than one infected individual?
Is it just that these adjustments for coherence are so small
and the percentage of infections, infection rate, and transmission
rates are so low that it doesn't affect the answer much?

I don't see how making the sum of contacts constant helps you
identify transmission rate; you just fix a point arbitrarily
for one parameter to identify the other. Does this get corrected
in some kind of posterior summary?

Also, I'm still not sure what the data is that's being measured.

- Bob

[David Manheim]

Hi Bob,

The model includes both the probability that an uninfected member gets infected, and the rate at which they recover. The reason symmetry is important is that infected people in group 1 infect susceptible people in group 2 at a rate proportional to their contacts - and infected people in group 2 infect susceptible people in group 1 at a rate proportional to the contact. These two processes are constrained, because contact is symmetric. Regarding double counting, this is a continuous diffeq model, so the issue is simplified away; fractional people are infected, and the overlap is incorporated.

Keeping the sum of contacts constant primarily helps with identifiability in the bayesian approximation model. The summary statistic we care about is the R0, the average number of newly infected individuals per infected individual, in a susceptible population - and this is a function of both the contact matrix and the rate.

The data I'm actually modeling isn't included here - but, in part, it is the number of infections observed in each group the subset of patients that show up at clinics that are part of the influenza disease surveillance network. This is an output of the differential equation model, which is based on Michael Betancourt's earlier model (which doesn't contain any subgroups, so it doesn't have these problems.)

Happy to answer any further questions - especially if they help me think about modeling this generatively, instead of my current approach. (But in the mean time, I'm still hoping for a more efficient way to model the contact matrix and its constraints.)

Cheers,

David

[andre.p...@googlemail.com]

Dear David,

considering the advises from Bob - Markov Chain - and Michael to model counts,

and assuming that your counts are integers you could model them via 

a discrete count distribution (poisson, GPD, negbin) the resulting lambda's 

could form a latent structure corresponding mu's, link: lambda = exp(mu). 

Since you want probabilities for sure, you could get them by la_i / sum(la)

Take a look here:

https://onlinecourses.science.psu.edu/stat504/node/48

Also worth to have a look at the alligator model:

Especially the link between the Poisson and multinomial:

http://www.openbugs.net/Examples/Aligators.html

If your counts are contineous have a look at the gamma distribution.
Andre

[David Manheim]

Hi Andre,

The counts are continuous, and so I'm unsure how modeling them in any of the ways you mention makes it computationally efficient to implement the constraints, or to sample.

It may be, however, that it would be more elegant to use the poisson/multinomial connection you noted to parameterize the matrix than my current method.

On the other hand, it might be more effcient to do something similar with a normal distribution instead of the Poisson distribution, and then normalize - but I'm unsure if this would be more efficient. Hence why I'm here asking questions - ;) 

Thanks,

David

[andre.p...@googlemail.com]

Hi David,

you could also consider the log-normal distribution, to start with something easy.

Constrains are to be avoided. They are translated in STAN anyway. See in the manual 

of STAN.

I would suggest to model the "counts"  via distribution function, check the model with plots.

Look at the the trace plots, pairs plot, if the model fits nicely. 

At the latent structure:

First start with a simple intercept. Then exchange the intercept by a Markov chain...

hierarchical model, etc. or what makes sense.

Then improve at the distribution function level. 

Deriving the probability is simple math.

Thats just my 2 cents,

Andre

[Bob Carpenter]

If they're continuous, they're not counts!

- Bob

> On Aug 4, 2016, at 12:34 PM, David Manheim <davidm...@gmail.com> wrote:
>
> Hi Andre,
>

> The counts are continuous, and so I'm unsure how modeling them in any of the ways you mention makes it computationally more efficient to implement the constraints, or to sample.

>
> Thanks,
> David
>
> On Thursday, August 4, 2016 at 6:22:09 AM UTC-4, andre.p...@googlemail.com wrote:
> Dear David,
>
> considering the advises from Bob - Markov Chain - and Michael to model counts,
> and assuming that your counts are integers you could model them via
> a discrete count distribution (poisson, GPD, negbin) the resulting lambda's
> could form a latent structure corresponding mu's, link: lambda = exp(mu).
>
> Since you want probabilities for sure, you could get them by la_i / sum(la)
> Take a look here:
>
> https://onlinecourses.science.psu.edu/stat504/node/48
>
> Also worth to have a look at the alligator model:
> Especially the link between the Poisson and multinomial:
>
> http://www.openbugs.net/Examples/Aligators.html
>
> If your counts are contineous have a look at the gamma distribution.
> Andre
>
>