What is the likelihood for this process?

A patient is admitted to the hospital. Their length of stay depends on 2 things: The severity of their injury, and how much their insurance is willing to pay to keep them in the hospital. Some patients will leave prematurely if their insurance decides to stop paying for their stay.

Assume the following:

1) The length of stay is poisson distributed (just assume this for now, it may or may not be a realistic assumption) with parameter $$\lambda$$.

2) Various insurance plans cover 7, 14, and 21 day stays. Many patients will leave after 7,14, or 21 day stays (because their insurance runs out and they must leave).

If I were to obtain data from this process, it might look as follows:

As you can see, there are spikes at the 7, 14, and 21 day mark. These are patients that leave when their insurance ends.

Clearly, the data can be modeled as a mixture. I'm having a tough time writing down the likelihood for this distribution. It's like a zero inflated poisson, but the inflation is at 7, 14, and 21.

What is the likelihood for this data? What is the thought process behind the likelihood?

• To start, you'd need to know the probabilities of 7, 14, and 21-day forced leaving times. – BruceET Jul 13 at 17:26
• To me this sounds like a mixture of a Poisson and three right-truncated (at 7, 14 and 21) Poisson distributions. Writing them down is another step altogether. – Carsten Jul 13 at 17:43
• @BruceET I'm going to do Bayesian inference on this model, so I'd like to write it down in the most general case. – Demetri Pananos Jul 13 at 18:10

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 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
"""



¹ see Section 1.2 here