comparison of automatic and explicit modeling of censoring

truncation_misfit.r

## A simple study in implementing univariate censoring
library(rstan)
 
## simulated data
m <- 1 # boundary
z <- rlnorm(100) # latent variable
z <- z[z > m] # only kept when above m
N <- length(z)
z_hat <- rlnorm(N, log(z), 0.1) # observed with noise

## simple way of modelling, relying on Stan to do the transformation
## and adjust the likelihood for censoring
model1 <- stan_model(model_code =
"
data {
  int<lower=0> N;
  real<lower=0> z_hat[N];
}
parameters {
  real<lower=0> m;
  real<lower=m> z[N];
}
model {
  vector[N] log_z;
  for (i in 1:N) {
    log_z[i] <- log(z[i]);
    z[i] ~ lognormal(0,1) T[m, ];
  }
  z_hat ~ lognormal(log_z, 0.1);
  m ~ uniform(0, 2);
}
")

## explicit transformation and adjustment for censoring
model2 <- stan_model(model_code =
"
data {
  int<lower=0> N;
  real<lower=0> z_hat[N];
}
parameters {
  real<lower=0> m;
  real z_aux[N];
}
model {
  vector[N] z;
  vector[N] log_z;
  for (i in 1:N) {
    z[i] <- exp(z_aux[i]) + m;
    increment_log_prob(z_aux[i]);
    log_z[i] <- log(z[i]);
  }
  increment_log_prob(-lognormal_ccdf_log(m, 0, 1) * N);
  z ~ lognormal(0, 1);
  z_hat ~ lognormal(log_z, 0.1);
  m ~ uniform(0, 2);
}
")
   
## fitting the model
fit1 <- sampling(model1,
                 data = list(N = N, z_hat = z_hat),
                 init = list(list(m = 1.6, z = rep(5,N))),
                 iter=10000,
                 chains=1)

fit2 <- sampling(model2,
                 data = list(N = N, z_hat = z_hat),
                 init = list(list(m = 1.6, z_aux = exp(rep(5,N)))),
                 iter=10000,
                 chains=1)


## comparing chains
## png("/tmp/m-draws.png")
layout(1:2)
plot(extract(fit1, pars = "m", permuted = FALSE, inc_warmup = TRUE ), type="l",
     ylab="m (truncated distribution)")
plot(extract(fit2, pars = "m", permuted = FALSE, inc_warmup = TRUE ), type="l",
     ylab="m (manual transformation)")
## dev.off()