Take-Home Midterm 3

library(survival)
library(flexsurv)

## Warning: package 'flexsurv' was built under R version 4.4.3

---

Setup Code

Kaplan-Meier

survKM2 <- function(dat, cen=NA, conf=0.95, eps=1e-9) {
   n <- length(dat); if (is.na(cen[1])) cen <- rep(T, n); cen <- as.logical(cen)
   di <- table(dat[cen])
   ti <- as.numeric(names(di)) - eps
    m <- length(ti)
   ni <- NULL; for(i in 1:m) ni[i] <- sum(dat >= ti[i])
    S <- cumprod((ni-di)/ni)
    V <- S^2*cumsum(di/ni/(ni-di))
   ti <- c(0, ti); di <- c(0, as.vector(di)); ni <- c(n, ni); S <- c(1, S); V <- c(0, V)
   se <- sqrt(V); z <- qnorm(1-(1-conf)/2)
   CI.lb <- S - z*se; CI.lb <- CI.lb - (CI.lb < 0)* CI.lb
   CI.ub <- S + z*se; CI.ub <- CI.ub - (CI.ub > 1)*(CI.ub-1)

     h <- di/ni
     H <- cumsum(h)
   V.H <- cumsum(di/ni^2)
   S.H <- exp(-H)

   cat("Kaplan-Meier estimate, a.k.a. the product-limit estimate of S(t) \n")
   cat("Nonparametric estimate of survival function S(t) \n\n")   
   tab <- data.frame(time=ti, di=di, ni=ni, S=round(S,3), Var.S=round(V,3), se=round(se,3),
                     CI.lb=round(CI.lb,3), CI.ub=round(CI.ub,3), XXXXX="",
                     h=round(h,3), H=round(H,3), Var.H=round(V.H,3), S.na=round(S.H,3))
   print(tab, row.names=F)
   cat("\nh is nonparametric estimate of the hazard function h(t).")
   cat("\nH is Nelson-Aallen estimate of the cumulative hazard function H(t).")
   cat("\nS.na is Nelson-Aallen estimate of S(t). \n\n")

   par(mar=c(5,5,4,1))
   plot(c(ti, 1.2*max(ti)), c(S, min(S)), type="s", lwd=2, 
      xlab="Time, t", ylab=expression("Estimated Survival Probability,   "*hat(S)(t)),
      main="Kaplan-Meier estimate, a.k.a. the product-limit estimate of S(t)")   
   points(ti[-1], S[-1], pch=19)
   points(ti    , CI.lb, type="s", lty=2)
   points(ti    , CI.ub, type="s", lty=2)

    m <- m + 1
   St <- c(S[-m]*diff(ti), 0)
   mu <- sum(St)
   Ai <- rev(cumsum(rev(St)))   # = mu - cumsum(c(0, St[-m]))
   Vm <- sum((Ai^2*di/ni/(ni-di))[-m])
   se <- sqrt(Vm)

   # if (!cen[which.max(dat)]), could replace w/ limit L 
   # for Mean survival time limited to a time L.

   cat("\nMean survival time based on Kaplan-Meier estimate of S(t) \n\n")   
   print(data.frame(Estimate=round(mu,3), Var=round(Vm,3), se=round(se,3),
      CI.lb=round(mu-z*se,3), CI.ub=round(mu+z*se,3)), row.names=F)

   invisible(tab) }

Res

quiet <- function(fn) {   # suppress output messages
    sink(tempfile())
    on.exit(sink())
    invisible(force(fn)) }


###########################
plotR <- function(r, cen) {
   cen <- as.logical(cen)
   r[!cen] <- r[!cen] + 1   # or add 0.693
   KMobj <- quiet(survKM2(r, cen))
     r <- KMobj$time
    Sr <- KMobj$S
   fit <- lm(-log(Sr) ~ r, subset=(Sr != 0))
   tab <- list(residual=r, Sr=Sr, R2=summary(fit)$r.sq)

   plot(r, -log(Sr), xlab="Residual, r", ylab=expression(-log(hat(S)(r))), 
        main="Cox-Snell residual plot with uncensored/censored data")   
   abline(fit, lty=2, col="red")
   abline(a=0, b=1, lwd=2) 

   invisible(tab) }

MLE

