standardize_level: Get standardized estimates using the g-formula with and separate models for each exposure level in the data

Description

Get standardized estimates using the g-formula with and separate models for each exposure level in the data

Usage

standardize_level(
  fitter_list,
  arguments,
  predict_fun_list,
  data,
  values,
  B = NULL,
  ci_level = 0.95,
  contrasts = NULL,
  reference = NULL,
  seed = NULL,
  times = NULL,
  transforms = NULL,
  progressbar = TRUE
)

Value

An object of class std_custom. Obtain numeric results using tidy.std_custom. This is a list with the following components:

res_contrast

An unnamed list with one element for each of the requested contrasts. Each element is itself a list with the elements:

B: The number of bootstrap replicates
estimates: Estimated counterfactual means and standard errors for each exposure level
fit_outcome: The estimated regression model for the outcome
estimates_boot: A list of estimates, one for each bootstrap resample
exposure_names: A character vector of the exposure variable names
times: The vector of times at which the calculation is done, if relevant
est_table: Data.frame of the estimates of the contrast with inference
transform: The transform argument used for this contrast
contrast: The requested contrast type
reference: The reference level of the exposure
ci_level: Confidence interval level

res

A named list with the elements:

B: The number of bootstrap replicates

estimates

Estimated counterfactual means and standard errors for each exposure level

fit_outcome

The estimated regression model for the outcome

estimates_boot

A list of estimates, one for each bootstrap resample

exposure_names

A character vector of the exposure variable names

times

The vector of times at which the calculation is done, if relevant

Arguments

fitter_list: The function to call to fit the data (as a list).
arguments: The arguments to be used in the fitter function as a list.
predict_fun_list: The function used to predict the means/probabilities for a new data set on the response level. For survival data, this should be a matrix where each column is the time, and each row the data (as a list).
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.
B: Number of nonparametric bootstrap resamples. Default is NULL (no bootstrap).
ci_level: Coverage probability of confidence intervals.
contrasts: A vector of contrasts in the following format: If set to "difference" or "ratio", then $ψ (x) - ψ (x_{0})$ or $ψ (x) / ψ (x_{0})$ are constructed, where $x_{0}$ is a reference level specified by the reference argument. Has to be NULL if no references are specified.
reference: A vector of reference levels in the following format: If contrasts is not NULL, the desired reference level(s). This must be a vector or list the same length as contrasts, and if not named, it is assumed that the order is as specified in contrasts.
seed: The seed to use with the nonparametric bootstrap.
times: For use with survival data. Set to NULL otherwise.
transforms: A vector of transforms in the following format: If set to "log", "logit", or "odds", the standardized mean $θ (x)$ is transformed into $ψ (x) = \log {θ (x)}$ , $ψ (x) = \log [θ (x) / {1 - θ (x)}]$ , or $ψ (x) = θ (x) / {1 - θ (x)}$ , respectively. If the vector is NULL, then $ψ (x) = θ (x)$ .
progressbar: Logical, if TRUE will print bootstrapping progress to the console

Details

See standardize. The difference is here that different models can be fitted for each value of x in values.

References

Rothman K.J., Greenland S., Lash T.L. (2008). Modern Epidemiology, 3rd edition. Lippincott, Williams & Wilkins.

Sjölander A. (2016). Regression standardization with the R-package stdReg. European Journal of Epidemiology 31(6), 563-574.

Sjölander A. (2016). Estimation of causal effect measures with the R-package stdReg. European Journal of Epidemiology 33(9), 847-858.

Examples

Run this code


require(survival)
prob_predict.coxph <- function(object, newdata, times) {
  fit.detail <- suppressWarnings(basehaz(object))
  cum.haz <- fit.detail$hazard[sapply(times, function(x) max(which(fit.detail$time <= x)))]
  predX <- predict(object = object, newdata = newdata, type = "risk")
  res <- matrix(NA, ncol = length(times), nrow = length(predX))
  for (ti in seq_len(length(times))) {
    res[, ti] <- exp(-predX * cum.haz[ti])
  }
  res
}

set.seed(68)
n <- 500
Z <- rnorm(n)
X <- rbinom(n, 1, prob = 0.5)
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)
x <- standardize_level(
  fitter_list = list("coxph", "coxph"),
  arguments = list(
    list(
      formula = Surv(U, D) ~ X + Z + X * Z,
      method = "breslow",
      x = TRUE,
      y = TRUE
    ),
    list(
      formula = Surv(U, D) ~ X,
      method = "breslow",
      x = TRUE,
      y = TRUE
    )
  ),
  predict_fun_list = list(prob_predict.coxph, prob_predict.coxph),
  data = dd,
  times = seq(1, 5, 0.1),
  values = list(X = c(0, 1)),
  B = 100,
  reference = 0,
  contrasts = "difference"
)
print(x)

Run the code above in your browser using DataLab