我从一次果蝇实验中获得了生存数据,该实验检查了各种基因型的衰老率。我可以通过几种布局获得数据,因此哪种选择最适合您。

一个数据帧(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/

10-12 04:04
查看更多