bivariate censoring

[Tamas Papp]

Hi,

I would like to ask for some help in formulating the following model in
Stan.

I have N individuals, indexed with i=1, ..., N. I have two _latent_
variables w and z, constrained to be > 0, so I would initially start my
analysis with lognormal, eg

for (i in 1:N) {
w[i] ~ lognormal(mu_w, sigma_w);
z[i] ~ lognormal(mu_z, sigma_z);
}

So far, nothing is observable. The two observables y and n have the
following data generating process: there is a parameter m > 0, and when

z[i]/w[i] > m, n[i]=0 and y[i]=0;

otherwise n[i] is a linear function of z[i]/w[i] and y[i] is a linear
function of z[i] and w[i] (I need the coefficients, and I can always add
a lognormal error etc).

My problem is with the part when z[i]/w[i] > m. It is clear that I need
to separate the observations into two groups, one where n[i]=y[i]=0, the
other is when it isn't.

But I don't know how to tell Stan that in the first group, observations
for which z/w < m have a 0 likelihood --- should I somehow set lp__ to
minus infinity conditional on this happening?

Best,

Tamas

[Bob Carpenter]

If you set the log probability to -infinity, it will do rejection.
It's better to actually implement the constraints. The problem is
that Stan doesn't let you separately implement constraints for
each variable.

So you have to enforce this constraint yourself:

> z[i]/w[i] > m

Rewriting, that's

z[i] > m / w[i]

So you need a lower-bound on z[i] equal to m / w[i].

The way to do this is just do what Stan does internally for
transformed parameters. Start with a raw, unconstrained parameter,
then transform it and apply the Jacobian adjustment:

parameters {
vector[N] z_raw;
}
transformed parameters {
vector[N] z;

for (i in 1:N)

z[i] <- exp(z_raw) + m / w[i];
increment_log_prob(z_raw); // adds complete Jacobian
}
model {
z ~ ...; // whatever prior you want here
}

Here's the relevant log Jacobian calculation:

log abs( d/dz_raw exp(z_raw) )
= log abs(exp(z_raw))
= z_raw

You can make the whole thing explicit in a big Jacobian
matrix which has identity transforms for m and w.

Then you can put whatever prior on z you want.

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

[Tamas Papp]

On Mon, Aug 18 2014, Bob Carpenter <ca...@alias-i.com> wrote:

> If you set the log probability to -infinity, it will do rejection.

Just out of curiosity (not that I indend to do it, the reparametrization
below is better), how would I do that? I thought of

lp__ <- negative_infinity()

but lp__ is deprecated, so I am not sure about this. Is there a function
that sets lp__ to -infinity, _and_ possibly exits early from the model
posterior calculation? If not, will just open an issue.

Some clarifying questions on the approach you suggested, below:

> It's better to actually implement the constraints. The problem is
> that Stan doesn't let you separately implement constraints for
> each variable.
>
> So you have to enforce this constraint yourself:
>
>> z[i]/w[i] > m
>
> Rewriting, that's
>
> z[i] > m / w[i]

For the archives, that is z[i] > m * w[i], but I get the point.

> So you need a lower-bound on z[i] equal to m / w[i].
>
> The way to do this is just do what Stan does internally for
> transformed parameters. Start with a raw, unconstrained parameter,
> then transform it and apply the Jacobian adjustment:
>
> parameters {
> vector[N] z_raw;
> }

Just to clarify: I should use an unconstrained variable (or vector)
because Stan will automatically assign a flat (\propto 1) prior to it,
so it will not mess up my likelihood. And if I wanted to use something
like

parameters {
real<lower=0> z_raw[N];

}
transformed parameters {
vector[N] z;
for (i in 1:N)

z[i] <- z_raw + m / w[i];
// no need to increment lp__ with Jacobian correction, transformation linear

}
model {
z ~ ...; // whatever prior you want here
}

then (maybe some future version of) Stan may use a prior that would add
transformed terms to my posterior? Or would this be innocuous?

Best,

Tamas

[Tamas Papp]

