# \[Turing.jl\] Censored Values

**URL:** <https://discourse.julialang.org/t/turing-jl-censored-values/48990>\
**Category:** Probabilistic Programming\
**Created:** [October 25, 2020, 6:28pm UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990 "2020-10-25T18:28:33Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![anon92994695](https://avatars.discourse-cdn.com/v4/letter/a/ce7236/32.png) [@anon92994695](https://discourse.julialang.org/u/anon92994695)\
**Post date:** [October 25, 2020, 6:28pm UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/1 "2020-10-25T18:28:33Z")

</div>

Hi there Turing.jl peoples. Still new as ever to Bayes stuff, but having a blast learning(when I can). I ran into a situation where I want to explore censored modelling.

I’ve scoured the documentation in Turing.jl and didn’t see any examples or obvious ways to go about doing it.

PyMC3 Syntax: [https://docs.pymc.io/notebooks/censored\_data.html](https://docs.pymc.io/notebooks/censored_data.html)  
STAN Syntax: [Stan User’s Guide](https://mc-stan.org/docs/2_18/stan-users-guide/truncated-or-censored-data.html)

Is there a preferred way to go about doing this in Turing? I have an idea - but, would rather hear from the pros :).

---

<div class="post-metadata">

**Author:** ![cpfiffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cpfiffer/32/208747_2.png) [@cpfiffer](https://discourse.julialang.org/u/cpfiffer)\
**Post date:** [October 27, 2020, 12:40am UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/2 "2020-10-27T00:40:01Z")

</div>

Hm, I think if you want to draw the censored data from a distribution, you can leave the datapoint as `missing` while sampling. For example, if you had data

```julia
y = [1.0, missing]

```

and a model

```julia
@model function f(y)
    m ~ Normal(0, 1)
    for i in 1:length(y)
        y[i] ~ Normal(m, 1)
    end
end

```

You’ll actually end up imputing the distribution of `y[2]` since it’s a missing value. It’s equivalent to it being a parameter in the model.

Not sure if that answered your question?

---

<div class="post-metadata">

**Author:** ![anon92994695](https://avatars.discourse-cdn.com/v4/letter/a/ce7236/32.png) [@anon92994695](https://discourse.julialang.org/u/anon92994695)\
**Post date:** [October 27, 2020, 1:07am UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/3 "2020-10-27T01:07:58Z")

</div>

First of all - really awesome Turing supports modifying posterior pdf’s via missing values. Had no idea, but that’s very useful and convenient.

so I guess if our data is left or right censored, we could do something like the following?

```julia
@model function f(y, l_censored, r_censored)
    m ~ Normal(0, 1)
    y ~ filldist(Normal(m, 1), length(y))
    for (n,censorand) in enumerate(l_censored)
        missing ~ truncated(Normal(m, 1), censorand, Inf)
    end
    for (n, censorand) in enumerate(r_censored)
        missing ~ truncated(Normal(m, 1), -Inf, censorand)
    end
end

```

(consider this pseudocode I didn’t run it)

Or is there a more preferred way of dealing with this?

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [October 27, 2020, 2:42am UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/4 "2020-10-27T02:42:41Z")

</div>

That’s really cool. Is there specific to handling missing values in a nice way, or does something about the technical approach get you that for free?

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [October 27, 2020, 3:23am UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/5 "2020-10-27T03:23:09Z")

</div>

> [@anon92994695](#):
>
> `missing ~ `

you probably want to have `l_censored[n]` here and similarly for the other one

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [October 27, 2020, 3:24am UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/6 "2020-10-27T03:24:05Z")

</div>

> [@cscherrer](#):
>
> Is there specific to handling missing values in a nice way, or does something about the technical approach get you that for free?

We check for missing values and treat them differently. The check happens at compile-time if the whole variable is missing or at runtime if it’s a `Vector{T, Missing}` for example.

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [October 27, 2020, 3:51am UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/7 "2020-10-27T03:51:46Z")

</div>

Nice. Soss should be able to do this in principle, but I haven’t worked out how I’d want to have the data structures to make imputation efficient.

Censored days problems like @anon92994695 's are [MNAR](https://en.wikipedia.org/wiki/Missing_data?wprov=sfla1). Representing the missingness process for this case is a nice PPL challenge.

---

<div class="post-metadata">

**Author:** ![anon92994695](https://avatars.discourse-cdn.com/v4/letter/a/ce7236/32.png) [@anon92994695](https://discourse.julialang.org/u/anon92994695)\
**Post date:** [October 27, 2020, 7:34pm UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/8 "2020-10-27T19:34:49Z")

</div>

@cscherrer - I agree this is a really good way to challenge a PPL. I’d love to see how Soss evolves to handle this case. It seems niche at first, but, as we’ve chit-chatted it does show up a lot. Having Julia support missing values as a first class citizen makes this very interesting. I REALLY like the design decision made in Turing, but I feel like left/right censoring may benefit from a convenience function? Also not sure how derivative based samplers would work here - my guess is they’d pop?

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [October 28, 2020, 3:12am UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/9 "2020-10-28T03:12:06Z")

</div>

> [@anon92994695](#):
>
> Also not sure how derivative based samplers would work here - my guess is they’d pop?

Guess again 🙂

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [October 28, 2020, 3:34am UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/10 "2020-10-28T03:34:17Z")

</div>

In a lot of applications, a sensor with a limited range will just clamp the result, literally like Julia’s `clamp` function. So instead of `missing`, an overflows might lead to lots of repetitions of whatever the max value is.

I wouldn’t think this would be that hard to implement, and I’m kind of surprised it’s not already in Distributions, especially given `Truncated` is there. That would probably be my approach. Implement it using `Truncated` as a guide, and then you could use the result in any PPL.

---

<div class="post-metadata">

**Author:** ![anon92994695](https://avatars.discourse-cdn.com/v4/letter/a/ce7236/32.png) [@anon92994695](https://discourse.julialang.org/u/anon92994695)\
**Post date:** [October 28, 2020, 7:08pm UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/11 "2020-10-28T19:08:52Z")

</div>

So depending on how you do the censoring you get vastly different solutions

```julia
x_l = rand(truncated(Normal(0,1), 1.0, Inf), 300)

@model function missing_trunc(l_censored)
    m ~ Normal(0, 1)
    for n in 1:length(l_censored)
        missing ~ truncated(Normal(m, 1), l_censored[n], Inf)
    end
end
Model = missing_trunc( x_l )
chain1 = sample(Model, PG(100), 300)

@model function l_censor_trunc(l_censored)
    m ~ Normal(0, 1)
    for n in 1:length(l_censored)
        l_censored[n] ~ truncated(Normal(m, 1), l_censored[n], Inf)
    end
end
Model = l_censor_trunc( x_l )
chain2 = sample(Model, PG(100), 300)

```

I think you’re right in that the correct way to do this is via `l_censored[n] ~`, but the predictions it offers don’t make sense to me. Using `Missing` gives a value I’d expect, but weirdly Turing appears to be treating the term missing as a variable?

Hm. Not so sure if what I’ve suggested is correct or not at this point…

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [October 28, 2020, 7:38pm UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/12 "2020-10-28T19:38:09Z")

</div>

You `x_l` doesn’t have any missing values, so it’s all treated as observations. The only parameter/random variable inferred in this case is `m`. The first model doesn’t make any sense in Turing, it just creates a random variable called `missing` and you are sampling it multiple times which should error in Turing. If it doesn’t, please open an issue. The second model is the correct one.

---

<div class="post-metadata">

**Author:** ![anon92994695](https://avatars.discourse-cdn.com/v4/letter/a/ce7236/32.png) [@anon92994695](https://discourse.julialang.org/u/anon92994695)\
**Post date:** [October 28, 2020, 9:44pm UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/13 "2020-10-28T21:44:41Z")

</div>

Turing doesn’t throw an error so I’ll open an issue - but, it’s okay.

I agree it is correct. Kind of hard to get my head around what a left/right censored value really means I guess. I feel like its incorrect to represent it this way, it almost seems like it needs a hyperprior? But I’ll come to terms with it :).

---

<div class="post-metadata">

**Author:** ![joshualeond](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joshualeond/32/230_2.png) [@joshualeond](https://discourse.julialang.org/u/joshualeond)\
**Post date:** [October 29, 2020, 1:38am UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/14 "2020-10-29T01:38:46Z")

</div>

I played around with this a few months ago and had something working with the `logpdf`, the `logcdf`, and the `logccdf`. In the model function I brought in a variable named ‘censor’ related to what type of censoring (left/right). Something like this:

```julia
    for i = 1:length(years)
        if censor[i] == 0 # observed
            Turing.acclogp!(_varinfo, logpdf(Weibull(beta, eta), years[i]) * weights[i])
        elseif censor[i] == 1 # right censored
            Turing.acclogp!(_varinfo, logccdf(Weibull(beta, eta), years[i]) * weights[i])
        elseif censor[i] == -1 # left censored
            Turing.acclogp!(_varinfo, logcdf(Weibull(beta, eta), years[i]) * weights[i])
        end
    end

```

I’d check this blog post out for something similar with Stan:

[https://ermeel86.github.io/case\_studies/surv\_stan\_example.html](https://ermeel86.github.io/case_studies/surv_stan_example.html)

---

<div class="post-metadata">

**Author:** ![anon92994695](https://avatars.discourse-cdn.com/v4/letter/a/ce7236/32.png) [@anon92994695](https://discourse.julialang.org/u/anon92994695)\
**Post date:** [October 29, 2020, 6:07pm UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/15 "2020-10-29T18:07:59Z")

</div>

That’s the syntax I was looking for. Okay neat- now I can compare the `missing` approach to directly adding to the logp. Once I do this I’ll probably write about the results in a brief post pick a “solution” and we can go on our way.

Can I ask why you opted to use `Turing.acclogp` instead of sampling from the weibull directly? Shouldn’t we only need to use that in cases where things are censored?

---

<div class="post-metadata">

**Author:** ![anon92994695](https://avatars.discourse-cdn.com/v4/letter/a/ce7236/32.png) [@anon92994695](https://discourse.julialang.org/u/anon92994695)\
**Post date:** [October 29, 2020, 10:03pm UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/16 "2020-10-29T22:03:26Z")

</div>

Alright so now I am back to thinking that `missing~` makes sense… @mohamed82008 can you talk me out of it? So the idea is, Turing handles `missing ~` or the index into an array containing `missing`s as a special case? So we can sample a value thats left/right censored, by using a truncated distribution - why not. But, we aren’t going to claim that the posterior distribution should map to the censored value? So Missing does make sense, and it is giving me - at least in this contrived case a very similar distribution for our parameter `m` as using `acclogp` directly.

```nohighlight
x_l = rand(truncated(Normal(0,1), 1.0, Inf), 300)

@model function missing_trunc(l_censored)
    m ~ Normal(0, 1)
    for n in 1:length(l_censored)
        missing ~ truncated(Normal(m, 1), l_censored[n], Inf)
    end
end
Model = missing_trunc( x_l )
chain1 = sample(Model, PG(100), 300)

@model function l_censor_trunc(l_censored)
    m ~ Normal(0, 1)
    for n in 1:length(l_censored)
        l_censored[n] ~ truncated(Normal(m, 1), l_censored[n], Inf)
    end
end
Model = l_censor_trunc( x_l )
chain2 = sample(Model, PG(100), 300)

@model function logpccdf(l_censored)
    m ~ Normal(0, 1)
    for n in 1:length(l_censored)
        Turing.acclogp!(_varinfo, logccdf(Normal(m, 1), l_censored[n] ) )
    end
end
Model = logpccdf( x_l )
chain3 = sample(Model, PG(100), 300)

"missing_distribution" => median( get_params( chain1 ).m )
"censored value" => median( get_params( chain2 ).m )
"logp_hacking" => median( get_params( chain3 ).m )

```

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [October 30, 2020, 6:51pm UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/17 "2020-10-30T18:51:47Z")

</div>

> [@anon92994695](#):
>
> @mohamed82008 can you talk me out of it?

It’s doing the wrong thing! I wrote the DSL code so I know 🙂

> [@anon92994695](#):
>
> the index into an array containing `missing` s as a special case

It’s not very special. A scalar input variable that’s missing will be treated the same way. But using `missing` like that actually defines a variable called `missing`. So think of it as `x` or something. The “wrong” part is that you are defining the variable multiple times in a loop. You should use an array to tell Turing that these are different variables. This array will be partly numbers (observations) and partly `Base.missing` which will be sampled.

---

<div class="post-metadata">

**Author:** ![anon92994695](https://avatars.discourse-cdn.com/v4/letter/a/ce7236/32.png) [@anon92994695](https://discourse.julialang.org/u/anon92994695)\
**Post date:** [October 30, 2020, 6:54pm UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/18 "2020-10-30T18:54:21Z")

</div>

Alright we’re good now - thanks for the low down :).

---

<div class="post-metadata">

**Author:** ![drbenvincent](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drbenvincent/32/3046_2.png) [@drbenvincent](https://discourse.julialang.org/u/drbenvincent)\
**Post date:** [January 23, 2021, 5:41pm UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/19 "2021-01-23T17:41:13Z")

</div>

Hi. I’m working on an example of censored regression. I’m trying to adapt the approach outlined by @joshualeond. I know the basics of what we are doing here… we are not treating censored data as missing, we are not imputing what their values are, instead we are manually incrementing the model log probability. But I am pretty new to Turing, so would appreciate some implementation help if possible?

> [@joshualeond](#):
>
> ```julia
> for i = 1:length(years)
> if censor[i] == 0 # observed
> Turing.acclogp!(_varinfo, logpdf(Weibull(beta, eta), years[i]) * weights[i])
> elseif censor[i] == 1 # right censored
> Turing.acclogp!(_varinfo, logccdf(Weibull(beta, eta), years[i]) * weights[i])
> elseif censor[i] == -1 # left censored
> Turing.acclogp!(_varinfo, logcdf(Weibull(beta, eta), years[i]) * weights[i])
> end
> end
> 
> ```

Here is a MWE of a censored linear regression that I have so far.

```julia
using Turing, StatsPlots

function censor_data(x, y, L, U)
    censor = fill(0, length(x))
    for i in 1:length(x)
        if y[i] > U
            censor[i] = 1
            y[i] = U
        elseif y[i] < L
            censor[i] = -1
            y[i] = L
        else
            censor[i] = 0
        end
    end
    return x, y, censor
end

@model function censored_regression(x, y, censor)
    m ~ Normal(0, 2)
    c ~ Normal(0, 2)
    σ ~ truncated(Normal(0, 20), 0, Inf)
    for i in 1:length(x)
        if censor[i] == 0 # observed
            Turing.acclogp!(_varinfo, logpdf(Normal(m*x[i]+c, σ), y[i]))
        elseif censor[i] == 1 # right/upper censored
            Turing.acclogp!(_varinfo, logccdf(Normal(m*x[i]+c, σ), y[i]))
        elseif censor[i] == -1 # left/lower censored
            Turing.acclogp!(_varinfo, logcdf(Normal(m*x[i]+c, σ), y[i]))
        end
    end
end

# generate data
m, c, σ, N = 1, 0, 2, 200
x = rand(Uniform(-10, 10), N)
y = rand.(Normal.(m.*x.+c, σ))
# Apply censoring
L, U = -5, 5
x, y, censor = censor_data(x, y, L, U)

chain = sample(censored_regression(x, y, censor), NUTS(), 2000)
plot(chain)

```

I think I’m definitely learning here, as the model does actually run! The chains look good.

 ![Screenshot 2021-01-23 at 17.33.08](https://global.discourse-cdn.com/julialang/original/3X/5/e/5ecd521cdf265338ab59ba092ad6c5e0fad00a09.png)

However, the parameter estimates are biased (ie `m` is quite far below it’s true value of 1). So I don’t think the model is doing what it is intended to do. It is as if it is just doing regular regression, so I expect this issue is with the model definition. I would appreciate any insights here.

```julia
function visualise(x, y, chain, subplot, true_params)
    m, c = true_params
    scatter!(x, y, color=:grey, markerstrokewidth=0, label=:none, subplot=subplot)
    xi = [minimum(x), maximum(x)]
    n_samples = length(chain[:m])
    samples_to_draw = 100
    samples_to_plot = sample([1:n_samples;], samples_to_draw; replace=false)
    for i in samples_to_plot
        plot!(xi, chain[:m][i].*xi.+chain[:c][i], color=:black, alpha=0.1, label=false, subplot=subplot)
    end
    plot!(xi, mean(chain[:m]).*xi.+mean(chain[:c]), color=:red, lw=2, label="posterior mean", ylim=[minimum(y)*1.5, maximum(y)*1.5], subplot=subplot)
    plot!(xi, m.*xi.+c, color=:black, lw=3, label="true", ylim=[minimum(y)*1.5, maximum(y)*1.5], subplot=subplot)
    plot!(tick_direction=:out)
end

visualise(x, y, chain, 1, (m, c))

```

 ![Screenshot 2021-01-23 at 17.32.54](https://global.discourse-cdn.com/julialang/original/3X/d/4/d4b68ddd91faff5102f750b14bf17e5fcd7d7e12.png)

---

<div class="post-metadata">

**Author:** ![drbenvincent](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drbenvincent/32/3046_2.png) [@drbenvincent](https://discourse.julialang.org/u/drbenvincent)\
**Post date:** [January 23, 2021, 5:46pm UTC](https://discourse.julialang.org/t/turing-jl-censored-values/48990/20 "2021-01-23T17:46:29Z")

</div>

Ha! My bad… it does actually work! I think it was just a bad run / dataset. Repeatedly running with different generated datasets seems to work fine 🙂

 ![Screenshot 2021-01-23 at 17.45.46](https://global.discourse-cdn.com/julialang/original/3X/7/c/7c6d9681e6131d2b3280a4e8dfb74d02cdb3b492.png)

[Next page](https://discourse.julialang.org/t/turing-jl-censored-values/48990.md?page=2)