nlL <- function(par, dist=c("exponential", "exp2par", "weibull", "lognormal", "gamma", "gengamma", "loglogistic", "gompertz"), 
                dat, cen=NA) {
   n <- length(dat); if (is.na(cen[1])) cen <- rep(T, n); cen <- as.logical(cen)

   if (dist[1]=="exponential") { 
      lambda <- par; if (!(lambda > 0)) return(Inf)
      pdf <- function(tt)   dexp(tt, rate=lambda)
       sf <- function(tt) 1-pexp(tt, rate=lambda) }
   if (dist[1]=="exp2par") {       
      lambda <- par[1]; if (!(lambda > 0   )) return(Inf)
           G <- par[2]; if (!(G <= min(dat))) return(Inf)
      pdf <- function(tt)   dexp(tt-G, rate=lambda)
       sf <- function(tt) 1-pexp(tt-G, rate=lambda) }
   if (dist[1]=="weibull") { 
      if (!all(par>0)) return(Inf)
      lambda <- par[1]; gam <- par[2]
      pdf <- function(tt)   dweibull(tt, shape=gam, scale=1/lambda)
       sf <- function(tt) 1-pweibull(tt, shape=gam, scale=1/lambda) }
   if (dist[1]=="lognormal") {       
         mu <- par[1]
      sigma <- par[2]; if (!(sigma > 0)) return(Inf)
      pdf <- function(tt)   dlnorm(tt, meanlog=mu, sdlog=sigma)
       sf <- function(tt) 1-plnorm(tt, meanlog=mu, sdlog=sigma) }
   if (dist[1]=="gamma") { 
      if (!all(par>0)) return(Inf)
      lambda <- par[1]; gam <- par[2]
      pdf <- function(tt)   dgamma(tt, shape=gam, scale=1/lambda)
       sf <- function(tt) 1-pgamma(tt, shape=gam, scale=1/lambda) }
   if (dist[1]=="gengamma") { 
      if (!all(par[-2]>0)) return(Inf)
      lambda <- par[1]; alpha <- par[2]; gam <- par[3]
      pdf <- function(tt) if (alpha < 0) 
           dgengamma(tt, mu=-log(lambda), sigma=-1/alpha/sqrt(gam), Q=-1/sqrt(gam)) else
           dgengamma(tt, mu=-log(lambda), sigma= 1/alpha/sqrt(gam), Q= 1/sqrt(gam)) 
       sf <- function(tt) if (alpha < 0) 
         1-pgengamma(tt, mu=-log(lambda), sigma=-1/alpha/sqrt(gam), Q=-1/sqrt(gam)) else
         1-pgengamma(tt, mu=-log(lambda), sigma= 1/alpha/sqrt(gam), Q= 1/sqrt(gam)) }
   if (dist[1]=="loglogistic") { 
      if (!all(par>0)) return(Inf)
      alpha <- par[1]; gam <- par[2]
      pdf <- function(tt)   dllogis(tt, shape=gam, scale=1/alpha^(1/gam))
       sf <- function(tt) 1-pllogis(tt, shape=gam, scale=1/alpha^(1/gam)) }
   if (dist[1]=="gompertz") { 
      if (!all(par>0)) return(Inf)
      lambda <- par[1]; gam <- par[2]
      pdf <- function(tt)   dgompertz(tt, shape=gam, rate=exp(lambda))
       sf <- function(tt) 1-pgompertz(tt, shape=gam, rate=exp(lambda)) }

   logL <- sum(log(pdf(dat[cen]))) + sum(log(sf(dat[!cen])))
   tt <- sort(dat)
   attr(logL, "H") <- -log(sf(tt))
   -logL }


###########################################################
MLE <- function(par0, dist, dat, cen=NA, se=T, conf=0.95) {
   suppressWarnings(
   MLobj <- optim(par0, nlL, dist=dist, dat=dat, cen=cen, hessian=se))
   mle <- MLobj$par
     V <- try(solve(MLobj$hessian), silent=T)
   if (se) se <- sqrt(diag(V)) else se <- NA
   z <- qnorm(1-(1-conf)/2)
   CI.lb <- mle - z*se
   CI.ub <- mle + z*se   # approximate confidence interval
 
   par <- switch(dist, 
      exponential=c("lambda"         ),
      exp2par    =c("lambda", "G"    ),
      weibull    =c("lambda", "gamma"),
      lognormal  =c("mu"    , "sigma"),
      gamma      =c("lambda", "gamma"),
      gengamma   =c("lambda", "alpha" , "gamma"),
      loglogistic=c("alpha" , "gamma"),
      gompertz   =c("lambda", "gamma"))
   tab <- data.frame(parameter=par, MLE=round(mle,3), se=round(se,3), 
                     CI.lb=round(CI.lb,3), CI.ub=round(CI.ub,3))
   logL <- -nlL(mle, dist, dat, cen)   # = -MLobj$val
   attr(tab, "max.logL") <- logL
   attr(tab, "residual") <- attr(logL, "H")   # Cox-Snell residuals 
   tab }


###############################################################
MLEexp <- function(dat, cen=NA, shift=F, conf=0.95, eps=1e-9) {
   n <- length(dat); if (is.na(cen[1])) cen <- rep(T, n); cen <- as.logical(cen)
   r <- n-sum(!cen) - shift; nc <- length(unique(dat[!cen]))
   dat[!cen] <- dat[!cen] + eps

   t1 <- min(dat)
   mle <- r/sum(dat - t1*shift); mle[2] <- mu <- 1/mle
   se <- mle/sqrt(r)
   z <- qnorm(1-(1-conf)/2)
   CI.lb <- mle - z*se
   CI.ub <- mle + z*se   # approximate confidence interval

   tab <- data.frame(parameter=c("lambda", "mu"), MLE=round(mle,3), se=round(se,3), 
                     approxCI.lb=round(CI.lb,3), approxCI.ub=round(CI.ub,3))

   if (shift) {   # dist="exp2par"
       G <- t1 - mu/n
      se <- mu/n*(1+1/r)
      CI.lb <- G - z*se
      CI.ub <- G + z*se   # approximate confidence interval

      tab <- rbind(tab, data.frame(parameter="G", MLE=round(G,3), se=round(se,3), 
                          approxCI.lb=round(CI.lb,3), approxCI.ub=round(CI.ub,3))) }

   if ((nc==0) || ((nc==1) && (!cen[which.max(dat)]))) {
      chiL <- qchisq(  (1-conf)/2, df=2*r)
      chiU <- qchisq(1-(1-conf)/2, df=2*r)
      CI <- 2*r*mu/c(chiU, chiL)
      CI <- rbind(rev(1/CI), CI)   # exact confidence interval

      if (shift) {   # dist=exp2par
         fU <- qf(conf, df1=2, df2=2*r)   # approx r*(1/(1-conf)^(1/r) - 1)
         CI <- rbind(CI, c(t1 - fU*mu/n, t1)) }

      colnames(CI) <- c("exactCI.lb", "exactCI.ub"); rownames(CI) <- NULL      
      tab <- cbind(tab, round(CI,3)) }

   if (shift) {
      cat("Note: For two-parameter exponential distribution (with a shift parameter), \n")
      cat("      the minimum variance unbiased estimates are provided rather than MLE. \n\n") }
   tab }