Hi Bob,

I tried something like you suggested on simulated data, and everything
works out except for the truncation parameter. In particular, my data
generating process is

--8<---------------cut here---------------start------------->8---
N <- 10000
z <- rlnorm(N)
w <- rlnorm(N)
d <- data.frame(z=z, w=w)
d[, "z_hat"] <- d$z * rlnorm(N, 0, 0.1)
d[, "w_hat"] <- d$w * rlnorm(N, 0, 0.1)
d <- d[d$z/d$w > 1, ] # m=1
--8<---------------cut here---------------end--------------->8---

ie z, w are unit lognormal, only a noisy version z_hat, w_hat is
observed, data is kept iff z/w > m = 1.

And the corresponding stan model is

--8<---------------cut here---------------start------------->8---
data {
int<lower=0> N;
real<lower=0> z_hat[N];
real<lower=0> w_hat[N];
}

parameters {
real w_aux[N];
real z_aux[N];
real m_aux;
}

transformed parameters {
real m;
real w[N];
real z[N];
m <- exp(m_aux);

for (i in 1:N) {

w[i] <- exp(w_aux[i]);
z[i] <- exp(z_aux[i]) + m*w[i];
}
}

model {
increment_log_prob(m_aux);

for (i in 1:N) {

z[i] ~ lognormal(0.0, 1.0);
w[i] ~ lognormal(0.0, 1.0);
z_hat[i]/z[i] ~ lognormal(0, 0.1);
w_hat[i]/w[i] ~ lognormal(0, 0.1);
increment_log_prob(z_aux[i]);
increment_log_prob(w_aux[i]);
}
}
--8<---------------cut here---------------end--------------->8---

I run and check it with

--8<---------------cut here---------------start------------->8---
model <- stan_model(file = "bivariate_censoring1.stan")
fit <- sampling(model,
data =
list(N = nrow(d),
z_hat = d$z_hat,
w_hat = d$w_hat),
iter=1000,
chains=1)
posterior <- extract(fit, permuted = TRUE)
plot(d$z, colMeans(posterior$z))
abline(a=0, b=1, col="gray50")
plot(d$w, colMeans(posterior$w))
abline(a=0, b=1, col="gray50")
hist(posterior$m)
abline(v=1, col="gray50")
--8<---------------cut here---------------end--------------->8---

and everything checks out except that the posterior mean of m is well
below 1, and also, it depends on the sample size. I know that censoring
boundaries are tricky to identify well, but I am still wondering if I
miscoded something.

Also, I get a couple of
--8<---------------cut here---------------start------------->8---
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Error in function stan::prob::lognormal_log(N4stan5agrad3varE): Random variable is -nan:0, but must not be nan!
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
--8<---------------cut here---------------end--------------->8---
messages but not more than 1-2 per chain, so I think these can be
ignored safely.

Best,

Tamas

On Mon, Aug 18 2014, Bob Carpenter <ca...@alias-i.com> wrote:

[Bob Carpenter]

I'm answering both e-mails inline below.

On Aug 18, 2014, at 7:25 AM, Tamas Papp <tkp...@gmail.com> wrote:

> On Mon, Aug 18 2014, Bob Carpenter <ca...@alias-i.com> wrote:
>
>> If you set the log probability to -infinity, it will do rejection.
>

> Just out of curiosity (not that I indend to do it, the reparametrization
> below is better), how would I do that? I thought of
>
> lp__ <- negative_infinity()
>
> but lp__ is deprecated, so I am not sure about this. Is there a function
> that sets lp__ to -infinity, _and_ possibly exits early from the model
> posterior calculation? If not, will just open an issue.

increment_log_prob(negative_infinity());

will do it. We're also going to add a raise_exception(...) [perhaps to be
renamed] statement in 2.5 that'll short-circuit to the same effect.

> Some clarifying questions on the approach you suggested, below:
>

