這個過程的可能性有多大？
一名患者被送往醫院。他們的住院時間取決於兩件事：他們受傷的嚴重程度，以及他們的保險願意支付多少讓他們留在醫院。如果他們的保險決定停止支付住院費用，一些患者會提前離開。
假設如下：
1）停留時間是泊松分佈的（現在只是假設，它可能是也可能不是一個現實的假設）與參數 $ \lambda $ .
- 各種保險計劃涵蓋 7、14 和 21 天的住宿。許多患者會在住院 7,14 或 21 天后離開（因為他們的保險用完了，他們必須離開）。
如果我要從這個過程中獲取數據，它可能如下所示：
如您所見，在第 7 天、第 14 天和第 21 天出現峰值。這些是在保險到期時離開的患者。
顯然，數據可以建模為混合。我很難寫下這種分佈的可能性。這就像一個零膨脹泊松，但膨脹是 7、14 和 21。
這個數據的可能性是多少？可能性背後的思考過程是什麼？
在這種情況下，如果我們戴上生存分析的帽子，我相信存在一條解決方案。請注意，即使該模型沒有審查對象（在傳統意義上），我們仍然可以使用生存分析並談論對象的危害。
我們需要按此順序對三件事進行建模：i) 累積風險，ii) 風險，iii) 對數似然。
i) 我們將逐步完成 i) 部分。什麼是累積風險， $ H(t) $ ，泊松隨機變量？對於離散分佈，有兩種定義方式¹，但我們將使用定義 $ H(t) = -\log{S(t)} $ . 所以累積風險為 $ T \sim Poi(\lambda) $ 是
$$ H_T(t) = -\log{(1 - Q(t, \lambda))} = -\log{P(t, \lambda)} $$
在哪裡 $ Q, P $ 分別是上、下正則化 gamma 函數。
現在我們要添加保險用完的“危害”。累積風險的好處是它們是相加的，所以我們只需要在時間 7、14、21 處添加“風險”：
$$ H_{T'}(t) = -\log{P(t, \lambda)} + a\cdot\mathbb{1}{(t>7)} + b\cdot\mathbb{1}{(t>14)} + c\cdot\mathbb{1}_{(t>21)} $$
啟發式地，患者受到背景“泊松”風險，然後是 7、14 和 21 的逐點風險。（因為這是累積風險，我們累積這些逐點風險，因此 $ > $ .) 我們不知道是什麼 $ a, b $ ， 和 $ c $ 是，但我們稍後會將它們與我們的保險用完的概率聯繫起來。
實際上，由於我們知道21 是上限，之後所有患者都被移除，我們可以設置 $ c $ 成為無窮大。
$$ H_{T'}(t) = -\log{P(t, \lambda)} + a\cdot\mathbb{1}{(t>7)} + b\cdot\mathbb{1}{(t>14)} + \infty \cdot\mathbb{1}_{(t>21)} $$
ii) 接下來我們使用累積風險來獲得風險， $ h(t) $ . 這個公式是：
$$ h(t) = 1 - \exp{(H(t) - H(t+1))} $$
插入我們的累積風險，並簡化：
$$ h_{T'}(t) = 1 - \frac{P(t+1, \lambda)}{P(t, \lambda)} \exp(-a\cdot\mathbb{1}{(t=7)} - b\cdot\mathbb{1}{(t=14)} - \infty \cdot\mathbb{1}_{(t=21)}) $$
iii）最後，一旦我們有了風險和累積風險，為生存模型（沒有審查）編寫對數可能性是非常容易的：
$$ ll(\lambda, a, b ;|; t) = \sum_{i=1}^N \left(\log h(t_i) - H(t_i)\right) $$
它就在那裡！
存在連接我們的逐點風險係數和保險長度概率的關係： $ a = -\log(1 - p_a), b = -\log(1 - p_a - p_b) - \log(1 - p_a), p_c = 1 - (p_a + p_b) $ .
證據就在布丁裡。讓我們使用生命線的自定義模型語義進行一些模擬和推理。
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:"InsuranceDischargeModel_estimate", fitted with 40000 total observations, 0 right-censored observations> number of observations = 40000 number of events observed = 40000 log-likelihood = -78754.92088 hypothesis = lbd != 1, a != 1, b != 1 --- coef se(coef) coef lower 95% coef upper 95% z p -log2(p) lbd 18.01220 0.03351 17.94652 18.07789 507.62368 <5e-06 inf a 0.51426 0.00411 0.50620 0.52232 -118.14024 <5e-06 inf b 0.40674 0.00557 0.39582 0.41767 -106.43953 <5e-06 inf --- """
¹ 請參閱此處的第 1.2 節