In this case, I believe a path to a solution exists if we put on our survival analysis hat. Note that even though this model has no censored subjects (in the traditional sense), we can still use survival analysis and talk about hazards of subjects.
We need to model three things in this order: i) the cumulative hazard, ii) the hazard, iii) the log likelihood.
i) We'll do part i) in steps. What is the cumulative hazard, $H(t)$, of a Poisson random variable? For a discrete distribution, there are two ways to define it¹, but we will use the definition $H(t) = -\log{S(t)}$. So the cumulative hazard for $T \sim Poi(\lambda)$ is
$$ H_T(t) = -\log{(1 - Q(t, \lambda))} = -\log{P(t, \lambda)} $$
where $Q, P$ is the upper, lower regularized gamma function respectively.
Now we want to add the "hazards" of the insurance running out. The nice thing about cumulative hazards is that they are additive, so we simply need to add "risks" at the times 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)} $$
Heuristically, a patient is subject to a background "Poisson" risk, and then point-wise risks at 7, 14, and 21. (Because this is a cumulative hazard, we accumulate those point-wise risks, hence the $>$.) We don't know what $a, b$, and $c$ are, but we will later connect them to our probabilities of insurance running out.
Actually, since we know 21 is the upper limit and all patients are removed after that, we can set $c$ to be infinity.
$$ 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) Next we use the cumulative hazard to get the hazard, $h(t)$. The formula for this is:
$$h(t) = 1 - \exp{(H(t) - H(t+1))}$$
Plugging in our cumulative hazard, and simplifying:
$$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) Finally, writing the log likelihood for survival models (without censoring) is super easy once we have the hazard and cumulative hazard:
$$ll(\lambda, a, b \;|\; t) = \sum_{i=1}^N \left(\log h(t_i) - H(t_i)\right)$$
And there it is!
There exists the relationships that connects our point-wise hazard coefficients and the probabilities of insurance lengths: $a = -\log(1 - p_a), b = -\log(1 - p_a - p_b) - \log(1 - p_a), p_c = 1 - (p_a + p_b)$.
The proof is in the pudding. Let's do some simulations and inference using lifelines' custom model semantics.
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
---
"""
¹ see Section 1.2 here
