Generazione/tracciatura di una funzione di sopravvivenza normale del log

Question

Ho un modello di tempo di guasto accelerato in SAS LIFEREG che vorrei tracciare. Dato che SAS è profondamente pessima per il graficing, mi piacerebbe in realtà rigenerare i dati per le curve in R e tracciarli lì. SAS emette una scala (nel caso della distribuzione esponenziale fissata a 1), un'intercetta e un coefficiente di regressione per essere nella popolazione esposta o non esposta.

Ci sono due curve, una per l'esposta e una per la popolazione non esposta. Uno dei modelli è una distribuzione esponenziale, e ho prodotto i dati ed il grafico in questo modo:

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

che mi dà questo:

Non il grafico più bella di sempre, ma io non so davvero come aggirare ggplot2 abbastanza da farlo sprizzare. Ma ancora più importante, ho un secondo set di dati che proviene da una distribuzione normale di log, piuttosto che esponenziale, e i miei tentativi di generare i dati per quello sono falliti completamente - l'incorporazione del cdf per la distribuzione normale e simili mette oltre le mie abilità R.

Chiunque è in grado di indicarmi la direzione giusta, usando gli stessi numeri e un parametro di scala di 1?

Answer 1

La funzione di sopravvivenza al tempo t per un modello log-normale può essere rappresentata in R con 1 - plnorm(), dove plnorm() è la funzione di distribuzione cumulativa log-normale. Per illustrare, avremo prima messo vostro diagramma in una funzione per comodità:

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

Successivamente, specificheremo i coefficienti, variabili e modelli, e quindi generare le probabilità di sopravvivenza per l'esponenziale e log-normale modelli:

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

Infine, siamo in grado di tracciare le probabilità di sopravvivenza:

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

e le cifre risultanti:

noti che

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

è uguale

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

che è la Φ nell'equazione canonica per la funzione di sopravvivenza log-normale: S (t) = 1 - Φ ((ln (t) - μ)/σ)

Come già saprete, esistono un numero di pacchetti R in grado di gestire modelli con tempi di guasto accelerati con censura a sinistra, a destra o a intervalli, come elencato nello survival task view, nel caso in cui si sviluppi una preferenza per R su SAS.