Appropriate Application of Survival Analysis

Question

I have an experiment that will produce observations of the time until an event occurs. Some basic properties are that

1. We count the number of events that have occurred at some point $t_1,...,t_n$.
2. Event times are interval censored, between $(t-1,t]$,
3. Individuals will not leave the trial between $t_1,...t_n$, that is an individual either experiences the event by $t_1$ or does not, at which point they are censored,
4. A large proportion of the individuals will not receive the event by $t_n$, when we terminate the experiment, and
5. I cannot assume any underlying parametric decay models.

It appears to be a natural application for survival analysis. However, it is complicated by the fact that it is trivial to repeat the experiment from the same initial configuration multiple times. In effect, we will have a set of $m_{1,...,s}$ event counts (where $s$ is the number of samples) for each observation time $t_i$. I'm relatively new to statistics, and I'm struggling to see how to apply survival analysis to this data (if it is even applicable and there aren't more appropriate methods to measure this kind of time-to-event data). My inclination is to construct the survival function around the mean number of observed events at each interval (i.e. $\bar{m}_{1,...,n}$), which should better approximate the expected number of events at each interval in the population, however I have no idea if this is appropriate or of its implications.

I've searched to no avail on Google Scholar, if anybody could point me to more material (or give me the correct nomenclature for what I'm trying to do), it would be appreciated.

Edit

Given that the $(t-1, t]$ intervals are uniform across all samples, say that I had the following matrix describing the cumulative number of individuals for which an event has occurred in each interval

$M = \left(\begin{array}{ccc} 0 & 24 & 35 & 52 & 60 & 71 \\ 0 & 22 & 38 & 57 & 64 & 75 \\ 0 & 26 & 34 & 55 & 62 & 72 \\ 0 & 21 & 32 & 52 & 61 & 73 \end{array}\right)$

where each row gives the event count for the same set of individuals at risk at $t=0$ across all samples (i.e. multiple instances of the experiment), and each column is an observation interval. I assume that by taking the mean number of events for each interval, I can get a better estimation of the expected population survival, so let $n$ denote the number of time intervals, $s$ denote the number of samples (experiment instances), then the vector

$\bar{M} = \left[ {{\sum_{i=1}^{s}M_{it}}\over{s}} \right]_{t=1...n}$

will be the mean number of observed events for each time interval.

My goal, then, is to use this as the input to the survival estimation. Let $f$ be the number of individuals at risk when $t=0$. Using the naive estimator (for now, given that the event intervals are uniform across all samples and there is no censorship until $t_n$), the survivor function could be estimated as:

$S(t) = {{f - \bar{M}_t}\over{f}}$

Which would (hopefully) be a better estimate of the population survival than any individual sample (a single row from $M$). To reformulate my question:

1. Is $\bar{M}$ an appropriate input to a survival function estimation? I haven't seen this approach in any of the material I've read.
2. As I'm really, painfully a novice at statistics, can somebody point me to some material (academic papers, textbooks, wikis etc. would be fine) on estimating the confidence interval and variance for this survival function estimation? I presume it will not be identical to standard formulations.

Apologies if my original question was confusing, I probably didn't include sufficient information.

Answer 1

I recently had a set of interval censored survival data, so I know exactly what you need. If you have ever used R, this should help.

If you don't want to assume a parametric form, how about an interval censored Cox proportional hazards model? The intcox package that would do this is no longer in the R repository. I would suggest imputing survival times and then using the coxph function from the survival library. Keep in mind that your standard errors will be too low using this method; you have not accounted for the uncertainty of not knowing the exact survival time. If you want interval censored survival estimates, use the icfit function from the interval package.

Another way analyze the effect of covariates on survival time is by using interval censored, nonparametric regression. See the R package ICE: http://cran.r-project.org/web/packages/ICE/ICE.pdf. You first need to impute the midpoints of survival time, then you do a local linear regression using the locpoly function from the np package. It's not as hard as it sounds.

Answer 2

Survival function is usually right-continuous since it is a distribution function, I will use $a_k:=[t_{k-1}, t_k),k=1,2,\cdots,n$ as the interval.