######################################
MLElnorm <- function(dat, conf=0.95) {
   n <- length(dat)
   mu <- mean(log(dat))
   s2 <-  var(log(dat)); sigma <- sqrt(s2*(n-1)/n)
   mle <- c(mu, sigma)
   se <- c(sigma/sqrt(n), sigma*sqrt(sqrt(2*(n-1))/n))
      tq <- qt(1-(1-conf)/2, df=n-1)
   CI.lb <- mu - tq*sqrt(s2/(n-1))
   CI.ub <- mu + tq*sqrt(s2/(n-1))         # exact confidence interval for mu
    chiL <- qchisq(  (1-conf)/2, df=n-1)
    chiU <- qchisq(1-(1-conf)/2, df=n-1)
   CI.lb <- c(CI.lb, sigma*sqrt(n/chiU))
   CI.ub <- c(CI.ub, sigma*sqrt(n/chiL))   # exact confidence interval for sigma

   tab <- data.frame(parameter=c("mu", "sigma"), MLE=round(mle,3), se=round(se,3), 
                     exactCI.lb=round(CI.lb,3), exactCI.ub=round(CI.ub,3))
   tab }

Plots

plotP <- function(dist=c("exponential", "weibull", "lognormal", "loglogistic"), dat) {
   n <- length(dat)
   mi <- table(dat)
   ti <- as.numeric(names(mi)); idx <- order(ti); ti <- ti[idx]; mi <- as.vector(mi[idx])
    F <- (cumsum(mi)-0.5)/n

   if (dist[1]=="exponential") { 
      x <- log(1/(1-F)); xlab <- expression(-log(1-hat(F)(t)))
      y <- ti          ; ylab <- "t"
      fit <- lm(y ~ -1 + x)
      tab <- data.frame(lambda=1/fit$coef) }
   if (dist[1]=="weibull") { 
      x <- log(log(1/(1-F))); xlab <- expression(log(-log(1-hat(F)(t))))
      y <- log(ti)          ; ylab <- "log(t)"
      fit <- lm(y ~ x)
      tab <- data.frame(lambda=exp(-fit$coef[1]), gam=1/fit$coef[2]) }
   if (dist[1]=="lognormal") {       
      x <- qnorm(F); xlab <- expression(Phi^-1*(hat(F)(t)))
      y <- log(ti) ; ylab <- "log(t)"
      fit <- lm(y ~ x)
      tab <- data.frame(mu=fit$coef[1], sigma=fit$coef[2]) }
   if (dist[1]=="loglogistic") { 
      x <- log(1/(1-F)-1); xlab <- expression(-log((1-hat(F)(t))/hat(F)(t)))
      y <- log(ti)       ; ylab <- "log(t)"
      fit <- lm(y ~ x)
      tab <- data.frame(alpha=exp(-fit$coef[1]/fit$coef[2]), gam=1/fit$coef[2]) }

   plot(x, y, xlab=xlab, ylab=ylab, main="Probability plot with uncensored data")   
   abline(fit, lwd=2)
   tab$R2 <- summary(fit)$r.sq; rownames(tab) <- NULL
   round(tab,3) }
 

##############################################################################################
plotH <- function(dist=c("exponential", "weibull", "lognormal", "loglogistic"), dat, cen=NA) {
   n <- length(dat); if (is.na(cen[1])) cen <- rep(T, n); cen <- as.logical(cen)
   idx <- order(dat)
   cen <- cen[idx]
   dat <- dat[idx]; dat <- dat[cen]

   ti <- unique(dat)
   h <- 1/(n:1); h.old <- h[cen]
   h <- NULL; for(tt in ti) h <- c(h, sum(h.old[tt==dat])) 
   H <- cumsum(h)

   if (dist[1]=="exponential") { 
      x <- H ; xlab <- expression(hat(H)(t))
      y <- ti; ylab <- "t"
      fit <- lm(y ~ -1 + x)
      tab <- data.frame(lambda=1/fit$coef) }
   if (dist[1]=="weibull") { 
      x <- log(H) ; xlab <- expression(log(hat(H)(t)))
      y <- log(ti); ylab <- "log(t)"
      fit <- lm(y ~ x)
      tab <- data.frame(lambda=exp(-fit$coef[1]), gam=1/fit$coef[2]) }
   if (dist[1]=="lognormal") {       
      x <- qnorm(1-exp(-H)); xlab <- expression(Phi^-1*(1-e^-hat(H)(t)))
      y <- log(ti)         ; ylab <- "log(t)"
      fit <- lm(y ~ x)
      tab <- data.frame(mu=fit$coef[1], sigma=fit$coef[2]) }
   if (dist[1]=="loglogistic") { 
      x <- log(exp(H)-1); xlab <- expression(log(e^hat(H)(t)*-1))
      y <- log(ti)      ; ylab <- "log(t)"
      fit <- lm(y ~ x)
      tab <- data.frame(alpha=exp(-fit$coef[1]/fit$coef[2]), gam=1/fit$coef[2]) }

   plot(x, y, xlab=xlab, ylab=ylab, main="Hazard plot with uncensored/censored data")   
   abline(fit, lwd=2)
   tab$R2 <- summary(fit)$r.sq; rownames(tab) <- NULL
   round(tab,3) }

Tests

sigcode <- function(p.val) ifelse(p.val < 0.001, "***", ifelse(p.val < 0.01, "**", ifelse(p.val < 0.05, "*", ifelse(p.val < 0.1, ".", ""))))

quiet <- function(fn) {   # suppress output messages
    sink(tempfile())
    on.exit(sink())
    invisible(force(fn)) }


