Какова вероятность этого процесса?

10

Пациент госпитализирован. Продолжительность их пребывания зависит от двух факторов: тяжести травмы и суммы, которую их страховка готова заплатить, чтобы удержать их в больнице. Некоторые пациенты уйдут преждевременно, если их страховка решит прекратить оплачивать свое пребывание.

Предположим следующее:

1) Продолжительность пребывания распределена по Пуассону (предположим, что на данный момент это может быть или не быть реалистичным предположением) с параметром λ .

2) Различные планы страхования охватывают 7, 14 и 21 день пребывания. Многие пациенты уходят после 7, 14 или 21 дня пребывания (потому что их страховка заканчивается, и они должны уйти).

Если бы я получил данные из этого процесса, он мог бы выглядеть следующим образом:

введите описание изображения здесь

Как вы можете видеть, на 7, 14 и 21 день есть всплески. Это пациенты, которые уходят, когда заканчивается их страховка.

Очевидно, что данные могут быть смоделированы как смесь. Мне тяжело записать вероятность этого распределения. Это похоже на пуассон с нулевым надуванием, но инфляция составляет 7, 14 и 21.

Какова вероятность для этих данных? Какой мыслительный процесс стоит за вероятностью?

Димитрий Пананос
источник
Для начала вам нужно знать вероятности 7, 14 и 21-дневного времени вынужденного отъезда.
BruceET
1
Для меня это звучит как смесь пуассоновского и трех усеченных справа (в 7, 14 и 21) пуассоновских распределений. Записать их - это еще один шаг.
Карстен
@BruceET Я собираюсь сделать байесовский вывод по этой модели, поэтому я хотел бы записать это в самом общем случае.
Деметрий Пананос

Ответы:

9

В этом случае, я полагаю, что путь к решению существует, если мы наденем шляпу для анализа выживания. Обратите внимание, что, хотя в этой модели нет цензурированных предметов (в традиционном смысле), мы все равно можем использовать анализ выживаемости и говорить об опасностях предметов.

Нам нужно смоделировать три вещи в следующем порядке: i) совокупная опасность, ii) опасность, iii) логарифмическая вероятность.

H(t)H(t)=logS(t)TPoi(λ)

HT(t)=log(1Q(t,λ))=logP(t,λ)

Q,P

Теперь мы хотим добавить «опасности» на исходе страховки. Хорошая вещь о совокупных опасностях состоит в том, что они являются аддитивными, поэтому нам просто нужно добавить «риски» в моменты времени 7, 14, 21:

HT(t)=logP(t,λ)+a1(t>7)+b1(t>14)+c1(t>21)

>a,bc

c

HT(t)=logP(t,λ)+a1(t>7)+b1(t>14)+1(t>21)

h(t)

h(t)=1exp(H(t)H(t+1))

Включение нашей кумулятивной опасности и упрощение:

hT(t)=1P(t+1,λ)P(t,λ)exp(a1(t=7)b1(t=14)1(t=21))

iii) Наконец, запись логарифмической вероятности для моделей выживания (без цензуры) становится очень простой, если мы имеем опасность и совокупную опасность:

ll(λ,a,b|t)=i=1N(logh(ti)H(ti))

И вот оно!

a=log(1pa),b=log(1papb)log(1pa),pc=1(pa+pb)


Доказательство в пудинге. Давайте делать некоторые моделирования и логического вывода с использованием Lifelines' пользовательские модели семантики.

from lifelines.fitters import ParametericUnivariateFitter
from autograd_gamma import gammaincln, gammainc
from autograd import numpy as np

MAX = 1e10

class InsuranceDischargeModel(ParametericUnivariateFitter):
    """
    parameters are related by
    a = -log(1 - p_a)
    b = -log(1 - p_a - p_b) - log(1 - p_a)
    p_c = 1 - (p_a + p_b)
    """
    _fitted_parameter_names = ["lbd", "a", "b"]
    _bounds = [(0, None), (0, None), (0, None)]

    def _hazard(self, params, t):
        # from (1.64c) in http://geb.uni-giessen.de/geb/volltexte/2014/10793/pdf/RinneHorst_hazardrate_2014.pdf
        return 1 - np.exp(self._cumulative_hazard(params, t) - self._cumulative_hazard(params, t+1))

    def _cumulative_hazard(self, params, t):
        lbd, a, b = params
        return -gammaincln(t, lbd) + a * (t > 7) + b * (t > 14) + MAX * (t > 21)


def gen_data():
    p_a, p_b = 0.4, 0.2
    p = [p_a, p_b, 1 - p_a - p_b]
    lambda_ = 18
    death_without_insurance = np.random.poisson(lambda_)
    insurance_covers_until = np.random.choice([7, 14, 21], p=p)
    if death_without_insurance < insurance_covers_until:
        return death_without_insurance
    else:
        return insurance_covers_until


durations = np.array([gen_data() for _ in range(40000)])
model = InsuranceDischargeModel()
model.fit(durations)
model.print_summary(5)
"""
<lifelines.InsuranceDischargeModel: fitted with 40000 observations, 0 censored>
number of subjects = 40000
  number of events = 40000
    log-likelihood = -78845.10392
        hypothesis = lbd != 1, a != 1, b != 1

---
        coef  se(coef)  lower 0.95  upper 0.95      p  -log2(p)
lbd 18.05026   0.03353    17.98455    18.11598 <5e-06       inf
a    0.50993   0.00409     0.50191     0.51794 <5e-06       inf
b    0.40777   0.00557     0.39686     0.41868 <5e-06       inf
"""

¹ см. Раздел 1.2 здесь

Cam.Davidson.Pilon
источник