Tạo/lập biểu đồ chức năng tồn tại log-normal

Question

Tôi có mô hình thời gian lỗi nhanh trong SAS LIFEREG mà tôi muốn vẽ. Bởi vì SAS là sâu sắc về đồ họa, tôi thực sự muốn tạo lại dữ liệu cho các đường cong trong R và vẽ chúng ở đó. SAS đặt ra một thang đo (trong trường hợp phân phối theo cấp số nhân cố định thành 1), một hệ thống đánh chặn và hệ số hồi quy để ở trong quần thể tiếp xúc hoặc không phơi sáng.

Có hai đường cong, một cho đường viền và một đường cho dân số chưa phơi sáng. Một trong những mô hình là một phân phối mũ, và tôi đã tạo ra các dữ liệu và đồ thị như sau:

intercept <- 5.00 
effect<- -0.500 
data<- data.frame(time=seq(0:180)-1) 
data$s_unexposed <- apply(data,1,function(row) exp(-(exp(-intercept))*row[1])) 
data$s_exposed <- apply(data,1,function(row) exp(-(exp(-(intercept+effect))*row[1]))) 

plot(data$time,data$s_unexposed, type="l", ylim=c(0,1) ,xaxt='n', 
    xlab="Days since Infection", ylab="Percent Surviving", lwd=2) 
axis(1, at=c(0, 20, 40, 60, 80, 100, 120, 140, 160, 180)) 
lines(data$time,data$s_exposed, col="red",lwd=2) 
legend("topright", c("ICU Patients", "Non-ICU Patients"), lwd=2, col=c("red","black")) 

nào mang lại cho tôi điều này:

Không đồ thị đẹp nhất bao giờ hết, nhưng tôi không thực sự biết cách của tôi xung quanh ggplot2 đủ để spruce nó lên. Nhưng quan trọng hơn, tôi có một bộ dữ liệu thứ hai xuất phát từ bản phân phối Log Normal, chứ không phải là số mũ và nỗ lực của tôi để tạo dữ liệu cho điều đó đã thất bại hoàn toàn - sự kết hợp của cdf cho phân phối bình thường và tương tự nó vượt ra ngoài kỹ năng R của tôi.

Bất kỳ ai cũng có thể chỉ cho tôi đúng hướng, sử dụng cùng số và tham số tỷ lệ 1?

Answer 1

Chức năng sống tại thời điểm t cho mô hình log-normal có thể được biểu diễn bằng R với 1 - plnorm(), trong đó plnorm() là hàm phân phối tích lũy log-normal. Để minh họa, đầu tiên chúng tôi sẽ đặt cốt truyện của bạn thành một chức năng để thuận tiện:

## Function to plot aft data 
plot.aft <- function(x, legend = c("ICU Patients", "Non-ICU Patients"), 
    xlab = "Days since Infection", ylab="Percent Surviving", lwd = 2, 
    col = c("red", "black"), at = c(0, 20, 40, 60, 80, 100, 120, 140, 160, 180), 
     ...) 
{ 
    plot(x[, 1], x[, 2], type = "l", ylim = c(0, 1), xaxt = "n", 
      xlab = xlab, ylab = ylab, col = col[2], lwd = 2, ...) 
    axis(1, at = at) 
    lines(x[, 1], x[, 3], col = col[1], lwd=2) 
    legend("topright", legend = legend, lwd = lwd, col = col) 
} 

Tiếp theo, chúng tôi sẽ xác định các hệ số, các biến, và các mô hình, và sau đó tạo ra khả năng sống sót cho mũ và đăng nhập bình thường mô hình:

## Specify coefficients, variables, and linear models 
beta0 <- 5.00 
beta1 <- -0.500 
icu <- c(0, 1) 
t <- seq(0, 180) 
linmod <- beta0 + (beta1 * icu) 
names(linmod) <- c("unexposed", "exposed") 

## Generate s(t) from exponential AFT model 
s0.exp <- dexp(exp(-linmod["unexposed"]) * t) 
s1.exp <- dexp(exp(-linmod["exposed"]) * t) 

## Generate s(t) from lognormal AFT model 
s0.lnorm <- 1 - plnorm(t, meanlog = linmod["unexposed"]) 
s1.lnorm <- 1 - plnorm(t, meanlog = linmod["exposed"]) 

Cuối cùng, chúng ta có thể vẽ đồ thị xác suất sống sót:

## Plot survival 
plot.aft(data.frame(t, s0.exp, s1.exp), main = "Exponential model") 
plot.aft(data.frame(t, s0.lnorm, s1.lnorm), main = "Log-normal model") 

Và con số kết quả:

Lưu ý rằng

plnorm(t, meanlog = linmod["exposed"]) 

cũng giống như

pnorm((log(t) - linmod["exposed"])/1) 

là Φ trong phương trình chính tắc cho chức năng sinh tồn log-bình thường: S (t) = 1 - Φ ((ln (t) - µ)/σ)

Như tôi chắc chắn bạn biết, có một số gói R có thể xử lý các mô hình thời gian lỗi tăng tốc với điều khiển trái, phải hoặc khoảng thời gian, như được liệt kê trong survival task view, trong trường hợp bạn phát triển tùy chọn cho R hơn SAS.