##############################################
LR.test <- function(par0, dist, dat, cen=NA) {
    res <- MLE(par0, dist, dat, cen, se=F)
   logL  <- attr(res, "max.logL")
   logL0 <- -nlL(par0, dist, dat, cen) 
     LRT <- -2*(logL0 - logL)
   p.val <- 1-pchisq(LRT, df=length(par0))

   cat("Likelihood ratio test for a specific distribution with known parameters \n")
   cat("   H0: S = S0 \n")
   cat(paste("       or equivalently, the underlying distribution is", dist, 
             "with specified parameter values. \n\n"))
   res <- res[,1:2]; res$null.value <- par0
   print(res, row.names=F); cat("\n")
   print(data.frame(logL=round(logL,3), logL0=round(logL0,3), LRT=round(LRT,3), 
                    p.value=round(p.val,4), signif=sigcode(p.val)), row.names=F) }


#################################################
HP.test <- function(sf0, dat, cen=NA, eps=1e-9) {
   n <- length(dat); if (is.na(cen[1])) cen <- rep(T, n); cen <- as.logical(cen)
   dat[!cen] <- dat[!cen] + eps
   idx <- order(dat)
   dat <- dat[idx]
   del <- cen <- cen[idx]

   ft <- 1/n*c(1, cumprod(((n:2)/((n-1):1))^(1-del[-n])))   # jump of Kaplan-Meier estimate
   S0 <- sf0(dat)
   s2 <- 1/16*sum(n/(n:1)*(c(1, S0[-n]^4) - S0^4))
   C <- sum(del*S0*ft)
   Z <- (C - 0.5)/sqrt(s2/n)
   p.val <- c(2*pnorm(-abs(Z)), 1-pnorm(Z), pnorm(Z))

   cat("Hollander and Proschan's test for a specific distribution with known parameters \n")
   cat("   H0: S = S0 \n")
   cat("       where S0 is the underlying survival function with specified parameter values \n\n")
   print(data.frame(C=round(C,3), sigma2=round(s2,3), C.star=round(Z,3)), row.names=F); cat("\n")
   print(data.frame(H1=c("S != S0", "S < S0", "S > S0"), 
      p.value=round(p.val,4), signif=sapply(p.val, sigcode)), row.names=F)

   quiet(survKM2(dat, cen))
   tt <- seq(0, 1.2*max(dat[cen]), len=100) 
   S0 <- sf0(tt)
   points(tt, S0, type="l", lwd=2, col="red") }   # overlay hypothesized survival function S0(t)

---

Problem 1

The following data are on remission times, in weeks, for a group of 30 leukemia patients in a certain type of therapy.

            Remission Times (in weeks)
 1  1  2  4  4  6  6   6  7   8   9  9 10  12  13
14 18 19 24 26 29 31+ 42 45+ 50+ 57 60 71+ 85+ 91

Assume that one-parameter exponential model holds for this data.

(a) Find the MLE of the rate parameter λ.

dat <- c(1, 1, 2, 4, 4, 6, 6, 6, 7, 8, 9, 9, 10, 12, 13,
         14, 18, 19, 24, 26, 29, 31, 42, 45, 50, 57, 60, 71, 85, 91)

cen <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
         1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1)   # 1 if observed, 0 if censored

res <- MLEexp(dat, cen)

print(res)

##   parameter    MLE    se approxCI.lb approxCI.ub
## 1    lambda  0.033 0.007       0.020       0.046
## 2        mu 30.400 6.080      18.483      42.317

So, we get a lambda of 0.033 as shown above.

(b) Find the MLE of the median remission time and S(26) = P(T > 26), where T is the remission time (in weeks).

lambda <- res$MLE[1]

# Median estimate
median_time <- log(2) / lambda

# Survival probability at 26 weeks
S26 <- 1-pexp(26, rate=lambda);

# Output the results
cat("Estimated median remission time:", round(median_time, 2), "weeks\n")

## Estimated median remission time: 21 weeks

cat("Estimated survival probability at 26 weeks:", round(S26, 4), "\n")

## Estimated survival probability at 26 weeks: 0.424

We got an estimated median remission time of 21 weeks as provided by the code, while the MLE of S(26) = P(T > 26) is 0.424 as shown above.

(c) Find a 95% confidence interval for λ.

r <- sum(cen)  # number of observed (uncensored) events
mu_hat <- 1 / res$MLE[1]  # estimated mean

chiL <- qchisq(0.025, df = 2 * r)
chiU <- qchisq(0.975, df = 2 * r)

CI_mu <- 2 * r * mu_hat / c(chiU, chiL)  # CI for μ
CI_lambda <- rev(1 / CI_mu)              # CI for λ (inverted)

cat("Manual exact 95% CI for lambda:", round(CI_lambda, 4), "\n")

## Manual exact 95% CI for lambda: 0.0214 0.0471

We got a 95% confidence interval of (0.0214, 0.0471).

(d) Find a 95% confidence interval for the mean parameter µ.

r <- sum(cen)  # number of uncensored events
mu_hat <- 1 / res$MLE[1]  # estimated mean

chiL <- qchisq(0.025, df = 2 * r)
chiU <- qchisq(0.975, df = 2 * r)

CI_mu <- 2 * r * mu_hat / c(chiU, chiL)  # CI for μ

cat("Exact 95% CI for mu:", round(CI_mu, 2), "\n")

## Exact 95% CI for mu: 21.21 46.83

The exact 95% confidence interval for the mean remission time is (21.21, 46.83) weeks.

---

Problem 2

Refer to the data in Question 1.

(a) Compute and plot the Kaplan-Meier estimate of the survival function S(t).

dat <- c(1, 1, 2, 4, 4, 6, 6, 6, 7, 8, 9, 9, 10, 12, 13,
         14, 18, 19, 24, 26, 29, 31, 42, 45, 50, 57, 60, 71, 85, 91)

