Survival models and efficiency (optional)

A survival-research extension of the score-projection argument

Before reading: Model restrictions and efficiency, likelihood scores, and basic survival concepts. This extension introduces risk-set and counting-process notation; it is not required for the missing-outcome mean construction.

1 The Cox model

Let T be an event time, C a censoring time, and X a baseline covariate. One observed record is W=(U,\Delta,X), where U=\min(T,C) and \Delta=1\{T\le C\}. We observe an event when \Delta=1; otherwise we know only that the event had not occurred by U. Assume independent people, continuous times, and conditionally independent censoring, T\perp C\mid X. The censoring law has nuisance parameters distinct from the event-model parameters. This permits censoring to depend on X, but excludes unrecorded dependence between T and C after conditioning on X.

The Cox model specifies the event hazard

\lambda(t\mid X)=\lambda_0(t)e^{\beta X}.

The hazard is an instantaneous event rate among people who have survived to time t. It is not the probability of an event at an exact time. The unknown baseline hazard \lambda_0 is unrestricted; e^\beta is the hazard ratio for a one-unit increase in X. The target is the common log hazard ratio \beta under this model.

1.1 Express the score as observed events minus predicted events

Define the at-risk indicator R(t)=1\{U\ge t\} and event-counting process N(t)=1\{U\le t,\Delta=1\}. For a single person, N jumps from zero to one if an observed event occurs. The notation dN(t) denotes that jump. Write

dM(t)=dN(t)-R(t)e^{\beta X}\lambda_0(t)\,dt.

The second term is the event count predicted over a small time interval by the model and observed history. Their difference is a martingale innovation: its conditional expectation given the past is zero. The integral \int f\,dM adds weighted event jumps and subtracts the corresponding integrated predictions. We use a stated covariance identity below; developing stochastic integration is beyond this optional extension.

To see where the score comes from, put \Lambda_0(u)=\int_0^u\lambda_0(t)\,dt. Under the censoring assumption, the event-model part of one person’s log likelihood is

\ell=\Delta\{\log\lambda_0(U)+\beta X\} -e^{\beta X}\Lambda_0(U).

The remaining factors describe censoring and the covariate law. Differentiating the displayed expression with respect to \beta gives

\begin{aligned} S_\beta&=X\Delta-Xe^{\beta X}\Lambda_0(U)\\ &=\int X\,dM(t). \end{aligned}

A bounded baseline direction a(t) gives the path \lambda_{0,u}(t)=\lambda_0(t)\{1+u a(t)\}, positive for sufficiently small |u| wherever the original hazard is positive. Differentiating its log likelihood gives the nuisance score S_a=\int a(t)\,dM(t). Unlike S_\beta, its integrand depends only on time.

1.2 Project away the baseline-hazard directions

For suitable predictable, square-integrable integrands f and a, the martingale covariance identity is

E\left[\left(\int f\,dM\right)\left(\int a\,dM\right)\right]=E\int f(t)a(t)R(t)e^{\beta X}\lambda_0(t)\,dt.

Here “predictable” means that the integrand uses information available just before the time in question. The identity replaces the ordinary product calculation used in the partially linear model. At each time with a nonempty population risk set, the best time-only predictor of X minimizes E[R(t)e^{\beta X}\{X-a(t)\}^2]. Setting its derivative to zero gives E[R(t)e^{\beta X}\{X-a(t)\}]=0, hence

\bar x(t) = \frac{E\big[X\,e^{\beta X}R(t)\big]}{E\big[e^{\beta X}R(t)\big]},

This is a hazard-weighted risk-set mean: people with larger e^{\beta X} receive more weight. It is an ordinary mean among people at risk only when those weights are constant. Subtracting it makes the score orthogonal to every baseline-hazard direction. Under the stated censoring model, event-martingale scores are also orthogonal to the censoring and covariate-law nuisance scores. The efficient score, information per person, and EIF are therefore

\begin{aligned} S_{\text{eff}}&=\int\{X-\bar x(t)\}\,dM(t),\\ I_{\text{eff}}&=E\int\{X-\bar x(t)\}^2 \times R(t)e^{\beta X}\lambda_0(t)\,dt,\\ \phi_{\text{eff}}&=I_{\text{eff}}^{-1}S_{\text{eff}}. \end{aligned}

Assume 0<I_{\text{eff}}<\infty. The influence-variance bound is I_{\text{eff}}^{-1}, the estimator-variance bound is 1/(NI_{\text{eff}}), and the corresponding SE bound is 1/\sqrt{NI_{\text{eff}}}. Risk-set covariate variation and observed events supply information; severe censoring can reduce both.

