9
$\begingroup$

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:

enter image description here

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?

$\endgroup$
  • $\begingroup$ To start, you'd need to know the probabilities of 7, 14, and 21-day forced leaving times. $\endgroup$ – BruceET Jul 13 at 17:26
  • 1
    $\begingroup$ 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. $\endgroup$ – Carsten Jul 13 at 17:43
  • $\begingroup$ @BruceET I'm going to do Bayesian inference on this model, so I'd like to write it down in the most general case. $\endgroup$ – Demetri Pananos Jul 13 at 18:10
8
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.