cen <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
         1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1)   # 1 if observed, 0 if censored

km_result <- survKM2(dat, cen)

## Kaplan-Meier estimate, a.k.a. the product-limit estimate of S(t) 
## Nonparametric estimate of survival function S(t) 
## 
##  time di ni     S Var.S    se CI.lb CI.ub XXXXX     h     H Var.H  S.na
##     0  0 30 1.000 0.000 0.000 1.000 1.000       0.000 0.000 0.000 1.000
##     1  2 30 0.933 0.002 0.046 0.844 1.000       0.067 0.067 0.002 0.936
##     2  1 28 0.900 0.003 0.055 0.793 1.000       0.036 0.102 0.003 0.903
##     4  2 27 0.833 0.005 0.068 0.700 0.967       0.074 0.176 0.006 0.838
##     6  3 25 0.733 0.007 0.081 0.575 0.892       0.120 0.296 0.011 0.743
##     7  1 22 0.700 0.007 0.084 0.536 0.864       0.045 0.342 0.013 0.710
##     8  1 21 0.667 0.007 0.086 0.498 0.835       0.048 0.390 0.015 0.677
##     9  2 20 0.600 0.008 0.089 0.425 0.775       0.100 0.490 0.020 0.613
##    10  1 18 0.567 0.008 0.090 0.389 0.744       0.056 0.545 0.023 0.580
##    12  1 17 0.533 0.008 0.091 0.355 0.712       0.059 0.604 0.027 0.547
##    13  1 16 0.500 0.008 0.091 0.321 0.679       0.062 0.666 0.031 0.514
##    14  1 15 0.467 0.008 0.091 0.288 0.645       0.067 0.733 0.035 0.480
##    18  1 14 0.433 0.008 0.090 0.256 0.611       0.071 0.805 0.040 0.447
##    19  1 13 0.400 0.008 0.089 0.225 0.575       0.077 0.881 0.046 0.414
##    24  1 12 0.367 0.008 0.088 0.194 0.539       0.083 0.965 0.053 0.381
##    26  1 11 0.333 0.007 0.086 0.165 0.502       0.091 1.056 0.062 0.348
##    29  1 10 0.300 0.007 0.084 0.136 0.464       0.100 1.156 0.072 0.315
##    42  1  8 0.262 0.007 0.081 0.103 0.422       0.125 1.281 0.087 0.278
##    57  1  5 0.210 0.006 0.080 0.053 0.367       0.200 1.481 0.127 0.227
##    60  1  4 0.158 0.006 0.075 0.010 0.305       0.250 1.731 0.190 0.177
##    91  1  1 0.000   NaN   NaN   NaN   NaN       1.000 2.731 1.190 0.065
## 
## h is nonparametric estimate of the hazard function h(t).
## H is Nelson-Aallen estimate of the cumulative hazard function H(t).
## S.na is Nelson-Aallen estimate of S(t).

## 
## Mean survival time based on Kaplan-Meier estimate of S(t) 
## 
##  Estimate   Var    se  CI.lb  CI.ub
##     29.65 34.05 5.835 18.213 41.087

plot(km_result$time, km_result$S, type = "s",
     xlab = "Time (weeks)", ylab = "Survival Probability",
     main = "Kaplan-Meier Estimate of S(t)")

(b) Obtain an estimate of the median remission time. Compare it with your answer obtained in Question 1(b).

median_KM <- km_result$time[min(which(km_result$S <= 0.5))]
cat("Kaplan-Meier Median Time:", median_KM, "weeks\n")

## Kaplan-Meier Median Time: 13 weeks

The median remission time estimated by the Kaplan-Meier method is 13 weeks, which is lower than the 21-week median estimated under the exponential model. This indicates that the exponential assumption may overestimate survival time, especially early in the follow-up period.

(c) Compute and plot the nonparametric estimate of the hazard function h(t) against time. Does it support the exponential assumption?

S <- km_result$S
time <- km_result$time
haz <- -diff(S) / head(S, -1)   # hazard increments
t_mid <- (head(time, -1) + tail(time, -1)) / 2  # midpoints for plotting

# Plot h(t)
plot(t_mid, haz, type = "b", pch = 16,
     xlab = "Time (weeks)", ylab = "Estimated Hazard h(t)",
     main = "Nonparametric Estimate of the Hazard Function")

The nonparametric hazard estimate is relatively flat early on but increases noticeably after 40 weeks, especially beyond 60 weeks. This suggests the hazard is not constant over time, indicating that the exponential model may not be appropriate across the full time range.

(d) Plot the MLE of the hazard rate obtained in Question 1(a) over the plot in (c). Comment on the assumption of exponentiality.

# Add the constant MLE hazard line (from Problem 1a)
plot(t_mid, haz, type = "b", pch = 16,
     xlab = "Time (weeks)", ylab = "Estimated Hazard h(t)",
     main = "Nonparametric Hazard Function with Exponential MLE")
abline(h = 0.033, col = "red", lwd = 2, lty = 2)
legend("topright", legend = c("Nonparametric", "MLE (Exponential)"),
       col = c("black", "red"), lty = c(1, 2), lwd = c(1, 2), pch = c(16, NA))

The MLE hazard rate under the exponential model is constant at 𝜆= 0.033, shown as a red dashed line. The nonparametric hazard rises over time, indicating that the assumption of constant hazard is violated. This suggests that the exponential model may not provide an adequate fit, particularly at later times.

---

Problem 3

Refer to the data in Question 1. Now assume that two-parameter Weibull distribution is appropriate for the data.

(a) Obtain the MLE of two parameters.