1.3 Connect the projection to partial likelihood

The sample counterpart of the weighted mean is

\bar X_N(\beta,t)= \frac{\sum_iR_i(t)e^{\beta X_i}X_i} {\sum_iR_i(t)e^{\beta X_i}}.

By construction, \sum_iR_i(t)e^{\beta X_i}\{X_i-\bar X_N(\beta,t)\}=0. Consequently the integrated prediction terms cancel when the individual residual scores are summed. What remains is

U_N(\beta)=\sum_{i:\Delta_i=1} \{X_i-\bar X_N(\beta,U_i)\}.

At each event, compare that person’s covariate with the hazard-weighted mean among those still at risk, then add the differences. This is Cox’s partial-likelihood score for untied event times. The cancellation explains why the score can be computed without estimating \lambda_0. Under the relevant identification, integrability, and large-sample conditions, its root attains the EIF and is efficient. Deriving the population projection suggests the estimator; the counting-process asymptotic argument establishes attainment.

Compare this with the partially linear example in Model restrictions and efficiency. Both remove the component explained by the nuisance directions, although their scores and inner-product weights differ. In the Cox model the projection is onto functions of time under risk-set weights. What remains supplies information about \beta that changing the baseline hazard cannot reproduce.

2 Checking the Cox model

The simulated Cox model has a covariate X\sim N(0,1), a Weibull baseline hazard \lambda_0(t)=1.5t^{0.5}, \beta=0.5, and independent censoring uniform on (0,2), which censors about 42% of observations. The continuous draws avoid ties. Four estimators are compared: coxph; a score with each event weighted by w(t)=e^{-2t}; Weibull regression; and exponential regression.

The time-weighted score still has population expectation zero at the true coefficient. With a unique root and appropriate regularity, it can consistently estimate \beta, but its nonconstant weighting generally loses efficiency. Weibull regression makes an additional, correct restriction on the baseline family while estimating its parameters. Exponential regression imposes a constant baseline hazard, incorrect in this experiment.

Code

library(survival)
beta <- 0.5
SimCox <- function(n) {
  x <- rnorm(n); u <- runif(n)
  t <- (-log(u) / exp(beta * x))^(1 / 1.5)        # hazard 1.5 * t^0.5 * exp(beta * x)
  cens <- runif(n, 0, 2)
  data.frame(time = pmin(t, cens), status = as.numeric(t <= cens), x = x)
}

# At each event time: the covariate of the person who fails minus the average covariate of the
# risk set, multiplied by a weight w(t). With w = 1 this is Cox's partial likelihood score.
WeightedCox <- function(d, w = function(t) rep(1, length(t))) {
  d <- d[order(d$time), ]; wt <- w(d$time)
  Score <- function(b) {
    r <- exp(b * d$x)
    s0 <- rev(cumsum(rev(r))); s1 <- rev(cumsum(rev(d$x * r)))   # sums over the risk set at each time
    sum(wt * d$status * (d$x - s1 / s0))
  }
  uniroot(Score, c(-5, 5))$root
}

set.seed(2027)
d <- SimCox(400)
c(censored = 1 - mean(d$status), coxph = unname(coef(coxph(Surv(time, status) ~ x, d))),
  unit_weight_score = WeightedCox(d))

         censored             coxph unit_weight_score 
        0.4175000         0.6932629         0.6932629

With every weight equal to 1, the custom score agrees numerically with coxph on this sample. This is a useful implementation check for the untied-time example. It does not itself establish efficiency; that follows from the model-specific argument. Now compare repeated samples.

Code

EstCox <- function(d) {
  cf <- coxph(Surv(time, status) ~ x, d)
  wb <- survreg(Surv(time, status) ~ x, d, dist = "weibull")
  ex <- survreg(Surv(time, status) ~ x, d, dist = "exponential")
  # survreg uses an accelerated failure time scale; dividing by the Weibull scale gives the hazard-ratio coefficient
  c(coxph = unname(coef(cf)), coxph_se = sqrt(vcov(cf)[1, 1]),
    weighted = WeightedCox(d, function(t) exp(-2 * t)),
    weibull = -unname(coef(wb)["x"]) / wb$scale,
    exponential = -unname(coef(ex)["x"]))
}
set.seed(2028)
R <- 1000; n <- 400
res <- t(replicate(R, EstCox(SimCox(n))))
est <- setdiff(colnames(res), "coxph_se")
data.frame(estimator = est, mean = round(colMeans(res[, est]), 3),
           sd = round(apply(res[, est], 2, sd), 4), row.names = NULL)

    estimator  mean     sd
