Approaching a multivariate ODE fit with missing data

[H. X.]

Hi All,

I having some trouble conceptualizing the pseudocode for a multivariate ODE fit with missing data. I read the univariate example in the documentation, but am not sure how to go about a case where data is missing in different times for each row of the matrix. Let's say, this simplified case, for example:

 _      _        _
| 2.1 1.2  x  3.5 | 
|_3.2  x  3.6 4.4_|

In pseudocode, how should I go about it? Anybody have some examples they can point me to?

Thanks

[Bob Carpenter]

ODEs aren't any different than other parts of Stan for
missing data.

The first thing to do is think about what the model would look
like if you had all the data. Then you take the missing data and
declare it as parameters. It's a strategy that's outlined in
the missing data chapter of the manual and will work for any
Bayesian application.

If you give us some more specifics of your problem, we might be
able to provide more concrete advice---I didn't know what to make
of two lines of data with some missing elements---where does it
fit into your ODE 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[H. X.]

Hi Bob,

So what I showed was not my actual data, but a toy example.
So right now, I have a soil carbon model that I want to fit to microbial biomass and carbon respiration data. Microbial biomass is a variable of the model. Carbon respiration is something I have to calculate for each time t from variable values. I am missing some microbial biomass data, and then some carbon respiration data, but not for the same times.

I read the missing data chapter of the manual, but did not have enough information to proceed from the given example to my ODE model since I am not sure how to create the matrix representing the missing data as parameters when the times for which carbon respiration and microbial biomass data are missing are different.

I have attached the model code, which is univariate right now.

Thanks

[H. X.]

Hi Bob,
Just hoping you could follow up when you had some extra time.

Thanks

[Daniel Lee]

Hi,

I'm taking a look, but I don't know how 

>  _      _        _ 
> | 2.1 1.2  x  3.5 | 
> |_3.2  x  3.6 4.4_| 

relates to your program.

Why does the first and second rows not line up? Are there observations at the same time? Why do you need to know the solution of the ode integration at the missing times?

Perhaps I'm not understanding your question. Can you provide a simpler working example?

Daniel

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Michael Betancourt]

There’s nothing special about the missing data here.  You just

solve the ODE for all measurement times, and then include a

measurement model only for those components that you do

measure at time point.  So sometimes you’d have a likelihood

contributes for each components, sometimes just a few.

If you want to impute those values then in generated quantities

you sample from the corresponding random number generator

(i.e. if your measurement model is Gaussian then use the

normal_rng) around the latent state.

[H. X.]

@Daniel,

The first example does not relate to my program. I was trying to draw a smaller example of a two row matrix of data I was trying to fit that was missing data in certain entries. For my actual example, which I will not post, there are many observations at the same time. The two rows of my actual data matrix consist of soil respiration and microbial biomass measurements. I am trying to fit soil respiration and microbial biomass solved by the model to those outputs. So in my integrate_ode line, I need to feed the function the list of times that I want output for. I imagine I feed it all my times -- the union of the set of times for soil respiration and microbial biomass. In the missing data example given in the Stan manual, it shows all missing times being gathered into one vector. However, I am not sure I would create a two row matrix representing missing times, considering that the missing times are different.

@Michael,
So I think I'm starting to get it. I'm going to put my thoughts in words pseudocode, and you can let me know if I am on the right track.

So I import some two-row matrix of data into Stan with columns corresponding to measurement times that is missing different entries on each row.

Let's say I am missing data in Data[1,3] and Data[2,5]. I will represent the missing as parameters. So in the parameters block, I will have

real<lower=0> missing1;
real<lower=0> missing2;

In the transformed parameters block, I would then state

Data[1,3] = missing1;  Data[2,5] = missing2;

And finally in my model block, I would have

missing1 ~ normal(0,1);
missing2 ~ normal(0,1);

Or whichever the ranges of values are and then finally have the line

Data ~ multi_normal(model_output_matrix, sigma);

[Bob Carpenter]

You need to do it the way it's outlined in the
chapter on missing data in the manual. So you have
to make that thing you assign to a transformed parameter.

But if it's missing states for the ODE, you just want
to do what Michael did, which is a much more efficient (both
computationally and statistically) way to code the same
quantities.

- Bob

