modelbound: Estimating qgcomp regression line confidence bounds

Description

Calculates: expected outcome (on the link scale), and upper and lower confidence intervals (both pointwise and simultaneous)

Usage

modelbound(x, emmval = NULL, alpha = 0.05, pwonly = FALSE, ...)

Value

A data frame containing

linpred:: The linear predictor from the marginal structural model
r/o/m:: The canonical measure (risk/odds/mean) for the marginal structural model link
se....:: the stndard error of linpred
ul..../ll....:: Confidence bounds for the effect measure, and bounds centered at the canonical measure (for plotting purposes)

The confidence bounds are either "pointwise" (pw) and "simultaneous" (simul) confidence intervals at each each quantized value of all exposures.

Arguments

x: "qgcompemmfit" object from qgcomp.emm.glm.boot or qgcomp.emm.glm.ee,
emmval: fixed value for effect measure modifier at which pointwise comparisons are calculated
alpha: alpha level for confidence intervals
pwonly: logical: return only pointwise estimates (suppress simultaneous estimates)
...: not used

Details

This method leverages the distribution of qgcomp model coefficients to estimate pointwise regression line confidence bounds. These are defined as the bounds that, for each value of the independent variable X (here, X is the joint exposure quantiles) the 95% bounds (for example) for the model estimate of the regression line E(Y|X) are expected to include the true value of E(Y|X) in 95% of studies. The "simultaneous" bounds are also calculated, and the 95% simultaneous bounds contain the true value of E(Y|X) for all values of X in 95% of studies. The latter are more conservative and account for the multiple testing implied by the former. Pointwise bounds are calculated via the standard error for the estimates of E(Y|X), while the simultaneous bounds are estimated using the method of Cheng (reference below). All bounds are large sample bounds that assume normality and thus will be underconservative in small samples. These bounds may also include illogical values (e.g. values less than 0 for a dichotomous outcome) and should be interpreted cautiously in small samples.

Following Cheng, this approach is possible on bootstrap fitted models (aside from Cox models). For estimating equation approaches, a normality assumption is used to draw values from the sampling distribution of the model parameters based on the covariance matrix of the parameters of the marginal structural model. Because those parameters are not separately estimated for the noboot approaches, this method is not available for noboot methods.

Reference:

Cheng, Russell CH. "Bootstrapping simultaneous confidence bands." Proceedings of the Winter Simulation Conference, 2005.. IEEE, 2005.

Examples

Run this code

if (FALSE) {
set.seed(50)
dat <- data.frame(y=runif(50), x1=runif(50), x2=runif(50),
                  z=rbinom(50, 1, 0.5), r=rbinom(50, 1, 0.5))
(qfit <- qgcomp.emm.glm.noboot(f=y ~ z + x1 + x2, emmvar="z",
                           expnms = c('x1', 'x2'), data=dat, q=4, family=gaussian()))
(qfit2 <- qgcomp.emm.glm.boot(f=y ~ z + x1 + x2, emmvar="z",
                          degree = 1,
                          expnms = c('x1', 'x2'), data=dat, q=4, family=gaussian()))
# modelbound(qfit) # this will error (only works with bootstrapped objects)
 modelbound(qfit2)
# logistic model
set.seed(200)
dat2 <- data.frame(y=rbinom(200, 1, 0.3), x1=runif(200), x2=runif(200),
                  z=rbinom(200, 1, 0.5))
(qfit3 <- qgcomp.emm.glm.boot(f=y ~ z + x1 + x2, emmvar="z",
                          degree = 1,
                          expnms = c('x1', 'x2'), data=dat2, q=4, rr = FALSE, family=binomial()))
modelbound(qfit3)
# risk ratios instead (check for upper bound > 1.0, indicating implausible risk)
(qfit3b <- qgcomp.emm.glm.boot(f=y ~ z + x1 + x2, emmvar="z",
                          degree = 1,
                          expnms = c('x1', 'x2'), data=dat2, q=4, rr = TRUE, family=binomial()))
modelbound(qfit3b)
# categorical modifier
set.seed(50)
dat3 <- data.frame(y=runif(50), x1=runif(50), x2=runif(50),
   z=sample(0:2, 50, replace=TRUE), r=rbinom(50, 1, 0.5))
dat3$z = as.factor(dat3$z)
(qfit4 <- qgcomp.emm.glm.boot(f=y ~ z + x1 + x2, emmvar="z",
                          degree = 1,
                          expnms = c('x1', 'x2'), data=dat3, q=4, family=gaussian()))
modelbound(qfit4, emmval=2)
}