MortalityTables LsqFit

gistfile1.txt

using LsqFit, MortalityTables, Plots, Distributions, Optim
#simple parametric survival function
data = [
1 0.99
2 0.98
3 0.95
4 0.9
5 0.8
6 0.65
7 0.5
8 0.38
9 0.25
10 0.2
11 0.1
12 0.05
13 0.02
14 0.01]

#1
plot(data[:,1], data[:,2])
#Модель на основе функции выживания
#Survival function model
@. model(x, p) = survival(MortalityTables.Weibull(;m = p[1],σ = p[2]), x)
fit = curve_fit(model, data[:,1], data[:,2], [1.0, 1.0])
plot!(data[:,1], model(data[:,1], fit.param))

#2
#ML оценка
#Генерируем 100 случайных наблюдений
#ML estimale with; generation of 100 datapoints
t = rand(Weibull(fit.param[2], fit.param[1]), 100)
#Нет цензурированных
#No censored data
fit_mle(Weibull, t)

#3
#Некоторые данные отмечены как цензурированные
#Вектор цензурированных данных
#Some datapoints marked as censored
c = collect(trues(100))
c[[1,3,7,9]] .= false
#ML function
survmle(x) = begin
    ml = 0.0
    for i = 1:length(t)
        if c[i]
            ml += logpdf(Weibull(x[2], x[1]), t[i]) #if not censored log(f(x))
        else
            ml += logccdf(Weibull(x[2], x[1]), t[i]) #if censored log(1-F)
        end
    end
    -ml
end
opt = Optim.optimize(survmle, [2.0,5.0]; method = Optim.Newton())
Optim.minimizer(opt)

#4
#Подгонка модели по эмпирической функции KM
#Fitting for survival function from Kaplan Meier
import LsqFit
using MortalityTables, Plots, Survival
#t- время;c - событие
#t- time vector;c - censored events vector
km =fit(KaplanMeier, t, c)
plot(km.times, km.survival; labels="Empirical")
@. model(x, p) = survival(MortalityTables.Weibull(;m = p[1],σ = p[2]), x)
mfit = LsqFit.curve_fit(model, km.times, km.survival, [2.0, 2.0])
plt = plot!(km.times, model(km.times, mfit.param), labels="Theoretical")