Let $T_{ij}$ and $C_{ij}$ be the true continuous survival and censoring time for subject $j$ in sample $i$, respectively. Both variables may not be observed directly, but only in one of the intervals $a_1,a_2,\cdots$. Furthermore, let $X_{ij}$ denote the interval within which $T_{ij}$ falls, essentially a discrete survival time, and similarly $\mathcal{C}_{ij}$ for $C_{ij}$. Then the censoring indicator is given by $\delta_{ij}=\mathbf{1}(X_{ij} \le \mathcal{C}_{ij})$.

The hazard function $h_{ij}(x)$ for the discrete survival time is defined as the conditional probability of the event occurring in the $x$th time interval given that it has not occurred prior to the $x-1$th interval, i.e

$$ h_{ij}(x)=P(X_{ij}=x | X_{ij} \ge x) $$

and the corresponding survival function $S_{ij}=P(X>x)$ can be written recursively using conditional product law:

$$ S_{ij}(x)=P(X_{ij}>x | X_{ij}\ge x)\cdots P(X_{ij}>1 | X_{ij} \ge 1)=\prod_{m=1}^x (1-h_{ij}(m)) $$

The likelihood function of the pair $(x_{ij},\delta_{ij})$ can be constructed as the product of two types of subjects, namely those who experienced an event at $x_{ij}$ ($X_{ij}=x_{ij},\delta_{ij}=1$) and those who were censored at $x_{ij}$ ($X_{ij}>x_{ij},\delta_{ij}=0$):

$$ \begin{split} \mathcal{L} &= \prod_{i=1}^{s}\prod_{j=1}^{n_{i}}[P(X_{ij}=x_{ij})]^{\delta_{ij}}[P(X_{ij}>x_{ij})]^{1-\delta_{ij}}\\ & = \prod_{i=1}^{s}\prod_{j=1}^{n_{i}}\left\{\left(h_{ij}(x_{ij})\prod_{m=1}^{x_{ij}-1}[1-h_{ij}(m)]\right)^{\delta_{ij}}\left( \prod_{m=1}^{x_{ij}}[1-h_{ij}(m)]\right)^{1-\delta_{ij}}\right\}\\ & = \prod_{i=1}^{s}\prod_{j=1}^{n_{i}}\left\{\left[\frac{h_{ij}(x_{ij})}{1-h_{ij}(x_{ij})}\right]^{\delta_{ij}}\prod_{m=1}^{x_{ij}}[1-h_{ij}(m)]\right\}. \end{split} $$ and corresponding log-likelihood function is:

$$ \ell=\sum_{i=1}^{M}\sum_{j=1}^{n_{i}}\left\{\delta_{ij}\log[\frac{h_{ij}(x_{ij})}{1-h_{ij}(x_{ij})}]+\sum_{m=1}^{x_{ij}}\log[1-h_{ij}(m)]\right\} $$

Now if we reconstruct our data into event history structure, that is recording at each interval, a event indicator variable $y_{ijk}$ for the $k$th interval of the $j$th subject from sample $i$, we can see that $\delta_{ij}\log[h_{ij}(x_{ij})/(1-h_{ij}(x_{ij})]$ in the above rewritten into $\sum_{k=1}^{x_{ij}}y_{ijk} \log[h_{ij}(k)/(1-h_{ij}(k))]$ (basically summing up all the 0's until the last observed interval of this subject, if he has event, it will be 1, if censored 0). Then we can rewrite our log-likelihood as

$$ \ell=\sum_{i=1}^s \sum_{j=1}^{n_i} \sum_{k=1}^{x_{ij}} \left\{ y_{ijk} \log h_{ij}(k) + (1 - y_{ijk}) \log [1 - h_{ij}(k)] \right\}. $$ This is identical to the log likelihood for a binary random variable $y_{ijk}$, but now with the proportion $p_{ijk}$ for an event in interval $k$ defined by $h_{ij}(k)$.

Now we finally can answer your question. If we would assume, that at same interval, $y_{ijk}$ is i.i.d. for different subject $j$ in sample $i$, and also across different sample, then $\bar{M_j}=(\sum_i n_i)^{-1}\sum_{i=1}^{s}\sum_{j=1}^{n_i}y_{ijk}$ is the appropriate estimator for $h_{ij}(k)=h(k)$.

And the appropriate estimator for $S(x)$ is therefore $\hat{S}(x)=\prod_{k=1}^{x}(1-\bar{M}_j)$.