Fitting models and displaying output
In ameras: Analyze Multiple Exposure Realizations in Association Studies

knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)

library(ameras)

Introduction

All functionality of the package is included in the function ameras. This vignette demonstrates the basic functionality and the generated output.

Example data

The included example dataset contains 3,000 individuals with 10 exposure replicates V1-V10, binary covariates X1 and X2 and risk 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)

There are exposure replicates are in columns V1-V10, so define dosevars as follows:

dosevars <- paste0("V", 1:10)

Linear regression & displaying output

Now we fit all methods to the data through one function call:

set.seed(12345)
fit.ameras.linreg <- ameras(Y="Y.gaussian", dosevars=dosevars, X=c("X1","X2"), data=data, 
                            family="gaussian", methods=c("RC", "ERC", "MCML", "FMA", "BMA"), 
                            niter.BMA = 5000, nburnin.BMA = 1000, CI=c("wald.orig","percentile"))

The output is a list object with a call component, and one component for each method, each being a list:

str(fit.ameras.linreg)

Access the results for e.g. regression calibration as follows:

fit.ameras.linreg$RC

For a summary of results, use summary (note that Rhat and n.eff only apply to BMA results):

summary(fit.ameras.linreg)

To extract model coefficients and compare between methods, use coef:

coef(fit.ameras.linreg)

To produce traceplots and density plots for BMA results, use traceplot:

traceplot(fit.ameras.linreg)

Logistic regression

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 doseRRmod="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="Y.binomial", dosevars=dosevars, X=c("X1","X2"), data=data, 
                            family="binomial", deg=2, doseRRmod = "EXP", 
                            methods=c("RC", "ERC", "MCML", "FMA", "BMA"), niter.BMA = 5000, 
                            nburnin.BMA = 1000, CI=c("wald.orig","percentile"))

summary(fit.ameras.logreg)

coef(fit.ameras.logreg)

traceplot(fit.ameras.logreg)

Without the quadratic term (deg=1):

set.seed(3521216)
fit.ameras.logreg.lin <- ameras(Y="Y.binomial", dosevars=dosevars, X=c("X1","X2"), data=data, 
                                family="binomial", deg=1, doseRRmod = "EXP", 
                                methods=c("RC", "ERC", "MCML", "FMA", "BMA"), niter.BMA = 5000, 
                                nburnin.BMA = 1000, CI=c("wald.orig","percentile"))

summary(fit.ameras.logreg.lin)

coef(fit.ameras.logreg.lin)

traceplot(fit.ameras.logreg.lin)

Poisson regression

Again, we first fit models including a quadratic exposure term by setting deg=2.

set.seed(332101)
fit.ameras.poisson <- ameras(Y="Y.poisson", dosevars=dosevars, X=c("X1","X2"), data=data, 
                             family="poisson", deg=2, doseRRmod = "EXP", 
                             methods=c("RC", "ERC", "MCML", "FMA", "BMA"), niter.BMA = 5000, 
                             nburnin.BMA = 1000, CI=c("wald.orig","percentile"))

summary(fit.ameras.poisson)

coef(fit.ameras.poisson)

traceplot(fit.ameras.poisson)

Without the quadratic term (deg=1):

set.seed(24252)
fit.ameras.poisson.lin <- ameras(Y="Y.poisson", dosevars=dosevars, X=c("X1","X2"), data=data, 
                                 family="poisson", deg=1, doseRRmod = "EXP", 
                                 methods=c("RC", "ERC", "MCML", "FMA", "BMA"), niter.BMA = 5000, 
                                 nburnin.BMA = 1000, CI=c("wald.orig","percentile"))

summary(fit.ameras.poisson.lin)

coef(fit.ameras.poisson.lin)

traceplot(fit.ameras.poisson.lin)

Proportional hazards regression

Proportional hazards regression uses the same syntax, with Y the 0-1 status variable and exit for the exit time variable. In case of left truncation, entry times can be specified through the entry argument. 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. This 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(Y="status", exit="time", dosevars=dosevars, X=c("X1","X2"), 
                             data=data, family="prophaz", deg=2, doseRRmod = "EXP", 
                             methods=c("RC", "ERC", "MCML", "FMA", "BMA"), niter.BMA = 5000, 
                             nburnin.BMA = 1000, CI=c("wald.orig","percentile"))

summary(fit.ameras.prophaz)

coef(fit.ameras.prophaz)

The BMA output now contains the intervals with piecewise constant baseline hazards, corresponding to the estimates h0:

fit.ameras.prophaz$BMA$prophaz.timepoints

traceplot(fit.ameras.prophaz)

Without the quadratic term (deg=1):

set.seed(24978252)
fit.ameras.prophaz.lin <- ameras(Y="status", exit="time", dosevars=dosevars, X=c("X1","X2"), 
                                 data=data, family="prophaz", deg=1, doseRRmod = "EXP", 
                                 methods=c("RC", "ERC", "MCML", "FMA", "BMA"), niter.BMA = 5000, 
                                 nburnin.BMA = 1000, CI=c("wald.orig","percentile"))

summary(fit.ameras.prophaz.lin)

coef(fit.ameras.prophaz.lin)

traceplot(fit.ameras.prophaz.lin)

Multinomial logistic regression

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="Y.multinomial", dosevars=dosevars, X=c("X1","X2"), data=data, 
                            family="multinomial", deg=2, doseRRmod = "EXP", 
                            methods=c("RC", "ERC", "MCML", "FMA", "BMA"), niter.BMA = 5000, 
                            nburnin.BMA = 1000, CI=c("wald.orig","percentile"))

summary(fit.ameras.multinomial)

coef(fit.ameras.multinomial)

traceplot(fit.ameras.multinomial)

Without the quadratic term (deg=1):

set.seed(44)
fit.ameras.multinomial.lin <- ameras(Y="Y.multinomial", dosevars=dosevars, X=c("X1","X2"), data=data, 
                            family="multinomial", deg=1, doseRRmod = "EXP", 
                            methods=c("RC","ERC", "MCML", "FMA", "BMA"), niter.BMA = 5000, 
                            nburnin.BMA = 1000, CI=c("wald.orig","percentile"))

summary(fit.ameras.multinomial.lin)

coef(fit.ameras.multinomial.lin)

traceplot(fit.ameras.multinomial.lin)

Conditional logistic regression

Again, we first fit models including a quadratic exposure term by setting deg=2.

set.seed(3301)
fit.ameras.clogit <- ameras(Y="Y.clogit", dosevars=dosevars, X=c("X1","X2"), data=data, 
                            family="clogit", deg=2, doseRRmod = "EXP", setnr="setnr",
                            methods=c("RC", "ERC", "MCML", "FMA", "BMA"), niter.BMA = 5000, 
                            nburnin.BMA = 1000, CI=c("wald.orig","percentile"))

summary(fit.ameras.clogit)

coef(fit.ameras.clogit)

traceplot(fit.ameras.clogit)

Without the quadratic term (deg=1):

set.seed(4401)
fit.ameras.clogit.lin <- ameras(Y="Y.clogit", dosevars=dosevars, X=c("X1","X2"), data=data, 
                            family="clogit", deg=1, doseRRmod = "EXP", setnr="setnr",
                            methods=c("RC","ERC", "MCML", "FMA", "BMA"), niter.BMA = 5000, 
                            nburnin.BMA = 1000, CI=c("wald.orig","percentile"))