stdReg2 (version 1.0.3)

standardize: Get standardized estimates using the g-formula with a custom model

Description

Get standardized estimates using the g-formula with a custom model

Usage

standardize(
  fitter,
  arguments,
  predict_fun,
  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

The function to call to fit the data.

arguments

The arguments to be used in the fitter function as a list.

predict_fun

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.

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 \(\psi(x)-\psi(x_0)\) or \(\psi(x) / \psi(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 \(\theta(x)\) is transformed into \(\psi(x)=\log\{\theta(x)\}\), \(\psi(x)=\log[\theta(x)/\{1-\theta(x)\}]\), or \(\psi(x)=\theta(x)/\{1-\theta(x)\}\), respectively. If the vector is NULL, then \(\psi(x)=\theta(x)\).

progressbar

Logical, if TRUE will print bootstrapping progress to the console

Details

Let \(Y\), \(X\), and \(Z\) be the outcome, the exposure, and a vector of covariates, respectively. standardize uses a model to estimate the standardized mean \(\theta(x)=E\{E(Y|X=x,Z)\}\), where \(x\) is a specific value of \(X\), and the outer expectation is over the marginal distribution of \(Z\). With survival data, \(Y=I(T > t)\), and a vector of different time points times (\(t\)) can be given, where \(T\) is the uncensored survival time.

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

set.seed(6)
n <- 100
Z <- rnorm(n)
X <- rnorm(n, mean = Z)
Y <- rbinom(n, 1, prob = (1 + exp(X + Z))^(-1))
dd <- data.frame(Z, X, Y)
prob_predict.glm <- function(...) predict.glm(..., type = "response")

x <- standardize(
  fitter = "glm",
  arguments = list(
    formula = Y ~ X * Z,
    family = "binomial"
  ),
  predict_fun = prob_predict.glm,
  data = dd,
  values = list(X = seq(-1, 1, 0.1)),
  B = 100,
  reference = 0,
  contrasts = "difference"
)
x

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 <- 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)
x <- standardize(
fitter = "coxph",
  arguments = list(
    formula = Surv(U, D) ~ X + Z + X * Z,
    method = "breslow",
    x = TRUE,
    y = TRUE
  ),
  predict_fun = prob_predict.coxph,
  data = dd,
  times = 1:5,
  values = list(X = c(-1, 0, 1)),
  B = 100,
  reference = 0,
  contrasts = "difference"
)
x