> On Nov 18, 2016, at 9:35 PM, H. X. <yaois...@gmail.com> wrote:
>
> @Daniel,
>
> The first example does not relate to my program. I was trying to draw a smaller example of a two row matrix of data I was trying to fit that was missing data in certain entries. For my actual example, which I will not post, there are many observations at the same time. The two rows of my actual data matrix consist of soil respiration and microbial biomass measurements. I am trying to fit soil respiration and microbial biomass solved by the model to those outputs. So in my integrate_ode line, I need to feed the function the list of times that I want output for. I imagine I feed it all my times -- the union of the set of times for soil respiration and microbial biomass. In the missing data example given in the Stan manual, it shows all missing times being gathered into one vector. However, I am not sure I would create a two row matrix representing missing times, considering that the missing times are different.
>
> @Michael,
> So I think I'm starting to get it. I'm going to put my thoughts in words pseudocode, and you can let me know if I am on the right track.
>
> So I import some two-row matrix of data into Stan with columns corresponding to measurement times that is missing different entries on each row.
>
> Let's say I am missing data in Data[1,3] and Data[2,5]. I will represent the missing as parameters. So in the parameters block, I will have
>
> real<lower=0> missing1;
> real<lower=0> missing2;
>
> In the transformed parameters block, I would then state
>
> Data[1,3] = missing1; Data[2,5] = missing2;
>
> And finally in my model block, I would have
>
> missing1 ~ normal(0,1);
> missing2 ~ normal(0,1);
>
> Or whichever the ranges of values are and then finally have the line
>

> Data ~ multi_normal(model_output, sigma);

[Daniel Lee]

On Fri, Nov 18, 2016 at 9:46 PM, Bob Carpenter <ca...@alias-i.com> wrote:

You need to do it the way it's outlined in the
chapter on missing data in the manual.  So you have
to make that thing you assign to a transformed parameter.

But if it's missing states for the ODE, you just want
to do what Michael did, which is a much more efficient (both
computationally and statistically) way to code the same
quantities.

+1. Thanks, Bob. I was just in the middle of writing that.

 

- Bob

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

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

--
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+unsubscribe@googlegroups.com.

[H. X.]

So the pseudocode I wrote is based on the code from the missing data chapter in the manual.

What Michael says makes sense verbally and I am looking up some terms right now, like measurement model and latent state. I'm just having trouble picturing the code due to my lack of familiarity with statistics (which I am seeking to rectify). What would the measurement model consist of in pseudocode?

[Michael Betancourt]

Let’s go to a simpler example.  Let’s say that we have a two-dimensional

Gaussian with known covariance and we measure only the first component.

The corresponding Stan model would look something like

data {

  int N;

  vector[2] y[N];

  matrix[2, 2] Sigma;

}

parameters {

  real mu;

}

model {

  mu ~ normal(0, 1);

  for (n in 1:N)

    y[1] ~ normal(mu, Sigma[1, 1]);

}

In other words, we marginalize out over the unobserved components which

in the Gaussian case reduces to dropping any terms that depend on the

second component.

The important statistical concept here is that you can’t learn from unobserved 

data, all you can do it predict it.  So you want to build a model without the

unobserved components and if you really care about trying to figure out what

went on then you predict those components in the generated quantities block.

[Bob Carpenter]

> On Nov 18, 2016, at 11:00 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
> ....

> In other words, we marginalize out over the unobserved components which
> in the Gaussian case reduces to dropping any terms that depend on the
> second component.
>
> The important statistical concept here is that you can’t learn from unobserved
> data, all you can do it predict it. So you want to build a model without the
> unobserved components and if you really care about trying to figure out what
> went on then you predict those components in the generated quantities block.

You can't learn from completely unobserved data, but you can
learn a bit from partially observed data (e.g., one or more of many
predictors is missing, or the outcome is missing). But like
Michael says, the trick is to provide a joint model for everything
as if it were all observed, then marginalize out what's missing. Sometimes
the easiest way to do that is to declare the missing quantity as a parameter.
Then you get a posterior that includes it, but if you ignore it, it amounts
to marginalizing it out. Sometimes the model is neatly formed and you
can marginalize analytically (this is common in multivariate normal models).

- Bob

[H. X.]

Thanks for the help, gentlemen.

So Googling the phrase "marginalizing out," I am led to the Wikipedia page about marginal distributions. So it appears that marginal distributions are univariate distributions of variables by themselves in a joint distribution. They are different from conditional distributions because there is no contingency on other variables. To get marginal distributions, we "marginalize out" the other variable(s) as unknown(s).

So in simple terms, it sounds like what I should do in my case is:
1) Generate 2-row matrix solution for all times.
2) "Marginalize out" the unknown entries in the data matrix as parameters following some distribution.
3) Fit as I would usually.

Alternatively, (and forgive me if this is an absolutely idiotic and heretical question), what would be the effect of setting the missing data matrix entries to 0 and setting the corresponding entries of the model output matrix to 0, so as to not fit them?

[Bob Carpenter]

- Bob

> On Nov 19, 2016, at 5:00 PM, H. X. <yaois...@gmail.com> wrote:
>
> Thanks for the help, gentlemen.
>
> So Googling the phrase "marginalizing out," I am led to the Wikipedia page about marginal distributions. So it appears that marginal distributions are univariate distributions of variables by themselves in a joint distribution. They are different from conditional distributions because there is no contingency on other variables. To get marginal distributions, we "marginalize out" the other variable(s) as unknown(s).

