# How to make this AFT regression stable?

**URL:** <https://discourse.julialang.org/t/how-to-make-this-aft-regression-stable/53685>\
**Category:** Statistics\
**Tags:** question, optimization, forwarddiff, distributions\
**Created:** [January 20, 2021, 8:37pm UTC](https://discourse.julialang.org/t/how-to-make-this-aft-regression-stable/53685 "2021-01-20T20:37:44Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![mkarikom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkarikom/32/7096_2.png) [@mkarikom](https://discourse.julialang.org/u/mkarikom)\
**Post date:** [January 20, 2021, 8:37pm UTC](https://discourse.julialang.org/t/how-to-make-this-aft-regression-stable/53685/1 "2021-01-20T20:37:45Z")

</div>

I have adapted [this excellent post by @aaowens](https://github.com/JuliaStats/Distributions.jl/issues/702#issuecomment-516979703) to fit a right-censored AFT regression model, but for some reason the likelihood is not compatible with forwarddiff.

```julia
using Distributions,Optim
function AFT(pvec,data)
    function l(theta)
        α = exp(theta[1])
        θ = exp(theta[2])
        wdist = Weibull(α,θ)
        -sum(logpdf(wdist,theta[3:end]' * d[3:end] - d[2])^(d[1]) + logccdf(wdist,theta[3:end]' * d[3:end] - d[2])^(1-d[1]) for d in data)
    end
    l(pvec)
    optimize(l,pvec,BFGS(),autodiff = :forward)
end

a = 2.0
b = 2.0
n = 100
beta = [0.1,0.3,.8]
wb = Weibull(a,b)
X = rand(n,3)
delta = rand([0,1],n) # delta[i] = 1 iff Y[i] censored
Y = X * beta
for i in 1:n
    if delta[i] == 1
        Y[i] = Y[i] - Y[i]*rand()
    end
end

paramInit = fill(0.1,5)
resp = AFT(paramInit, collect(eachrow([delta Y X])))

julia> resp = AFT(paramInit, collect(eachrow([delta Y X])))
 * Status: failure (line search failed)

 * Candidate solution
    Final objective value: Inf

 * Found with
    Algorithm: BFGS

 * Convergence measures
    |x - x'| = 0.00e+00 ≤ 0.0e+00
    |x - x'|/|x'| = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')| = NaN ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = NaN ≰ 0.0e+00
    |g(x)| = 1.15e+02 ≰ 1.0e-08

 * Work counters
    Seconds run: 0 (vs limit Inf)
    Iterations: 1
    f(x) calls: 1
    ∇f(x) calls: 1

```

---

<div class="post-metadata">

**Author:** ![Paul\_Soderlind](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paul_soderlind/32/1753_2.png) [@Paul\_Soderlind](https://discourse.julialang.org/u/Paul_Soderlind)\
**Post date:** [January 20, 2021, 8:57pm UTC](https://discourse.julialang.org/t/how-to-make-this-aft-regression-stable/53685/2 "2021-01-20T20:57:50Z")

</div>

I am not really sure, but my guess is that your `l` function cannot handle the dual numbers of ForwardDiff.jl. In short, `theta` will be a dual number and the output of the function needs to be one too.

As an example, see cell 4 or this [notebook](https://github.com/PaulSoderlind/FinancialEconometrics/blob/master/Ch15_Garch.ipynb). In particular, see how `s2` is set up so it can host dual numbers.

---

<div class="post-metadata">

**Author:** ![mkarikom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkarikom/32/7096_2.png) [@mkarikom](https://discourse.julialang.org/u/mkarikom)\
**Post date:** [January 20, 2021, 9:35pm UTC](https://discourse.julialang.org/t/how-to-make-this-aft-regression-stable/53685/3 "2021-01-20T21:35:36Z")

</div>

Thanks @Paul_Soderlind, can you explain what it is about the explicit allocation of `s2` or the subsequent updates that ensures `LL_t` will be a dual?

```julia
    s2 = zeros(typeof(alpha),T) #works with ForwardDiff 
    s2[1] = omega + alpha*s2_0 + beta1*s2_0 #simple, but slow approach
    for t = 2:T #using filter() is perhaps quicker
      s2[t] = omega + alpha*u[t-1]^2 + beta1*s2[t-1]
    end

```

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [January 20, 2021, 10:21pm UTC](https://discourse.julialang.org/t/how-to-make-this-aft-regression-stable/53685/4 "2021-01-20T22:21:33Z")

</div>

This line doesn’t make sense.

```julia
-sum(logpdf(wdist,d[3:end]' * d[2:end])^(d[1]) * logccdf(wdist,d[3:end]' * d[2:end])^(1-d[1]) for d in data)

```

First, you are iterating over `data` which is a `Vector{Float64}`, so each `d` is an individual `Float64`, not a vector (this is the source of the indexing error you see). I suspect you meant for `data` to be a `Vector{Vector{Float64}}`, something like `[rand(wb, 10) for _ in 1:10]`? Or maybe some reshape of `Y`?

When you fix the shape of the data, `d[3:end]` is a different length than `d[2:end]`, so that product will not work. Perhaps that should be `d[3:end]' * d[2:end-1]`?

---

<div class="post-metadata">

**Author:** ![Paul\_Soderlind](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paul_soderlind/32/1753_2.png) [@Paul\_Soderlind](https://discourse.julialang.org/u/Paul_Soderlind)\
**Post date:** [January 21, 2021, 10:00am UTC](https://discourse.julialang.org/t/how-to-make-this-aft-regression-stable/53685/5 "2021-01-21T10:00:58Z")

</div>

`s2 = zeros(typeof(alpha),T)` means that when alpha is a dual number, then `s2` will be as well.

However, it is not the explicit allocation per se that is important. What matters is that the output (likelihood value) can be a dual number when the parameters are. You can easily check this by trying `ForwardDiff.gradient(l,paramInit)`.

As far as I can see, the sum, pdf and cdf of the Weibull should survive this, so the problem might be elsewhere.

---

<div class="post-metadata">

**Author:** ![mkarikom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkarikom/32/7096_2.png) [@mkarikom](https://discourse.julialang.org/u/mkarikom)\
**Post date:** [January 22, 2021, 1:01am UTC](https://discourse.julialang.org/t/how-to-make-this-aft-regression-stable/53685/6 "2021-01-22T01:01:00Z")

</div>

Fixed a bunch of bugs, thanks for all the suggestions.  
The following seems to run fine, but has identifiability issues. While the relative magnitude of positive beta seems roughly correct, it can’t get the sign right.  
I tried adding a ridge regularizer only on the linear coefficients (commented) but that led to stability issues as well

```julia
using Revise,Distributions,Optim,DistributedTopicModels,Survival,LineSearches,Plots,Random,LinearAlgebra

n = 10000
# beta = [1,10]
beta = [10,1]
# beta = [1,-1]
# beta = [-1,1]
s = Random.seed!(1)
X = rand(s,n,2)
delta = rand([0,1],n) # delta[i] = 1 iff Y[i] censored
μ = X * beta

histogram(μ)

λ = 1.1
α = .1
wd = Weibull(λ,α)
s = Random.seed!(1)
ϵ = rand(s,wd,n)

Y = log.(0.1(μ) + ϵ)
# Y = log.(0.1(μ) + ϵ)*0.1

histogram(Y)

Ycen = [delta[i] == 1 ? Y[i]-rand()*Y[i] : Y[i] for i in 1:n]
paramInit = fill(0.1,3)

function tAFT(pvec,data)
    function f(w)
        exp(w - exp(w))
    end
    function F(w)
        exp(-exp(w))
    end
    function l(theta)
        -sum(log(f((d[2]-d[3:end]'*theta[2:end])/theta[1])^d[1]) + log(F((d[2]-d[3:end]'*theta[2:end])/theta[1])^(1-d[1])) for d in data)
        # -sum(log(f((d[2]-d[3:end]'*theta[2:end])/theta[1])^d[1]) + log(F((d[2]-d[3:end]'*theta[2:end])/theta[1])^(1-d[1])) + 0.1*norm(theta[3:end])^2 for d in data)
    end
    optimize(l,pvec,BFGS(linesearch=LineSearches.BackTracking()),autodiff = :forward)
    # optimize(l,pvec,LBFGS(),autodiff = :forward)
end

resp = tAFT(paramInit, collect(eachrow([delta Y X])))
resp.minimizer

julia> resp = tAFT(paramInit, collect(eachrow([delta Y X])))
 * Status: success

 * Candidate solution
    Final objective value: 8.489688e+03

 * Found with
    Algorithm: BFGS

 * Convergence measures
    |x - x'| = 4.88e-15 ≰ 0.0e+00
    |x - x'|/|x'| = 1.87e-15 ≰ 0.0e+00
    |f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)| = 7.78e-08 ≰ 1.0e-08

 * Work counters
    Seconds run: 1 (vs limit Inf)
    Iterations: 76
    f(x) calls: 150
    ∇f(x) calls: 77

julia> resp.minimizer
3-element Array{Float64,1}:
 2.6089785118820537
 2.0258941072666534
 0.362426844935407

```