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 binary effect modification variables. Finally,
X1 and X2 above represent optional additional
covariates, which can include factor variables and interactions such as
X1*X2. Note that interactions in the modifier term are not
allowed, e.g. M1*M2. When deg,
modifier, and model are not supplied, the
defaults are deg=1, no effect modifiers, and
model="ERR".
The included example dataset contains 3,000 individuals with 10
exposure realizations V1-V10, binary
covariates X1 and X2 and effect modifiers
M1 and M2, and outcomes of all types
(Y.gaussian, Y.binomial,
Y.poisson, status, time,
Y.multinomial, Y.clogit, and
setnr).
data(data, package="ameras")
head(data)
#> 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=data,
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=data, 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=data, 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=data, 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=data, 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=data, family="clogit",
methods=c("RC", "ERC", "MCML", "FMA", "BMA"),
niter.BMA = 5000, nburnin.BMA = 1000)Without the quadratic term (deg=1):