qgcomp: Quantile g-computation for continuous, binary, count, and...

Description Usage Arguments Value See Also Examples

View source: R/base.R

Description

This function automatically selects between qgcomp.noboot, qgcomp.boot, qgcomp.cox.noboot, and qgcomp.cox.boot for the most efficient approach to estimate the average expected change in the (log) outcome per quantile increase in the joint exposure to all exposures in expnms', given the underlying model. For example, if the underlying model (specified by the formula f) is a linear model with all linear terms for exposure, then qgcomp.noboot“ will be called to fit the model. Non-linear terms or requesting the risk ratio for binomial outcomes will result in the qgcomp.boot function being called. For a given linear model, boot and noboot versions will give identical inference, though when using survival outcomes, the 'boot' version uses simulation based inference, which can vary from the 'noboot' version due to simulation error (which can be minimized via setting the MCsize parameter very large - see qgcomp.cox.boot for details).

Usage

1
qgcomp(f, data = data, family = gaussian(), rr = TRUE, ...)

Arguments

f

R style formula (may include survival outcome via Surv)

data

data frame

family

gaussian(), binomial(), cox(), poisson() (works as argument to 'family' parameter in glm‘ or ’dist' parameter in zeroinfl)

rr

logical: if using binary outcome and rr=TRUE, qgcomp.boot will estimate risk ratio rather than odds ratio. Note, to get population average effect estimates for a binary outcome, set rr=TRUE (default: ORs are generally not of interest as population average effects, so if rr=FALSE then a conditional OR will be estimated, which cannot be interpreted as a population average effect

...

arguments to qgcomp.noboot or qgcomp.boot (e.g. q) or glm

Value

a qgcompfit object, which contains information about the effect measure of interest (psi) and associated variance (var.psi), as well as information on the model fit (fit) and possibly information on the marginal structural model (msmfit) used to estimate the final effect estimates (qgcomp.boot, qgcomp.cox.boot only). If appropriate, weights are also reported, which represent the proportion of a directional (positive/negative) effect that is accounted for by each exposure.

See Also

qgcomp.noboot, qgcomp.boot, qgcomp.cox.noboot, qgcomp.cox.boot qgcomp.zi.noboot and qgcomp.zi.boot (qgcomp is just a wrapper for these functions)

Examples

 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
set.seed(50)
dat <- data.frame(y=runif(50), x1=runif(50), x2=runif(50), z=runif(50))
qgcomp.noboot(y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=2)
qgcomp.boot(y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=2, B=10, seed=125)
# automatically selects appropriate method
qgcomp(y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=2)
# note for binary outcome this will choose the risk ratio (and bootstrap methods) by default
dat <- data.frame(y=rbinom(100, 1, 0.5), x1=runif(100), x2=runif(100), z=runif(100))
## Not run: 
qgcomp.noboot(y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=2, family=binomial())
set.seed(1231)
qgcomp.boot(y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=2, family=binomial())
set.seed(1231)
qgcomp(y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=2, family=binomial())

# automatically selects appropriate method when specifying rr or degree explicitly
qgcomp(y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=2, family=binomial(), rr=FALSE)
qgcomp(y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=2, family=binomial(), rr=TRUE)
qgcomp(y ~ z + factor(x1) + factor(x2), degree=2, expnms = c('x1', 'x2'), data=dat, q=4,
family=binomial())

#survival objects
set.seed(50)
N=200
dat <- data.frame(time=(tmg <- pmin(.1,rweibull(N, 10, 0.1))),
                d=1.0*(tmg<0.1), x1=runif(N), x2=runif(N), z=runif(N))
expnms=paste0("x", 1:2)
f = survival::Surv(time, d)~x1 + x2
qgcomp(f, expnms = expnms, data = dat)
# note if B or MCsize are set but the model is linear, an error will result
try(qgcomp(f, expnms = expnms, data = dat, B1=, MCsize))
# note that in the survival models, MCsize should be set to a large number
#  such that results are repeatable (within an error tolerance such as 2 significant digits)
# if you run them under different  seed values
f = survival::Surv(time, d)~x1 + x2 + x1:x2
qgcomp(f, expnms = expnms, data = dat, B=10, MCsize=100)

## End(Not run)

Example output

Scaled effect size (positive direction, sum of positive coefficients = 0.021)
x2 
 1 

Scaled effect size (negative direction, sum of negative coefficients = -0.0641)
x1 
 1 

Mixture slope parameters (Delta method CI):

             Estimate Std. Error Lower CI Upper CI t value  Pr(>|t|)
(Intercept)  0.618002   0.100256  0.42150   0.8145  6.1642 1.518e-07
psi1        -0.043103   0.108728 -0.25621   0.1700 -0.3964    0.6936
Mixture slope parameters (bootstrap CI):

             Estimate Std. Error Lower CI Upper CI t value  Pr(>|t|)
(Intercept)  0.469103   0.050758  0.36962 0.568586  9.2420 5.787e-12
psi1        -0.043103   0.066853 -0.17413 0.087925 -0.6448    0.5224
Scaled effect size (positive direction, sum of positive coefficients = 0.021)
x2 
 1 

Scaled effect size (negative direction, sum of negative coefficients = -0.0641)
x1 
 1 

Mixture slope parameters (Delta method CI):

             Estimate Std. Error Lower CI Upper CI t value  Pr(>|t|)
(Intercept)  0.618002   0.100256  0.42150   0.8145  6.1642 1.518e-07
psi1        -0.043103   0.108728 -0.25621   0.1700 -0.3964    0.6936
Scaled effect size (positive direction, sum of positive coefficients = 0.357)
    x1     x2 
0.9226 0.0774 

Scaled effect size (negative direction, sum of negative coefficients = 0)
None

Mixture log(OR) (Delta method CI):

            Estimate Std. Error Lower CI Upper CI Z value Pr(>|z|)
(Intercept)  0.23283    0.53619 -0.81808   1.2837  0.4342   0.6641
psi1         0.35686    0.54916 -0.71947   1.4332  0.6498   0.5158
Mixture log(RR) (bootstrap CI):

            Estimate Std. Error Lower CI Upper CI Z value  Pr(>|z|)
(Intercept) -0.92267    0.21955 -1.35298 -0.49237 -4.2027 2.638e-05
psi1         0.19417    0.32065 -0.43428  0.82263  0.6056    0.5448
Mixture log(RR) (bootstrap CI):

            Estimate Std. Error Lower CI Upper CI Z value  Pr(>|z|)
(Intercept) -0.92267    0.21955 -1.35298 -0.49237 -4.2027 2.638e-05
psi1         0.19417    0.32065 -0.43428  0.82263  0.6056    0.5448
Scaled effect size (positive direction, sum of positive coefficients = 0.357)
    x1     x2 
0.9226 0.0774 

Scaled effect size (negative direction, sum of negative coefficients = 0)
None

Mixture log(OR) (Delta method CI):

            Estimate Std. Error Lower CI Upper CI Z value Pr(>|z|)
(Intercept)  0.23283    0.53619 -0.81808   1.2837  0.4342   0.6641
psi1         0.35686    0.54916 -0.71947   1.4332  0.6498   0.5158
Mixture log(RR) (bootstrap CI):

            Estimate Std. Error Lower CI Upper CI Z value  Pr(>|z|)
(Intercept) -0.92267    0.23145 -1.37631 -0.46904 -3.9865 6.705e-05
psi1         0.19417    0.35355 -0.49878  0.88713  0.5492    0.5829
Mixture log(RR) (bootstrap CI):

              Estimate Std. Error Lower CI Upper CI Z value Pr(>|z|)
(Intercept) -0.9731701  0.3489822 -1.65716 -0.28918 -2.7886 0.005294
psi1        -0.0057746  0.5643864 -1.11195  1.10040 -0.0102 0.991837
psi2         0.0428528  0.1736242 -0.29744  0.38315  0.2468 0.805053
Scaled effect size (positive direction, sum of positive coefficients = 0)
None

Scaled effect size (negative direction, sum of negative coefficients = -0.244)
   x1    x2 
0.634 0.366 

Mixture log(hazard ratio) (Delta method CI):

     Estimate Std. Error Lower CI  Upper CI Z value Pr(>|z|)
psi1 -0.24398    0.12855 -0.49593 0.0079704  -1.898   0.0577
Error in qgcomp(f, expnms = expnms, data = dat, B1 = , MCsize) : 
  object 'MCsize' not found
Mixture log(hazard ratio) (bootstrap CI):

     Estimate Std. Error Lower CI Upper CI Z value Pr(>|z|)
psi1 -0.25411    0.14869 -0.54554  0.03732  -1.709  0.08746

qgcomp documentation built on Jan. 24, 2022, 5:08 p.m.