5

我有来自苍蝇实验的生存数据，该实验检查了各种基因型的衰老率。这些数据以多种布局提供给我，因此您可以选择哪种布局，以最适合答案的为准。

一个数据框（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可能包含答案的软件包，但我未能尝试找到并实施它。

4

1 回答 1

3

这应该让你开始......

首先，要使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，形状。

于 2013-09-11T22:50:06.033 回答