1       coxph 0.501 0.0729
2    weighted 0.503 0.0828
3     weibull 0.503 0.0724
4 exponential 0.408 0.0574

Code

round(mean(res[, "coxph_se"]), 4)    # the average model-based standard error of coxph

[1] 0.072

Compare each mean with the true value 0.5 and the standard deviations with one another. The Monte Carlo standard error is about 0.002 for each mean and 0.0016 for each standard deviation.

Cox’s estimator is approximately centered at the target in this simulation, and its standard deviation, 0.073, is close to its average model-based SE, 0.072. Its squared SE estimates 1/(NI_{\text{eff}}), the estimator-variance bound; the SE itself estimates the square root of that bound.
The weighted estimator is also approximately centered in this simulation. Its standard deviation, 0.083, is about 14% above Cox’s, which is a variance about 29% higher.
Weibull regression uses the correct baseline family. Its standard deviation, 0.072, is close to Cox’s within Monte Carlo error. This experiment detects little precision gain from that additional restriction; it does not prove equal efficiency bounds for the two models. The cost of nuisance freedom depends on the model, the population law, and which restrictions are compared.
The exponential regression has a smaller standard deviation, 0.057, and the wrong center, 0.41. It assumes a constant hazard, which is a working model in the sense of What is your estimator actually estimating?, and when the hazard is not constant its limit differs from \beta.

3 Reading

Andersen and Gill (1982), Cox’s regression model for counting processes: a large sample study, develops the counting-process large-sample argument. Begun, Hall, Huang, and Wellner (1983), Information and asymptotic efficiency in parametric–nonparametric models, supplies the efficiency framework.

Nan, Emond, and Wellner (2004), §3.1, gives the full-data Cox event, censoring, and covariate-law scores used in the projection. The rest of that paper treats missing covariates and goes beyond this extension.

4 Additional practice (optional)

Replace the time weight by an indicator for events before the median event time. Predict the variance change, then check it.
Simulate a time-varying coefficient and investigate what a constant-coefficient score estimates. Keep the working model, target, and censoring assumptions distinct.