Right---that's the marginalization we're talking about.
Though it's really more about removing one of the variables
than getting something down to univariate---that's just the
textbook examples. It's a pretty fundamental thing to understand
for applied and computational stats.

There's a chapter in the manual on latent discrete parameters
that goes through how to marginalize those out.

> So in simple terms, it sounds like what I should do in my case is:
> 1) Generate 2-row matrix solution for all times.
> 2) "Marginalize out" the unknown entries in the data matrix as parameters following some distribution.
> 3) Fit as I would usually.
>
> Alternatively, (and forgive me if this is an absolutely idiotic and heretical question), what would be the effect of setting the missing data matrix entries to 0 and setting the corresponding entries of the model output matrix to 0, so as to not fit them?

What happens is that you get inferences as if the values
were 0---they won't be ignored.

- Bob

[xi...@uci.edu]

Hi Bob,

So I want to try both approaches in terms of declaring the missing data elements as parameters or ignoring them. When I ignore the missing data, I need to get rid of the corresponding the matrix entries of my ODE solution. This is easy in a row vector, when I can just eliminate elements and shrink the length of the ODE solution vector. However, for a bivariate matrix in which different elements and different numbers of elements are missing on each row, I can't just delete elements, since it is obviously impossible to have a matrix with different sized rows. What would be the protocol for ignoring elements in this multivariate case? This is a question I was trying to get at earlier, but could not mentally phrase.

[Bob Carpenter]

Just allocate a dense N x 2 matrix and ignore
the entries you don't need. You can encode
the indexes to ignore, or even create an N x 2
matrix of 1 and 0 values to indicate which ones
are included.

If the bivariate distribution is easily marginalized
into univariate distributions, like say a normal,
you'd use those on the observed components. So
that'd be:

vector[2] y[N];
int<lower=0, upper=1> y_observed[N, 2];

for (n in 1:N) {
if (y_observed[n, 1] && y_observed[n, 2])
y[n] ~ multi_normal(mu, Sigma);
else if (y_observed[n, 1])
y[n, 1] ~ normal(mu[1], sqrt(Sigma[1, 1]));
else if (y_observed[n, 2])
y[n, 2] ~ normal(mu[2], sqrt(Sigma[2, 2]));
}

But... that's not very efficient. It's much better
to vectorize all of this. So, in transformed data,
build up three vectors of indices, for the three cases
above:

transformed data {
int ns12[observed_12(y_observed)];
int ns1[observed_1(y_observed)];
int ns2[observed_2(y_observed)];

you need to do this with functions to count because you need to declare
sizes. Then the rest of transformed data just fills (I'm going to
use the syntax about to release in 2.13):

int n12 = 1;
int n1 = 1;
int n2 = 1;
for (n in 1:N) {
if (y_observed[n, 1] && y_observed[n, 2]) {
ns12[n12] = n;
n12 = n12 + 1;
} else if (y_observed[n, 1]) {
ns1[n1] = n;
n1 = n1 + 1;
} else if (y_observed[n, 2]) {
ns2[n2] = n;
n2 = n2 + 1;
}
}

and then, in the model block, everything's nice and vectorizable
using those indexes constructed once in transformed data:

y[ns12] ~ multi_normal(mu, Sigma);
y[ns1] ~ normal(mu[1], sqrt(Sigma[1, 1]));
y[ns2] ~ normal(mu[2], sqrt(Sigma[2, 2]));

This is much more efficient than using latent variables for the
missing data, but requires the multivariate distribution to be marginalized
analytically.

I'll add this discussion to the missing data chapter:

https://github.com/stan-dev/stan/issues/2122#issuecomment-263031536

- Bob

[H. X.]

Thanks Bob, this was insanely helpful!

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

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

--
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/2eRJKICL81c/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[H. X.]

Since I'm going to do this in 2.12, what are the 2.13 syntax updates so I know not to include them in my code?

[H. X.]

Also Bob, in your example code, why did you instantiate y in your example code as a multi-dimensional vector rather than as a matrix?

[H. X.]

Actually, disregard this thread. Enough issues are cropping up and this one is getting too long I'll create a new one.

--
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/2eRJKICL81c/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

I'll collect up a few answers inline below.

> On Nov 26, 2016, at 12:42 AM, H. X. <yaois...@gmail.com> wrote:
>
> Since I'm going to do this in 2.12, what are the 2.13 syntax updates so I know not to include them in my code?

The compound declare-define is not in 2.12

int n12 = 1;
int n1 = 1;
int n2 = 1;

> On Nov 26, 2016, at 6:22 PM, H. X. <yaois...@gmail.com> wrote:
>
> Also Bob, in your example code, why did you instantiate y in your example code as a multi-dimensional vector rather than as a matrix?

That'll allow vectorization most easy. Then
y[ns12] and so on become arrays of vectors, which
is what you need to vectorize the multivariate normal.

> Actually, disregard this thread. Enough issues are cropping up and this one is getting too long I'll create a new one.

Too late!

- Bob