I did an experiment where I have ten different treatments, with six replicates each. I noted the survival of my animals from Day 4 to Day 12. Each replicate had 12 animals so I noted how many survived each day. I don't have individual data for each animal. So my data looks something like this:

Treatment      Replicate     Survival_D4  ... Survival_D12

Control        1             12              12
Control        2             11              11
Control        3             12              4
Control        4             9               6
Control        5             10              8
Control        6             11              12
Treatment 1    1             12              12
Treatment 9    6             12              8

My supervisor recommended a Kaplan Meier Survival Analysis, however, I don't have the 0-1 data for status.

I calculated the cumulative survival probability for each replicate and was able to make a plot.

I did this by calculating the survival probability each day by dividing the number of survivors by 12 (since I started with 12 animals). Then for the cumulative survival probability, I multiplied the survival probability with the cumulative survival probability of the day before.

But I want to know whether there are significant differences between my treatments. I want to see whether f.e. a certain treatment leads to faster death.

So my question is, can I still use Kaplan Meier and how would I do so? If I just use my own cumulative survival probability, how would I check for a significant difference?

Also, I averaged the cumulative survival probabilities of my replicates so that I only have one value per day for each treatment. Is this okay to do?


2 Answers 2


You do have the individuals' data. It's just not formatted in a nice way at the moment, and it will be some work to manipulate the data to extract it.

Let's say Treament X, Replicate Y has this data:

Survivors D4 = 10
Survivors D5 = 10
Survivors D6 = 10
Survivors D7 = 9
Survivors D8 = 9
Survivors D9 = 9
Survivors D10 = 4
Survivors D11= 4
Survivors D12 = 4

You know you started with 12, and Survivors D4 = 10. Therefore, 2 died between the start and D4. So you have two records with the event in the interval (0,4).

The next event happens between D6 and D7. One more individual died. So you have one record with the event in the interval (6,7).

5 events happened between D9 & D10. So you have five records with the event in the interval (9,10).

No more events happened after this - the rest survived, so their event is in the interval (12, infinity) (4 records).

So you individual-level data set for this dummy example is

(12, infinity) 
(12, infinity) 
(12, infinity) 
(12, infinity) 

Now you have to program that logic for all your groups. This should give data that you can input into a survival analysis software (at least, I know you should be able to do it for the R survival package).


As the answer from @AlexJ says (+1), you do have data that can be reformatted in a way that can take advantage of well-documented statistical survival functions. What can be surprising is that, when there's at most one event possible per individual, the individual identities often don't matter.

Recall the formula for the Kaplan-Meier estimator. It's calculated at each event time $t_i$ based on the number that had the event at that time, $d_i$, and the number that were at risk of the event, $n_i$. The estimated survival at time $t$,$\widehat S(t)$, is then:

$$\widehat S(t) = \prod\limits_{i:\ t_i\le t} \left(1 - \frac{d_i}{n_i}\right).$$

There's nothing there that requires knowing individual identities, just how many were at risk and how many experienced the event at each time. It's only in more complicated cases (e.g., more than one event possible per individual, some parametric survival models with time-varying covariates) that you need to keep track of individual identities.


In your experimental design, however, you probably do need to keep track of the "Replicate" associated with each individual. That will allow you to take potential correlations of outcomes within a Replicate into account. I suggest that you re-name the Replicate values, now all numbers 1 to 6, to include both the Treatment and the Replicate number, so there's no confusion that Replicate 4 in Control is the same as Replicate 4 in Treatment 5: e.g., use names like "R4C" and "R4T5" to keep them clearly distinguished.

"Interval-censored" data

The format that @AlexJ recommends represents the data as "interval censored"; that is, you know that the event occurred sometime between the two limits of the time interval, but not exactly when. In that format, the 0/1 status data you are thinking of is effectively replaced by the time value at the end of each interval: if it's finite, there was an event during the interval (status = 1), if it's infinite (Inf in R) the time to the event is right-censored (status = 0).

Need for regression model

The problem with a Kaplan-Meier analysis in this case is that it will be hard to distinguish treatment differences very well when you have so many treatments. The usual "log-rank" test done as part of such analysis can tell you if there are any differences among treatments, but not which ones are significantly different from control or which differ from each other. To get those types of results it's best to build a regression model of some type.

A Cox proportional hazards model is a frequent choice. In your situation that would not require keeping track of individuals, either. There are tools available in the R icenReg package for Cox models on interval-censored data. I'm not sure, however, that those tools can take the "Replicate" values into account in an efficient way.

One choice, often made in practice even though it's not theoretically exact, would be to perform a standard Cox model where the event time is recorded as the end of the time interval even though you don't know exactly when the event occurred. You can think of that as a model of when you observed the death rather than when it occurred. A cluster term in the Cox model could then take "Replicate" values into account.

With your current data setup, you can generate a new data table for that pretty simply. For each time point, for each row in your current data, create a number of rows in the new table equal to the difference in survival since the previous time point,* copying over all the information about Treatment/Replicate, setting the Time for those rows to be that time point, and the Event marker to be 1. At the last time point, make copies with the Event marker of 0 instead for the number that survived at last observation. That way, you end up with a number of rows equal to the original number of animals, with the last observation time and the status at that time (even though this doesn't distinguish which animals are which).

Discrete-time survival

As discussed on this page, however, with so few time points you might be better off with a discrete-time survival model. That's just a binomial generalized linear model with data in a long "person-period" format: one row for each individual at risk during each time interval, with indicators of the time interval, the treatment, the replicate, and 1/0 for event/no-event during the time interval.

If you have a new data table as recommended above, with one row per individual including the time point and status at the last observation, then you can use the dataLong() function in the R discSurv package to generate the (very long) person-period data format. A choice of timeAsFactor=TRUE will be most similar to the Kaplan-Meier and Cox models, in that there is no assumption about a numeric form of baseline survival over time. That generates the long-form data with a "timeInt" value representing the time interval and a new outcome variable "y" for the status of the individual at the end of that interval.

With such person-period data, you could use a binomial generalized mixed model to take the Replicate values into account, treating Replicate as a random intercept. That's implemented by the glmer() function in the R lme4 package. To be closest to a Cox model you could specify the "cloglog" link instead of the default "logit" link used for binomial logistic regression. The model could look something like this:

discreteModel <- glmer(y ~ timeInt + Treatment + 
                         data = dataInLongFormat,
                         family = binomial(link = "cloglog"))

Then you can use post-modeling tools like those in the R emmeans package to evaluate specific comparisons of interest. For example, you can evaluate the difference of each treatment from control if you specify "trt.vs.ctrl" contrasts in that package.

*Before you do this, double-check your data so that you don't have any "resurrections" like you show for Replicate 6 of Control in your sample data.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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