Generar/trazar una función de supervivencia de log-normal

Question

Tengo un modelo de tiempo de falla acelerado en SAS LIFEREG que me gustaría trazar. Como SAS es profundamente malo en la creación de gráficos, me gustaría volver a generar los datos para las curvas en R y trazarlos allí. SAS saca una escala (en el caso de la distribución exponencial fijada a 1), una intersección y un coeficiente de regresión para estar en la población expuesta o no expuesta.

Hay dos curvas, una para la población expuesta y otra para la población no expuesta. Uno de los modelos es una distribución exponencial, y yo he producido los datos y el gráfico de esta manera:

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

que me da esto:

No es la más bonita gráfica nunca, pero yo realmente no conozco mi camino alrededor de ggplot2 lo suficiente como para arreglarlo. Pero lo más importante es que tengo un segundo conjunto de datos que proviene de una distribución Log Normal, en lugar de una exponencial, y mis intentos de generar los datos para eso han fallado por completo: la incorporación del cdf para la distribución normal y similares. más allá de mis habilidades de R.

¿Alguien capaz de señalarme en la dirección correcta, usando los mismos números, y un parámetro de escala de 1?

Answer 1

La función de supervivencia en el tiempo t para un modelo log-normal se puede representar en R con 1 - plnorm(), donde plnorm() es la función de distribución acumulativa logarítmica normal. Para ilustrar esto, vamos a poner su primera trama en una función por conveniencia:

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

A continuación, especificaremos los coeficientes, variables y modelos, y luego generar las probabilidades de supervivencia para la exponencial y log-normal modelos:

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

Por último, podemos trazar las probabilidades de supervivencia:

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

y las cifras resultantes:

Nota que

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

es la misma que

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

que es el Φ en la ecuación canónica para la función de supervivencia log-normal: S (t) = 1 - Φ ((ln (t) - μ)/σ)

Como estoy seguro, hay una cantidad de paquetes R que pueden manejar modelos de tiempo de falla acelerada con censor de izquierda, derecha o intervalo, como se detalla en el survival task view, en caso de que desarrolle una preferencia por R sobre SAS.