dat <- c(1, 1, 2, 4, 4, 6, 6, 6, 7, 8, 9, 9, 10, 12, 13,
         14, 18, 19, 24, 26, 29, 31, 42, 45, 50, 57, 60, 71, 85, 91)
cen <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
         1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1)   # 1 if observed, 0 if censored

MLE(c(0.04,0.831), "weibull", dat, cen);

##   parameter   MLE    se CI.lb CI.ub
## 1    lambda 0.034 0.008 0.018 0.051
## 2     gamma 0.835 0.135 0.570 1.099

We are given 0.034 for the scale parameter, and 0.835 for the shape parameter.

(b) Obtain an estimate of S(26), and compare it with your answer obtained in Question 1(b). Comment.

1-pweibull(26, shape=0.835, scale=1/0.034)

## [1] 0.405689

We get a Weibull estimate of 0.405689 for S(26), which is slightly below the exponential and KM values (~0.424). This suggests a modestly quicker early decline, but the difference is minimal—indicating a comparable fit at 26 weeks.

(c) Make an appropriate plot to see if the data follows a Weibull distribution.

plotH("weibull", dat, cen)

##   lambda   gam    R2
## 1   0.04 0.831 0.978

The log-log hazard plot shows that the points align closely with the reference line, suggesting that the Weibull distribution provides a reasonable fit for the data.

---

Problem 4

7.2 In order to computerize patients’ records, a data clerk is hired to transcribe medical data from the patients’ charts to computer coding forms. The number of correct entries between errors is listed in chronological order of occurrence over a period of five days as follows: 73, 12, 40, 65, 100, 15, 70, 40, 110, 64, 200, 6, 90, 102, 20, 102, 90, and 34. Assuming that the correct entries between errors follow the two-parameter exponential distribution, obtain:

(a) An estimate of G.

dat <- c(73, 12, 40, 65, 100, 15, 70, 40, 110, 64, 200, 6, 90, 102, 20, 102, 90, 34)
res <- MLE(c(1,1), "exp2par", dat, se=F)

G <- res$MLE[2]
cat("Estimate of G:", G, "\n")

## Estimate of G: 6

The estimate of G is 6, which corresponds to the minimum value in the dataset.

(b) The MLE of l.

lambda <- res$MLE[1]
cat("MLE of lambda:", lambda, "\n")

## MLE of lambda: 0.016

The MLE of λ is 0.016, based on the two-parameter exponential fit.

(c) The MLE of m.

1/lambda + G

## [1] 68.5

The MLE of μ is approximately 68.5

(d) The probability of at least 100 correct entries between two errors.

S100 <- 1-pexp(100-G, rate=lambda)

print(S100)

## [1] 0.2222394

The probability of at least 100 correct entries between two errors is approximately 0.2222.

---

Problem 5

7.3 In a clinical study, 28 patients with cancer of the head and neck did not respond to chemotherapy. Their survival times in weeks are given in the following.

1.7  8.3   14.0  22.7   6.0+  13.1+
5.1  9.6   15.9  33.0   7.4+  13.4+
5.3  11.3  16.7   3.7+  8.0+  16.1+
6.0  12.1  17.0   5.0+  8.3+
8.3  12.3  21.0   5.9+  9.1+

Obtain the MLE of the parameter(s) and mean survival times, assuming:

(a) A one-parameter exponential distribution.

dat <- c(1.7, 5.1, 5.3, 6.0, 8.3, 8.3, 9.6, 11.3, 12.1, 12.3, 14.0, 15.9, 16.7, 17.0, 21.0, 22.7, 33.0, 3.7, 5.0, 5.9, 6.0, 7.4, 8.0, 8.3, 9.1, 13.1, 13.4, 16.1)

cen <- c(rep(1, 17), rep(0, 11))   # 1 if observed, 0 if censored

res_exp <- MLEexp(dat, cen)
res_exp

##   parameter    MLE    se approxCI.lb approxCI.ub
## 1    lambda  0.054 0.013       0.028       0.079
## 2        mu 18.606 4.513       9.761      27.450

Based on the results above, under the one-parameter exponential model, the estimated mean survival time is approximately 18.6 weeks, with a rate of λ = 0.054.

(b) A Weibull distribution.

res_wei <- MLE(c(0.05,1), "weibull", dat, cen)
res_wei

##   parameter   MLE    se CI.lb CI.ub
## 1    lambda 0.059 0.007 0.046 0.073
## 2     gamma 2.042 0.351 1.353 2.731

l_hat <- res_wei$MLE[1]
g_hat <- res_wei$MLE[2]

mu_weib <- (1 / l_hat) * gamma(1 + 1 / g_hat)

The MLEs for the Weibull model are λ = 0.059 and γ = 2.042, both with a mean survival time of approximately 17.45 weeks

---

Problem 6

7.4 In a study of deep venous thrombosis, the following blood clot lysis times in hours were recorded from 20 patients: 2, 3, 4, 5.5, 9, 13, 16.5, 17.5, 12.5, 7, 6, 17.5, 11.5, 6, 14, 25, 49, 37.5, 49, and 28. Assume that the blood clot lysis times follow the lognormal distribution.

(a) Obtain MLEs of the parameters μ and σ^2.

dat <- c(2, 3, 4, 5.5, 9, 13, 16.5, 17.5, 12.5, 7, 6, 17.5, 11.5, 6, 14, 25, 49, 37.5, 49, 28)

res_6 <- MLE(c(1,1), "lognormal", dat, se=T)
res_6

##   parameter   MLE    se CI.lb CI.ub
## 1        mu 2.464 0.195 2.081 2.846
## 2     sigma 0.872 0.138 0.602 1.142

sigma2_hat <- res_6$MLE[2]^2
sigma2_hat

## [1] 0.760384

