5

さまざまな遺伝子型の老化率を調べたハエの実験からの生存データがあります。データはいくつかのレイアウトで利用できるので、どれを選択するかはあなた次第です。

1 つのデータフレーム (wide.df) は次のようになります。各遺伝子型 (Exp、うち 640 まであります) には行があり、日は 4 日目から 98 日目まで水平方向に連続して実行され、2 日ごとに新しい死亡者数がカウントされます。

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")

私がやりたいことは、ゴンペルツ分析を実行することです(ここで「生命表」の 2 番目の段落を参照してください)。方程式は次のとおりです。

μx = α*e β*x

ここで、 μxは特定の時点での死亡確率、αは初期死亡率、βは老化率です。

後でさらに分析するために、約 640 の遺伝子型のそれぞれについてαとβの推定値を持つデータフレームを取得できるようにしたいと考えています。

上記のデータフレームから、R の遺伝子型ごとにこれらの値を出力するための助けが必要です。

答えが含まれている可能性のあるパッケージflexsurvを調べましたが、それを見つけて実装する試みに失敗しました。

4

1 に答える 1

3

これで始められるはずです...

まず、flexsurvreg関数を機能させるには、入力データをオブジェクトとして指定する必要がありますSurv(からpackage:survival)。これは、観測ごとに 1 行を意味します。

最初に、提供した要約テーブルから「生」データを再作成します。（効率的ではないことはわかっていますが、大きなセットの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を持つゴンペルツ分布にはハザード関数があります

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

したがって、アルファはb、レートであり、ベータはa、形状であるように見えます。

于 2013-09-11T22:50:06.033 に答える