stdParfrailty {stdReg}

Regression standardization in shared frailty gamma-Weibull models

Description

stdParfrailty performs regression standardization in shared frailty gamma-Weibull 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. stdParfrailty 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

stdParfrailty(fit, data, X, x, t, clusterid, subsetnew)

Arguments

fit

an object of class "parfrailty", as returned by the parfrailty function in the stdReg package.

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 the mean of X. If x is set to NA, then X is not altered. This produces an estimate of the marginal survival function S(t)=E\{S(t|X,Z)\}.

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.

clusterid

a string containing the name of the cluster identification variable.

subsetnew

an optional logical statement specifying a subset of observations to be used in the standardization. This set is assumed to be a subset of the subset (if any) that was used to fit the regression model.

Details

stdParfrailty assumes that a shared frailty gamma-Weibull model

\lambda(t_{ij}|X_{ij},Z_{ij})=\lambda(t_{ij};\alpha,\eta)U_iexp\{h(X_{ij},Z_{ij};\beta)\}

has been fitted, with parametrization as descibed in the help section for parfrailty. Integrating out the gamma frailty gives the survival function

S(t|X,Z)=[1+\phi\Lambda_0(t;\alpha,\eta)exp\{h(X,Z;\beta)\}]^{-1/\phi},

where \Lambda_0(t;\alpha,\eta) is the cumulative baseline hazard

(t/\alpha)^{\eta}.

The ML estimates of (\alpha,\eta,\phi,\beta) are used to obtain estimates of the survival function S(t|X=x,Z):

\hat{S}(t|X=x,Z)=[1+\hat{\phi}\Lambda_0(t;\hat{\alpha},\hat{\eta})exp\{h(X,Z;\hat{\beta})\}]^{-1/\hat{\phi}}.

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.

Value

An object of class "stdParfrailty" is a list containing

call

the matched call.

input

input is a list containing all input arguments.

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]).

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]).

Note

Standardized survival functions are sometimes referred to as (direct) adjusted survival functions in the literature.

stdParfrailty does not currently handle time-varying exposures or covariates.

stdParfrailty internally loops over all values in the t argument. Therefore, the function will usually be considerably faster if length(t) is small.

The variance calculation performed by stdParfrailty does not condition on the observed covariates \bar{Z}=(Z_1,...,Z_n). To see how this matters, note that

var\{\hat{\theta}(t,x)\}=E[var\{\hat{\theta}(t,x)|\bar{Z}\}]+var[E\{\hat{\theta}(t,x)|\bar{Z}\}].

The usual parameter \beta in a Cox proportional hazards model does not depend on \bar{Z}. Thus, E(\hat{\beta}|\bar{Z}) is independent of \bar{Z} as well (since E(\hat{\beta}|\bar{Z})=\beta), so that the term var[E\{\hat{\beta}|\bar{Z}\}] in the corresponding variance decomposition for var(\hat{\beta}) becomes equal to 0. However, \theta(t,x) depends on \bar{Z} through the average over the sample distribution for Z, and thus the term var[E\{\hat{\theta}(t,x)|\bar{Z}\}] is not 0, unless one conditions on \bar{Z}. The variance calculation by Gail and Byar (1986) ignores this term, and thus effectively conditions on \bar{Z}.

Author(s)

Arvid Sjolander

References

Chang I.M., Gelman G., Pagano M. (1982). Corrected group prognostic curves and summary statistics. Journal of Chronic Diseases 35, 669-674.

Dahlqwist E., Pawitan Y., Sjolander A. (2019). Regression standardization and attributable fraction estimation with between-within frailty models for clustered survival data. Statistical Methods in Medical Research 28(2), 462-485.

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.

Makuch R.W. (1982). Adjusted survival curve estimation using covariates. Journal of Chronic Diseases 35, 437-443.

Examples

## Not run: 

require(survival)

#simulate data
n <- 1000
m <- 3
alpha <- 1.5
eta <- 1
phi <- 0.5
beta <- 1
id <- rep(1:n, each=m)
U <- rep(rgamma(n, shape=1/phi, scale=phi), each=m)
X <- rnorm(n*m)
#reparametrize scale as in rweibull function
weibull.scale <- alpha/(U*exp(beta*X))^(1/eta)
T <- rweibull(n*m, shape=eta, scale=weibull.scale)

#right censoring
C <- runif(n*m, 0, 10)
D <- as.numeric(T<C)
T <- pmin(T, C)
  
#strong left-truncation
L <- runif(n*m, 0, 2)
incl <- T>L
incl <- ave(x=incl, id, FUN=sum)==m
dd <- data.frame(L, T, D, X, id)
dd <- dd[incl, ]  
 
fit <- parfrailty(formula=Surv(L, T, D)~X, data=dd, clusterid="id")
fit.std <- stdParfrailty(fit=fit, data=dd, X="X", x=seq(-1,1,0.5), t=1:5, clusterid="id")
print(summary(fit.std, t=3))
plot(fit.std)


## End(Not run)