Generation / construction of a log-normal survival function

Question

I have an accelerated failure time model in SAS LIFEREG that I would like to build. Since SAS has a deep bad line when plotting, I would like to actually generate data for the curves in R and build them there. SAS sets the scale (in the case of an exponential distribution fixed at 1), interception and regression coefficient to be in the exposed or untreated population.

There are two curves: one for the exposed and one for the unexposed population. One of the models is the exponential distribution, and I created the data and the graph as follows:

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

What gives me this:

Not the most beautiful graph in history, but I really don’t know my way around ggplot2 to decorate it. But more importantly, I have a second dataset that comes from the log-normal distribution, and not exponential, and my attempts to generate data for this completely failed - turning on cdf for normal distribution, etc. it's above my R. skills

Anyone who can point me in the right direction using the same numbers and scale option 1?

Answer 1

The survival function at time t for the log-normal model can be represented in R using 1 - plnorm() , where plnorm() is the logarithmically normal cumulative distribution function. To illustrate this, we first put your graph in a function:

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

Then we specify the coefficients, variables and models, and then we generate the survival probabilities for exponential and logarithmically normal models:

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

Finally, we can build the probability of survival:

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

And the resulting numbers:

note that

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

coincides with

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

which is Φ in the canonical equation for the logarithmic survival function: S (t) = 1 - Φ ((ln (t) - μ) / σ)

As I’m sure, you know, there are many R packages that can process models with accelerated failure times with left, right or interval censorship, listed in the form of survival , in case you have to develop an R preference over SAS.