我从一次果蝇实验中获得了生存数据,该实验检查了各种基因型的衰老率。我可以通过几种布局获得数据,因此哪种选择最适合您。
一个数据帧(wide.df)看起来像这样,其中每个基因型(Exp,其中有〜640个)都有一行,并且各天从第4天到第98天按水平顺序运行,每两天有新的死亡计数。
Exp Day4 Day6 Day8 Day10 Day12 Day14 ...
A 0 0 0 2 3 1 ...
我使用以下示例:
wide.df2<-data.frame("A",0,0,0,2,3,1,3,4,5,3,4,7,8,2,10,1,2)
colnames(wide.df2)<-c("Exp","Day4","Day6","Day8","Day10","Day12","Day14","Day16","Day18","Day20","Day22","Day24","Day26","Day28","Day30","Day32","Day34","Day36")
另一个版本是这样的,其中每天每个“Exp”都有一行,并记录该天的死亡人数。
Exp Deaths Day
A 0 4
A 0 6
A 0 8
A 2 10
A 3 12
.. .. ..
举个例子:
df2<-data.frame(c("A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A","A"),c(0,0,0,2,3,1,3,4,5,3,4,7,8,2,10,1,2),c(4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36))
colnames(df2)<-c("Exp","Deaths","Day")
我想做的是执行 Gompertz分析(See second paragraph of "the life table" here)。等式是:
μx=α* eβ* x
其中μx是给定时间的死亡概率,α是初始死亡率,β是衰老率。
我希望能够获得一个大约640个基因型具有α和β估计的数据框,以便以后进行进一步分析。
我需要帮助,从上述数据框到我在R中每个基因型的这些值的输出。
我已经浏览了
flexsurv
包,该包可能包含答案,但我尝试找到并实现它失败。 最佳答案
这应该让您开始...
首先,为了使flexsurvreg
函数起作用,您需要将输入数据指定为Surv
对象(来自package:survival
)。这意味着每个观察结果一行。
第一件事是从您提供的摘要表中重新创建“原始”数据。
(我知道rbind
效率不高,但是对于大集合,您始终可以切换到data.table
)。
### get rows with >1 death
df3 <- df2[df2$Deaths>1, 2:3]
### expand to give one row per death per time
df3 <- sapply(df3, FUN=function(x) rep(df3[, 2], df3[, 1]))
### each death is 1 (occurs once)
df3[, 1] <- 1
### add this to the rows with <=1 death
df3 <- rbind(df3, df2[!df2$Deaths>1, 2:3])
### convert to Surv object
library(survival)
s1 <- with(df3, Surv(Day, Deaths))
### get parameters for Gompertz distribution
library(flexsurv)
f1 <- flexsurvreg(s1 ~ 1, dist="gompertz")
给予
> f1$res
est L95% U95%
shape 0.165351912 0.1281016481 0.202602176
rate 0.001767956 0.0006902161 0.004528537
请注意,这是一个仅拦截的模型,因为您的所有基因型均为
A
。一旦如上所述重新创建了每个观测数据,就可以在多个生存对象上循环执行此操作。
从
flexsurv
文档:这样看来,您的alpha就是 b 的比率,而beta就是和的形状。
关于r - R中的Gompertz老化分析,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/17718566/