研究蜜蜂常常需要做生存分析，就是比较两组或者多组处理之间的生存曲线有没有差别。 也可以用来比较行为， 比如两个不同处理的蜜蜂变成采集蜂的时间。老的办法是，每一个蜜蜂都需要一行。 累人， 而且容易出错。

在R里面， 数据是这样的结构。 第一列是处理，第二是日龄或者其他（这里是小时），第三是活着（0） 或者死了（1）。

下面这个数据实际上是大蜂螨被注射不同基因的双链RNA以后的成活情况。

这样很费时间， 而且数据记录本来就不是这样，是48小时死了2个，96小时死了3个，等等。 你还得在EXCEL里面一个一个去数， 数据对不对。

数据如果是下面这样多好！设下面的数据在名为sa-short-data.csv 的文件中。

N就是那个时间段的死亡数（1）， 如果实验结束时还有20个没有死，就是20，0（第四行数据）。 这样原来要20行的数据， 你一行搞定。

trt age N 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

写一个简单的code， 就可以使用了。

test < - read.csv(file="sa-short-data.csv",head=TRUE, sep=",")

test2=NULL

for (i in 1:nrow(test))
{
for (j in 1:test[i,3]) test2=rbind(test2, test[i,c(1,2,4)]) 
#loops for "N" times (element 3) in test and keeps rows of 1, 2 and 4 
#in the new matrix test2.
}

使用Package survival 里面的survdiff，得到下面结果。

rho=0是要的Log-Rank test的卡方和P值。

rho=1 会是出来 Wilcoxon test的卡方和P值。

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

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

Chisq= 28.3  on 1 degrees of freedom, p= 1.01e-07

这是R画出的图（昨天的另外一组数据）：

这个方法不一定只用在蜜蜂身上。 啥动物（或者植物，细菌）都可以用的。

SAS 的简易快速code：

DATA survival;  input treatment day dead censor;
Do bee=1 to dead;
output;
end;
DATALINES;
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
;
RUN;
proc lifetest plot=(s) graphics;
time day*censor(0); strata treatment;
run;

SAS会自动出来三个不同Tests的结果。 但是SAS每年要交钱， 我比较穷。 所以自己学习R了。 虽然研究生开始时学的SAS，用了几十年了。

Test of Equality over Strata
Test Chi-Square DF Chi-Square Pr

Log-Rank 28.3464 1 〈.0001
Wilcoxon 24.3984 1 〈.0001
-2Log(LR) 18.5757 1 〈.0001