Phân tích Agper Gompertz trong R

Question

Tôi có dữ liệu sống sót từ một thí nghiệm trên ruồi kiểm tra tốc độ lão hóa trong các kiểu gen khác nhau. Các dữ liệu có sẵn cho tôi trong một số bố cục để lựa chọn trong đó là tùy thuộc vào bạn, tùy theo điều kiện nào phù hợp nhất với câu trả lời.

Một khung dữ liệu (wide.df) trông như thế này, trong đó mỗi kiểu gen (Exp, trong đó có ~ 640) có một hàng và ngày chạy theo thứ tự theo chiều ngang từ ngày 4 đến ngày 98 với số lượng tử vong mới mỗi hai ngày.

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

tôi làm ví dụ sử dụng này:

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

phiên bản khác là như thế này, nơi mỗi ngày có hàng cho mỗi 'Exp' và số ca tử vong vào ngày hôm đó được ghi lại.

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

Để làm ví dụ này:

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

Những gì tôi muốn làm là thực hiện một Gompertz Phân tích (See second paragraph of "the life table" here). Phương trình là:

μx = α * e ^{β * x}

đâu μx là xác suất tử vong tại một thời điểm nhất định, α là tỷ lệ tử vong ban đầu, và β là tỷ lệ lão hóa.

Tôi muốn để có thể có được một dataframe trong đó có α và β ước tính cho mỗi tôi ~ 640 kiểu gen để phân tích sau đó.

Tôi cần giúp đi từ dataframes trên để một đầu ra của các giá trị cho mỗi kiểu gen của tôi trong R.

tôi đã xem xét thông qua các gói flexsurv có thể chứa câu trả lời nhưng tôi đã thất bại trong nỗ lực để tìm và thực hiện nó.

Answer 1

này sẽ giúp bạn bắt đầu ...

Thứ nhất, cho flexsurvreg chức năng để làm việc, bạn cần phải xác định dữ liệu đầu vào của bạn như là một đối tượng Surv (từ package:survival). Điều này có nghĩa là một hàng cho mỗi quan sát.

Điều đầu tiên là tạo lại dữ liệu 'thô' từ bảng tóm tắt mà bạn cung cấp. (Tôi biết rbind không hiệu quả, nhưng bạn luôn có thể chuyển sang data.table cho các bộ lớn).

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

cho

> f1$res 
       est   L95%  U95% 
shape 0.165351912 0.1281016481 0.202602176 
rate 0.001767956 0.0006902161 0.004528537 

Lưu ý rằng đây là một mô hình đánh chặn chỉ như tất cả các kiểu gen của bạn là A. Bạn có thể lặp vòng điều này qua nhiều đối tượng tồn tại khi bạn đã tạo lại dữ liệu mỗi quan sát như trên.

Từ flexsurv tài liệu:

phân phối với hình dạng tham số một tham số và tỷ lệ Gompertz b có chức năng nguy hiểm

H (x: a, b) = được^{rìu }

Vì vậy, có vẻ như alpha của bạn là b, t anh ấy xếp hạng và beta là a, hình dạng.