>> It's better to actually implement the constraints. The problem is
>> that Stan doesn't let you separately implement constraints for
>> each variable.
>>
>> So you have to enforce this constraint yourself:
>>
>>> z[i]/w[i] > m
>>
>> Rewriting, that's
>>
>> z[i] > m / w[i]
>

> For the archives, that is z[i] > m * w[i], but I get the point.

Thanks for fixing my algebra!

>
>> So you need a lower-bound on z[i] equal to m / w[i].
>>
>> The way to do this is just do what Stan does internally for
>> transformed parameters. Start with a raw, unconstrained parameter,
>> then transform it and apply the Jacobian adjustment:
>>
>> parameters {
>> vector[N] z_raw;
>> }
>

> Just to clarify: I should use an unconstrained variable (or vector)
> because Stan will automatically assign a flat (\propto 1) prior to it,
> so it will not mess up my likelihood.

Correct. Stan assigns flat priors over the valid values
of each variable.

> And if I wanted to use something
> like
>
> parameters {
> real<lower=0> z_raw[N];
> }

> transformed parameters {
> vector[N] z;
> for (i in 1:N)

> z[i] <- z_raw + m / w[i];
> // no need to increment lp__ with Jacobian correction, transformation linear
> }

> model {
> z ~ ...; // whatever prior you want here
> }
>

> then (maybe some future version of) Stan may use a prior that would add
> transformed terms to my posterior? Or would this be innocuous?

After fixing the algebra, then you can use whatever prior you want on z.

Or you could put a prior on z_raw, and then you don't need the Jacobian
adjustment. See below.

Is w[i] constrained to be positive? It's not necessary
for the model here. If it is, then you can just
define the parameter with a constraint lower=0 and then you
won't need another term. Same for m if that's constrained
to be positive.

> model {
> increment_log_prob(m_aux);
> for (i in 1:N) {
> z[i] ~ lognormal(0.0, 1.0);
> w[i] ~ lognormal(0.0, 1.0);
> z_hat[i]/z[i] ~ lognormal(0, 0.1);
> w_hat[i]/w[i] ~ lognormal(0, 0.1);
> increment_log_prob(z_aux[i]);
> increment_log_prob(w_aux[i]);
> }
> }

You can't use

z_hat[i] / z[i] ~ ...

without applying another Jacobian, this time for the u -> const/u
transform.

Every time you transform a variable on the left side of a sampling
statement, you need to add the log Jacobian determinant.

> --8<---------------cut here---------------end--------------->8---
>
> I run and check it with
>
> --8<---------------cut here---------------start------------->8---
> model <- stan_model(file = "bivariate_censoring1.stan")
> fit <- sampling(model,
> data =
> list(N = nrow(d),
> z_hat = d$z_hat,
> w_hat = d$w_hat),
> iter=1000,
> chains=1)
> posterior <- extract(fit, permuted = TRUE)
> plot(d$z, colMeans(posterior$z))
> abline(a=0, b=1, col="gray50")
> plot(d$w, colMeans(posterior$w))
> abline(a=0, b=1, col="gray50")
> hist(posterior$m)
> abline(v=1, col="gray50")
> --8<---------------cut here---------------end--------------->8---
>
> and everything checks out except that the posterior mean of m is well
> below 1, and also, it depends on the sample size. I know that censoring
> boundaries are tricky to identify well, but I am still wondering if I
> miscoded something.

I'm not sure what you were trying to code, but I think you missed
a Jacobian.

>
> Also, I get a couple of
> --8<---------------cut here---------------start------------->8---
> Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
> Error in function stan::prob::lognormal_log(N4stan5agrad3varE): Random variable is -nan:0, but must not be nan!
> If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
> but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
> --8<---------------cut here---------------end--------------->8---
> messages but not more than 1-2 per chain, so I think these can be
> ignored safely.

Correct --- this is usually just an early adaptation problem until
the variable scales and step sizes are adjusted properly. Even then,
you can run over boundaries in some cases, but the rejection will do
the appropriate Metropolis adjustment.

- Bob

[Tamas Papp]

Hi Bob and others,

