Empirical Survival Function with Doubly Interval Censored and Right Truncated Data

DICRT.jl

using Distributions
using Random
using Plots

function p_up!(p,n,m,j1,j2,dj)
    bj = zeros(m)
    for i in 1:n
        if j1[i]>0 && j2[i]>0
        q = sum(p[j1[i]:j2[i]])
        if q > 0.0
        b = (1.0-q)/q
        ind = setdiff(1:m,j1[i]:j2[i])
        r = vec(p[ind])
        r ./= sum(r)
        bj[ind] += b*r
    end
    end
    end
    num = dj + bj
    copy!(p, num ./ sum(num))
end

function d_up(p,j1,j2,n,m)
    dj = zeros(m)
    for i in 1:n
        if j1[i]>0 && j2[i]>0
            q = p[j1[i]:j2[i]]
            q ./= sum(q)
            dj[j1[i]:j2[i]] += q
        end
    end
    return dj
end

function acount(L,R,breaks,n,m)
    afirst = zeros(Int,n)
    alast = zeros(Int,n)
    for i in 1:n
        alpha = L[i] .<= breaks[1:(m-1)] .&& breaks[2:m] .<= R[i]
        if any(alpha)
            afirst[i] = findfirst(alpha)
            alast[i] =  findlast(alpha)
        end
    end
    return afirst, alast
end

function bcount(U,breaks,n,m)
    bfirst = zeros(Int,n)
    blast = zeros(Int,n)
    for i in 1:n
        beta = breaks .<= U[i]
        if any(beta)
            bfirst[i] = findfirst(beta)
            blast[i] =  findlast(beta)
        end
    end
    return bfirst, blast
end

function setbreaks(trow, O=0.0)
    tj = sort(unique(max.(O, trow)))
    return tj
end

function make_dic(rng, y, at, tau_e, tau_s)
    S = y + at
    ue = rand(rng)
    E_R = at + tau_e * ue
    E_L = max(0.0, at - tau_e * (1.0 - ue))
    us = rand(rng)
    S_R = S + tau_s * us
    S_L = max(0.0, S - tau_s * (1.0 - us))
    S = y + at
    return E_L, E_R, S_L, S_R
end

function make_dicrt(rng, dist, Tmax,  N)
    EL = zeros(N)
    ER = zeros(N)
    SL = zeros(N)
    SR = zeros(N)
    for i in 1:N
        y = rand(rng, dist)
        at = rand(rng, Uniform(0.0, Tmax))
        tau_e = rand(rng, Exponential(1.0))
        tau_s = rand(rng, Exponential(1.0))
        EL[i], ER[i], SL[i], SR[i] = make_dic(rng, y, at, tau_e, tau_s)
    end
    d = SR .<= Tmax
    return EL[d], ER[d], SL[d], SR[d]
end

#EL, ER, SL, SR, = dat
function SurvDICRT(EL, ER, SL, SR, Tmax, iter=100, tol=1e-8)
    S = 0.5*(SL+SR)
    tj = setbreaks([S-EL;S-ER])
    n = length(EL)
    m = length(tj)
    aind = acount(S-ER, S-EL,tj,n,m)
    p = (1.0/m)*ones(m)
    bind = bcount(Tmax.-SR,tj,n,m)
    dj = d_up(p,aind[1],aind[2],n,m)
    p_up!(p, n, m, bind[1], bind[2],dj)
    con = false
    count = 0
    p2 = p #こういうとこcopy!()とか使ったほうがいいですか？
    for it in 1:iter
        count += 1
        dj = d_up(p,aind[1],aind[2],n,m)
        p_up!(p, n, m, bind[1], bind[2], dj)
        con = all(abs.(p2-p) .< tol)
        p = p2
    if con || any(isnan.(p))
        break
    end
    end
    return tj, 1.0 .- cumsum(p), con, count
end

function simit(rng, N, it, dist)
    dat = make_dicrt(rng, dist, 100.0, N)
    fit0 = SurvDICRT(dat[1], dat[2], dat[3], dat[4], 100.0)
    tim = fit0[1]
    surv = fit0[2]
    for i in 2:it
        dat = make_dicrt(rng, dist, 100.0, N)
        fit0 = SurvDICRT(dat[1], dat[2], dat[3], dat[4], 100.0)
        append!(tim, fit0[1])
        append!(surv, fit0[2])
    end
    return tim, surv
end

rng = MersenneTwister()
mlnorm = MixtureModel(LogNormal[LogNormal(1.0, 1.0), LogNormal(2.0, 0.1)], [0.3, 0.7])
res = simit(rng, 200, 100, mlnorm)
plot(res[1], res[2], legend=false, linetype=:steppost, color=:gray, alpha=0.3)
plot!(x -> ccdf(mlnorm,x),color=:royalblue, lw=2)
png("./Desktop/test.png")

"Maximum Likelihood or Bayesian Estimation for Parametric and Non-Parametric Survival Function of Incubation Period with Doubly Interval Censored and Right Truncated Data" (in preparation)