--- title: "Survival models and efficiency (optional)" subtitle: "A survival-research extension of the score-projection argument" --- **Before reading:** [Model restrictions and efficiency](models-and-efficiency.qmd), likelihood scores, and basic survival concepts. This extension introduces risk-set and counting-process notation; it is not required for the missing-outcome mean construction. ## The Cox model Let $T$ be an event time, $C$ a censoring time, and $X$ a baseline covariate. One observed record is $W=(U,\Delta,X)$, where $U=\min(T,C)$ and $\Delta=1\{T\le C\}$. We observe an event when $\Delta=1$; otherwise we know only that the event had not occurred by $U$. Assume independent people, continuous times, and **conditionally independent censoring**, $T\perp C\mid X$. The censoring law has nuisance parameters distinct from the event-model parameters. This permits censoring to depend on $X$, but excludes unrecorded dependence between $T$ and $C$ after conditioning on $X$. The Cox model specifies the event hazard $$ \lambda(t\mid X)=\lambda_0(t)e^{\beta X}. $$ The hazard is an instantaneous event rate among people who have survived to time $t$. It is not the probability of an event at an exact time. The unknown baseline hazard $\lambda_0$ is unrestricted; $e^\beta$ is the hazard ratio for a one-unit increase in $X$. The target is the common log hazard ratio $\beta$ under this model. ### Express the score as observed events minus predicted events Define the at-risk indicator $R(t)=1\{U\ge t\}$ and event-counting process $N(t)=1\{U\le t,\Delta=1\}$. For a single person, $N$ jumps from zero to one if an observed event occurs. The notation $dN(t)$ denotes that jump. Write $$ dM(t)=dN(t)-R(t)e^{\beta X}\lambda_0(t)\,dt. $$ The second term is the event count predicted over a small time interval by the model and observed history. Their difference is a martingale innovation: its conditional expectation given the past is zero. The integral $\int f\,dM$ adds weighted event jumps and subtracts the corresponding integrated predictions. We use a stated covariance identity below; developing stochastic integration is beyond this optional extension. To see where the score comes from, put $\Lambda_0(u)=\int_0^u\lambda_0(t)\,dt$. Under the censoring assumption, the event-model part of one person's log likelihood is $$ \ell=\Delta\{\log\lambda_0(U)+\beta X\} -e^{\beta X}\Lambda_0(U). $$ The remaining factors describe censoring and the covariate law. Differentiating the displayed expression with respect to $\beta$ gives $$ \begin{aligned} S_\beta&=X\Delta-Xe^{\beta X}\Lambda_0(U)\\ &=\int X\,dM(t). \end{aligned} $$ A bounded baseline direction $a(t)$ gives the path $\lambda_{0,u}(t)=\lambda_0(t)\{1+u a(t)\}$, positive for sufficiently small $|u|$ wherever the original hazard is positive. Differentiating its log likelihood gives the nuisance score $S_a=\int a(t)\,dM(t)$. Unlike $S_\beta$, its integrand depends only on time. ### Project away the baseline-hazard directions For suitable predictable, square-integrable integrands $f$ and $a$, the martingale covariance identity is $$ E\left[\left(\int f\,dM\right)\left(\int a\,dM\right)\right]=E\int f(t)a(t)R(t)e^{\beta X}\lambda_0(t)\,dt. $$ Here “predictable” means that the integrand uses information available just before the time in question. The identity replaces the ordinary product calculation used in the partially linear model. At each time with a nonempty population risk set, the best time-only predictor of $X$ minimizes $E[R(t)e^{\beta X}\{X-a(t)\}^2]$. Setting its derivative to zero gives $E[R(t)e^{\beta X}\{X-a(t)\}]=0$, hence $$ \bar x(t) = \frac{E\big[X\,e^{\beta X}R(t)\big]}{E\big[e^{\beta X}R(t)\big]}, $$ This is a **hazard-weighted** risk-set mean: people with larger $e^{\beta X}$ receive more weight. It is an ordinary mean among people at risk only when those weights are constant. Subtracting it makes the score orthogonal to every baseline-hazard direction. Under the stated censoring model, event-martingale scores are also orthogonal to the censoring and covariate-law nuisance scores. The efficient score, information per person, and EIF are therefore $$ \begin{aligned} S_{\text{eff}}&=\int\{X-\bar x(t)\}\,dM(t),\\ I_{\text{eff}}&=E\int\{X-\bar x(t)\}^2 \times R(t)e^{\beta X}\lambda_0(t)\,dt,\\ \phi_{\text{eff}}&=I_{\text{eff}}^{-1}S_{\text{eff}}. \end{aligned} $$ Assume $0<I_{\text{eff}}<\infty$. The influence-variance bound is $I_{\text{eff}}^{-1}$, the estimator-variance bound is $1/(NI_{\text{eff}})$, and the corresponding SE bound is $1/\sqrt{NI_{\text{eff}}}$. Risk-set covariate variation and observed events supply information; severe censoring can reduce both. ### Connect the projection to partial likelihood The sample counterpart of the weighted mean is $$ \bar X_N(\beta,t)= \frac{\sum_iR_i(t)e^{\beta X_i}X_i} {\sum_iR_i(t)e^{\beta X_i}}. $$ By construction, $\sum_iR_i(t)e^{\beta X_i}\{X_i-\bar X_N(\beta,t)\}=0$. Consequently the integrated prediction terms cancel when the individual residual scores are summed. What remains is $$ U_N(\beta)=\sum_{i:\Delta_i=1} \{X_i-\bar X_N(\beta,U_i)\}. $$ At each event, compare that person's covariate with the hazard-weighted mean among those still at risk, then add the differences. This is Cox's partial-likelihood score for untied event times. The cancellation explains why the score can be computed without estimating $\lambda_0$. Under the relevant identification, integrability, and large-sample conditions, its root attains the EIF and is efficient. Deriving the population projection suggests the estimator; the counting-process asymptotic argument establishes attainment. Compare this with the partially linear example in [Model restrictions and efficiency](models-and-efficiency.qmd). Both remove the component explained by the nuisance directions, although their scores and inner-product weights differ. In the Cox model the projection is onto functions of time under risk-set weights. What remains supplies information about $\beta$ that changing the baseline hazard cannot reproduce. ## Checking the Cox model The simulated Cox model has a covariate $X\sim N(0,1)$, a Weibull baseline hazard $\lambda_0(t)=1.5t^{0.5}$, $\beta=0.5$, and independent censoring uniform on $(0,2)$, which censors about 42% of observations. The continuous draws avoid ties. Four estimators are compared: `coxph`; a score with each event weighted by $w(t)=e^{-2t}$; Weibull regression; and exponential regression. The time-weighted score still has population expectation zero at the true coefficient. With a unique root and appropriate regularity, it can consistently estimate $\beta$, but its nonconstant weighting generally loses efficiency. Weibull regression makes an additional, correct restriction on the baseline family while estimating its parameters. Exponential regression imposes a constant baseline hazard, incorrect in this experiment. ```{r} library(survival) beta <- 0.5 SimCox <- function(n) { x <- rnorm(n); u <- runif(n) t <- (-log(u) / exp(beta * x))^(1 / 1.5) # hazard 1.5 * t^0.5 * exp(beta * x) cens <- runif(n, 0, 2) data.frame(time = pmin(t, cens), status = as.numeric(t <= cens), x = x) } # At each event time: the covariate of the person who fails minus the average covariate of the # risk set, multiplied by a weight w(t). With w = 1 this is Cox's partial likelihood score. WeightedCox <- function(d, w = function(t) rep(1, length(t))) { d <- d[order(d$time), ]; wt <- w(d$time) Score <- function(b) { r <- exp(b * d$x) s0 <- rev(cumsum(rev(r))); s1 <- rev(cumsum(rev(d$x * r))) # sums over the risk set at each time sum(wt * d$status * (d$x - s1 / s0)) } uniroot(Score, c(-5, 5))$root } set.seed(2027) d <- SimCox(400) c(censored = 1 - mean(d$status), coxph = unname(coef(coxph(Surv(time, status) ~ x, d))), unit_weight_score = WeightedCox(d)) ``` With every weight equal to 1, the custom score agrees numerically with `coxph` on this sample. This is a useful implementation check for the untied-time example. It does not itself establish efficiency; that follows from the model-specific argument. Now compare repeated samples. ```{r} EstCox <- function(d) { cf <- coxph(Surv(time, status) ~ x, d) wb <- survreg(Surv(time, status) ~ x, d, dist = "weibull") ex <- survreg(Surv(time, status) ~ x, d, dist = "exponential") # survreg uses an accelerated failure time scale; dividing by the Weibull scale gives the hazard-ratio coefficient c(coxph = unname(coef(cf)), coxph_se = sqrt(vcov(cf)[1, 1]), weighted = WeightedCox(d, function(t) exp(-2 * t)), weibull = -unname(coef(wb)["x"]) / wb$scale, exponential = -unname(coef(ex)["x"])) } set.seed(2028) R <- 1000; n <- 400 res <- t(replicate(R, EstCox(SimCox(n)))) est <- setdiff(colnames(res), "coxph_se") data.frame(estimator = est, mean = round(colMeans(res[, est]), 3), sd = round(apply(res[, est], 2, sd), 4), row.names = NULL) round(mean(res[, "coxph_se"]), 4) # the average model-based standard error of coxph ``` Compare each mean with the true value 0.5 and the standard deviations with one another. The Monte Carlo standard error is about 0.002 for each mean and 0.0016 for each standard deviation. - Cox's estimator is approximately centered at the target in this simulation, and its standard deviation, 0.073, is close to its average model-based SE, 0.072. Its squared SE estimates $1/(NI_{\text{eff}})$, the estimator-variance bound; the SE itself estimates the square root of that bound. - The weighted estimator is also approximately centered in this simulation. Its standard deviation, 0.083, is about 14% above Cox's, which is a variance about 29% higher. - Weibull regression uses the correct baseline family. Its standard deviation, 0.072, is close to Cox's within Monte Carlo error. This experiment detects little precision gain from that additional restriction; it does not prove equal efficiency bounds for the two models. The cost of nuisance freedom depends on the model, the population law, and which restrictions are compared. - The exponential regression has a smaller standard deviation, 0.057, and the wrong center, 0.41. It assumes a constant hazard, which is a working model in the sense of [What is your estimator actually estimating?](../questions-and-design/estimands-and-misspecification.qmd), and when the hazard is not constant its limit differs from $\beta$. ## Reading Andersen and Gill (1982), *Cox’s regression model for counting processes: a large sample study*, develops the counting-process large-sample argument. Begun, Hall, Huang, and Wellner (1983), *Information and asymptotic efficiency in parametric–nonparametric models*, supplies the efficiency framework. [Nan, Emond, and Wellner (2004), §3.1](https://faculty.washington.edu/jonw/JAW-papers/jaw-NanEmond-04AS.pdf), gives the full-data Cox event, censoring, and covariate-law scores used in the projection. The rest of that paper treats missing covariates and goes beyond this extension. ## Additional practice (optional) 1. Replace the time weight by an indicator for events before the median event time. Predict the variance change, then check it. 2. Simulate a time-varying coefficient and investigate what a constant-coefficient score estimates. Keep the working model, target, and censoring assumptions distinct.