having learned from the univariate example, I fixed this one too. I
missed not the Jacobian, but the truncation adjustment. A complete
example with simulated data is below. The transformation (z->x) is
needed because AFAIK Stan does not allow expressions after lower= to be
a vector (correct me if I am wrong, also, this is not a complaint, I am
fine with this), and using

real<lower=0> x[N];

and then a truncated distribution (from y[i]*m) instead of what I have
below does not play nicely with the sampler and leads to bad mixing (lot
of rejections). I added comments for people unfamiliar with Jacobians
and transformations, so I am posting this in the hope that some other
newbies to Stan will may it useful.

Best,

Tamas

--8<---------------cut here---------------start------------->8---
## A simple study in implementing univariate censoring
library(rstan)

## simulated data
m <- 1 # boundary
x <- rlnorm(100) # latent variable
y <- rlnorm(length(x)) # latent variable
keep_p <- x/y > m # only keep when ratio > m
x <- x[keep_p]
y <- y[keep_p]
N <- length(keep_p)
x_hat <- rlnorm(N, log(x), 0.1) # observed with noise
y_hat <- rlnorm(N, log(y), 0.1) # observed with noise

## transformation: z[i] = log(x[i]-y[i]*m)
## inverse transformation x[i] = exp(z[i])+y[i]*m
## if we had N=1: f: (m, y, z) -> (m, y, x)
## the Jacobian would look like:
## /dm /dy /dz
## dm 1 0 0
## dy 0 1 0
## dx y m exp(z)
##
## triangular, so determinant is exp(z). this generalizes to N > 0.

model <- stan_model(model_code =
"
data {
int<lower=0> N;
real<lower=0> x_hat[N];
real<lower=0> y_hat[N];
}
parameters {
real<lower=0> m;
real<lower=0> y[N];
real z[N];
}
model {
vector[N] x;
vector[N] log_x;
vector[N] log_y;
for (i in 1:N) {
x[i] <- exp(z[i]) + m*y[i];
log_x[i] <- log(x[i]);
log_y[i] <- log(y[i]);
}
y ~ lognormal(0, 1);
x ~ lognormal(0, 1);
increment_log_prob(N*log(m)); # truncation adjustment
increment_log_prob(z); # Jacobian adjustment
x_hat ~ lognormal(log_x, 0.1);
y_hat ~ lognormal(log_y, 0.1);
m ~ uniform(0, 2);
}
")


## fitting the model
fit <- sampling(model,
data = list(N = N, x_hat = x_hat, y_hat = y_hat),
init = function(...) list(m = 1, y = rep(2,N), z = rep(0, N)),
iter=1000,
chains=3)

## very nice chains for m
plot(extract(fit, pars = "m", permuted = FALSE, inc_warmup = TRUE ),
type="l", ylab="m (chain)")
hist(extract(fit, pars = "m", permuted = TRUE)$m)
--8<---------------cut here---------------end--------------->8---

[Bob Carpenter]

On Aug 22, 2014, at 11:32 AM, Tamas Papp <tkp...@gmail.com> wrote:

> Hi Bob and others,
>
> having learned from the univariate example, I fixed this one too. I
> missed not the Jacobian, but the truncation adjustment. A complete
> example with simulated data is below. The transformation (z->x) is
> needed because AFAIK Stan does not allow expressions after lower= to be
> a vector (correct me if I am wrong, also, this is not a complaint, I am
> fine with this),

So far, it doesn't. I'm going to generalize the constraints so that
they can accept vectors, and hopefully being able to use functions
will make it easy to code the vectors. Don't know how soon I'll
be able to get to it.

> and using
>
> real<lower=0> x[N];
>
> and then a truncated distribution (from y[i]*m) instead of what I have
> below does not play nicely with the sampler and leads to bad mixing (lot
> of rejections).

The problem is that you want the support in the model
to match that of the constraints. And even then, if
there are parameters involved, you'll also need the truncation.

> I added comments for people unfamiliar with Jacobians
> and transformations, so I am posting this in the hope that some other
> newbies to Stan will may it useful.

Thanks!

- Bob