Generating/plotting a log-normal survival function

Question

I have an accelerated failure time model in SAS LIFEREG that I'd like to plot. Because SAS is to profoundly bad at graphing, I'd like to actually re-generate the data for the curves in R and plot them there. SAS puts out a scale (in the case of the exponential distribution fixed to 1), an intercept, and a regression coefficient for being in the exposed or unexposed population. There's two curves, one for the exposed and one for the unexposed population. One of the models is an exponential distribution, and I've produced the data and graph like so: ``` 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,dat…

Accepted Answer

The survival function at time t for a log-normal model can be represented in R with `1 - plnorm()`, where `plnorm()` is the log-normal cumulative distribution function. To illustrate, we'll first put your plot into a function for convenience:

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

Next, we'll specify the coefficients, variables, and models, and then generate the survival probabilities for the exponential and log-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 plot the survival probabilities:

```
## 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 figures:

Note that

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

is the same as

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

which is the Φ in the canonical equation for the log-normal survival function: S(t) = 1 − Φ((ln(t) − µ) / σ)

As I'm sure you know, there are a number of R packages that can handle accelerated failure time models with left, right, or interval censoring, as listed in the survival task view, in case you happen to develop a preference for R over SAS.

Generating/plotting a log-normal survival function

Solution

Problem