View source: R/parfrailty_methods.R
standardize_parfrailty | R Documentation |
standardize_parfrailty
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.
standardize_parfrailty
fits a parametric frailty 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
.
standardize_parfrailty(
formula,
data,
values,
times,
clusterid,
ci_level = 0.95,
ci_type = "plain",
contrasts = NULL,
family = "gaussian",
reference = NULL,
transforms = NULL
)
formula |
The formula which is used to fit the model for the outcome. |
data |
The data. |
values |
A named list or data.frame specifying the variables and values at which marginal means of the outcome will be estimated. |
times |
A vector containing the specific values of |
clusterid |
An optional string containing the name of a cluster identification variable when data are clustered. |
ci_level |
Coverage probability of confidence intervals. |
ci_type |
A string, indicating the type of confidence intervals. Either "plain", which gives untransformed intervals, or "log", which gives log-transformed intervals. |
contrasts |
A vector of contrasts in the following format:
If set to |
family |
The family argument which is used to fit the glm model for the outcome. |
reference |
A vector of reference levels in the following format:
If |
transforms |
A vector of transforms in the following format:
If set to |
standardize_parfrailty
fits 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)\}
, with parameterization as described 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.
An object of class std_surv
. Obtain numeric results by using tidy.std_surv.
This is a list with the following components:
An unnamed list with one element for each of the requested contrasts. Each element is itself a list with the elements:
The function call
A list with components used in the estimation
Either "survival" or "rmean"
Estimated counterfactual means and standard errors for each exposure level
Estimated covariance matrix of counterfactual means for each time
Data.frame of the estimates of the contrast with inference
The vector of times used in the calculation
The transform argument used for this contrast
The requested contrast type
The reference level of the exposure
Confidence interval type
Confidence interval level
A named list with the elements:
The function call
A list with components used in the estimation
Either "survival" or "rmean"
Estimated counterfactual means and standard errors for each exposure level
Estimated covariance matrix of counterfactual means for each time
Standardized survival functions are sometimes referred to as (direct) adjusted survival functions in the literature.
standardize_coxph/standardize_parfrailty
does not currently handle time-varying exposures or
covariates.
standardize_coxph/standardize_parfrailty
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 standardize_coxph
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}
.
Arvid Sjölander
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., Sjölander 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 treatment effect. Biometrical Journal 28(5), 587-599.
Makuch R.W. (1982). Adjusted survival curve estimation using covariates. Journal of Chronic Diseases 35, 437-443.
require(survival)
# simulate data
set.seed(6)
n <- 300
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)
# reparameterize 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.std <- standardize_parfrailty(
formula = Surv(L, T, D) ~ X,
data = dd,
values = list(X = seq(-1, 1, 0.5)),
times = 1:5,
clusterid = "id"
)
print(fit.std)
plot(fit.std)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.