# Censored Weibull

**URL:** <https://discourse.pymc.io/t/censored-weibull/14173>\
**Category:** General\
**Created:** [April 3, 2024, 3:30pm UTC](https://discourse.pymc.io/t/censored-weibull/14173 "2024-04-03T15:30:11Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Simon](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/simon/32/5566_2.png) [@Simon](https://discourse.pymc.io/u/Simon)\
**Post date:** [April 3, 2024, 3:30pm UTC](https://discourse.pymc.io/t/censored-weibull/14173/1 "2024-04-03T15:30:11Z")

</div>

Hi all. I’m trying to adapt the approach to estimate incubation periods from this paper: [https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4757028/#pone.0148506.ref002](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4757028/#pone.0148506.ref002)  
R code and data here: [Dryad | Data -- Association between the severity of influenza A(H7N9) virus infections and length of the incubation period](https://datadryad.org/stash/dataset/doi:10.5061/dryad.g2j62) .  
However, I cannot get similar results to those reported in the paper, and I cannot figure out why. Here’s my code:

```auto
# -*- coding: utf-8 -*-
import pandas as pd
import numpy as np
import pymc as pm
import arviz as az
import pytensor.tensor as pt
from scipy.special import gammainc as gammai

np.random.seed(27)

def weibull_cdf(x, kappa, theta):
    return 1 - pt.exp(-(x/theta)**kappa)

def censored(name, kappa, theta, lower, upper):
    L = weibull_cdf(lower, kappa, theta)
    U = weibull_cdf(upper, kappa, theta)
    return U - L

data = pd.read_csv("./original_data/data_h7n9_severity.csv")

lower = data.IncP_min.values #incubation periods lower boundary
upper = data.IncP_max.values #incubation periods upper boundary

fatal = data.death_status.values   

with pm.Model() as mod:
    # kappa = pm.HalfNormal('kappa', 5)
    # theta = pm.HalfNormal('theta', 5)
    kappa = pm.Uniform('kappa', 0, 100, shape=2) #as used in the paper
    theta = pm.Uniform('theta', 0, 100, shape=2) #as used in the paper
    mu = pm.Deterministic('mu', theta*pt.gamma(1 + 1 / kappa)) #incubation period mean
    y = pm.Potential('y', censored('censored', kappa[fatal], theta[fatal], lower, upper))
    
with mod:
    idata = pm.sample(2000, tune=2000, nuts_sampler='numpyro', random_seed=27)

az.summary(idata)

Out[27]: 
            mean sd hdi_3% hdi_97% ... mcse_sd ess_bulk ess_tail r_hat
kappa[0] 51.804 27.256 7.450 96.496 ... 0.643 831.0 917.0 1.01
kappa[1] 55.424 26.378 14.275 99.970 ... 0.574 969.0 829.0 1.00
theta[0] 3.569 1.385 1.252 6.081 ... 0.043 530.0 627.0 1.01
theta[1] 3.825 1.130 1.371 4.991 ... 0.048 423.0 681.0 1.01
f[0] 3.503 1.379 1.180 5.957 ... 0.042 546.0 738.0 1.01
f[1] 3.771 1.131 1.309 4.986 ... 0.046 427.0 721.0 1.01
mu[0] 3.508 1.363 1.260 5.991 ... 0.042 531.0 644.0 1.01
mu[1] 3.769 1.115 1.365 4.936 ... 0.047 420.0 698.0 1.01

[8 rows x 9 columns]

```

Any help will be really appreciated.

---

<div class="post-metadata">

**Author:** ![iavicenna](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/iavicenna/32/7923_2.png) [@iavicenna](https://discourse.pymc.io/u/iavicenna)\
**Post date:** [April 3, 2024, 3:41pm UTC](https://discourse.pymc.io/t/censored-weibull/14173/2 "2024-04-03T15:41:28Z")

</div>

There is a wrapper for creating censored distributions

[https://www.pymc.io/projects/docs/en/stable/api/distributions/censored.html](https://www.pymc.io/projects/docs/en/stable/api/distributions/censored.html)

Also for censoring what happens is that you have a bunch of data say [10, 10, 10 ,180, 180, 1240]  
with an upper an lower boundary (say 10, 1240). You supply this data as the observed input to your censored distribution and 10 and 1240 as lower and upper boundaries so that censored distribution decides how to tackle each data point. Probability of 10 and 1240 is then automatically calculated with cdf where as 180’s are calculated with the pdf.

Since I cant see what lower and upper looks like in your case without the data, I am not sure what you are doing here but it does not seem like it. Can you post an example of what lower and upper looks like?

---

<div class="post-metadata">

**Author:** ![Simon](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/simon/32/5566_2.png) [@Simon](https://discourse.pymc.io/u/Simon)\
**Post date:** [April 4, 2024, 8:24am UTC](https://discourse.pymc.io/t/censored-weibull/14173/3 "2024-04-04T08:24:38Z")

</div>

Hi, thank you. Yes, sorry about that. When reading the paper I felt a similar confusion, as it does not seem like a censored approach ‘proper’. But it is (partially, I guess) a standard in epidemiology to estimate incubation periods in this way. The main issue is that there are not unique lower and upper boundaries. Instead, every observation has its own boundaries. The common case is that each patient has a possible lower and upper boundary of an hypothetical/estimated exposure period, where incubation is time of symptoms onset minus time of exposure. So we have, for instance, minimum exposure moment (minExp) and maximum exposure moment (maxExp) and symptoms onset (onset) So the point is estimating a ‘true’ value (often called ‘exact’) of the incubation period (i.e. that lies between onset - minExp and onset - maxExp).

The approach in the paper is by sampling from the CDFs of Weibull distributions. You can see that in the linked paper and R code. You can see the data and code in the second link of my initial post. But the key bit of their code is this one (I guess):

```auto
# prior function
logprior = function(current){
  current$logprior=dunif(current$k,0,100,log=TRUE)+dunif(current$theta,0,100,log=T)
  current
}

# The likelihood function
loglikelihood = function(current,data){
  current$loglikelihood <- log(pweibull(data$IncP_max,current$k,current$theta)-pweibull(data$IncP_min,current$k,current$theta))
  current
}

```

Which is what I’m trying to translate to PyMC. They sample it with Metropolis, I’d prefer to use HMC.

I have checked the implementation of the PyMC Censored and multiple posts on the topic before posting. But it seems that the data I’m working with is not matching the requirements. PyMC censored requires access to the censored data (e.g. a list of values) and two absolute boundaries (upper and lower). While in the data I have corresponds to two lists of variable boundaries (395 lower boundaries, 395 upper boundaries).

---

<div class="post-metadata">

**Author:** ![iavicenna](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/iavicenna/32/7923_2.png) [@iavicenna](https://discourse.pymc.io/u/iavicenna)\
**Post date:** [April 4, 2024, 9:14am UTC](https://discourse.pymc.io/t/censored-weibull/14173/4 "2024-04-04T09:14:45Z")

</div>

> [@Simon](#):
>
> The main issue is that there are not unique lower and upper boundaries. Instead, every observation has its own boundaries.

> > PyMC censored requires access to the censored data (e.g. a list of values) and two absolute boundaries (upper and lower)

If you want to achieved this, it is easy with the Censored class. Lets say your data is

```auto
[10,10,20,180,180,1200]

```

If you supply

```auto
lower = [10, 10, 20, -np.inf, -np.inf, -np.inf] 
upper = [np.inf, np.inf, np.inf, np.inf, 180, 1200]

```

it will apply lower (left) censoring to 10,10,20 and upper (right) censoring to 180 and 1200.

Regarding the code you posted

```auto
pweibull(data$IncP_max,current$k,current$theta)

```

I assume this is the cumulative distribution function. If so, perhaps you do not really want censoring. Do you just have some measurements that you know falls within an interval (but do not know the exact measurement)? In that case difference between CDFs is indeed the right way to go and I would also do it via Potential terms. I have looked at the data provided by the paper and there are some instances where the lower is equal to upper so for those cases they might have used pdf for weibul (that is what I would have done atleast). So for instance something like this:

```auto
import pandas as pd
import numpy as np
import pymc as pm
import arviz as az
import pytensor.tensor as pt
from scipy.special import gammainc as gammai

np.random.seed(27)

def weibull_cdf(x, kappa, theta):
    return 1 - pt.exp(-(x/theta)**kappa)

def censored(name, kappa, theta, lower, upper):
    L = weibull_cdf(lower, kappa, theta)
    U = weibull_cdf(upper, kappa, theta)
    return U - L

data = pd.read_csv("./data_h7n9_severity.csv")

lower = data.IncP_min.values #incubation periods lower boundary
upper = data.IncP_max.values #incubation periods upper boundary

fatal = data.death_status.values

interval = np.array([[x,y,f] for x,y,f in zip(lower,upper,fatal) if x!=y])
exact = np.array([[x,f] for x,y,f in zip(lower, upper, fatal) if x==y])

with pm.Model() as mod:
    # kappa = pm.HalfNormal('kappa', 5)
    # theta = pm.HalfNormal('theta', 5)
    kappa = pm.Uniform('kappa', 0, 5, shape=2) #as used in the paper
    theta = pm.Uniform('theta', 0, 5, shape=2) #as used in the paper
    mu = pm.Deterministic('mu', theta*pt.gamma(1 + 1 / kappa)) #incubation period mean
    y = pm.Potential('y', censored('censored', kappa[interval[:,2]],
                                   theta[interval[:,2]], interval[:,0],
                                   interval[:,1]))
    pm.Weibull("exact", kappa[exact[:,1]], theta[exact[:,1]],
               observed=exact[:,0])

with mod:
    idata = pm.sample(1000, tune=1000)#, nuts_sampler='numpyro')#, random_seed=27)

az.summary(idata)

az.plot_posterior(idata)

```

from which I get

 ![posterior](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/3/3162d71882d6d73c710d0bb2248033bf3f99247c.png)

---

<div class="post-metadata">

**Author:** ![iavicenna](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/iavicenna/32/7923_2.png) [@iavicenna](https://discourse.pymc.io/u/iavicenna)\
**Post date:** [April 4, 2024, 9:21am UTC](https://discourse.pymc.io/t/censored-weibull/14173/5 "2024-04-04T09:21:29Z")

</div>

There is one thing else missing in your code which might help too. You are not taking log of the difference as far as I can see. You should return

```auto
return pm.math.log(U - L)

```

however for this to work you have to either discard data where U=L or do like I did (which I think is more sensible if I did not misunderstand your data). Nevertheless if you throw away U=L and return log you get

 ![posterior](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/6/6dae66e35c11c00fc7cd63ba83deaca29a8b422a.png)

---

<div class="post-metadata">

**Author:** ![Simon](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/simon/32/5566_2.png) [@Simon](https://discourse.pymc.io/u/Simon)\
**Post date:** [April 4, 2024, 9:28am UTC](https://discourse.pymc.io/t/censored-weibull/14173/6 "2024-04-04T09:28:03Z")

</div>

Thank you very much for the quick and helpful reply. This solved the issue.

Indeed, I had removed that log as in my code it was not working, because I didn’t have the sampling distribution for ‘exact’ values, I couldn’t guess that from their code, but they mention it on the paper a bit implicitly. I should learn to look a bit more carefully to the data as well 😅