# Integrating out censored data in negative binomial model

**URL:** <https://discourse.mc-stan.org/t/integrating-out-censored-data-in-negative-binomial-model/567>\
**Category:** Modeling\
**Tags:** specification\
**Created:** [May 16, 2017, 4:02pm UTC](https://discourse.mc-stan.org/t/integrating-out-censored-data-in-negative-binomial-model/567 "2017-05-16T16:02:42Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Blaine\_Mooers](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/blaine_mooers/32/174_2.png) [@Blaine\_Mooers](https://discourse.mc-stan.org/u/Blaine_Mooers)\
**Post date:** [May 16, 2017, 4:02pm UTC](https://discourse.mc-stan.org/t/integrating-out-censored-data-in-negative-binomial-model/567/1 "2017-05-16T16:02:43Z")

</div>

[http://discourse.mc-stan.org/t/right-censored-data-in-a-negative-binomial-model/562?source\_topic\_id=566&source\_topic\_id=567](http://discourse.mc-stan.org/t/right-censored-data-in-a-negative-binomial-model/562?source_topic_id=566&source_topic_id=567)

By accident, my initial post did not go to a new topic. Please click on above the link.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [May 18, 2017, 10:21pm UTC](https://discourse.mc-stan.org/t/integrating-out-censored-data-in-negative-binomial-model/567/2 "2017-05-18T22:21:02Z")

</div>

I’m afraid that link’s broken.

---

<div class="post-metadata">

**Author:** ![Blaine\_Mooers](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/blaine_mooers/32/174_2.png) [@Blaine\_Mooers](https://discourse.mc-stan.org/u/Blaine_Mooers)\
**Post date:** [May 19, 2017, 3:18am UTC](https://discourse.mc-stan.org/t/integrating-out-censored-data-in-negative-binomial-model/567/3 "2017-05-19T03:18:39Z")

</div>

Thank you Bob for letting me know that the link is broken. I can still access it, so it may be permission issue.  
I have reformatted my post below.

I suspect that I am missing something really basic. This is my first Stan model. I am trying to apply the two procedures for right censored data in section 11.3 in the manual. I adapted the example for continuous data to my count data and have run into trouble.

The first procedure imputes the censored data while

Contents of SimpleNBrightCensored.stan

```
data{
    int<lower=0> N;
    int<lower=0> N_obs;
    int<lower=0> N_cens;
    int y_obs[N_obs];
    int y_cens[N_cens];
    int<lower=max(y_obs)> U;
}
parameters{
  real mu;
  real<lower=0.001> scale;
}
model{
  scale ~ exponential( 2 );
  y_obs ~ neg_binomial_2( mu, scale );
  y_cens ~ neg_binomial_2( mu, scale );
}
generated quantities{
  real dev;
  dev = 0;
  dev = dev + (-2)*neg_binomial_2_lpmf( y_obs | mu , scale );
}

```

Rcode for Model1

```
model_file = 'simpleNBrightCensored.stan'
iterations = 10000

# inits
N = 54
N_obs = 44
N_cens = 10
U = 601
mu = 256
scale = 1.7

# data 
y_obs=c(7,23,54,59,61,76,79,91,99,127,129,133,141,147,147,150,
151,155,176,184,187,189,207,225,238,239,251,276,280,300,305,332,
369,391,425,432,462,481,483,506,518,544,569,581)
y_cens=c(601,601,601,601,601,601,601,601,601,601)

stan_data = list(N=N, N_obs=N_obs, y_obs=y_obs, N_cens=N_cens, y_cens=y_cens) # data passed to stan
# set up the model
stan_model = stan(model_file, data = stan_data, chains = 1)

# fit model
stanfit2 = stan(fit = stan_model, data = stan_data,
                iter=iterations) # run the model

print(stanfit2,digits=3)

```

The second procedure involves integrating out the censored data.

The model simpleNBrightCensored2.stan:

```
data{
    int<lower=0> N;
    int<lower=0> N_obs;
    int<lower=0> N_cens;
    int y_obs[N_obs];
    int y_cens[N_cens];
    int<lower=max(y_obs)> U;
}
parameters{
  real mu;
  real<lower=0.1> scale;
}
model{
  scale ~ exponential( 2 );
  y_obs ~ neg_binomial_2( mu , scale );
  target += N_cens * neg_binomial_2_lccdf(U | mu, scale);
}
generated quantities{
  real dev;
  dev = 0;
  dev = dev + (-2)*neg_binomial_2_lpmf( y_obs | mu , scale );
}

```

The Rcode to run this model:

```
model_file = 'simpleNBrightCensored2.stan'
iterations = 10000
# inits
N = 54
N_obs = 44
N_cens = 10
#U = 601
#mu = 256
#scale = 1.7
# data 
y_obs=c(7,23,54,59,61,76,79,91,99,127,129,133,141,147,147,150,151,155,176,184,187,189,207,225,238,239,251,276,280,300,305,332,369,391,425,432,462,481,483,506,518,544,569,581)
y_cens=c(601,601,601,601,601,601,601,601,601,601)

stan_data = list(N=N, N_obs=N_obs, y_obs=y_obs, N_cens=N_cens, y_cens=y_cens) # data passed to stan
# set up the model
stan_model2 = stan(model_file, data = stan_data, chains = 1)

# fit model
stanfit3 = stan(fit = stan_model2, data = stan_data,
                iter=iterations) # run the model
print(stanfit3,digits=3)

```

I get the following error message:

```
"Rejecting initial value:
Log probability evaluates to log(0), i.e. negative infinity.
Stan can't start sampling from this initial value."
```

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [May 19, 2017, 1:56pm UTC](https://discourse.mc-stan.org/t/integrating-out-censored-data-in-negative-binomial-model/567/4 "2017-05-19T13:56:57Z")

</div>

When I tried to run the second model I got an error about U missing. I added a value for U in the stan\_data and got the error you’re describing.

I added bounds on mu and scale as specified in the back of the manual and things seem to work:

```
real<lower=0.0> mu;
real<lower=0.0> scale;

```

Hope that helps!

---

<div class="post-metadata">

**Author:** ![Blaine\_Mooers](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/blaine_mooers/32/174_2.png) [@Blaine\_Mooers](https://discourse.mc-stan.org/u/Blaine_Mooers)\
**Post date:** [May 19, 2017, 5:02pm UTC](https://discourse.mc-stan.org/t/integrating-out-censored-data-in-negative-binomial-model/567/5 "2017-05-19T17:02:31Z")

</div>

Thank you very much bbbales2!

In addition to updating the model and uncommenting U, I had to uncomment mu and scale in the R code and it then ran to completion.

The corrected model code is here

```
data{
  int<lower=0> N;
  int<lower=0> N_obs;
  int<lower=0> N_cens;
  int y_obs[N_obs];
  int y_cens[N_cens];
  int<lower=max(y_obs)> U;
  }
  parameters{
    real<lower=0.0> mu;
    real<lower=0.0> scale;
  }
  model{
    scale ~ exponential( 2 );
    y_obs ~ neg_binomial_2( mu , scale );
    target += N_cens * neg_binomial_2_lccdf(U | mu, scale);
  }
  generated quantities{
    real dev;
   dev = 0;
   dev = dev + (-2)*neg_binomial_2_lpmf( y_obs | mu , scale );
}

```

and the corrected R code is here:

```
rm(list=ls())
library(rstan)
model_file = 'simpleNBrightCensored2.stan'
iterations = 10000
# inits
N = 54
N_obs = 44
N_cens = 10
U = 601
mu = 256
scale = 1.7

# data 
y_obs=c(7,23,54,59,61,76,79,91,99,127,129,133,141,147,147,150,151,155,176,184,187,189,207,225,238,239,251,276,280,300,305,332,369,391,425,432,462,481,483,506,518,544,569,581)
y_cens=c(601,601,601,601,601,601,601,601,601,601)

stan_data = list(N=N, N_obs=N_obs, y_obs=y_obs, N_cens=N_cens, y_cens=y_cens) # data passed to stan
# set up the model
stan_model2 = stan(model_file, data = stan_data, chains = 1)

# fit the model
stanfit3 = stan(fit = stan_model2, data = stan_data,
            iter=iterations)
```

---

<div class="post-metadata">

**Author:** ![ABUJARAD](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/abujarad/32/2876_2.png) [@ABUJARAD](https://discourse.mc-stan.org/u/ABUJARAD)\
**Post date:** [August 1, 2018, 10:14am UTC](https://discourse.mc-stan.org/t/integrating-out-censored-data-in-negative-binomial-model/567/6 "2018-08-01T10:14:21Z")

</div>

But how can I get WAIC and AIC for censored data??

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [August 1, 2018, 10:26am UTC](https://discourse.mc-stan.org/t/integrating-out-censored-data-in-negative-binomial-model/567/7 "2018-08-01T10:26:57Z")

</div>

@ABUJARAD Your question deserves it’s own thread I think. People might not notice it here.

---

<div class="post-metadata">

**Author:** ![ABUJARAD](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/abujarad/32/2876_2.png) [@ABUJARAD](https://discourse.mc-stan.org/u/ABUJARAD)\
**Post date:** [August 1, 2018, 10:29am UTC](https://discourse.mc-stan.org/t/integrating-out-censored-data-in-negative-binomial-model/567/8 "2018-08-01T10:29:48Z")

</div>

If possible help solve this problem

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [August 1, 2018, 10:34am UTC](https://discourse.mc-stan.org/t/integrating-out-censored-data-in-negative-binomial-model/567/9 "2018-08-01T10:34:27Z")

</div>

@ABUJARAD I don’t know a good answer or where to point you to for one. I can at least tell you the answer probably isn’t going to be simple, haha. Starting a new thread is your best bet.