# Turing: censored regression example

**URL:** <https://discourse.julialang.org/t/turing-censored-regression-example/53837>\
**Category:** Probabilistic Programming\
**Created:** [January 23, 2021, 6:22pm UTC](https://discourse.julialang.org/t/turing-censored-regression-example/53837 "2021-01-23T18:22:42Z")\
**Posts on this page:** 10\
**Page:** 1

<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, 6:22pm UTC](https://discourse.julialang.org/t/turing-censored-regression-example/53837/1 "2021-01-23T18:22:42Z")

</div>

Building upon my previous example of [truncated regression](https://discourse.julialang.org/t/turing-truncated-regression-example/53832), I thought I’d do a worked example of censored regression. To clarify, censoring is where observations below a lower bound and/or observations above an upper bound are recorded at being at those bounds. An example would be weighting something heavier than the maximum weight that a measuring device can handle, and recording it as the max weight on the dial.

If we think about that in the context of a simple linear regression model, we could be concrete about censoring, defining it in a function:

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

```

The vector `censor` has values `0` for uncensored data, `-1` for left/lower censored data, and `1` for right/upper censored data. This is necessary for the censored regression model that we’ll construct below.

We can then generate some fake data, censor it, and visualise the process of censoring

```julia
# 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
xc, yc, censor = censor_data(x, y, L, U)

```

 ![Screenshot 2021-01-23 at 18.10.18](https://global.discourse-cdn.com/julialang/original/3X/a/9/a926c7232a48606fc235da68b5da56aa1e070522.png)

Our censored regression model can be defined as below. I borrowed from the suggestion by @joshualeond in [this post](https://discourse.julialang.org/t/turing-jl-censored-values/48990/14).

```julia
using Turing

@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.@addlogprob! logpdf(Normal(m*x[i]+c, σ), y[i])
        elseif censor[i] == 1 # right/upper censored
            Turing.@addlogprob! logccdf(Normal(m*x[i]+c, σ), y[i])
        elseif censor[i] == -1 # left/lower censored
            Turing.@addlogprob! logcdf(Normal(m*x[i]+c, σ), y[i])
        end
    end
end

```

We can then sample and check the chains

```julia
chain = sample(censored_regression(xc, yc, censor), NUTS(), 2000)
plot(chain)

```

 ![Screenshot 2021-01-23 at 18.12.00](https://global.discourse-cdn.com/julialang/original/3X/f/3/f3a417ff4114cb8f7a76588fbe7ff85a8c061a52.png)  
Everything looks good. The posterior for m is pretty close, and the posterior for c is slightly low, but nothing overly problematic.

Visualising the posterior predictive confirms that the model is a good fit to the data. By using the censored regression model, we have avoided the bias in estimating the slope that would have occurred with regular regression alone.

 ![Screenshot 2021-01-23 at 18.15.45](https://global.discourse-cdn.com/julialang/original/3X/5/7/57bf3162747f4a88cb7dc9796bb6cdd7d0b045ce.png)

---

<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:** [January 23, 2021, 6:31pm UTC](https://discourse.julialang.org/t/turing-censored-regression-example/53837/2 "2021-01-23T18:31:05Z")

</div>

Very nice!

It’s interesting, I usually think of PPLs as a way to combine existing distributions (or meaures). So I’d probably approach this by first defining a new `Censored` measure, and then using that. It’s especially fun that `Censored` would be singular wrt Lebesgue measure, so we’d need to represent it on top of as a superposition of Lebesgue and Dirac measures.

I’ll need to see about getting this set up in MeasureTheory.jl, then we can compare that with your Turing implementation.

---

<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, 6:35pm UTC](https://discourse.julialang.org/t/turing-censored-regression-example/53837/3 "2021-01-23T18:35:10Z")

</div>

Thanks @cscherrer. I’m less of a mathematician/statistician than you, so I might not be following you exactly. But are you suggesting to replicate `truncated()` into a new `censored` function to work with distributions? From the user / model builder perspective, being able to write something like

```julia
y[i] ~ censored(Normal(m*x[i]+c, σ, L, U)

```

would be pretty amazing.

---

<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:** [January 23, 2021, 7:04pm UTC](https://discourse.julialang.org/t/turing-censored-regression-example/53837/4 "2021-01-23T19:04:38Z")

</div>

Sorry if I was a little cryptic, but it sounds like you got the idea.

On the “singular wrt Lebesgue measure” thing… Say you took billions of samples from a censored distribution, and tried to build histograms. The samples at the limits of the distribution would tend to overwhelm the rest. The distribution is concentrated there, so in a sense the units don’t work out properly, and it’s a mess.

But this is the kind of use case we’re targeting for MeasureTheory, so it makes for a nice test.

---

<div class="post-metadata">

**Author:** ![grahamstark](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/grahamstark/32/5200_2.png) [@grahamstark](https://discourse.julialang.org/u/grahamstark)\
**Post date:** [January 23, 2021, 8:05pm UTC](https://discourse.julialang.org/t/turing-censored-regression-example/53837/5 "2021-01-23T20:05:06Z")

</div>

The problem you’re describing is modelled with [Tobit Regression](https://web.ics.purdue.edu/~jltobias/674/lectures/tobit.pdf) in my world. That would be an interesting thing to add to one of the econometrics packages.

---

<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, 8:09pm UTC](https://discourse.julialang.org/t/turing-censored-regression-example/53837/6 "2021-01-23T20:09:52Z")

</div>

Yep, I believe (but don’t know for sure) that this is doing exactly the same as Tobit regression, but with priors over `m`, `c`, and `σ`.

The lines

```julia
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

```

modify the likelihood in the same way as the Tobit model, I believe. But happy to be corrected on that.

---

<div class="post-metadata">

**Author:** ![devmotion](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devmotion/32/9201_2.png) [@devmotion](https://discourse.julialang.org/u/devmotion)\
**Post date:** [January 24, 2021, 10:15am UTC](https://discourse.julialang.org/t/turing-censored-regression-example/53837/7 "2021-01-24T10:15:55Z")

</div>

You can use `Turing.@addlogprob!` instead of `Turing.acclogp!`:

```julia
Turing.@addlogprob! logpdf(Normal(m*x[i]+c, σ), y[i])

```

The macro just rewrites it with `Turing.acclogp!` but you don’t have to specify (and hence are not warned about) the internal variable `_varinfo`.

---

<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 24, 2021, 11:30am UTC](https://discourse.julialang.org/t/turing-censored-regression-example/53837/8 "2021-01-24T11:30:00Z")

</div>

> [@devmotion](#):
>
> `Turing.@addlogprob! logpdf(Normal(m*x[i]+c, σ), y[i])`

That’s pretty nice. I’ll edit the post to use that approach.

---

<div class="post-metadata">

**Author:** ![EvoArt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evoart/32/25357_2.png) [@EvoArt](https://discourse.julialang.org/u/EvoArt)\
**Post date:** [April 25, 2022, 2:14pm UTC](https://discourse.julialang.org/t/turing-censored-regression-example/53837/9 "2022-04-25T14:14:20Z")

</div>

Distributions.jl now supports censored distributions, which you can use in the likelihood of your Turing model [https://juliastats.org/Distributions.jl/stable/censored/](https://juliastats.org/Distributions.jl/stable/censored/)

---

<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:** [April 30, 2022, 3:07pm UTC](https://discourse.julialang.org/t/turing-censored-regression-example/53837/10 "2022-04-30T15:07:28Z")

</div>

Very nice

```julia
@model function censored_regression(x, y, L, U)
    m ~ Normal(0, 2)
    c ~ Normal(0, 2)
    σ ~ truncated(Normal(0, 20), lower=0)
    y .~ censored.(Normal.(m .* x .+ c, σ), lower=L, upper=U)
end

```