Skip to content

Instantly share code, notes, and snippets.

@infotroph
Created November 17, 2015 20:07
Show Gist options
  • Save infotroph/8632db0e8d485bb20e2b to your computer and use it in GitHub Desktop.
Save infotroph/8632db0e8d485bb20e2b to your computer and use it in GitHub Desktop.
Is "soft" censoring a thing?
# Testing a soft-left-censored Stan model:
# Values above some minimum are detected normally,
# values less than minimum go detected with some probability > 0.
set.seed(2345767)
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = 7)
sim_mu = 3
sim_sigma = 1.5
detection_limit = 4
# Simulated data: lognormal, 1/3 of values below detection limit are missing
lcens_log = rlnorm(1000, sim_mu, sim_sigma)
lcens_ptrunc_log = lcens_log
lcens_ptrunc_log[lcens_ptrunc_log < detection_limit] = sapply(
lcens_ptrunc_log[lcens_ptrunc_log < detection_limit],
function(x)sample(c(x,x,NA),size=1))
print(summary(lcens_log))
print(summary(lcens_ptrunc_log))
stanstr = "
data{
real<lower=0> L;
int<lower=0> N_obs;
int<lower=0> N_cens;
real<lower=0> y_obs[N_obs];
}
parameters{
real mu;
real<lower=0> sigma;
}
model{
sigma ~ normal(0,100); // Weak prior on sigma to keep sampler away from Inf
y_obs ~ lognormal(mu, sigma);
increment_log_prob(N_cens * lognormal_cdf_log(L, mu, sigma));
}"
fit = stan(
model_code=stanstr,
data=list(
L = detection_limit,
N_obs=length(which(!is.na(lcens_ptrunc_log))),
N_cens=length(which(is.na(lcens_ptrunc_log))),
y_obs=lcens_ptrunc_log[!is.na(lcens_ptrunc_log)]),
chains=7,
iter=10000,
warmup=1000)
print(fit)
plot(fit, pars=c("mu", "sigma"))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment