stdReg (version 1.0)

stdCoxph: Regression standardization in Cox proportional hazards models

Description

stdCoxph performs regression standardization in Cox proportional hazards models, at specified values of the exposure, over the sample covariate distribution. Let $T$, $X$, and $Z$ be the survival outcome, the exposure, and a vector of covariates, respectively. stdCoxph uses a fitted Cox proportional hazards model to estimate the standardized survival function $\theta(t,x)=E{S(t|X=x,Z)}$, where $t$ is a specific value of $T$, $x$ is a specific value of $X$, and the expectation is over the marginal distribution of $Z$.

Usage

stdCoxph(fit, data, X, x, t, clusters)

Arguments

fit

an object of class "coxph", as returned by the coxph function in the survival package, but without special terms strata, cluster, or tt.

data

a data frame containing the variables in the model. This should be the same data frame as was used to fit the model in fit.

X

a string containing the name of the exposure variable $X$ in data.

x

an optional vector containing the specific values of $X$ at which to estimate the standardized survival function. If $X$ is binary (0/1) or a factor, then x defaults to all values of $X$. If $X$ is numeric, then x defaults to

t

an optional vector containing the specific values of $T$ at which to estimate the standardized survival function. It defaults to all the observed event times in data.

clusters

an optional vector of cluster identification numbers to be used in the fitting process when data are clustered. Should be NULL or a numeric vector.

Value

• An object of class "stdCoxph" is a list containing all input arguments to the stdCoxph function, except data. In addition the list contains:
• fit
  the fitted coxph object.
• surv.est
  a matrix with length(t) rows and length(x) columns, where the element on row i and column j is equal to $\hat{\theta}$(t[i],x[j]).
• surv.vcov
  a list with length(t) elements. Each element is a square matrix with length(x) rows. In the k:th matrix, the element on row i and column j is the (estimated) covariance of $\hat{\theta}$(t[k],x[i]) and $\hat{\theta}$(t[k],x[j]).

Details

stdCoxph assumes that a Cox proportional hazards model $$\lambda(t|X,Z)=\lambda_0(t)exp{h(X,Z;\beta)}$$ has been fitted. Breslow's estimator of the cumulative baseline hazard $\Lambda_0(t)=\int_0^t\lambda_0(u)du$ is used together with the partial likelihood estimate of $\beta$ to obtain estimates of the survival function $S(t|X=x,Z)$: $$\hat{S}(t|X=x,Z)=exp[-\hat{\Lambda}_0(t)exp{h(X=x,Z;\hat{\beta})}].$$ For each $t$ in the t argument and for each $x$ in the x argument, these estimates are averaged across all subjects (i.e. all observed values of $Z$) to produce estimates $$\hat{\theta}(t,x)=\sum_{i=1}^n \hat{S}(t|X=x,Z_i)/n.$$ The variance for $\hat{\theta}(t,x)$ is obtained by the sandwich formula.

References

Makuch R.W. (1982). Adjusted survival curve estimation using covariates. Journal of Chronic Diseases 35, 437-443. Chang I.M., Gelman G., Pagano M. (1982). Corrected group prognostic curves and summary statistics. Journal of Chronic Diseases 35, 669-674. Gail M.H. and Byar D.P. (1986). Variance calculations for direct adjusted survival curves, with applications to testing for no treatement effect. Biometrical Journal 28(5), 587-599.

Examples

require(survival)

n <- 1000
Z <- rnorm(n)
X <- rnorm(n, mean = Z)
T <- rexp(n, rate = exp(X + Z + X*Z)) #survival time
C <- rexp(n, rate = exp(X + Z + X*Z)) #censoring time
U <- pmin(T, C) #time at risk
D <- as.numeric(T < C) #event indicator
dd <- data.frame(Z,X,U,D)
fit <- coxph(formula = Surv(U, D) ~ X + Z + X*Z, data = dd)
fit.std <- stdCoxph(fit = fit, data = dd, X = "X", x = seq(-1,1,0.5), t = 1:5)
print(summary(fit.std, t = 3))
plot(fit.std)