standardize {stdReg2}R Documentation

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
)

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.

Value

An object of class std_custom. This is a list with components estimates and fit for the outcome model.

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

[Package stdReg2 version 1.0.1 Index]