Based on the results above, we have μ = 2.464 and σ^2 = 0.760384.

(b) Obtain 95% confidence intervals for μ and σ^2.

As shown on the results above, we have a 95% confidence intervals of (2.081, 2.846) for μ, and (0.602, 1.142) for σ^2.

---

Problem 7

7.5 Consider the following tumor-free times in days of 10 animals: 2, 3.5, 5, 7, 9, 10, 15, 20, 30, and 40. Assume that the tumor-free times follow the log-logistic distribution. Estimate the parameters α and γ.

dat <- c(2, 3.5, 5, 7, 9, 10, 15, 20, 30, 40)
cen <- rep(1, length(dat))

res_7 <- MLE(c(1,1), "loglogistic", dat, cen, se=T)
res_7

##   parameter   MLE    se  CI.lb CI.ub
## 1     alpha 0.014 0.017 -0.019 0.047
## 2     gamma 1.868 0.470  0.947 2.789

Based on the results, by using the log-logistic model, the estimated parameters are: the scale of α = 0.014 (so approximately 71.43 days when inverted), and the shape of γ = 1.868. With this, the wide confidence interval for α suggests some uncertainty. The shape estimate indicates a moderately increasing hazard over time.

---

Problem 8

7.9 Twenty-five rats were injected with a give tumor inoculum. Their times, in days, to the development of a tumor of a certain size are given in the following.

30 53 77  91 118
38 54 78  95 120
45 58 81 101 125
46 66 84 108 134
50 69 85 115 135

Which of the distributions discussed in this chapter provide a reasonable fit to the data? Estimate graphically the parameters of the distribution chosen.

dat <- c(30, 38, 45, 46, 50, 53, 54, 58, 66, 69, 77, 78, 81, 84, 85, 91, 95, 101, 108, 115, 118, 120, 125, 134, 135)

plotP("exponential", dat)

##   lambda    R2
## 1  0.017 0.808

plotP("weibull", dat)

##   lambda   gam    R2
## 1  0.011 3.035 0.973

plotP("lognormal", dat)

##      mu sigma    R2
## 1 4.332 0.411 0.962

plotP("loglogistic", dat)

##   alpha   gam    R2
## 1     0 4.343 0.947

After comparing all four probability plots, we can see that the Weibull distribution fits the best. The points align closely to the straight line with the highest R-square value of 0.973, and with the given parameters of λ ≈ 0.011 and γ ≈ 3.035.

---

Problem 9

7.10 Consider the data in Exercise 7.3.

(a) Make a hazard plot for each of the following distributions: exponential, Weibull, lognormal, and log-logistic.

dat <- c(1.7, 5.1, 5.3, 6.0, 8.3, 8.3, 9.6, 11.3, 12.1, 12.3, 14.0, 15.9, 16.7, 17.0, 21.0, 22.7, 33.0, 3.7, 5.0, 5.9, 6.0, 7.4, 8.0, 8.3, 9.1, 13.1, 13.4, 16.1)
cen <- c(rep(1, 17), rep(0, 11))

plotH("exponential", dat, cen)

##   lambda    R2
## 1  0.081 0.917

plotH("weibull", dat, cen)

##   lambda   gam    R2
## 1   0.06 1.761 0.962

plotH("lognormal", dat, cen)

##      mu sigma    R2
## 1 2.493 0.702 0.903

plotH("loglogistic", dat, cen)

##   alpha   gam    R2
## 1 0.002 2.471 0.916

(b) Which distribution provides a reasonable fit to the data? Estimate graphically the parameters of the distribution chosen.

Based on the plots above, the exponential hazard is constant but does not match the trend in the data, as the log-normal and log-logistic appear to be non-monotonic. The Weibull plot is the best fit with the R-square of 0.962 with the parameters of λ ≈ 0.06 and γ ≈ 1.761.

---

Problem 10

Redo Problem 7.10 by plotting the Cox-Snell residuals.

dat <- c(1.7, 5.1, 5.3, 6.0, 8.3, 8.3, 9.6, 11.3, 12.1, 12.3, 14.0, 15.9, 16.7, 17.0, 21.0, 22.7, 33.0, 3.7, 5.0, 5.9, 6.0, 7.4, 8.0, 8.3, 9.1, 13.1, 13.4, 16.1)
cen <- c(rep(1, 17), rep(0, 11))

res_exp <- MLE(c(1), "exponential", dat, cen)
r_exp <- attr(res_exp, "residual")
plotR(r_exp, cen)

res_wei <- MLE(c(0.05, 1), "weibull", dat, cen)
r_wei <- attr(res_wei, "residual")
plotR(r_wei, cen)

res_norm <- MLE(c(1,1), "lognormal", dat, cen)
r_lnorm <- attr(res_norm, "residual")
plotR(r_lnorm, cen)

res_log <- MLE(c(1,1), "loglogistic", dat, cen, se=F)
r_llog <- attr(res_log, "residual")
plotR(r_llog, cen)

Each Cox-Snell residual plot shows noticeable deviation from the 45-degree reference line, suggesting none of the distributions fit the data perfectly, but the lognormal distribution appears to come closest.

---

Problem 11

8.5 Consider the survival time of 28 cancer patients in Exercise 7.3

(a) Obtain the log-likelihoods for the exponential, Weibull, log-normal, and the extended generalized gamma distributions. Perform the likelihood ratio test and select the best distribution among these four distributions.

dat <- c(1.7, 5.1, 5.3, 6.0, 8.3, 8.3, 9.6, 11.3, 12.1, 12.3, 14.0, 15.9, 16.7, 17.0, 21.0, 22.7, 33.0, 3.7, 5.0, 5.9, 6.0, 7.4, 8.0, 8.3, 9.1, 13.1, 13.4, 16.1)
cen <- c(rep(1, 17), rep(0, 11))

