Survival analysis in R: time-saving method

In my last blog, we successfully reproduced the SAS method for survival analysis, but it was a pain (and error prone) to do one line each for one testing subject. In other words, if on day 10, you had 50 dead bees, you will have to enter 50 lines of data!

With the help of a R guru from Penn State, the following code works for our example:

  1. assume the external file has the following data:

>test
trt age dead censor
gfp  24    1      1
gfp  48    2      1
gfp  96    3      1
gfp  96   20      0
rpl8  24    5      1
rpl8  48   12      1
rpl8  96   25      1
rpl8  96    3      0

the code below will translate each line into multiple lines so you do not have to do the work!

for example, line 4 will produce “gfp 96 0” 20 times so survival analysis package will interpret the data correctly.

d3=NULL d4=NULL
for (i in 1:nrow(test))
{
  d4=rep(c(test$trt[i],test$age[i],test$censor[i]), times=test$dead[i])
 dim(d4)=c(3,test$dead[i])
 d4=t(d4)
 d3=rbind(d3,d4)
}

dd=data.frame(d3)
>names(dd)=c("trt","age","censor")

dd now contains the correct data for survival analysis (same as outputed by SAS, except that the trt has 1 and 2 instead of gfp and rpl8).

dd
trt age censor
1    1  24      1
2    1  48      1
3    1  48      1
4    1  96      1
5    1  96      1
6    1  96      1
7    1  96      0
8    1  96      0
9    1  96      0
10   1  96      0
11   1  96      0
12   1  96      0
13   1  96      0
14   1  96      0
15   1  96      0
16   1  96      0
17   1  96      0
18   1  96      0
19   1  96      0
20   1  96      0
21   1  96      0
22   1  96      0
23   1  96      0
24   1  96      0
25   1  96      0
26   1  96      0
27   2  24      1
28   2  24      1
29   2  24      1
30   2  24      1
31   2  24      1
32   2  48      1
33   2  48      1
34   2  48      1
35   2  48      1
36   2  48      1
37   2  48      1
38   2  48      1
39   2  48      1
40   2  48      1
41   2  48      1
42   2  48      1
43   2  48      1
44   2  96      1
45   2  96      1
46   2  96      1
47   2  96      1
48   2  96      1
49   2  96      1
50   2  96      1
51   2  96      1
52   2  96      1
53   2  96      1
54   2  96      1
55   2  96      1
56   2  96      1
57   2  96      1
58   2  96      1
59   2  96      1
60   2  96      1
61   2  96      1
62   2  96      1
63   2  96      1
64   2  96      1
65   2  96      1
66   2  96      1
67   2  96      1
68   2  96      1
69   2  96      0
70   2  96      0
71   2  96      0

> library(survival)
> survdiff(Surv(age,censor)~trt, data=dd, rho=0)
Call:
survdiff(formula = Surv(age, censor) ~ trt, data = dd, rho = 0)

N Observed Expected (O-E)^2/E (O-E)^2/V
trt=1 26        6     20.2      9.99      28.3
trt=2 45       42     27.8      7.27      28.3

Chisq= 28.3  on 1 degrees of freedom, p= 1.01e-07
Author: Zachary Huang

Leave a Reply

Your email address will not be published. Required fields are marked *