R中的Gompertz老化分析

Question

我有来自苍蝇实验的生存数据,该实验检查各种基因型的衰老率.我可以通过多种布局获取数据,因此您可以选择最适合您的答案.

一个数据帧(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分析(参见此处"生命表"的第二段).等式是:

μx=α*eβ ^*x

其中μx是给定时间的死亡概率,α是初始死亡率,β是衰老率.

我希望能够获得一个数据框,其中包含我的~640个基因型中的每一个的α和β估计值,以供稍后进一步分析.

我需要帮助从上面的数据帧到R中每个基因型的这些值的输出.

我查看了flexsurv可能包含答案的软件包,但我没有尝试查找和实现它.

Answer 1

这应该让你开始...

首先，为了使该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文档：

具有形状参数a和速率参数 b的 Gompertz 分布具有危险函数

H(x: a, b) = be^{ax}

所以看起来你的 alpha 是b，即速率，而 beta 是a，即形状。