stdParfrailty: Regression standardization in shared frailty gamma-Weibull...
In stdReg: Regression Standardization

Description Usage Arguments Details Value Note Author(s) References Examples

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 θ(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.

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

`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.

stdParfrailty assumes that a shared frailty gamma-Weibull model

λ(t_{ij}|X_{ij},Z_{ij})=λ(t_{ij};α,η)U_iexp\{h(X_{ij},Z_{ij};β)\}

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+φΛ_0(t;α,η)exp\{h(X,Z;β)\}]^{-1/φ},

where Λ_0(t;α,η) is the cumulative baseline hazard

(t/α)^{η}.

The ML estimates of (α,η,φ,β) are used to obtain estimates of the survival function S(t|X=x,Z):

\hat{S}(t|X=x,Z)=[1+\hat{φ}Λ_0(t;\hat{α},\hat{η})exp\{h(X,Z;\hat{β})\}]^{-1/\hat{φ}}.

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{θ}(t,x)=∑_{i=1}^n \hat{S}(t|X=x,Z_i)/n.

The variance for \hat{θ}(t,x) is obtained by the sandwich formula.

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{θ}(`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{θ}(`t[k]`,`x[i]`) and \hat{θ}(`t[k]`,`x[j]`).

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{θ}(t,x)\}=E[var\{\hat{θ}(t,x)|\bar{Z}\}]+var[E\{\hat{θ}(t,x)|\bar{Z}\}].

The usual parameter β in a Cox proportional hazards model does not depend on \bar{Z}. Thus, E(\hat{β}|\bar{Z}) is independent of \bar{Z} as well (since E(\hat{β}|\bar{Z})=β), so that the term var[E\{\hat{β}|\bar{Z}\}] in the corresponding variance decomposition for var(\hat{β}) becomes equal to 0. However, θ(t,x) depends on \bar{Z} through the average over the sample distribution for Z, and thus the term var[E\{\hat{θ}(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}.

Arvid Sjolander

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.

## 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)