Gompertz Aging Analyse in R

Question

Ich habe Überlebensdaten aus einem Experiment in Fliegen, die Raten des Alterns in verschiedenen Genotypen untersucht. Die Daten stehen mir in mehreren Layouts zur Verfügung, so dass die Auswahl Ihnen überlassen bleibt, je nachdem, welche Antwort am besten passt.

Ein Datenrahmen (wide.df) sieht so aus, wobei jeder Genotyp (Exp, von dem es ~ 640 gibt) eine Zeile hat und die Tage horizontal von Tag 4 bis Tag 98 mit Zählungen neuer Todesfälle verlaufen jeden zweiten Tag.

Exp  Day4 Day6 Day8 Day10 Day12 Day14 ... 
A  0  0  0  2  3  1  ... 

Ich mache das Beispiel dies mit:

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

Einem anderen Version wie das ist, wo jeder Tag eine Zeile für jeden ‚Exp‘ hat und die Zahl der Todesfälle an diesem Tag aufgezeichnet werden.

Exp  Deaths Day  
A  0  4  
A  0  6 
A  0  8 
A  2  10 
A  3  12 
..  ..  .. 

Um dieses Beispiel zu machen:

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

Was würde ich tun, wie eine Gompertz Analyse (See second paragraph of "the life table" here) durchzuführen. Die Gleichung lautet:

uX = α * e ^{β * x}

Wo uX zu einer gegebenen Zeit Wahrscheinlichkeit des Todes ist, α anfängliche Mortalitätsrate ist, und β ist die Rate des Alterns.

Ich möchte in der Lage sein, einen Datenrahmen zu erhalten, die später für jedes meiner ~ 640 Genotypen für die weitere Analyse α und β Schätzungen hat.

Ich brauche Hilfe von den oben genannten Datenrahmen an einen Ausgang dieser Werte für jeden meiner Genotypen in R. geht

ich durch das Paket ausgesehen haben flexsurv, die die Antwort Haus kann, aber ich habe in Versuche schlugen fehl um es zu finden und umzusetzen.

Answer 1

Dies sollte Ihnen den Einstieg ...

Zum einen für die flexsurvreg Funktion, Sie arbeiten müssen Ihre Eingangsdaten als Surv Objekt bestimmen (von package:survival). Dies bedeutet eine Reihe pro Beobachtung.

Als Erstes müssen die Rohdaten aus den von Ihnen bereitgestellten Übersichtstabellen neu erstellt werden. (Ich weiß, rbind ist nicht effizient, aber Sie können immer auf data.table für große Sätze wechseln).

### 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 

Hinweis geben, dass dies ein Intercept-Modell nur als alle Ihre Genotypen A sind. Sie können diese Schleife auf mehrere Überlebensobjekte anwenden, nachdem Sie die Daten pro Beobachtung wie oben beschrieben neu erstellt haben.

Aus dem flexsurv docs:

Gompertz Verteilung mit Formparametern a und Ratenparameter b hat Gefahrenfunktion

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

So scheint es, Ihr Alpha ist b, t er Rate und Beta ist ein, die Form.