Get standardized estimates using the g-formula with a custom model
Source:R/standardize_custom.R
standardize.RdGet 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 thereferenceargument. Has to beNULLif no references are specified.- reference
A vector of reference levels in the following format: If
contrastsis notNULL, the desired reference level(s). This must be a vector or list the same length ascontrasts, 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
NULLotherwise.- 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 isNULL, then \(\psi(x)=\theta(x)\).- progressbar
Logical, if TRUE will print bootstrapping progress to the console
Value
An object of class std_custom.
This is a list with components estimates and fit for the outcome model.
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"
)
#>
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 5%
|
|=== | 6%
|
|==== | 7%
|
|==== | 8%
|
|===== | 9%
|
|===== | 10%
|
|====== | 11%
|
|====== | 12%
|
|======= | 13%
|
|======= | 14%
|
|======== | 15%
|
|======== | 16%
|
|========= | 17%
|
|========= | 18%
|
|========== | 19%
|
|========== | 20%
|
|=========== | 21%
|
|=========== | 22%
|
|============ | 23%
|
|============ | 24%
|
|============= | 25%
|
|============= | 26%
|
|============== | 27%
|
|============== | 28%
|
|=============== | 29%
|
|=============== | 30%
|
|================ | 31%
|
|================ | 32%
|
|================= | 33%
|
|================= | 34%
|
|================== | 35%
|
|================== | 36%
|
|=================== | 37%
|
|=================== | 38%
|
|==================== | 39%
|
|==================== | 40%
|
|===================== | 41%
|
|===================== | 42%
|
|====================== | 43%
|
|====================== | 44%
|
|======================= | 45%
|
|======================= | 46%
|
|======================== | 47%
|
|======================== | 48%
|
|========================= | 49%
|
|========================= | 51%
|
|========================== | 52%
|
|========================== | 53%
|
|=========================== | 54%
|
|=========================== | 55%
|
|============================ | 56%
|
|============================ | 57%
|
|============================= | 58%
|
|============================= | 59%
|
|============================== | 60%
|
|============================== | 61%
|
|=============================== | 62%
|
|=============================== | 63%
|
|================================ | 64%
|
|================================ | 65%
|
|================================= | 66%
|
|================================= | 67%
|
|================================== | 68%
|
|================================== | 69%
|
|=================================== | 70%
|
|=================================== | 71%
|
|==================================== | 72%
|
|==================================== | 73%
|
|===================================== | 74%
|
|===================================== | 75%
|
|====================================== | 76%
|
|====================================== | 77%
|
|======================================= | 78%
|
|======================================= | 79%
|
|======================================== | 80%
|
|======================================== | 81%
|
|========================================= | 82%
|
|========================================= | 83%
|
|========================================== | 84%
|
|========================================== | 85%
|
|=========================================== | 86%
|
|=========================================== | 87%
|
|============================================ | 88%
|
|============================================ | 89%
|
|============================================= | 90%
|
|============================================= | 91%
|
|============================================== | 92%
|
|============================================== | 93%
|
|=============================================== | 94%
|
|=============================================== | 95%
|
|================================================ | 96%
|
|================================================ | 97%
|
|================================================= | 98%
|
|================================================= | 99%
|
|==================================================| 100%
x
#> Number of bootstraps: 100
#> Confidence intervals are based on percentile bootstrap confidence intervals
#>
#> Exposure: X
#> Tables:
#>
#> X Estimate lower.0.95 upper.0.95
#> 1 -1.0 0.689 0.537 0.882
#> 2 -0.9 0.671 0.526 0.865
#> 3 -0.8 0.653 0.521 0.845
#> 4 -0.7 0.635 0.515 0.822
#> 5 -0.6 0.617 0.508 0.797
#> 6 -0.5 0.600 0.493 0.770
#> 7 -0.4 0.583 0.476 0.745
#> 8 -0.3 0.566 0.466 0.720
#> 9 -0.2 0.550 0.450 0.699
#> 10 -0.1 0.534 0.429 0.679
#> 11 0.0 0.519 0.409 0.658
#> 12 0.1 0.504 0.389 0.638
#> 13 0.2 0.490 0.369 0.621
#> 14 0.3 0.476 0.348 0.608
#> 15 0.4 0.462 0.322 0.595
#> 16 0.5 0.449 0.298 0.587
#> 17 0.6 0.437 0.275 0.583
#> 18 0.7 0.425 0.254 0.577
#> 19 0.8 0.413 0.235 0.571
#> 20 0.9 0.401 0.222 0.561
#> 21 1.0 0.391 0.207 0.551
#>
#> Reference level: = 0
#> Contrast: difference
#> Exposure Estimate lower.0.95 upper.0.95
#> 1 -1.0 0.1697 0.04037 0.30879
#> 2 -0.9 0.1516 0.03562 0.28505
#> 3 -0.8 0.1335 0.03104 0.25348
#> 4 -0.7 0.1156 0.02663 0.22134
#> 5 -0.6 0.0980 0.02238 0.19201
#> 6 -0.5 0.0807 0.01829 0.16211
#> 7 -0.4 0.0637 0.01435 0.13054
#> 8 -0.3 0.0472 0.01056 0.09810
#> 9 -0.2 0.0310 0.00691 0.06525
#> 10 -0.1 0.0153 0.00339 0.03243
#> 11 0.0 0.0000 0.00000 0.00000
#> 12 0.1 -0.0148 -0.03177 -0.00326
#> 13 0.2 -0.0292 -0.06286 -0.00641
#> 14 0.3 -0.0432 -0.09365 -0.00944
#> 15 0.4 -0.0566 -0.12370 -0.01237
#> 16 0.5 -0.0697 -0.15192 -0.01519
#> 17 0.6 -0.0823 -0.17663 -0.01791
#> 18 0.7 -0.0945 -0.19954 -0.02054
#> 19 0.8 -0.1062 -0.22074 -0.02308
#> 20 0.9 -0.1176 -0.24033 -0.02553
#> 21 1.0 -0.1285 -0.25891 -0.02790
#>
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"
)
#>
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 5%
|
|=== | 6%
|
|==== | 7%
|
|==== | 8%
|
|===== | 9%
|
|===== | 10%
|
|====== | 11%
|
|====== | 12%
|
|======= | 13%
|
|======= | 14%
|
|======== | 15%
|
|======== | 16%
|
|========= | 17%
|
|========= | 18%
|
|========== | 19%
|
|========== | 20%
|
|=========== | 21%
|
|=========== | 22%
|
|============ | 23%
|
|============ | 24%
|
|============= | 25%
|
|============= | 26%
|
|============== | 27%
|
|============== | 28%
|
|=============== | 29%
|
|=============== | 30%
|
|================ | 31%
|
|================ | 32%
|
|================= | 33%
|
|================= | 34%
|
|================== | 35%
|
|================== | 36%
|
|=================== | 37%
|
|=================== | 38%
|
|==================== | 39%
|
|==================== | 40%
|
|===================== | 41%
|
|===================== | 42%
|
|====================== | 43%
|
|====================== | 44%
|
|======================= | 45%
|
|======================= | 46%
|
|======================== | 47%
|
|======================== | 48%
|
|========================= | 49%
|
|========================= | 51%
|
|========================== | 52%
|
|========================== | 53%
|
|=========================== | 54%
|
|=========================== | 55%
|
|============================ | 56%
|
|============================ | 57%
|
|============================= | 58%
|
|============================= | 59%
|
|============================== | 60%
|
|============================== | 61%
|
|=============================== | 62%
|
|=============================== | 63%
|
|================================ | 64%
|
|================================ | 65%
|
|================================= | 66%
|
|================================= | 67%
|
|================================== | 68%
|
|================================== | 69%
|
|=================================== | 70%
|
|=================================== | 71%
|
|==================================== | 72%
|
|==================================== | 73%
|
|===================================== | 74%
|
|===================================== | 75%
|
|====================================== | 76%
|
|====================================== | 77%
|
|======================================= | 78%
|
|======================================= | 79%
|
|======================================== | 80%
|
|======================================== | 81%
|
|========================================= | 82%
|
|========================================= | 83%
|
|========================================== | 84%
|
|========================================== | 85%
|
|=========================================== | 86%
|
|=========================================== | 87%
|
|============================================ | 88%
|
|============================================ | 89%
|
|============================================= | 90%
|
|============================================= | 91%
|
|============================================== | 92%
|
|============================================== | 93%
|
|=============================================== | 94%
|
|=============================================== | 95%
|
|================================================ | 96%
|
|================================================ | 97%
|
|================================================= | 98%
|
|================================================= | 99%
|
|==================================================| 100%
x
#> Number of bootstraps: 100
#> Confidence intervals are based on percentile bootstrap confidence intervals
#>
#> Exposure: X
#> Tables:
#>
#> Time: 1
#> X estimate lower upper
#> 1 -1 0.6691328 0.6078557 0.7171407
#> 2 0 0.3632410 0.3011138 0.4060458
#> 3 1 0.2466151 0.1878811 0.2902875
#>
#> Time: 2
#> X estimate lower upper
#> 1 -1 0.4676722 0.3847559 0.5519054
#> 2 0 0.2135942 0.1659279 0.2611208
#> 3 1 0.1656470 0.1168408 0.2079257
#>
#> Time: 3
#> X estimate lower upper
#> 1 -1 0.3505275 0.26905828 0.4459009
#> 2 0 0.1520739 0.11902194 0.2044295
#> 3 1 0.1315275 0.09380255 0.1761028
#>
#> Time: 4
#> X estimate lower upper
#> 1 -1 0.2546911 0.12386946 0.3508293
#> 2 0 0.1101689 0.06701637 0.1505347
#> 3 1 0.1069771 0.06241746 0.1430167
#>
#> Time: 5
#> X estimate lower upper
#> 1 -1 0.16309819 0.07709911 0.2613375
#> 2 0 0.07480425 0.03816440 0.1127717
#> 3 1 0.08456621 0.04640860 0.1207358
#>
#>
#> Reference level: = 0
#> Contrast: difference
#> Time: 1
#> Exposure estimate lower upper
#> 1 -1 0.3058918 0.2515877 0.35267464
#> 2 0 0.0000000 0.0000000 0.00000000
#> 3 1 -0.1166259 -0.1534746 -0.08540127
#>
#> Time: 2
#> Exposure estimate lower upper
#> 1 -1 0.25407801 0.18502799 0.33076959
#> 2 0 0.00000000 0.00000000 0.00000000
#> 3 1 -0.04794719 -0.07480939 -0.02350433
#>
#> Time: 3
#> Exposure estimate lower upper
#> 1 -1 0.19845364 0.12283854 0.284032243
#> 2 0 0.00000000 0.00000000 0.000000000
#> 3 1 -0.02054638 -0.04719901 0.001297897
#>
#> Time: 4
#> Exposure estimate lower upper
#> 1 -1 0.144522281 0.04976681 0.22983546
#> 2 0 0.000000000 0.00000000 0.00000000
#> 3 1 -0.003191768 -0.02345454 0.01911037
#>
#> Time: 5
#> Exposure estimate lower upper
#> 1 -1 0.088293936 0.018981386 0.17857030
#> 2 0 0.000000000 0.000000000 0.00000000
#> 3 1 0.009761965 -0.009490497 0.02585953
#>
#>