View source: R/standardize_custom.R
standardize | R Documentation |
Get standardized estimates using the g-formula with a custom model
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
)
fitter |
The function to call to fit the data. |
arguments |
The arguments to be used in the fitter function as a |
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 |
ci_level |
Coverage probability of confidence intervals. |
contrasts |
A vector of contrasts in the following format:
If set to |
reference |
A vector of reference levels in the following format:
If |
seed |
The seed to use with the nonparametric bootstrap. |
times |
For use with survival data. Set to |
transforms |
A vector of transforms in the following format:
If set to |
progressbar |
Logical, if TRUE will print bootstrapping progress to the console |
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.
An object of class std_custom
. Obtain numeric results using tidy.std_custom.
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 number of bootstrap replicates
Estimated counterfactual means and standard errors for each exposure level
The estimated regression model for the outcome
A list of estimates, one for each bootstrap resample
A character vector of the exposure variable names
The vector of times at which the calculation is done, if relevant
Data.frame of the estimates of the contrast with inference
The transform argument used for this contrast
The requested contrast type
The reference level of the exposure
Confidence interval level
A named list with the elements:
The number of bootstrap replicates
Estimated counterfactual means and standard errors for each exposure level
The estimated regression model for the outcome
A list of estimates, one for each bootstrap resample
A character vector of the exposure variable names
The vector of times at which the calculation is done, if relevant
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.
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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.