--- title: "Fitting models and displaying output" output: rmarkdown::html_vignette: toc: true vignette: > %\VignetteIndexEntry{Fitting models and displaying output} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", message = FALSE ) ``` ```{r setup} library(ameras) ``` # Introduction 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](confidenceintervals.html). ## Model specification 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"`. ## Example data The included example dataset contains 3,000 individuals with 10 exposure replicates `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`). ```{r} data(data, package="ameras") head(data) ``` # Linear regression & displaying output Now we fit all methods to the data through one function call: ```{r modelfit.linreg, cache=TRUE, eval = identical(Sys.getenv("NOT_CRAN"), "true")} 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: ```{r, eval = identical(Sys.getenv("NOT_CRAN"), "true")} str(fit.ameras.linreg) ``` Access the results for e.g. regression calibration as follows: ```{r, eval = identical(Sys.getenv("NOT_CRAN"), "true")} fit.ameras.linreg$RC ``` For a summary of results, use `summary` (note that `Rhat` and `n.eff` only apply to BMA results): ```{r, eval = identical(Sys.getenv("NOT_CRAN"), "true")} summary(fit.ameras.linreg) ``` To extract model coefficients and compare between methods, use `coef`: ```{r, eval = identical(Sys.getenv("NOT_CRAN"), "true")} coef(fit.ameras.linreg) ``` To produce traceplots and density plots for BMA results, use `traceplot`: ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} 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 `model="EXP"`. For more information, see the vignette [Relative risk models](relativeriskmodels.html). First, we fit models including a quadratic exposure term by setting `deg=2`. ```{r modelfit.logreg, cache=TRUE, eval = identical(Sys.getenv("NOT_CRAN"), "true")} 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) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} summary(fit.ameras.logreg) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} coef(fit.ameras.logreg) ``` Without the quadratic term (`deg=1`): ```{r modelfit.logreg.lin, cache=TRUE, eval = identical(Sys.getenv("NOT_CRAN"), "true")} set.seed(3521216) fit.ameras.logreg.lin <- ameras(Y.binomial~dose(V1:V10, deg=1, model="EXP")+X1+X2, data=data, family="binomial", methods=c("RC", "ERC", "MCML", "FMA", "BMA"), niter.BMA = 5000, nburnin.BMA = 1000) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} summary(fit.ameras.logreg.lin) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} coef(fit.ameras.logreg.lin) ``` # Poisson regression 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: ```{r modelfit.poisson, cache=TRUE, eval = identical(Sys.getenv("NOT_CRAN"), "true")} 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) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} summary(fit.ameras.poisson) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} coef(fit.ameras.poisson) ``` Without the quadratic term (`deg=1`): ```{r modelfit.poisson.lin, cache=TRUE, eval = identical(Sys.getenv("NOT_CRAN"), "true")} set.seed(24252) fit.ameras.poisson.lin <- ameras(Y.poisson~dose(V1:V10, deg=1, model="EXP")+X1+X2, data=data, family="poisson", methods=c("RC", "ERC", "MCML", "FMA", "BMA"), niter.BMA = 5000, nburnin.BMA = 1000) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} summary(fit.ameras.poisson.lin) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} coef(fit.ameras.poisson.lin) ``` # Proportional hazards regression 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. 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`. ```{r modelfit.prophaz, cache=TRUE, eval = identical(Sys.getenv("NOT_CRAN"), "true")} 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) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} summary(fit.ameras.prophaz) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} coef(fit.ameras.prophaz) ``` The BMA output now contains the intervals with piecewise constant baseline hazards, corresponding to the estimates `h0`: ```{r, eval = identical(Sys.getenv("NOT_CRAN"), "true")} fit.ameras.prophaz$BMA$prophaz.timepoints ``` Without the quadratic term (`deg=1`): ```{r modelfit.prophaz.lin, cache=TRUE, eval = identical(Sys.getenv("NOT_CRAN"), "true")} set.seed(24978252) fit.ameras.prophaz.lin <- ameras(Surv(time, status)~ dose(V1:V10, deg=1, model="EXP")+X1+X2, data=data, family="prophaz", methods=c("RC", "ERC", "MCML", "FMA", "BMA"), niter.BMA = 5000, nburnin.BMA = 1000) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} summary(fit.ameras.prophaz.lin) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} coef(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`. ```{r modelfit.multinomial, cache=TRUE, eval = identical(Sys.getenv("NOT_CRAN"), "true")} 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) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} summary(fit.ameras.multinomial) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} coef(fit.ameras.multinomial) ``` Without the quadratic term (`deg=1`): ```{r modelfit.multinomial.lin, cache=TRUE, eval = identical(Sys.getenv("NOT_CRAN"), "true")} set.seed(44) fit.ameras.multinomial.lin <- ameras(Y.multinomial~ dose(V1:V10, deg=1, model="EXP")+X1+X2, data=data, family="multinomial", methods=c("RC", "ERC", "MCML", "FMA", "BMA"), niter.BMA = 5000, nburnin.BMA = 1000) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} summary(fit.ameras.multinomial.lin) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} coef(fit.ameras.multinomial.lin) ``` # Conditional logistic regression 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`. ```{r modelfit.clogit, cache=TRUE, eval = identical(Sys.getenv("NOT_CRAN"), "true")} 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) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} summary(fit.ameras.clogit) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} coef(fit.ameras.clogit) ``` Without the quadratic term (`deg=1`): ```{r modelfit.clogit.lin, cache=TRUE, eval = identical(Sys.getenv("NOT_CRAN"), "true")} set.seed(4401) fit.ameras.clogit.lin <- 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) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} summary(fit.ameras.clogit.lin) ``` ```{r, fig.fullwidth=TRUE, fig.show="hold", out.width='100%', fig.width=6, fig.height=8, eval = identical(Sys.getenv("NOT_CRAN"), "true")} coef(fit.ameras.clogit.lin) ```