sigcode <- function(p.val) ifelse(p.val < 0.001, "***", ifelse(p.val < 0.01, "**", ifelse(p.val < 0.05, "*", ifelse(p.val < 0.1, ".", ""))))

res <- list()   ### Example 8.2 & 8.4 in the textbook 4E
res$exponential <- MLE(  1     , "exponential", dat, cen, se=F)
res$weibull     <- MLE(c(1,1)  , "weibull"    , dat, cen, se=F)
res$lognormal   <- MLE(c(1,1)  , "lognormal"  , dat, cen, se=F)
res$gengamma    <- MLE(c(1,1,1), "gengamma"   , dat, cen, se=F)

p <- unlist(lapply(res, function(tab) length(tab$par)))
logL <- unlist(lapply(res, function(tab) attr(tab, "max.logL")))
tab <- data.frame(dist=names(logL), p=p, logL=logL); rownames(tab) <- NULL
tab

##          dist p      logL
## 1 exponential 1 -66.69912
## 2     weibull 2 -60.88686
## 3   lognormal 2 -62.48707
## 4    gengamma 3 -60.88234

ll <- rbind(cbind(tab[1,], tab[2,]),
            cbind(tab[1,], tab[4,]),
            cbind(tab[2,], tab[4,]),
            cbind(tab[3,], tab[4,]))
colnames(ll)[1:3] <- c("dist0", "p0", "logL0")
LRT <- -2*(ll$logL0 - ll$logL)
p.val <- 1-pchisq(LRT, df=(ll$p-ll$p0))
GOF <- ll[,c(1,4)]; rownames(GOF) <- NULL
GOF <- cbind(GOF, LRT=round(LRT,3), p.value=round(p.val,4), signif=sapply(p.val, sigcode))
GOF

##         dist0     dist    LRT p.value signif
## 1 exponential  weibull 11.625  0.0007    ***
## 2 exponential gengamma 11.634  0.0030     **
## 3     weibull gengamma  0.009  0.9242       
## 4   lognormal gengamma  3.209  0.0732      .

The likelihood ratio test suggests that the generalized gamma distribution fits the data best, as it significantly improves the fit over the exponential and log-normal models.

(b) Use the BIC and AIC procedures to select the best distribution among the four distributions in (a) plus the log-logistic distribution.

n <- length(dat)
tab$AIC <- logL-2*p
tab$BIC <- logL-p/2*log(n)
tab

##          dist p      logL       AIC       BIC
## 1 exponential 1 -66.69912 -68.69912 -68.36523
## 2     weibull 2 -60.88686 -64.88686 -64.21907
## 3   lognormal 2 -62.48707 -66.48707 -65.81928
## 4    gengamma 3 -60.88234 -66.88234 -65.88065

Based on both AIC and BIC values, the Weibull distribution provides the best fit among the four considered models.

(c) Compare the results obtained in (a) and (b) earlier and those obtained in Exercise 7.10.

The Weibull distribution, which was favored in Exercises 7.10 and 8.5(b), consistently shows the best fit across both the Cox-Snell residual analysis and model selection criteria, reinforcing its suitability for modeling the survival data.

---

Problem 12

Consider the following survival time in weeks of 10 mice with injection of tumor cells: 5, 16, 18+, 20, 22+, 24+, 25, 30+, 35, 40+. Do the data follow the exponential distribution with l = 0.02?

(a) Use the likelihood ratio test.

dat <- c(5, 16, 18, 20, 22, 24, 25, 30, 35, 40)
cen <- c(1,  1,  0,  1,  0,  0,  1,  0,  1,  0)

lambda0 <- 0.02
LR.test(lambda0, "exponential", dat, cen)

## Likelihood ratio test for a specific distribution with known parameters 
##    H0: S = S0 
##        or equivalently, the underlying distribution is exponential with specified parameter values. 
## 
##  parameter   MLE null.value
##     lambda 0.021       0.02
## 
##     logL  logL0   LRT p.value signif
##  -24.251 -24.26 0.019  0.8911

Based on the output, we get a p-value of 0.8911, which is quite high So, we do not reject the null hypothesis. This suggests that the exponential distribution with λ = 0.02 provides a reasonable fit for the data.

(b) Use the Hollander and Proschan’s test statistic.

sf0 <- function(tt) 1-pexp(tt, rate=lambda0)
HP.test(sf0, dat, cen)

## Hollander and Proschan's test for a specific distribution with known parameters 
##    H0: S = S0 
##        where S0 is the underlying survival function with specified parameter values 
## 
##      C sigma2 C.star
##  0.471  0.093 -0.296
## 
##       H1 p.value signif
##  S != S0  0.7671       
##   S < S0  0.6164       
##   S > S0  0.3836

We can see from the output that the p-value is large, so we fail to reject the null hypothesis, which means that there is no significant deviation from the exponential distribution.

(c) Plot the Kaplan–Meier estimator of S(t) and the hypothesized distribution.

Since the Kaplan-Meier estimate is provided in part b, we can see that the red line closely follows the steps of the KM curve within confidence bounds, supporting a good fit.

---

Problem 13

Suppose that T has an exponential distribution with parameter λ. However, the hazard function λ varies across the individuals (i.e., population heterogeneity). Specifically, suppose that the distribution of T given λ has the probability density function expressed as

π(λ) = λ^(k−1) e^(−λ/α) / α^(k) Γ(k) ,    λ > 0.

(a) Find the unconditional probability density function f(t).

(b) Find the unconditional survival function S(t).

(c) Find the unconditional hazard function h(t). Does the unconditional distribution of T have the increasing failure rate (IFR) property? or the decreasing failure rate (DFR) property? Justify.

…