All functionality of the package is included in the function
ameras. This vignette demonstrates the basic functionality
and the generated output. For more information on obtaining
confidence/credible intervals, see Confidence intervals.
Models are specified through formulas of the form
Y~dose(dose_expression, model="ERR", deg=1, modifier=~M1+M2)+X1+X2.
Here dose_expression specifies the dose realization columns
and is parsed by eval_select from the
tidyselect package. Useful examples are
D1:D1000 if the doses are in a sequence of columns with
sequential names such as D1-D1000, and
all_of(dosevars) where dosevars is a vector
with the names of all dose columns. Further, model
specifies, for non-Gaussian families, whether to use the exponential
dose-response model (model="EXP"), the linear-exponential
model (model="LINEXP") or the linear ERR model
(model="ERR"). Next, deg is used to specify
whether a quadratic dose term should (deg=2) or should not
(deg=1) be estimated for the exponential or linear ERR
dose-response model. The modifier term is optional and used
to specify effect modification variables through a formula. Use
modifier=~M for reference-plus-contrast coding, and
modifier=~0+M or modifier=~M-1 for
subgroup-specific dose coefficients. Formula modifiers may be binary
numeric/logical variables or factors; factor modifiers use the existing
factor level order, with the first level as the reference. Interactions,
transformed modifier terms, continuous numeric modifiers, character
modifiers, and no-intercept coding with multiple modifier variables are
not currently supported. For ERR and linear-exponential models,
subgroup-specific modifier coding can be more numerically robust because
default bounds are applied directly to each subgroup dose-effect
parameter. See Effect modification
for more detail. When deg, modifier, and
model are not supplied, the defaults are
deg=1, no effect modifiers, and
model="ERR".
Additional right-hand-side covariates are expanded using
model.matrix(), so standard formula features such as factor
contrasts, interactions, I() terms, and basis functions can
be used. For example, formulas can include terms such as
X1*X2, I(age^2),
splines::ns(age, df = 3), or
splines::bs(age, df = 4). These terms must preserve one
value per row of the input data; terms that drop rows, such as
stats::na.omit(X), are not supported. When
keep.data = FALSE, later calls that use supplied data apply
the same expanded covariate columns as the original fit.
The included ameras_simulated dataset contains 3,000
simulated individuals with 10 exposure realizations
V1-V10, binary covariates X1 and
X2, effect modifiers M1 and M2,
and outcomes of all supported types (Y.gaussian,
Y.binomial, Y.poisson, status,
time, Y.multinomial, Y.clogit,
and setnr). The outcomes were generated under linear excess
relative risk models; see ?ameras_simulated for
details.
head(ameras_simulated)
#> Y.gaussian Y.binomial Y.poisson time status setnr Y.clogit
#> 1 -0.32647093 0 0 0.30276565 0 1 0
#> 2 -0.18734036 1 0 0.19735142 1 1 1
#> 3 0.08404044 0 2 0.30276565 0 1 0
#> 4 0.22432504 0 0 0.23602584 1 1 0
#> 5 -0.46317255 0 0 0.30276565 0 2 0
#> 6 -1.41036573 0 0 0.07838133 1 2 0
#> Y.multinomial X1 X2 M1 M2 V1 V2 V3 V4
#> 1 3 0 0 0 1 0.42868043 0.61542487 0.41960219 0.49265549
#> 2 2 1 0 1 0 0.73321154 0.35512449 0.41876478 0.49235658
#> 3 2 0 0 1 0 0.70369712 0.43407408 1.04115924 0.79882088
#> 4 3 0 0 1 0 0.01845324 0.01373367 0.02733303 0.01912686
#> 5 3 1 0 0 0 0.39389441 0.40087181 0.61932032 0.51715526
#> 6 2 1 1 1 0 0.01493158 0.02335143 0.01828983 0.02705350
#> V5 V6 V7 V8 V9 V10
#> 1 0.31363762 0.42218455 0.42464021 0.29630858 0.38211182 0.45751570
#> 2 0.49515815 0.56837639 0.61126842 0.67723449 0.53361810 0.49393510
#> 3 0.66613754 0.72346942 0.64077434 0.79894534 0.98278177 1.06068250
#> 4 0.01917956 0.03056413 0.01536966 0.02135999 0.01548655 0.01596626
#> 5 0.36440322 0.60255525 0.47512525 0.52567606 0.53391825 0.56026531
#> 6 0.02298922 0.02399258 0.01890339 0.02094013 0.02303085 0.02091743Now we fit all methods to the data through one function call:
set.seed(12345)
fit.ameras.linreg <- ameras(Y.gaussian~dose(V1:V10)+X1+X2, data = ameras_simulated,
family="gaussian", niter.BMA=5000, nburnin.BMA=1000,
methods=c("RC", "ERC", "MCML", "FMA", "BMA"))The output is a list object with a call component, and
one component for each method, each being a list:
Access the results for e.g. regression calibration as follows:
For a summary of results, use summary (note that
Rhat and n.eff only apply to BMA results):
To extract model coefficients and compare between methods, use
coef:
To produce traceplots and density plots for BMA results, use
traceplot:
To fit a logistic regression model, the syntax is very similar.
However, ameras offers three functional forms for modeling
the exposure-outcome relationship. Here, we will illustrate the standard
exponential relationship model="EXP". For more information,
see the vignette Relative risk
models.
First, we fit models including a quadratic exposure term by setting
deg=2.
set.seed(33521)
fit.ameras.logreg <- ameras(Y.binomial~dose(V1:V10, deg=2, model="EXP")+X1+X2,
data = ameras_simulated, family="binomial",
methods=c("RC", "ERC", "MCML", "FMA", "BMA"),
niter.BMA = 5000, nburnin.BMA = 1000)Without the quadratic term (deg=1):
For Poisson regression, an offset can be optionally used by including
an offset term in the formula, e.g.,
Y~dose(V1:V10)+X1+offset(PYR). Here, we show models without
using an offset, first including a quadratic exposure term:
set.seed(332101)
fit.ameras.poisson <- ameras(Y.poisson~dose(V1:V10, deg=2, model="EXP")+X1+X2,
data = ameras_simulated, family="poisson",
methods=c("RC", "ERC", "MCML", "FMA", "BMA"),
niter.BMA = 5000, nburnin.BMA = 1000)Without the quadratic term (deg=1):
Proportional hazards regression uses syntax similar to the
survival package, with models specified using formulas with
Surv(exit, status) or
Surv(entry, exit, status) on the left-hand side. Note that
BMA fits a piecewise constant baseline hazard h0 as the
proportional hazards model is not directly supported. By default, the
observed time interval is divided into 10 intervals using quantiles of
the observed event times among cases. The number of such intervals can
be specified through the prophaz.numints.BMA argument. The
BMA output contains the prophaz.numints.BMA+1 cutpoints
defining the intervals in addition to h0.
Again, we first fit models including a quadratic exposure term by
setting deg=2.
set.seed(332120000)
fit.ameras.prophaz <- ameras(Surv(time, status)~
dose(V1:V10, deg=2, model="EXP")+X1+X2,
data = ameras_simulated, family="prophaz",
methods=c("RC", "ERC", "MCML", "FMA", "BMA"),
niter.BMA = 5000, nburnin.BMA = 1000)The BMA output now contains the intervals with piecewise constant
baseline hazards, corresponding to the estimates h0:
Without the quadratic term (deg=1):
For multinomial logistic regression, the last category (in the case
of the example data, Y.multinomial='3') is used as the
referent category.
Again, we first fit models including a quadratic exposure term by
setting deg=2.
set.seed(33)
fit.ameras.multinomial <- ameras(Y.multinomial~
dose(V1:V10, deg=2, model="EXP")+X1+X2,
data = ameras_simulated, family="multinomial",
methods=c("RC", "ERC", "MCML", "FMA", "BMA"),
niter.BMA = 5000, nburnin.BMA = 1000)Without the quadratic term (deg=1):
For conditional logistic regression, the set number variable is
specified through a strata term in the formula. Again, we
first fit models including a quadratic exposure term by setting
deg=2.
set.seed(3301)
fit.ameras.clogit <- ameras(Y.clogit~dose(V1:V10, deg=2, model="EXP")+X1+X2+
strata(setnr), data = ameras_simulated, family="clogit",
methods=c("RC", "ERC", "MCML", "FMA", "BMA"),
niter.BMA = 5000, nburnin.BMA = 1000)Without the quadratic term (deg=1):