CODEWITHSUNDEEP
rdistributionmeansurvival-analysis
user2781406
user2781406

Reputation: 31

Finding the mean of the log-normal distribution in survival analysis in R

I am a novice with R. Currently I am fitting a log-normal distribution to some survival data I have, however I have become stuck when trying to calculate statistics such as the median and the mean. This is the code I have used so far, can anyone tell me what I should type next to find the mean?

# rm(list=ls(all=TRUE))
library(survival)
data<-read.table("M:\\w2k\\Diss\\Hoyle And Henley True IPD with number at risk known.txt",header=T)
attach(data)
data
times_start <-c(  rep(start_time_censor, n_censors), rep(start_time_event, n_events) )
times_end <-c(  rep(end_time_censor, n_censors), rep(end_time_event, n_events)  )
model <- survreg(Surv(times_start, times_end, type="interval2")~1, dist="lognormal")
intercept <- summary(model)$table[1]   
log_scale <- summary(model)$table[2]

this is where I got stuck, I have tried:

meantime<-exp(intercept+log_scale/2)

but this does not seem to give a realistic mean.

Upvotes: 2

Views: 2105

Answers (1)

IRTFM
IRTFM

Reputation: 263332

The place to look for a worked example is ?predict.survreg. (In general, using the help system for predict methods is a productive strategy for any regression method.)

Running the last example should give you enough basis to proceed. In particular you should see that the regression coefficients are not estimates of survival times or quantiles.

lfit <- survreg(Surv(time, status) ~ ph.ecog, data=lung)
pct <- 1:98/100   # The 100th percentile of predicted survival is at +infinity
ptime <- predict(lfit, newdata=data.frame(ph.ecog=2), type='quantile',
                 p=pct, se=TRUE)
matplot(cbind(ptime$fit, ptime$fit + 2*ptime$se.fit,
                          ptime$fit - 2*ptime$se.fit)/30.5, 1-pct,
         xlab="Months", ylab="Survival", type='l', lty=c(1,2,2), col=1)
 # The plot should be examined since you asked for a median survival time
 abline(h= 0.5)
 # You can  drop a vertical from the intersection to get that graphically

.... or ... 

 str(ptime)
List of 2
 $ fit   : num [1:98] 9.77 16.35 22.13 27.46 32.49 ...
 $ se.fit: num [1:98] 2.39 3.53 4.42 5.16 5.82 ...

You can extract the 50th percentile from that sequence of survival times with:

 ptime$fit[which((1-pct)==0.5)]
# [1] 221.6023

Measured in days which was why Therneau divided by 30.5 to display months

Upvotes: 2

Related Questions