# Question about Bayesian Survival Analysis with Cox (proportional hazards) regression method in Julia with Turing

**URL:** https://discourse.julialang.org/t/question-about-bayesian-survival-analysis-with-cox-proportional-hazards-regression-method-in-julia-with-turing/80672
**Category:** Biology, Health, and Medicine
**Tags:** question
**Created:** [May 7, 2022, 2:46pm UTC](https://discourse.julialang.org/t/question-about-bayesian-survival-analysis-with-cox-proportional-hazards-regression-method-in-julia-with-turing/80672 "2022-05-07T14:46:27Z")
**Posts on this page:** 7
**Page:** 1

<div class="post-metadata">

### Author: ![RyanKuo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryankuo/32/36050_2.png) [@RyanKuo](https://discourse.julialang.org/u/RyanKuo)
#### Post date: [May 7, 2022, 2:46pm UTC](https://discourse.julialang.org/t/question-about-bayesian-survival-analysis-with-cox-proportional-hazards-regression-method-in-julia-with-turing/80672/1 "2022-05-07T14:46:27Z")

</div>

Dear All,  
I am trying to implement Bayesian Survival Analysis with Cox (proportional hazards) regression method in Julia based on the PyMC3 code written by Austin Rochford [https://austinrochford.com/posts/2015-10-05-bayes-survival.html](https://austinrochford.com/posts/2015-10-05-bayes-survival.html).  
Here are my codes:

```julia
using RDatasets,DataFrames,Plots,StatsPlots,Turing,Distributions,FillArrays
HSAUR = dataset("HSAUR","mastectomy");
HSAUR[!,:CMetastized]=levelcode.(HSAUR[!,:Metastized]).-1;

interval_length = 3
interval_bounds =0:interval_length:maximum(HSAUR[!,:Time])+interval_length+1
n_intervals =length(interval_bounds)-1
n_patients=nrow(HSAUR)
patients=1:n_intervals
death =zeros(Int64, n_patients, n_intervals)
exposure =zeros(Int64, n_patients, n_intervals)
last_period =floor.(Int64,HSAUR[!,:Time] ./ interval_length);

```

I don’t know if Julia has a similar function like greater\_equal.outer in NumPy, so I just simply write the following code to covert:

```julia
for i in 1:n_patients
    death[i,last_period[i]]=HSAUR[i,:Event]
    for j in 1:n_intervals
       if HSAUR[i,:Time] >= interval_bounds[j]
          exposure[i,j]=1*interval_length
        else
          exposure[i,j]=0
        end   
    end
    exposure[i,last_period[i]]=HSAUR[i,:Time]-interval_bounds[last_period[i]]
end

```

Turing model and sampling code:

```julia
@model function cox_regression(death,exposure,x)
    λ0 ~Gamma(0.01, 0.01)
   # σ ~Uniform(0.0, 10.0)
   # τ = σ^2
   # μβ ~ Normal(0.0, 100)
   # β ~ Normal(μβ,τ)
     β ~ Normal(0.0,100)
    λ = exp.(β .* x).* fill(λ0,n_intervals)'
    μ = exposure .*λ    
    for i in 1:length(death)
      death[i] ~ Poisson(μ[i])
    end
end

chain = sample(cox_regression(death,exposure,HSAUR[!,:CMetastized]), NUTS(0.65),1000)

```

I am new to Julia and Turing, I am not sure what is wrong that causes very long sampling time (seems never stop)… Could experts here advise me to solve the issue?

Thanks,

Ryan

---

<div class="post-metadata">

### Author: ![kevbonham](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kevbonham/32/216165_2.png) [@kevbonham](https://discourse.julialang.org/u/kevbonham)
#### Post date: [May 10, 2022, 1:07am UTC](https://discourse.julialang.org/t/question-about-bayesian-survival-analysis-with-cox-proportional-hazards-regression-method-in-julia-with-turing/80672/2 "2022-05-10T01:07:33Z")

</div>

Sorry it’s taken a while to get a response, I’m guessing it’s because the overlap of folks in biology and probabilistic programming in this community is not (yet!) Very large - I’m going to post the link in some channels on slack and zulip to see if there’s anyone that can help. Hang in there!

---

<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: [May 10, 2022, 6:29am UTC](https://discourse.julialang.org/t/question-about-bayesian-survival-analysis-with-cox-proportional-hazards-regression-method-in-julia-with-turing/80672/3 "2022-05-10T06:29:38Z")

</div>

> [@RyanKuo](#):
>
> ```julia
> @model function cox_regression(death,exposure,x)
> λ0 ~Gamma(0.01, 0.01)
> # σ ~Uniform(0.0, 10.0)
> # τ = σ^2
> # μβ ~ Normal(0.0, 100)
> # β ~ Normal(μβ,τ)
> β ~ Normal(0.0,100)
> λ = exp.(β .* x).* fill(λ0,n_intervals)'
> μ = exposure .*λ    
> for i in 1:length(death)
> death[i] ~ Poisson(μ[i])
> end
> end
> 
> chain = sample(cox_regression(death,exposure,HSAUR[!,:CMetastized]), NUTS(0.65),1000)
> 
> ```

Hi Ryan

Is this the exact code you tried to run?

I get an error

```julia
ERROR: ArgumentError: column name :CMetastized not found in the data frame; existing most similar names are: :Metastized

```

and then after correcting that

```julia
MethodError: no method matching *(::Float64, ::CategoricalArrays.CategoricalValue{String, UInt8})

```

These are all quick little fixes, but it would be much easier if you have a working MWE. one quick thing I noticed is that you use global variables in the model like `n_intervals`. In Julia you will tend to get much better performance if you pass those variables into the actual function.

---

<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: [May 10, 2022, 6:34am UTC](https://discourse.julialang.org/t/question-about-bayesian-survival-analysis-with-cox-proportional-hazards-regression-method-in-julia-with-turing/80672/4 "2022-05-10T06:34:43Z")

</div>

in terms of the second error, try

```julia
x = HSAUR[!,:Metastized] .== "yes"

```

This will creat a Bool vector wher yes = 1 and no = 0  
then pass in `x` instead of `HSAUR[!,:Metastized]`

---

<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: [May 10, 2022, 9:28am UTC](https://discourse.julialang.org/t/question-about-bayesian-survival-analysis-with-cox-proportional-hazards-regression-method-in-julia-with-turing/80672/5 "2022-05-10T09:28:26Z")

</div>

this atleast runs

```julia
using RDatasets,DataFrames,Plots,StatsPlots,Turing,Distributions,FillArrays
HSAUR = dataset("HSAUR","mastectomy");
HSAUR[!,:CMetastized]=levelcode.(HSAUR[!,:Metastized]).-1;

interval_length = 3
interval_bounds =0:interval_length:maximum(HSAUR[!,:Time])+interval_length+1
n_intervals =length(interval_bounds)-1
n_patients=nrow(HSAUR)
patients=1:n_intervals
death =zeros(Int64, n_patients, n_intervals)
exposure =zeros(Int64, n_patients, n_intervals)
last_period =floor.(Int64,HSAUR[!,:Time] ./ interval_length);

for i in 1:n_patients
    death[i,last_period[i]]=HSAUR[i,:Event]
    for j in 1:n_intervals
       if HSAUR[i,:Time] >= interval_bounds[j]
          exposure[i,j]=1*interval_length
        else
          exposure[i,j]=0
        end   
    end
    exposure[i,last_period[i]]=HSAUR[i,:Time]-interval_bounds[last_period[i]]
end

x = HSAUR[!,:Metastized] .== "yes"
@model function cox_regression(death,exposure,x, n_intervals,nm = size(death))
    n, m = nm
    λ0 ~Gamma(0.01, 0.01)
     β ~ Normal(0.0,100)
    λ = exp.(β .* x) .* fill(λ0,n_intervals)'
    μ = exposure .*λ    
    for j in 1:m
        for i in 1:n
            death[i,j] ~ Poisson(μ[i,j])
        end
    end
end

chain = sample(cox_regression(death,exposure .+0.1,x,n_intervals), NUTS(0.65),1000)

```

I had to add 0.1 to your `exposure` array, (which I’m sure is not what you want, but I did it just to check that the model will run with valid inputs). The issue seems to be that there are instances with 0 exposure, which means 0 μ as well. and `logpdf(Poisson(0),any_number_other_than_zero) = -Inf` i.e. given your inputs, the posterior probability of any parameter values at all will always be zero and you’ll get flooded with warnings.

```julia
Warning: The current proposal will be rejected due to numerical error(s).
│ isfinite.((θ, r, ℓπ, ℓκ)) = (true, true, false, true)

```

I’d check that the calculations for exposure are equivalent to the python code. I’ve never used PyMC3, but essentially this is not a bug, you should not be able to sample from a model where the likelihood encounters `y ~ Poisson(0)` where y ≠ 0 on each iteration.

---

<div class="post-metadata">

### Author: ![RyanKuo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryankuo/32/36050_2.png) [@RyanKuo](https://discourse.julialang.org/u/RyanKuo)
#### Post date: [May 10, 2022, 1:00pm UTC](https://discourse.julialang.org/t/question-about-bayesian-survival-analysis-with-cox-proportional-hazards-regression-method-in-julia-with-turing/80672/6 "2022-05-10T13:00:13Z")

</div>

Hi Kevbonham,

Thanks for your help.

Best Regards,

Ryan

---

<div class="post-metadata">

### Author: ![RyanKuo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryankuo/32/36050_2.png) [@RyanKuo](https://discourse.julialang.org/u/RyanKuo)
#### Post date: [May 10, 2022, 1:17pm UTC](https://discourse.julialang.org/t/question-about-bayesian-survival-analysis-with-cox-proportional-hazards-regression-method-in-julia-with-turing/80672/7 "2022-05-10T13:17:11Z")

</div>

Hi EvoArt:

Thank you very much for looking into it. I did have the code to categorized CMetastized:  
HSAUR[!,:CMetastized]=levelcode.(HSAUR[!,:Metastized]).-1;

I forgot to post that I used CategoricalArrays.jl . But your method works easier and faster.  
And thanks again for your explanation. I agree with you, since the exposure has many 0, that sampler would be struggled with… the formula is: μ = exposure_λ = exposure_λ0 _exp.(β ._ x) that’s why posterior is Poisson. shift by 0.1 should not affect the result. I don’t know how PyMC3 samples such a scenario and R also has a package, I will see if I can find what posterior other people used from paper. I tried to run the PyMC3 code, but it has many package dependency issues.

Best Regards,

Ryan