Software 3.5.1 (The exponential fit by hand). Example 3.2.2: the number of deaths and total follow-up, the estimated rate and its standard error per year, the \(95\%\) interval on the log scale, and the fitted mean and median in days:

> library(survival)
> t <- lung$time; d <- as.numeric(lung$status == 2)
> c(deaths = sum(d), total_days = sum(t))
    deaths total_days
       165      69593
> lam <- sum(d) / sum(t); se <- lam / sqrt(sum(d))     # MLE and its SE
> round(c(rate_per_year = 365.25 * lam, se = 365.25 * se), 4)
rate_per_year            se
       0.8660        0.0674
> ci <- exp(log(lam) + c(-1, 1) * 1.96 / sqrt(sum(d)))  # CI on the log scale
> round(365.25 * ci, 4)
[1] 0.7434 1.0087
> round(c(mean_days = 1 / lam, median_days = log(2) / lam), 1)
  mean_days median_days
      421.8       292.4

Reading the output. The rate is \(0.866\) deaths per person-year (s.e. \(0.067\)), with interval \((0.743, 1.009)\); the fitted median survival is \(292\) days.

Python: the same arithmetic; lifelines.ExponentialFitter().fit(t, d) gives the fitted mean, \(1/\hat \lambda \), as lambda_.

Software 3.5.2 (Four families by survreg). The exponential, Weibull, log-normal and log-logistic models fitted with survreg, their log-likelihoods and AICs; the Weibull’s location and scale and their translation to \(k\) and \(\lambda \); and the likelihood ratio test of \(k = 1\):

> fits <- lapply(c("exponential", "weibull", "lognormal", "loglogistic"),
+                function(dist) survreg(Surv(time, status) ~ 1, data = lung,
+                                       dist = dist))
> names(fits) <- c("exponential", "weibull", "lognormal", "loglogistic")
> round(sapply(fits, function(f) c(logLik = f$loglik[1], AIC = AIC(f))), 2)
       exponential  weibull lognormal loglogistic
logLik    -1162.34 -1153.85  -1169.27    -1160.93
AIC        2326.68  2311.70   2342.54     2325.86
> w <- fits$weibull
> c(intercept = unname(coef(w)), scale = w$scale)      # log T = mu + sigma W
intercept     scale
6.0349039 0.7593936
> round(c(shape_k = 1 / w$scale, lambda_days = exp(coef(w))[[1]]), 3)
    shape_k lambda_days
      1.317     417.759
> lr <- 2 * (w$loglik[1] - fits$exponential$loglik[1])
> round(c(LR = lr, p = pchisq(lr, 1, lower.tail = FALSE)), 5)  # is k = 1?
      LR        p
16.97398  0.00004

Reading the output. The Weibull has the largest log-likelihood among the two-parameter models and the smallest AIC. Its scale \(\hat \sigma = 0.759\) gives \(\hat k = 1/0.759 = 1.32\), and its intercept \(6.035\) gives \(\hat \lambda = e^{6.035} = 418\) days. The test of \(k = 1\) gives \(\chi ^2_1 = 17.0\), \(p = 0.00004\).

Python: lifelines.WeibullFitter, LogNormalFitter, LogLogisticFitter, each with an AIC_ attribute.

Software 3.5.3 (The Weibull likelihood by hand). The censored Weibull log-likelihood of the subsection on the Weibull family, maximised with optim over \((\log k, \log \lambda )\) so that both stay positive; the standard errors from the inverse Hessian; and the fitted median \(\lambda (\log 2)^{1/k}\):

> nll <- function(p) {              # Weibull by hand, p = (log k, log lam)
+   k <- exp(p[1]); l <- exp(p[2])
+   -sum(d * (log(k) - log(l) + (k - 1) * log(t / l)) - (t / l)^k)
+ }
> o <- optim(c(0, log(300)), nll, hessian = TRUE)
> round(c(k = exp(o$par[1]), lambda = exp(o$par[2]), logLik = -o$value), 3)
        k    lambda    logLik
    1.316   417.723 -1153.851
> round(sqrt(diag(solve(o$hessian))), 4)   # SEs of log k and log lambda
[1] 0.0624 0.0591
> med <- exp(o$par[2]) * log(2)^(1 / exp(o$par[1])); round(med, 1)
[1] 316.2

Reading the output. The hand-written likelihood reproduces survreg: \(\hat k = 1.316\), \(\hat \lambda = 418\) days, \(\ell = -1153.85\). The standard error of \(\log \hat k\) is \(0.062\), so \(\hat k\) is about \(1.32 \pm 0.16\), clearly above \(1\). The median is \(316\) days.

Python: scipy.optimize.minimize with the same function.