Prediction: Discrete-Time Event History (Survival) Model in R

Question

I want to perform predictions of a logistic discrete time hazard model, given by:

$\lambda (t)= P(T = t| T \geq t)=\frac{exp(\alpha_{0t} + X_i^T \beta)}{1+ exp(\alpha_{0t} + X_i^T \beta)}$

How to fit this model is outligned in a previous post: Discrete-Time Event History (Survival) Model in R.

What I do not understand is the "jump" or equivalence of the Survival Analysis World and the Binary Model World.

Let's start with the likelihood for the discrete hazard model:

$L=\prod_{i=1}^{N}[(\frac{exp(\alpha_{0t} + X_i^T \beta)}{1+ exp(\alpha_{0t} + X_i^T \beta)})^{\sigma_i} (\frac{exp(\alpha_{0t} + X_i^T \beta)}{1+ exp(\alpha_{0t} + X_i^T \beta)})^{1-\sigma_i} \prod_{j=1}^{t_i-1}(1- \frac{exp(\alpha_{0j} + X_i^T \beta)}{1+ exp(\alpha_{0j} + X_i^T \beta)})]$

where $\sigma_i=1$ if the observation $i$ is censored at $t_i$. (Survival Analysis World)

The well known trick is now to introduce as many binary variables $y_{it}$ as subject $i$ survived. The underlying model is now $P(Y_{it}=1)=\frac{exp(\alpha_{0t} + X_i^T \beta)}{1+ exp(\alpha_{0t} + X_i^T \beta)}$ and one uses a modified model matrix such that the likelihood of this binary model is:(Binary Model World)

$L=\prod_{i=1}^{N} \prod_{k=1}^{t_i}(\frac{exp(\alpha_{0t} + X_i^T \beta)}{1+ exp(\alpha_{0t} + X_i^T \beta)})^{y_{it}} (1-\frac{exp(\alpha_{0t} + X_i^T \beta)}{1+ exp(\alpha_{0t} + X_i^T \beta)})^{1-y_{it}} $

This procedure is precisely described in: http://www.statisticalhorizons.com/wp-content/uploads/Allison.SM82.pdf

Now my questions are:

1. I understand that the likelihood functions are the same, but in the Survival Analysis World I specified a conditional probability, whereas in the Binary Model World I have an unconditional probability. So according to these specifications: $\lambda (t)= P(T = t| T \geq t)=P(Y_{it})$ right?

2. So to obtain the unconditional prob. of the event happening at time $t$ I could either use $P(T=t)=\lambda (t) \prod_{t=1}^{t_i-1} (1-\lambda (t))$ or simply $P(Y_{it}=1)$ -right? So for the latter approach -$P(Y_{it}=1)$ - I only need 1 predict statement in R whereas in the first case, I need multiple predict statements.

Since these two prediction approaches will yild different prob. I'm stuck on how to align these two approaches.

My final aim is to get the unconditional prob. of the event as outligned in 2.

Many thanks in advance!

Answer 1

There is a subtle mistake right here:

The well known trick is now to introduce as many binary variables $y_{it}$ as subject $i$ survived. The underlying model is now $P(Y_{it}=1) = \frac{exp(\alpha_{0t}+X^T_i\beta)}{1+exp(\alpha_{0t}+X^T_i\beta)}$

The underlying model did NOT change when you used this trick. The underlying model is still:

$P(T=t | T \ge t) = \frac{exp(\alpha_{0t}+X^T_i\beta)}{1+exp(\alpha_{0t}+X^T_i\beta)}$

The trick is to notice that this (one and only) underlying model just happens to have the exact same likelihood function as a logit in what you are calling Binary Model World. Binary model world is wholly hypothetical, non-existent, unreal, and unmodelled, however. It's just a numerical trick. $P(Y_{it}=1)$ corresponds to nothing in the real world. There is no random variable $Y_{it}$ unless $Y_{i1}=Y_{i2}=\cdots=Y_{it-1}=0$. $Y_{it}$ is kind of a conditional object. It does not exist (or, to finesse this issue, we don't see it) unless a bunch of other random variables take on the value 0.

To belabor the point, if you want to try to calculate $P(T=t)$ in Binary Model World, you have to calculate $P(Y_{it}=1 \text{ and we see }Y_{it})$. This is the same as $P(Y_{it}=1, Y_{i1}=Y_{i2}=\cdots=Y_{it-1}=0)$. And that is given by:

\begin{align} P(T=t) &= P(Y_{it}=1, Y_{i1}=Y_{i2}=\cdots=Y_{it-1}=0)\\ &= P(Y_{i1}=0) P(Y_{it}=1, Y_{i2}=\cdots=Y_{it-1}=0|Y_{i1}=0)\\ &= \ldots\\ &= P(Y_{i1}=0) \cdot P(Y_{i2}=0 | Y_{i1}=0) \cdot P(Y_{i3}=0 | Y_{i1}=0, Y_{i2}=0) \cdots P(Y_{it}=1 | Y_{i1}=Y_{i2}=\cdots=Y_{it-1}=0) \end{align}

Finally, that is equal to:

$P(T=t) = \frac{exp(\alpha_{0t}+X^T_i\beta)}{1+exp(\alpha_{0t}+X^T_i\beta)} \cdot \Pi_{\tau=1}^{t-1} \frac{1}{1+exp(\alpha_{0t}+X^T_i\beta)}$

This is the same as your expression calculating $P(T=t)$ using $\lambda(t)$.