Data Description

Diabetic Data is a dataset with 394 observations with the following 2 variables.

The 197 patients in this dataset were a 50% random sample of the patients with “high-risk” diabetic retinopathy as defined by the Diabetic Retinopathy Study (DRS). Each patient had one eye randomized to laser treatment and the other eye received no treatment. For each eye, the event of interest was the time from initiation of treatment to the time when visual acuity dropped below 5/200 two visits in a row. Thus there is a built-in lag time of approximately 6 months (visits were every 3 months). Survival times in this dataset are therefore the actual time to blindness in months, minus the minimum possible time to event (6.5 months). Censoring was caused by death, dropout, or end.

Load the required libraries

library(survival)
library(flexsurv)
library(ggplot2)
library(EnvStats)

Load the diabetic dataset and create the survival object

diabetic = diabetic[c("time", "status")]
head(diabetic)
##    time status
## 1 46.23      0
## 2 46.23      0
## 3 42.50      0
## 4 31.30      1
## 5 42.27      0
## 6 42.27      0
s <- with(diabetic,Surv(time,status))

Histogram of time variable

diabetic["status"] = as.factor(diabetic$status)
g_hist <- ggplot(diabetic, aes(x = time, color = status, fill = status))
g_hist + geom_histogram(position = "identity")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Fitting Weibull, Exponential, LogNormal and Gamma distributions

sWei  <- flexsurvreg(s ~ 1,dist='weibull',data= diabetic)
sLno  <- flexsurvreg(s ~ 1,dist='lnorm',data= diabetic)   
sExp <- flexsurvreg(s ~ 1, dist = "exponential", data = diabetic) 
sGam <- flexsurvreg(s ~ 1, dist = "gamma", data = diabetic)

Plotting S(t) for fitted distributions

Kaplan Meier curve

km = survfit(s ~ 1, data = diabetic)


kmi <- rep("KM",length(km$time))
km_df <- data.frame(km$time,km$surv,kmi)
names(km_df) <- c("Time","Surv","Model")


Weibull <- rep("Weibull",length(km$time))
Weibull_df <- data.frame(as.data.frame(summary(sWei))$time,as.data.frame(summary(sWei))$est, Weibull)
names(Weibull_df) <- c("Time","Surv","Model")

Exp <- rep("Exponential",length(km$time))
Exp_df <- data.frame(as.data.frame(summary(sExp))$time,as.data.frame(summary(sExp))$est,Exp)
names(Exp_df) <- c("Time","Surv","Model")

LNo <- rep("LogNormal",length(km$time))
LNo_df <- data.frame(as.data.frame(summary(sLno))$time,as.data.frame(summary(sLno))$est, LNo)
names(LNo_df) <- c("Time","Surv","Model")

Gamma <- rep("Gamma",length(km$time))
Gamma_df <- data.frame(as.data.frame(summary(sGam))$time,as.data.frame(summary(sGam))$est, Gamma)
names(Gamma_df) <- c("Time","Surv","Model")

plot_df <- rbind(km_df, Weibull_df, Exp_df, LNo_df, Gamma_df)

library(ggplot2)
p <- ggplot(plot_df, aes(x = Time, y = Surv, color = Model))
p + geom_line(size = 1)

QQPLOT for the fitted distributions

censored <- diabetic$status == "0"

par(mfrow = c(2, 2))

with(diabetic, qqPlotCensored(time, censored, prob.method = "kaplan-meier",
 points.col = "blue", add.line = TRUE,
 dist = "weibull", param.list = list(shape = 0.7974, scale = 109.2914), main="Q-Q Plot for\nCensored Data Set"))

with(diabetic, qqPlotCensored(time, censored, prob.method = "kaplan-meier",
 points.col = "blue", add.line = TRUE,
 dist = "lnorm", param.list = list(meanlog = 4.334, sdlog = 1.948), main="Q-Q Plot for\nCensored Data Set"))

with(diabetic, qqPlotCensored(time, censored, prob.method = "kaplan-meier",
 points.col = "blue", add.line = TRUE,
 distribution = "exp", param.list = list(rate = 0.011057), main="Q-Q Plot for\nCensored Data Set"))

with(diabetic, qqPlotCensored(time, censored, prob.method = "kaplan-meier",
 points.col = "blue", add.line = TRUE,
 distribution = "gamma", param.list = list(shape = 0.77490, scale = 0.00703), main="Q-Q Plot for\nCensored Data Set"))

Log Likelihood Test

exp_ll = sExp$loglik
wei_ll = sWei$loglik
lno_ll = sLno$loglik
gam_ll = sGam$loglik

GOF = 2 * (lno_ll - exp_ll)
GOF
## [1] 20.70881
GOF1 = 2 * (lno_ll - gam_ll)
GOF1
## [1] 11.93076
qchisq(0.95, 2)
## [1] 5.991465