# Fitting models and displaying output

    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)
    #>    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.02091743

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"))
    #> Note: BMA may require extensive computation time in the order of multiple hours
    #> Fitting RC
    #> Fitting ERC
    #> Fitting MCML
    #> Fitting FMA
    #> Fitting BMA
    #> Defining model
    #> Building model
    #> Setting data and initial values
    #> Running calculate on model
    #>   [Note] Any error reports that follow may simply reflect missing values in model variables.
    #> Checking model sizes and dimensions
    #>   [Note] This model is not fully initialized. This is not an error.
    #>          To see which variables are not initialized, use model$initializeInfo().
    #>          For more information on model initialization, see help(modelInitialization).
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> running chain 1...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|
    #> running chain 2...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|

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

    str(fit.ameras.linreg)
    #> List of 6
    #>  $ call: language ameras(data = data, family = "gaussian", Y = "Y.gaussian", dosevars = dosevars,      X = c("X1", "X2"), methods =| __truncated__ ...
    #>  $ RC  :List of 8
    #>   ..$ coefficients     : Named num [1:5] -1.362 0.481 -0.519 1.166 1.108
    #>   .. ..- attr(*, "names")= chr [1:5] "(Intercept)" "X1" "X2" "dose" ...
    #>   ..$ sd               : Named num [1:5] 0.0367 0.0405 0.0497 0.0204 0.0143
    #>   .. ..- attr(*, "names")= chr [1:5] "(Intercept)" "X1" "X2" "dose" ...
    #>   ..$ vcov             : num [1:5, 1:5] 1.34e-03 -8.36e-04 -4.86e-04 -4.15e-04 1.55e-09 ...
    #>   .. ..- attr(*, "dimnames")=List of 2
    #>   .. .. ..$ : chr [1:5] "(Intercept)" "X1" "X2" "dose" ...
    #>   .. .. ..$ : chr [1:5] "(Intercept)" "X1" "X2" "dose" ...
    #>   ..$ CI               :'data.frame':    5 obs. of  2 variables:
    #>   .. ..$ lower: num [1:5] -1.434 0.401 -0.616 1.126 1.08
    #>   .. ..$ upper: num [1:5] -1.29 0.56 -0.421 1.206 1.136
    #>   ..$ convergence.optim: int 0
    #>   ..$ counts.optim     : Named num [1:2] 509 1
    #>   .. ..- attr(*, "names")= chr [1:2] "function" "gradient"
    #>   ..$ loglik           : num -4563
    #>   ..$ runtime          : chr "0.1 seconds"
    #>  $ ERC :List of 8
    #>   ..$ coefficients     : Named num [1:5] -1.361 0.48 -0.519 1.165 1.106
    #>   .. ..- attr(*, "names")= chr [1:5] "(Intercept)" "X1" "X2" "dose" ...
    #>   ..$ sd               : Named num [1:5] 0.0366 0.0404 0.0497 0.0204 0.0142
    #>   .. ..- attr(*, "names")= chr [1:5] "(Intercept)" "X1" "X2" "dose" ...
    #>   ..$ vcov             : num [1:5, 1:5] 1.34e-03 -8.34e-04 -4.85e-04 -4.14e-04 1.53e-06 ...
    #>   .. ..- attr(*, "dimnames")=List of 2
    #>   .. .. ..$ : chr [1:5] "(Intercept)" "X1" "X2" "dose" ...
    #>   .. .. ..$ : chr [1:5] "(Intercept)" "X1" "X2" "dose" ...
    #>   ..$ CI               :'data.frame':    5 obs. of  2 variables:
    #>   .. ..$ lower: num [1:5] -1.433 0.4 -0.616 1.125 1.078
    #>   .. ..$ upper: num [1:5] -1.289 0.559 -0.421 1.205 1.134
    #>   ..$ convergence.optim: int 0
    #>   ..$ counts.optim     : Named num [1:2] 546 9
    #>   .. ..- attr(*, "names")= chr [1:2] "function" "gradient"
    #>   ..$ loglik           : num -4559
    #>   ..$ runtime          : chr "67.2 seconds"
    #>  $ MCML:List of 8
    #>   ..$ coefficients     : Named num [1:5] -1.28 0.484 -0.517 1.079 1.138
    #>   .. ..- attr(*, "names")= chr [1:5] "(Intercept)" "X1" "X2" "dose" ...
    #>   ..$ sd               : Named num [1:5] 0.0371 0.0416 0.0511 0.0199 0.0147
    #>   .. ..- attr(*, "names")= chr [1:5] "(Intercept)" "X1" "X2" "dose" ...
    #>   ..$ CI               :'data.frame':    5 obs. of  2 variables:
    #>   .. ..$ lower: num [1:5] -1.353 0.402 -0.617 1.04 1.109
    #>   .. ..$ upper: num [1:5] -1.207 0.565 -0.417 1.118 1.167
    #>   ..$ vcov             : num [1:5, 1:5] 1.38e-03 -8.83e-04 -5.14e-04 -3.98e-04 -1.91e-08 ...
    #>   .. ..- attr(*, "dimnames")=List of 2
    #>   .. .. ..$ : chr [1:5] "(Intercept)" "X1" "X2" "dose" ...
    #>   .. .. ..$ : chr [1:5] "(Intercept)" "X1" "X2" "dose" ...
    #>   ..$ convergence.optim: int 0
    #>   ..$ counts.optim     : Named num [1:2] 555 7
    #>   .. ..- attr(*, "names")= chr [1:2] "function" "gradient"
    #>   ..$ loglik           : num -4646
    #>   ..$ runtime          : chr "0.4 seconds"
    #>  $ FMA :List of 6
    #>   ..$ coefficients       : Named num [1:5] -1.28 0.484 -0.517 1.079 1.138
    #>   .. ..- attr(*, "names")= chr [1:5] "(Intercept)" "X1" "X2" "dose" ...
    #>   ..$ sd                 : Named num [1:5] 0.0373 0.0418 0.0512 0.0202 0.0147
    #>   .. ..- attr(*, "names")= chr [1:5] "(Intercept)" "X1" "X2" "dose" ...
    #>   ..$ CI                 :'data.frame':  5 obs. of  2 variables:
    #>   .. ..$ lower: num [1:5] -1.353 0.402 -0.618 1.04 1.109
    #>   .. ..$ upper: num [1:5] -1.206 0.565 -0.417 1.119 1.166
    #>   ..$ included.replicates: int [1:10] 1 2 3 4 5 6 7 8 9 10
    #>   ..$ included.samples   : int 100006
    #>   ..$ runtime            : chr "0.9 seconds"
    #>  $ BMA :List of 7
    #>   ..$ coefficients       : Named num [1:5] -1.28 0.483 -0.517 1.079 1.139
    #>   .. ..- attr(*, "names")= chr [1:5] "(Intercept)" "X1" "X2" "dose" ...
    #>   ..$ sd                 : Named num [1:5] 0.038 0.0424 0.0525 0.0202 0.015
    #>   .. ..- attr(*, "names")= chr [1:5] "(Intercept)" "X1" "X2" "dose" ...
    #>   ..$ CI                 :'data.frame':  5 obs. of  2 variables:
    #>   .. ..$ lower: num [1:5] -1.351 0.402 -0.623 1.039 1.111
    #>   .. ..$ upper: num [1:5] -1.208 0.567 -0.417 1.118 1.168
    #>   ..$ Rhat               :'data.frame':  5 obs. of  2 variables:
    #>   .. ..$ Rhat : num [1:5] 1 1.01 1 1 1.02
    #>   .. ..$ n.eff: num [1:5] 871 822 800 800 744
    #>   ..$ samples            :List of 2
    #>   .. ..$ chain1: num [1:400, 1:6] -1.32 -1.27 -1.33 -1.29 -1.33 ...
    #>   .. .. ..- attr(*, "dimnames")=List of 2
    #>   .. .. .. ..$ : NULL
    #>   .. .. .. ..$ : chr [1:6] "(Intercept)" "X1" "X2" "dose" ...
    #>   .. ..$ chain2: num [1:400, 1:6] -1.3 -1.29 -1.23 -1.3 -1.26 ...
    #>   .. .. ..- attr(*, "dimnames")=List of 2
    #>   .. .. .. ..$ : NULL
    #>   .. .. .. ..$ : chr [1:6] "(Intercept)" "X1" "X2" "dose" ...
    #>   ..$ included.replicates: int [1:10] 1 2 3 4 5 6 7 8 9 10
    #>   ..$ runtime            : chr "39.9 seconds"
    #>  - attr(*, "class")= chr "amerasfit"

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

    fit.ameras.linreg$RC
    #> $coefficients
    #> (Intercept)          X1          X2        dose       sigma
    #>  -1.3622239   0.4807463  -0.5187877   1.1660368   1.1075696
    #>
    #> $sd
    #> (Intercept)          X1          X2        dose       sigma
    #>  0.03665711  0.04045137  0.04971981  0.02039612  0.01429864
    #>
    #> $vcov
    #>               (Intercept)            X1            X2          dose
    #> (Intercept)  1.343744e-03 -8.358763e-04 -4.863064e-04 -4.154402e-04
    #> X1          -8.358763e-04  1.636314e-03 -1.376614e-05  1.526165e-05
    #> X2          -4.863064e-04 -1.376614e-05  2.472060e-03 -2.443289e-05
    #> dose        -4.154402e-04  1.526165e-05 -2.443289e-05  4.160019e-04
    #> sigma        1.547492e-09 -1.113331e-08  8.801420e-09  2.066813e-09
    #>                     sigma
    #> (Intercept)  1.547492e-09
    #> X1          -1.113331e-08
    #> X2           8.801420e-09
    #> dose         2.066813e-09
    #> sigma        2.044510e-04
    #>
    #> $CI
    #>                  lower      upper
    #> (Intercept) -1.4340719 -1.2903760
    #> X1           0.4014616  0.5600310
    #> X2          -0.6162385 -0.4213368
    #> dose         1.1260604  1.2060132
    #> sigma        1.0795443  1.1355949
    #>
    #> $convergence.optim
    #> [1] 0
    #>
    #> $counts.optim
    #> function gradient
    #>      509        1
    #>
    #> $loglik
    #> [1] -4563.325
    #>
    #> $runtime
    #> [1] "0.1 seconds"

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

    summary(fit.ameras.linreg)
    #> Call:
    #> ameras(data = data, family = "gaussian", Y = "Y.gaussian", dosevars = dosevars,
    #>     X = c("X1", "X2"), methods = c("RC", "ERC", "MCML", "FMA",
    #>         "BMA"), CI = c("wald.orig", "percentile"), nburnin.BMA = 1000,
    #>     niter.BMA = 5000)
    #>
    #> Total run time: 108.5 seconds
    #>
    #> Runtime in seconds by method:
    #>
    #>  Method Runtime
    #>      RC     0.1
    #>     ERC    67.2
    #>    MCML     0.4
    #>     FMA     0.9
    #>     BMA    39.9
    #>
    #> Summary of coefficients by method:
    #>
    #>  Method        Term Estimate      SE CI.lowerbound CI.upperbound Rhat  n.eff
    #>      RC (Intercept)  -1.3622 0.03666       -1.4341       -1.2904   NA     NA
    #>      RC          X1   0.4807 0.04045        0.4015        0.5600   NA     NA
    #>      RC          X2  -0.5188 0.04972       -0.6162       -0.4213   NA     NA
    #>      RC        dose   1.1660 0.02040        1.1261        1.2060   NA     NA
    #>      RC       sigma   1.1076 0.01430        1.0795        1.1356   NA     NA
    #>     ERC (Intercept)  -1.3611 0.03660       -1.4328       -1.2893   NA     NA
    #>     ERC          X1   0.4795 0.04041        0.4003        0.5587   NA     NA
    #>     ERC          X2  -0.5188 0.04969       -0.6162       -0.4214   NA     NA
    #>     ERC        dose   1.1649 0.02038        1.1250        1.2049   NA     NA
    #>     ERC       sigma   1.1062 0.01425        1.0783        1.1342   NA     NA
    #>    MCML (Intercept)  -1.2803 0.03714       -1.3530       -1.2075   NA     NA
    #>    MCML          X1   0.4837 0.04156        0.4022        0.5651   NA     NA
    #>    MCML          X2  -0.5171 0.05108       -0.6173       -0.4170   NA     NA
    #>    MCML        dose   1.0790 0.01993        1.0400        1.1181   NA     NA
    #>    MCML       sigma   1.1378 0.01469        1.1090        1.1666   NA     NA
    #>     FMA (Intercept)  -1.2803 0.03732       -1.3534       -1.2065   NA     NA
    #>     FMA          X1   0.4837 0.04181        0.4017        0.5652   NA     NA
    #>     FMA          X2  -0.5173 0.05124       -0.6183       -0.4172   NA     NA
    #>     FMA        dose   1.0792 0.02015        1.0397        1.1187   NA     NA
    #>     FMA       sigma   1.1377 0.01467        1.1091        1.1664   NA     NA
    #>     BMA (Intercept)  -1.2796 0.03797       -1.3511       -1.2077 1.00 871.00
    #>     BMA          X1   0.4831 0.04240        0.4024        0.5668 1.01 822.00
    #>     BMA          X2  -0.5169 0.05254       -0.6226       -0.4172 1.00 800.00
    #>     BMA        dose   1.0788 0.02024        1.0386        1.1176 1.00 800.00
    #>     BMA       sigma   1.1391 0.01503        1.1108        1.1685 1.02 744.00

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

    coef(fit.ameras.linreg)
    #>                     RC        ERC       MCML        FMA        BMA
    #> (Intercept) -1.3622239 -1.3610557 -1.2802532 -1.2803273 -1.2796163
    #> X1           0.4807463  0.4795308  0.4836506  0.4836869  0.4830670
    #> X2          -0.5187877 -0.5188299 -0.5171452 -0.5173188 -0.5169488
    #> dose         1.1660368  1.1649227  1.0790214  1.0791951  1.0787855
    #> sigma        1.1075696  1.1062486  1.1377732  1.1377321  1.1391166

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"))
    #> Note: BMA may require extensive computation time in the order of multiple hours
    #> Fitting RC
    #> Fitting ERC
    #> Fitting MCML
    #> Fitting FMA
    #> Fitting BMA
    #> Defining model
    #> Building model
    #> Setting data and initial values
    #> Running calculate on model
    #>   [Note] Any error reports that follow may simply reflect missing values in model variables.
    #> Checking model sizes and dimensions
    #>   [Note] This model is not fully initialized. This is not an error.
    #>          To see which variables are not initialized, use model$initializeInfo().
    #>          For more information on model initialization, see help(modelInitialization).
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> running chain 1...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|
    #> running chain 2...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|
    #> Warning in ameras.bma(family = family, dosevars = dosevars, data = data, :
    #> WARNING: Potential problems with MCMC convergence, consider using longer chains

    summary(fit.ameras.logreg)
    #> Call:
    #> ameras(data = data, family = "binomial", Y = "Y.binomial", dosevars = dosevars,
    #>     X = c("X1", "X2"), methods = c("RC", "ERC", "MCML", "FMA",
    #>         "BMA"), deg = 2, doseRRmod = "EXP", CI = c("wald.orig",
    #>         "percentile"), nburnin.BMA = 1000, niter.BMA = 5000)
    #>
    #> Total run time: 87.1 seconds
    #>
    #> Runtime in seconds by method:
    #>
    #>  Method Runtime
    #>      RC     0.2
    #>     ERC    50.6
    #>    MCML     0.8
    #>     FMA     2.0
    #>     BMA    33.5
    #>
    #> Summary of coefficients by method:
    #>
    #>  Method         Term Estimate      SE CI.lowerbound CI.upperbound Rhat   n.eff
    #>      RC  (Intercept) -0.94461 0.08409     -1.109439      -0.77979   NA      NA
    #>      RC           X1  0.44552 0.07667      0.295254       0.59579   NA      NA
    #>      RC           X2 -0.33376 0.09601     -0.521945      -0.14557   NA      NA
    #>      RC         dose  0.37904 0.10388      0.175441       0.58265   NA      NA
    #>      RC dose_squared  0.01943 0.02750     -0.034468       0.07334   NA      NA
    #>     ERC  (Intercept) -0.93189 0.08533     -1.099136      -0.76464   NA      NA
    #>     ERC           X1  0.44998 0.07678      0.299496       0.60046   NA      NA
    #>     ERC           X2 -0.33924 0.09614     -0.527664      -0.15081   NA      NA
    #>     ERC         dose  0.33858 0.10745      0.127985       0.54917   NA      NA
    #>     ERC dose_squared  0.03528 0.02841     -0.020409       0.09097   NA      NA
    #>    MCML  (Intercept) -0.91147 0.08356     -1.075241      -0.74770   NA      NA
    #>    MCML           X1  0.44619 0.07674      0.295784       0.59660   NA      NA
    #>    MCML           X2 -0.34431 0.09625     -0.532951      -0.15567   NA      NA
    #>    MCML         dose  0.31654 0.10412      0.112472       0.52061   NA      NA
    #>    MCML dose_squared  0.03800 0.02774     -0.016375       0.09237   NA      NA
    #>     FMA  (Intercept) -0.91244 0.08386     -1.076771      -0.74789   NA      NA
    #>     FMA           X1  0.44619 0.07686      0.295520       0.59764   NA      NA
    #>     FMA           X2 -0.34400 0.09627     -0.533249      -0.15610   NA      NA
    #>     FMA         dose  0.31871 0.10523      0.114880       0.52709   NA      NA
    #>     FMA dose_squared  0.03735 0.02816     -0.018766       0.09172   NA      NA
    #>     BMA  (Intercept) -0.90779 0.08101     -1.056223      -0.74279 1.03  214.00
    #>     BMA           X1  0.44552 0.07775      0.286128       0.59552 1.01  370.00
    #>     BMA           X2 -0.34361 0.09652     -0.529598      -0.15370 1.00 1242.00
    #>     BMA         dose  0.30658 0.09894      0.102541       0.49441 1.08   92.00
    #>     BMA dose_squared  0.04134 0.02564     -0.007435       0.09729 1.11  105.00

    coef(fit.ameras.logreg)
    #>                       RC         ERC        MCML         FMA         BMA
    #> (Intercept)  -0.94461327 -0.93188863 -0.91146831 -0.91243544 -0.90778810
    #> X1            0.44552273  0.44997904  0.44619209  0.44619391  0.44552058
    #> X2           -0.33375991 -0.33923955 -0.34430935 -0.34400234 -0.34360534
    #> dose          0.37904346  0.33857980  0.31654277  0.31870705  0.30657532
    #> dose_squared  0.01943381  0.03528262  0.03799845  0.03735107  0.04134337

    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"))
    #> Note: BMA may require extensive computation time in the order of multiple hours
    #> Fitting RC
    #> Fitting ERC
    #> Fitting MCML
    #> Fitting FMA
    #> Fitting BMA
    #> Defining model
    #> Building model
    #> Setting data and initial values
    #> Running calculate on model
    #>   [Note] Any error reports that follow may simply reflect missing values in model variables.
    #> Checking model sizes and dimensions
    #>   [Note] This model is not fully initialized. This is not an error.
    #>          To see which variables are not initialized, use model$initializeInfo().
    #>          For more information on model initialization, see help(modelInitialization).
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> running chain 1...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|
    #> running chain 2...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|

    summary(fit.ameras.logreg.lin)
    #> Call:
    #> ameras(data = data, family = "binomial", Y = "Y.binomial", dosevars = dosevars,
    #>     X = c("X1", "X2"), methods = c("RC", "ERC", "MCML", "FMA",
    #>         "BMA"), deg = 1, doseRRmod = "EXP", CI = c("wald.orig",
    #>         "percentile"), nburnin.BMA = 1000, niter.BMA = 5000)
    #>
    #> Total run time: 67.6 seconds
    #>
    #> Runtime in seconds by method:
    #>
    #>  Method Runtime
    #>      RC     0.1
    #>     ERC    33.9
    #>    MCML     0.5
    #>     FMA     1.0
    #>     BMA    32.1
    #>
    #> Summary of coefficients by method:
    #>
    #>  Method        Term Estimate      SE CI.lowerbound CI.upperbound Rhat  n.eff
    #>      RC (Intercept)  -0.9760 0.07190       -1.1169       -0.8351   NA     NA
    #>      RC          X1   0.4460 0.07667        0.2957        0.5963   NA     NA
    #>      RC          X2  -0.3359 0.09596       -0.5240       -0.1478   NA     NA
    #>      RC        dose   0.4471 0.04101        0.3668        0.5275   NA     NA
    #>     ERC (Intercept)  -0.9898 0.07216       -1.1312       -0.8483   NA     NA
    #>     ERC          X1   0.4533 0.07667        0.3030        0.6036   NA     NA
    #>     ERC          X2  -0.3437 0.09609       -0.5320       -0.1554   NA     NA
    #>     ERC        dose   0.4632 0.04086        0.3831        0.5433   NA     NA
    #>    MCML (Intercept)  -0.9725 0.07154       -1.1127       -0.8323   NA     NA
    #>    MCML          X1   0.4469 0.07674        0.2965        0.5973   NA     NA
    #>    MCML          X2  -0.3467 0.09618       -0.5352       -0.1582   NA     NA
    #>    MCML        dose   0.4498 0.04105        0.3693        0.5302   NA     NA
    #>     FMA (Intercept)  -0.9725 0.07138       -1.1131       -0.8333   NA     NA
    #>     FMA          X1   0.4472 0.07664        0.2975        0.5981   NA     NA
    #>     FMA          X2  -0.3465 0.09602       -0.5347       -0.1580   NA     NA
    #>     FMA        dose   0.4497 0.04101        0.3697        0.5302   NA     NA
    #>     BMA (Intercept)  -0.9791 0.06869       -1.1235       -0.8492 1.03 386.00
    #>     BMA          X1   0.4555 0.07669        0.3098        0.6074 1.00 482.00
    #>     BMA          X2  -0.3486 0.09589       -0.5382       -0.1523 1.02 882.00
    #>     BMA        dose   0.4523 0.04042        0.3743        0.5292 1.03 518.00

    coef(fit.ameras.logreg.lin)
    #>                     RC        ERC       MCML        FMA        BMA
    #> (Intercept) -0.9759817 -0.9897555 -0.9724854 -0.9725462 -0.9791394
    #> X1           0.4459778  0.4533083  0.4468997  0.4471557  0.4555060
    #> X2          -0.3358992 -0.3436867 -0.3467298 -0.3465319 -0.3485601
    #> dose         0.4471346  0.4632116  0.4497526  0.4497072  0.4522620

    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"))
    #> Note: BMA may require extensive computation time in the order of multiple hours
    #> Fitting RC
    #> Fitting ERC
    #> Fitting MCML
    #> Fitting FMA
    #> Fitting BMA
    #> Defining model
    #> Building model
    #> Setting data and initial values
    #> Running calculate on model
    #>   [Note] Any error reports that follow may simply reflect missing values in model variables.
    #> Checking model sizes and dimensions
    #>   [Note] This model is not fully initialized. This is not an error.
    #>          To see which variables are not initialized, use model$initializeInfo().
    #>          For more information on model initialization, see help(modelInitialization).
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> running chain 1...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|
    #> running chain 2...
    #> warning: logProb of data node Y[247]: logProb less than -1e12.
    #> warning: logProb of data node Y[676]: logProb less than -1e12.
    #> warning: logProb of data node Y[833]: logProb less than -1e12.
    #> warning: logProb of data node Y[1635]: logProb less than -1e12.
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|

    summary(fit.ameras.poisson)
    #> Call:
    #> ameras(data = data, family = "poisson", Y = "Y.poisson", dosevars = dosevars,
    #>     X = c("X1", "X2"), methods = c("RC", "ERC", "MCML", "FMA",
    #>         "BMA"), deg = 2, doseRRmod = "EXP", CI = c("wald.orig",
    #>         "percentile"), nburnin.BMA = 1000, niter.BMA = 5000)
    #>
    #> Total run time: 95.3 seconds
    #>
    #> Runtime in seconds by method:
    #>
    #>  Method Runtime
    #>      RC     0.2
    #>     ERC    55.8
    #>    MCML     1.1
    #>     FMA     2.4
    #>     BMA    35.8
    #>
    #> Summary of coefficients by method:
    #>
    #>  Method         Term Estimate       SE CI.lowerbound CI.upperbound Rhat  n.eff
    #>      RC  (Intercept) -1.09456 0.048655      -1.18992      -0.99920   NA     NA
    #>      RC           X1  0.49070 0.041922       0.40853       0.57287   NA     NA
    #>      RC           X2 -0.37625 0.055638      -0.48530      -0.26719   NA     NA
    #>      RC         dose  0.61976 0.040375       0.54062       0.69889   NA     NA
    #>      RC dose_squared -0.03849 0.007566      -0.05332      -0.02366   NA     NA
    #>     ERC  (Intercept) -1.09068 0.048474      -1.18569      -0.99567   NA     NA
    #>     ERC           X1  0.49180 0.042008       0.40946       0.57414   NA     NA
    #>     ERC           X2 -0.37855 0.055639      -0.48760      -0.26950   NA     NA
    #>     ERC         dose  0.61138 0.039328       0.53430       0.68847   NA     NA
    #>     ERC dose_squared -0.03626 0.007177      -0.05033      -0.02220   NA     NA
    #>    MCML  (Intercept) -1.07519 0.046954      -1.16722      -0.98316   NA     NA
    #>    MCML           X1  0.49897 0.041918       0.41681       0.58113   NA     NA
    #>    MCML           X2 -0.37711 0.055643      -0.48617      -0.26805   NA     NA
    #>    MCML         dose  0.60089 0.034218       0.53383       0.66796   NA     NA
    #>    MCML dose_squared -0.03644 0.005719      -0.04764      -0.02523   NA     NA
    #>     FMA  (Intercept) -1.07542 0.047056      -1.16736      -0.98339   NA     NA
    #>     FMA           X1  0.49896 0.041800       0.41696       0.58065   NA     NA
    #>     FMA           X2 -0.37716 0.055479      -0.48601      -0.26864   NA     NA
    #>     FMA         dose  0.60138 0.035163       0.53336       0.67123   NA     NA
    #>     FMA dose_squared -0.03656 0.005998      -0.04836      -0.02519   NA     NA
    #>     BMA  (Intercept) -1.07944 0.044741      -1.16649      -0.99362 1.02 215.00
    #>     BMA           X1  0.49801 0.039956       0.42219       0.57163 1.03 435.00
    #>     BMA           X2 -0.38167 0.055527      -0.49101      -0.27220 1.00 800.00
    #>     BMA         dose  0.60805 0.035663       0.54553       0.68040 1.03  87.00
    #>     BMA dose_squared -0.03798 0.006215      -0.05048      -0.02715 1.03  79.00

    coef(fit.ameras.poisson)
    #>                       RC         ERC        MCML         FMA        BMA
    #> (Intercept)  -1.09455980 -1.09068092 -1.07519346 -1.07541780 -1.0794372
    #> X1            0.49070108  0.49179875  0.49897425  0.49895993  0.4980120
    #> X2           -0.37624508 -0.37854805 -0.37711401 -0.37715581 -0.3816682
    #> dose          0.61975742  0.61138406  0.60089474  0.60137955  0.6080481
    #> dose_squared -0.03849039 -0.03626348 -0.03643505 -0.03656317 -0.0379750

    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"))
    #> Note: BMA may require extensive computation time in the order of multiple hours
    #> Fitting RC
    #> Fitting ERC
    #> Fitting MCML
    #> Fitting FMA
    #> Fitting BMA
    #> Defining model
    #> Building model
    #> Setting data and initial values
    #> Running calculate on model
    #>   [Note] Any error reports that follow may simply reflect missing values in model variables.
    #> Checking model sizes and dimensions
    #>   [Note] This model is not fully initialized. This is not an error.
    #>          To see which variables are not initialized, use model$initializeInfo().
    #>          For more information on model initialization, see help(modelInitialization).
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> running chain 1...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|
    #> running chain 2...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|

    summary(fit.ameras.poisson.lin)
    #> Call:
    #> ameras(data = data, family = "poisson", Y = "Y.poisson", dosevars = dosevars,
    #>     X = c("X1", "X2"), methods = c("RC", "ERC", "MCML", "FMA",
    #>         "BMA"), deg = 1, doseRRmod = "EXP", CI = c("wald.orig",
    #>         "percentile"), nburnin.BMA = 1000, niter.BMA = 5000)
    #>
    #> Total run time: 71 seconds
    #>
    #> Runtime in seconds by method:
    #>
    #>  Method Runtime
    #>      RC     0.1
    #>     ERC    36.2
    #>    MCML     0.5
    #>     FMA     1.3
    #>     BMA    32.9
    #>
    #> Summary of coefficients by method:
    #>
    #>  Method        Term Estimate      SE CI.lowerbound CI.upperbound Rhat  n.eff
    #>      RC (Intercept)  -0.9650 0.04127       -1.0458       -0.8841   NA     NA
    #>      RC          X1   0.5054 0.04195        0.4232        0.5876   NA     NA
    #>      RC          X2  -0.3640 0.05560       -0.4730       -0.2550   NA     NA
    #>      RC        dose   0.4204 0.01346        0.3940        0.4468   NA     NA
    #>     ERC (Intercept)  -0.9667 0.04134       -1.0477       -0.8856   NA     NA
    #>     ERC          X1   0.5049 0.04202        0.4226        0.5873   NA     NA
    #>     ERC          X2  -0.3667 0.05559       -0.4756       -0.2577   NA     NA
    #>     ERC        dose   0.4222 0.01361        0.3955        0.4489   NA     NA
    #>    MCML (Intercept)  -0.9173 0.04048       -0.9967       -0.8380   NA     NA
    #>    MCML          X1   0.5129 0.04231        0.4300        0.5958   NA     NA
    #>    MCML          X2  -0.3579 0.05582       -0.4674       -0.2485   NA     NA
    #>    MCML        dose   0.3823 0.01231        0.3581        0.4064   NA     NA
    #>     FMA (Intercept)  -0.9174 0.04044       -0.9969       -0.8384   NA     NA
    #>     FMA          X1   0.5121 0.04257        0.4284        0.5950   NA     NA
    #>     FMA          X2  -0.3584 0.05587       -0.4680       -0.2496   NA     NA
    #>     FMA        dose   0.3826 0.01255        0.3585        0.4077   NA     NA
    #>     BMA (Intercept)  -0.9180 0.04041       -0.9935       -0.8378 1.00 301.00
    #>     BMA          X1   0.5100 0.04121        0.4248        0.5896 1.01 364.00
    #>     BMA          X2  -0.3581 0.05584       -0.4661       -0.2460 1.00 800.00
    #>     BMA        dose   0.3833 0.01343        0.3561        0.4087 1.01 371.00

    coef(fit.ameras.poisson.lin)
    #>                     RC        ERC       MCML        FMA        BMA
    #> (Intercept) -0.9649529 -0.9666735 -0.9173445 -0.9173785 -0.9179930
    #> X1           0.5054270  0.5049256  0.5129102  0.5120519  0.5099994
    #> X2          -0.3639707 -0.3666592 -0.3579466 -0.3583942 -0.3581445
    #> dose         0.4204226  0.4222084  0.3822723  0.3826304  0.3833150

    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"))
    #> Note: BMA may require extensive computation time in the order of multiple hours
    #> Fitting RC
    #> Fitting ERC
    #> Warning in ameras.rc(family = family, dosevars = dosevars, data = data, :
    #> WARNING: Hessian was not invertible or inverse was not positive definite,
    #> variance matrix could not be obtained
    #> Fitting MCML
    #> Fitting FMA
    #> Fitting BMA
    #> Defining model
    #> Building model
    #> Setting data and initial values
    #> Running calculate on model
    #>   [Note] Any error reports that follow may simply reflect missing values in model variables.
    #> Checking model sizes and dimensions
    #>   [Note] This model is not fully initialized. This is not an error.
    #>          To see which variables are not initialized, use model$initializeInfo().
    #>          For more information on model initialization, see help(modelInitialization).
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> running chain 1...
    #> warning: logProb of data node zeros[7]: logProb less than -1e12.
    #> warning: logProb of data node zeros[17]: logProb less than -1e12.
    #> warning: logProb of data node zeros[41]: logProb less than -1e12.
    #> warning: logProb of data node zeros[46]: logProb less than -1e12.
    #> warning: logProb of data node zeros[52]: logProb less than -1e12.
    #> warning: logProb of data node zeros[58]: logProb less than -1e12.
    #> warning: logProb of data node zeros[71]: logProb less than -1e12.
    #> warning: logProb of data node zeros[157]: logProb less than -1e12.
    #> warning: logProb of data node zeros[166]: logProb less than -1e12.
    #> warning: logProb of data node zeros[175]: logProb less than -1e12.
    #> warning: logProb of data node zeros[178]: logProb less than -1e12.
    #> warning: logProb of data node zeros[212]: logProb less than -1e12.
    #> warning: logProb of data node zeros[238]: logProb less than -1e12.
    #> warning: logProb of data node zeros[243]: logProb less than -1e12.
    #> warning: logProb of data node zeros[246]: logProb less than -1e12.
    #> warning: logProb of data node zeros[247]: logProb less than -1e12.
    #> warning: logProb of data node zeros[260]: logProb less than -1e12.
    #> warning: logProb of data node zeros[261]: logProb less than -1e12.
    #> warning: logProb of data node zeros[263]: logProb less than -1e12.
    #> warning: logProb of data node zeros[267]: logProb less than -1e12.
    #> warning: logProb of data node zeros[270]: logProb less than -1e12.
    #> warning: logProb of data node zeros[278]: logProb less than -1e12.
    #> warning: logProb of data node zeros[311]: logProb less than -1e12.
    #> warning: logProb of data node zeros[335]: logProb less than -1e12.
    #> warning: logProb of data node zeros[341]: logProb less than -1e12.
    #> warning: logProb of data node zeros[348]: logProb less than -1e12.
    #> warning: logProb of data node zeros[352]: logProb less than -1e12.
    #> warning: logProb of data node zeros[370]: logProb less than -1e12.
    #> warning: logProb of data node zeros[372]: logProb less than -1e12.
    #> warning: logProb of data node zeros[378]: logProb less than -1e12.
    #> warning: logProb of data node zeros[386]: logProb less than -1e12.
    #> warning: logProb of data node zeros[420]: logProb less than -1e12.
    #> warning: logProb of data node zeros[432]: logProb less than -1e12.
    #> warning: logProb of data node zeros[433]: logProb less than -1e12.
    #> warning: logProb of data node zeros[438]: logProb less than -1e12.
    #> warning: logProb of data node zeros[440]: logProb less than -1e12.
    #> warning: logProb of data node zeros[449]: logProb less than -1e12.
    #> warning: logProb of data node zeros[450]: logProb less than -1e12.
    #> warning: logProb of data node zeros[458]: logProb less than -1e12.
    #> warning: logProb of data node zeros[508]: logProb less than -1e12.
    #> warning: logProb of data node zeros[509]: logProb less than -1e12.
    #> warning: logProb of data node zeros[528]: logProb less than -1e12.
    #> warning: logProb of data node zeros[535]: logProb less than -1e12.
    #> warning: logProb of data node zeros[546]: logProb less than -1e12.
    #> warning: logProb of data node zeros[547]: logProb less than -1e12.
    #> warning: logProb of data node zeros[552]: logProb less than -1e12.
    #> warning: logProb of data node zeros[554]: logProb less than -1e12.
    #> warning: logProb of data node zeros[577]: logProb less than -1e12.
    #> warning: logProb of data node zeros[608]: logProb less than -1e12.
    #> warning: logProb of data node zeros[611]: logProb less than -1e12.
    #> warning: logProb of data node zeros[620]: logProb less than -1e12.
    #> warning: logProb of data node zeros[627]: logProb less than -1e12.
    #> warning: logProb of data node zeros[631]: logProb less than -1e12.
    #> warning: logProb of data node zeros[656]: logProb less than -1e12.
    #> warning: logProb of data node zeros[661]: logProb less than -1e12.
    #> warning: logProb of data node zeros[676]: logProb less than -1e12.
    #> warning: logProb of data node zeros[684]: logProb less than -1e12.
    #> warning: logProb of data node zeros[698]: logProb less than -1e12.
    #> warning: logProb of data node zeros[714]: logProb less than -1e12.
    #> warning: logProb of data node zeros[716]: logProb less than -1e12.
    #> warning: logProb of data node zeros[727]: logProb less than -1e12.
    #> warning: logProb of data node zeros[749]: logProb less than -1e12.
    #> warning: logProb of data node zeros[806]: logProb less than -1e12.
    #> warning: logProb of data node zeros[809]: logProb less than -1e12.
    #> warning: logProb of data node zeros[821]: logProb less than -1e12.
    #> warning: logProb of data node zeros[833]: logProb less than -1e12.
    #> warning: logProb of data node zeros[839]: logProb less than -1e12.
    #> warning: logProb of data node zeros[864]: logProb less than -1e12.
    #> warning: logProb of data node zeros[873]: logProb less than -1e12.
    #> warning: logProb of data node zeros[910]: logProb less than -1e12.
    #> warning: logProb of data node zeros[962]: logProb less than -1e12.
    #> warning: logProb of data node zeros[971]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1002]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1004]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1009]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1014]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1036]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1040]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1052]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1070]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1074]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1076]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1081]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1116]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1122]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1137]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1151]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1179]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1187]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1208]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1214]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1245]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1247]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1282]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1296]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1300]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1319]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1322]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1327]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1340]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1354]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1357]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1368]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1373]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1375]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1406]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1425]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1449]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1464]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1470]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1475]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1480]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1492]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1496]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1517]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1524]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1527]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1541]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1553]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1566]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1569]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1572]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1588]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1603]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1617]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1635]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1661]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1663]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1664]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1693]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1716]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1755]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1768]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1783]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1796]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1807]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1818]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1836]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1838]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1851]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1888]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1900]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1921]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1930]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1948]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1968]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1971]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1972]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1979]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1983]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1991]: logProb less than -1e12.
    #> warning: logProb of data node zeros[1997]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2021]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2035]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2054]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2087]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2089]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2101]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2112]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2129]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2131]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2134]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2156]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2198]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2208]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2216]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2227]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2237]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2239]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2255]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2294]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2312]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2338]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2339]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2360]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2368]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2373]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2395]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2398]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2407]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2409]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2413]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2422]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2444]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2452]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2463]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2475]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2478]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2490]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2495]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2524]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2530]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2559]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2562]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2573]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2605]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2617]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2655]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2665]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2671]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2684]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2686]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2688]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2694]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2715]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2723]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2728]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2733]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2743]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2747]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2757]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2771]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2774]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2787]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2807]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2813]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2824]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2838]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2839]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2907]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2909]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2915]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2937]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2945]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2968]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2971]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2993]: logProb less than -1e12.
    #> warning: logProb of data node zeros[2998]: logProb less than -1e12.
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|
    #> running chain 2...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|

    summary(fit.ameras.prophaz)
    #> Call:
    #> ameras(data = data, family = "prophaz", Y = "status", dosevars = dosevars,
    #>     X = c("X1", "X2"), exit = "time", methods = c("RC", "ERC",
    #>         "MCML", "FMA", "BMA"), deg = 2, doseRRmod = "EXP", CI = c("wald.orig",
    #>         "percentile"), nburnin.BMA = 1000, niter.BMA = 5000)
    #>
    #> Total run time: 363.8 seconds
    #>
    #> Runtime in seconds by method:
    #>
    #>  Method Runtime
    #>      RC     0.2
    #>     ERC   296.1
    #>    MCML     0.6
    #>     FMA     1.8
    #>     BMA    65.1
    #>
    #> Summary of coefficients by method:
    #>
    #>  Method         Term  Estimate      SE CI.lowerbound CI.upperbound Rhat  n.eff
    #>      RC           X1  0.629674 0.08485       0.46338      0.795971   NA     NA
    #>      RC           X2 -0.423116 0.11150      -0.64165     -0.204584   NA     NA
    #>      RC         dose  0.587614 0.08915       0.41289      0.762339   NA     NA
    #>      RC dose_squared -0.033682 0.01807      -0.06910      0.001741   NA     NA
    #>     ERC           X1  0.636755      NA            NA            NA   NA     NA
    #>     ERC           X2 -0.422042      NA            NA            NA   NA     NA
    #>     ERC         dose  0.295488      NA            NA            NA   NA     NA
    #>     ERC dose_squared -0.003357      NA            NA            NA   NA     NA
    #>    MCML           X1  0.629162 0.08518       0.46221      0.796110   NA     NA
    #>    MCML           X2 -0.425156 0.11196      -0.64460     -0.205715   NA     NA
    #>    MCML         dose  0.590217 0.08046       0.43252      0.747910   NA     NA
    #>    MCML dose_squared -0.038038 0.01523      -0.06789     -0.008184   NA     NA
    #>     FMA           X1  0.629250 0.08557       0.46066      0.796701   NA     NA
    #>     FMA           X2 -0.425507 0.11190      -0.64526     -0.206481   NA     NA
    #>     FMA         dose  0.589537 0.08127       0.42964      0.749311   NA     NA
    #>     FMA dose_squared -0.037872 0.01546      -0.06813     -0.007464   NA     NA
    #>     BMA           X1  0.635157 0.08523       0.46388      0.802937 1.00 397.00
    #>     BMA           X2 -0.423897 0.11063      -0.63473     -0.209823 1.01 747.00
    #>     BMA         dose  0.617073 0.08181       0.46308      0.778760 1.01 101.00
    #>     BMA dose_squared -0.042851 0.01554      -0.07371     -0.013833 1.01 116.00
    #>     BMA        h0[1]  0.316373 0.05147       0.22764      0.437066 1.00 384.00
    #>     BMA        h0[2]  0.343831 0.05275       0.25174      0.455786 1.00 394.00
    #>     BMA        h0[3]  0.262707 0.04153       0.19213      0.348833 1.00 281.00
    #>     BMA        h0[4]  0.294011 0.04621       0.21334      0.388598 1.00 410.00
    #>     BMA        h0[5]  0.314660 0.04963       0.23387      0.423030 1.01 362.00
    #>     BMA        h0[6]  0.396503 0.06335       0.28802      0.529583 1.00 411.00
    #>     BMA        h0[7]  0.269226 0.03952       0.19758      0.351930 1.00 373.00
    #>     BMA        h0[8]  0.309864 0.04606       0.22259      0.405984 1.00 326.00
    #>     BMA        h0[9]  0.268182 0.04019       0.19841      0.354162 1.00 398.00
    #>     BMA       h0[10]  0.304145 0.04624       0.22255      0.396432 1.01 367.00

    coef(fit.ameras.prophaz)
    #>                       RC          ERC        MCML         FMA         BMA
    #> X1            0.62967434  0.636754979  0.62916199  0.62925032  0.63515662
    #> X2           -0.42311583 -0.422041734 -0.42515576 -0.42550683 -0.42389658
    #> dose          0.58761408  0.295487972  0.59021691  0.58953694  0.61707317
    #> dose_squared -0.03368205 -0.003356734 -0.03803829 -0.03787174 -0.04285062
    #> h0[1]                 NA           NA          NA          NA  0.31637277
    #> h0[2]                 NA           NA          NA          NA  0.34383087
    #> h0[3]                 NA           NA          NA          NA  0.26270690
    #> h0[4]                 NA           NA          NA          NA  0.29401059
    #> h0[5]                 NA           NA          NA          NA  0.31466013
    #> h0[6]                 NA           NA          NA          NA  0.39650271
    #> h0[7]                 NA           NA          NA          NA  0.26922630
    #> h0[8]                 NA           NA          NA          NA  0.30986447
    #> h0[9]                 NA           NA          NA          NA  0.26818217
    #> h0[10]                NA           NA          NA          NA  0.30414532

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

    fit.ameras.prophaz$BMA$prophaz.timepoints
    #> NULL

    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"))
    #> Note: BMA may require extensive computation time in the order of multiple hours
    #> Fitting RC
    #> Fitting ERC
    #> Warning in ameras.rc(family = family, dosevars = dosevars, data = data, :
    #> WARNING: Hessian was not invertible or inverse was not positive definite,
    #> variance matrix could not be obtained
    #> Fitting MCML
    #> Fitting FMA
    #> Fitting BMA
    #> Defining model
    #> Building model
    #> Setting data and initial values
    #> Running calculate on model
    #>   [Note] Any error reports that follow may simply reflect missing values in model variables.
    #> Checking model sizes and dimensions
    #>   [Note] This model is not fully initialized. This is not an error.
    #>          To see which variables are not initialized, use model$initializeInfo().
    #>          For more information on model initialization, see help(modelInitialization).
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> running chain 1...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|
    #> running chain 2...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|

    summary(fit.ameras.prophaz.lin)
    #> Call:
    #> ameras(data = data, family = "prophaz", Y = "status", dosevars = dosevars,
    #>     X = c("X1", "X2"), exit = "time", methods = c("RC", "ERC",
    #>         "MCML", "FMA", "BMA"), deg = 1, doseRRmod = "EXP", CI = c("wald.orig",
    #>         "percentile"), nburnin.BMA = 1000, niter.BMA = 5000)
    #>
    #> Total run time: 219.1 seconds
    #>
    #> Runtime in seconds by method:
    #>
    #>  Method Runtime
    #>      RC     0.1
    #>     ERC   154.0
    #>    MCML     0.3
    #>     FMA     1.0
    #>     BMA    63.7
    #>
    #> Summary of coefficients by method:
    #>
    #>  Method   Term Estimate      SE CI.lowerbound CI.upperbound Rhat  n.eff
    #>      RC     X1   0.6358 0.08488        0.4695        0.8022   NA     NA
    #>      RC     X2  -0.4161 0.11144       -0.6345       -0.1977   NA     NA
    #>      RC   dose   0.4284 0.03004        0.3695        0.4873   NA     NA
    #>     ERC     X1   0.6416      NA            NA            NA   NA     NA
    #>     ERC     X2  -0.4220      NA            NA            NA   NA     NA
    #>     ERC   dose   0.2812      NA            NA            NA   NA     NA
    #>    MCML     X1   0.6462 0.08565        0.4783        0.8141   NA     NA
    #>    MCML     X2  -0.4125 0.11253       -0.6330       -0.1919   NA     NA
    #>    MCML   dose   0.4002 0.02908        0.3432        0.4572   NA     NA
    #>     FMA     X1   0.6457 0.08567        0.4774        0.8142   NA     NA
    #>     FMA     X2  -0.4125 0.11263       -0.6344       -0.1925   NA     NA
    #>     FMA   dose   0.4001 0.02901        0.3431        0.4567   NA     NA
    #>     BMA     X1   0.6474 0.08153        0.4888        0.8140 1.00 388.00
    #>     BMA     X2  -0.4248 0.11473       -0.6443       -0.2092 1.00 800.00
    #>     BMA   dose   0.3987 0.02884        0.3427        0.4572 1.01 478.00
    #>     BMA  h0[1]   0.3646 0.05506        0.2675        0.4718 1.00 644.00
    #>     BMA  h0[2]   0.3926 0.05694        0.2857        0.5086 1.01 540.00
    #>     BMA  h0[3]   0.3036 0.04541        0.2239        0.3993 1.00 549.00
    #>     BMA  h0[4]   0.3378 0.04899        0.2479        0.4480 1.00 692.00
    #>     BMA  h0[5]   0.3615 0.05340        0.2688        0.4773 1.00 720.00
    #>     BMA  h0[6]   0.4547 0.06635        0.3319        0.5874 1.00 546.00
    #>     BMA  h0[7]   0.3080 0.04573        0.2288        0.4042 1.01 690.00
    #>     BMA  h0[8]   0.3551 0.05105        0.2575        0.4580 1.00 589.00
    #>     BMA  h0[9]   0.3091 0.04455        0.2277        0.4003 1.00 907.00
    #>     BMA h0[10]   0.3496 0.05211        0.2535        0.4542 1.01 542.00

    coef(fit.ameras.prophaz.lin)
    #>                RC        ERC       MCML        FMA        BMA
    #> X1      0.6358422  0.6415799  0.6462073  0.6457466  0.6474031
    #> X2     -0.4160762 -0.4219654 -0.4124721 -0.4125357 -0.4247814
    #> dose    0.4283893  0.2811645  0.4001775  0.4000975  0.3987353
    #> h0[1]          NA         NA         NA         NA  0.3645755
    #> h0[2]          NA         NA         NA         NA  0.3925686
    #> h0[3]          NA         NA         NA         NA  0.3035865
    #> h0[4]          NA         NA         NA         NA  0.3378431
    #> h0[5]          NA         NA         NA         NA  0.3615119
    #> h0[6]          NA         NA         NA         NA  0.4547288
    #> h0[7]          NA         NA         NA         NA  0.3079510
    #> h0[8]          NA         NA         NA         NA  0.3550875
    #> h0[9]          NA         NA         NA         NA  0.3091293
    #> h0[10]         NA         NA         NA         NA  0.3495984

    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"))
    #> Note: BMA may require extensive computation time in the order of multiple hours
    #> Fitting RC
    #> Fitting ERC
    #> Fitting MCML
    #> Fitting FMA
    #> Fitting BMA
    #> Defining model
    #> Building model
    #> Setting data and initial values
    #> Running calculate on model
    #>   [Note] Any error reports that follow may simply reflect missing values in model variables.
    #> Checking model sizes and dimensions
    #>   [Note] This model is not fully initialized. This is not an error.
    #>          To see which variables are not initialized, use model$initializeInfo().
    #>          For more information on model initialization, see help(modelInitialization).
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> running chain 1...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|
    #> running chain 2...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|

    summary(fit.ameras.multinomial)
    #> Call:
    #> ameras(data = data, family = "multinomial", Y = "Y.multinomial",
    #>     dosevars = dosevars, X = c("X1", "X2"), methods = c("RC",
    #>         "ERC", "MCML", "FMA", "BMA"), deg = 2, doseRRmod = "EXP",
    #>     CI = c("wald.orig", "percentile"), nburnin.BMA = 1000, niter.BMA = 5000)
    #>
    #> Total run time: 306.1 seconds
    #>
    #> Runtime in seconds by method:
    #>
    #>  Method Runtime
    #>      RC     0.9
    #>     ERC   116.5
    #>    MCML     6.4
    #>     FMA     7.8
    #>     BMA   174.5
    #>
    #> Summary of coefficients by method:
    #>
    #>  Method             Term  Estimate      SE CI.lowerbound CI.upperbound Rhat
    #>      RC  (1)_(Intercept) -1.134879 0.11395      -1.35823      -0.91153   NA
    #>      RC           (1)_X1  0.521885 0.10631       0.31352       0.73025   NA
    #>      RC           (1)_X2 -0.347821 0.14935      -0.64054      -0.05510   NA
    #>      RC         (1)_dose  0.541114 0.13250       0.28142       0.80081   NA
    #>      RC (1)_dose_squared -0.028978 0.03377      -0.09517       0.03721   NA
    #>      RC  (2)_(Intercept) -0.007547 0.08795      -0.17992       0.16483   NA
    #>      RC           (2)_X1 -0.513675 0.08620      -0.68264      -0.34471   NA
    #>      RC           (2)_X2  0.705782 0.10738       0.49532       0.91625   NA
    #>      RC         (2)_dose  0.555006 0.11504       0.32952       0.78049   NA
    #>      RC (2)_dose_squared -0.037842 0.03053      -0.09769       0.02201   NA
    #>     ERC  (1)_(Intercept) -1.139984 0.11278      -1.36104      -0.91893   NA
    #>     ERC           (1)_X1  0.524098 0.10651       0.31534       0.73286   NA
    #>     ERC           (1)_X2 -0.351142 0.14952      -0.64420      -0.05808   NA
    #>     ERC         (1)_dose  0.561788 0.12424       0.31827       0.80531   NA
    #>     ERC (1)_dose_squared -0.031891 0.02990      -0.09050       0.02672   NA
    #>     ERC  (2)_(Intercept) -0.008660 0.08626      -0.17772       0.16040   NA
    #>     ERC           (2)_X1 -0.511536 0.08638      -0.68084      -0.34223   NA
    #>     ERC           (2)_X2  0.703688 0.10757       0.49285       0.91452   NA
    #>     ERC         (2)_dose  0.563142 0.10587       0.35564       0.77065   NA
    #>     ERC (2)_dose_squared -0.036409 0.02606      -0.08748       0.01466   NA
    #>    MCML  (1)_(Intercept) -1.121541 0.11214      -1.34134      -0.90174   NA
    #>    MCML           (1)_X1  0.522503 0.10638       0.31400       0.73101   NA
    #>    MCML           (1)_X2 -0.349944 0.14943      -0.64282      -0.05707   NA
    #>    MCML         (1)_dose  0.526846 0.12548       0.28090       0.77280   NA
    #>    MCML (1)_dose_squared -0.027194 0.03020      -0.08639       0.03200   NA
    #>    MCML  (2)_(Intercept) -0.001584 0.08613      -0.17041       0.16724   NA
    #>    MCML           (2)_X1 -0.513787 0.08620      -0.68275      -0.34483   NA
    #>    MCML           (2)_X2  0.703779 0.10746       0.49316       0.91440   NA
    #>    MCML         (2)_dose  0.559807 0.10680       0.35049       0.76913   NA
    #>    MCML (2)_dose_squared -0.041118 0.02628      -0.09264       0.01040   NA
    #>     FMA  (1)_(Intercept) -1.125692 0.11268      -1.34744      -0.90539   NA
    #>     FMA           (1)_X1  0.522203 0.10653       0.31231       0.72993   NA
    #>     FMA           (1)_X2 -0.349524 0.14951      -0.64083      -0.05569   NA
    #>     FMA         (1)_dose  0.536663 0.12823       0.28843       0.79027   NA
    #>     FMA (1)_dose_squared -0.030264 0.03136      -0.09280       0.03040   NA
    #>     FMA  (2)_(Intercept) -0.001039 0.08607      -0.16976       0.16766   NA
    #>     FMA           (2)_X1 -0.513380 0.08602      -0.68183      -0.34437   NA
    #>     FMA           (2)_X2  0.704601 0.10763       0.49377       0.91608   NA
    #>     FMA         (2)_dose  0.557269 0.10747       0.34688       0.76959   NA
    #>     FMA (2)_dose_squared -0.040745 0.02652      -0.09316       0.01108   NA
    #>     BMA  (1)_(Intercept) -1.120533 0.11510      -1.35169      -0.90697 1.02
    #>     BMA           (1)_X1  0.518102 0.11473       0.30333       0.74778 1.00
    #>     BMA           (1)_X2 -0.342415 0.14116      -0.61963      -0.07669 1.01
    #>     BMA         (1)_dose  0.522811 0.12793       0.27114       0.79384 1.04
    #>     BMA (1)_dose_squared -0.024286 0.03270      -0.08275       0.04419 1.01
    #>     BMA  (2)_(Intercept)  0.003372 0.09478      -0.17065       0.19664 1.00
    #>     BMA           (2)_X1 -0.522597 0.08709      -0.68900      -0.35222 1.00
    #>     BMA           (2)_X2  0.708877 0.10811       0.49078       0.91896 1.00
    #>     BMA         (2)_dose  0.554817 0.12798       0.23771       0.76811 1.01
    #>     BMA (2)_dose_squared -0.037946 0.03326      -0.09258       0.04645 1.01
    #>   n.eff
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>  173.00
    #>  330.00
    #>  663.00
    #>   73.00
    #>   73.00
    #>  101.00
    #>  364.00
    #>  462.00
    #>   47.00
    #>   45.00

    coef(fit.ameras.multinomial)
    #>                            RC          ERC        MCML          FMA
    #> (1)_(Intercept)  -1.134879146 -1.139983536 -1.12154130 -1.125692409
    #> (1)_X1            0.521885281  0.524097727  0.52250322  0.522202729
    #> (1)_X2           -0.347820786 -0.351142281 -0.34994376 -0.349524398
    #> (1)_dose          0.541114276  0.561787528  0.52684585  0.536663279
    #> (1)_dose_squared -0.028978278 -0.031891085 -0.02719420 -0.030263819
    #> (2)_(Intercept)  -0.007546883 -0.008659617 -0.00158420 -0.001039459
    #> (2)_X1           -0.513675252 -0.511536376 -0.51378683 -0.513379778
    #> (2)_X2            0.705781705  0.703687593  0.70377938  0.704600521
    #> (2)_dose          0.555005580  0.563142229  0.55980729  0.557268589
    #> (2)_dose_squared -0.037841810 -0.036408828 -0.04111845 -0.040744704
    #>                           BMA
    #> (1)_(Intercept)  -1.120532838
    #> (1)_X1            0.518102027
    #> (1)_X2           -0.342415237
    #> (1)_dose          0.522811057
    #> (1)_dose_squared -0.024285677
    #> (2)_(Intercept)   0.003371856
    #> (2)_X1           -0.522596807
    #> (2)_X2            0.708876654
    #> (2)_dose          0.554817132
    #> (2)_dose_squared -0.037945537

    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"))
    #> Note: BMA may require extensive computation time in the order of multiple hours
    #> Fitting RC
    #> Fitting ERC
    #> Fitting MCML
    #> Fitting FMA
    #> Fitting BMA
    #> Defining model
    #> Building model
    #> Setting data and initial values
    #> Running calculate on model
    #>   [Note] Any error reports that follow may simply reflect missing values in model variables.
    #> Checking model sizes and dimensions
    #>   [Note] This model is not fully initialized. This is not an error.
    #>          To see which variables are not initialized, use model$initializeInfo().
    #>          For more information on model initialization, see help(modelInitialization).
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> running chain 1...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|
    #> running chain 2...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|

    summary(fit.ameras.multinomial.lin)
    #> Call:
    #> ameras(data = data, family = "multinomial", Y = "Y.multinomial",
    #>     dosevars = dosevars, X = c("X1", "X2"), methods = c("RC",
    #>         "ERC", "MCML", "FMA", "BMA"), deg = 1, doseRRmod = "EXP",
    #>     CI = c("wald.orig", "percentile"), nburnin.BMA = 1000, niter.BMA = 5000)
    #>
    #> Total run time: 246.6 seconds
    #>
    #> Runtime in seconds by method:
    #>
    #>  Method Runtime
    #>      RC     0.7
    #>     ERC    83.3
    #>    MCML     5.0
    #>     FMA     5.9
    #>     BMA   151.7
    #>
    #> Summary of coefficients by method:
    #>
    #>  Method            Term Estimate      SE CI.lowerbound CI.upperbound Rhat
    #>      RC (1)_(Intercept) -1.10156 0.10093      -1.29939      -0.90373   NA
    #>      RC          (1)_X1  0.52203 0.10628       0.31372       0.73034   NA
    #>      RC          (1)_X2 -0.34710 0.14929      -0.63971      -0.05448   NA
    #>      RC        (1)_dose  0.45264 0.05805       0.33885       0.56643   NA
    #>      RC (2)_(Intercept)  0.04417 0.07667      -0.10611       0.19445   NA
    #>      RC          (2)_X1 -0.51312 0.08614      -0.68196      -0.34428   NA
    #>      RC          (2)_X2  0.70793 0.10728       0.49767       0.91819   NA
    #>      RC        (2)_dose  0.43049 0.05170       0.32915       0.53183   NA
    #>     ERC (1)_(Intercept) -1.10285 0.10195      -1.30266      -0.90304   NA
    #>     ERC          (1)_X1  0.52577 0.10650       0.31704       0.73451   NA
    #>     ERC          (1)_X2 -0.35176 0.14953      -0.64483      -0.05868   NA
    #>     ERC        (1)_dose  0.46441 0.06099       0.34486       0.58396   NA
    #>     ERC (2)_(Intercept)  0.03947 0.07756      -0.11255       0.19148   NA
    #>     ERC          (2)_X1 -0.51003 0.08634      -0.67926      -0.34081   NA
    #>     ERC          (2)_X2  0.70435 0.10749       0.49368       0.91502   NA
    #>     ERC        (2)_dose  0.44543 0.05364       0.34031       0.55056   NA
    #>    MCML (1)_(Intercept) -1.09147 0.10115      -1.28972      -0.89322   NA
    #>    MCML          (1)_X1  0.52239 0.10633       0.31398       0.73081   NA
    #>    MCML          (1)_X2 -0.34934 0.14941      -0.64218      -0.05650   NA
    #>    MCML        (1)_dose  0.44410 0.05912       0.32822       0.55998   NA
    #>    MCML (2)_(Intercept)  0.05917 0.07624      -0.09027       0.20860   NA
    #>    MCML          (2)_X1 -0.51290 0.08614      -0.68175      -0.34406   NA
    #>    MCML          (2)_X2  0.70604 0.10735       0.49564       0.91645   NA
    #>    MCML        (2)_dose  0.41753 0.05156       0.31647       0.51860   NA
    #>     FMA (1)_(Intercept) -1.09114 0.10101      -1.28880      -0.89308   NA
    #>     FMA          (1)_X1  0.52245 0.10612       0.31541       0.73147   NA
    #>     FMA          (1)_X2 -0.34937 0.14960      -0.64255      -0.05564   NA
    #>     FMA        (1)_dose  0.44381 0.05893       0.32835       0.56009   NA
    #>     FMA (2)_(Intercept)  0.05869 0.07627      -0.09165       0.20799   NA
    #>     FMA          (2)_X1 -0.51272 0.08608      -0.68150      -0.34349   NA
    #>     FMA          (2)_X2  0.70611 0.10708       0.49662       0.91574   NA
    #>     FMA        (2)_dose  0.41769 0.05137       0.31727       0.51840   NA
    #>     BMA (1)_(Intercept) -1.09319 0.09280      -1.26893      -0.92000 1.00
    #>     BMA          (1)_X1  0.52119 0.10307       0.32635       0.72192 1.02
    #>     BMA          (1)_X2 -0.34994 0.15277      -0.63763      -0.04922 1.01
    #>     BMA        (1)_dose  0.44801 0.05670       0.33765       0.56194 1.00
    #>     BMA (2)_(Intercept)  0.05267 0.07548      -0.09026       0.19230 1.00
    #>     BMA          (2)_X1 -0.51037 0.08339      -0.67078      -0.34844 1.01
    #>     BMA          (2)_X2  0.70371 0.11185       0.49766       0.92805 1.01
    #>     BMA        (2)_dose  0.42386 0.05084       0.32103       0.52702 1.00
    #>   n.eff
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>      NA
    #>  245.00
    #>  390.00
    #>  447.00
    #>  254.00
    #>  324.00
    #>  373.00
    #>  650.00
    #>  313.00

    coef(fit.ameras.multinomial.lin)
    #>                          RC         ERC        MCML         FMA         BMA
    #> (1)_(Intercept) -1.10156112 -1.10285018 -1.09147154 -1.09114141 -1.09318983
    #> (1)_X1           0.52203186  0.52577389  0.52239330  0.52245186  0.52118788
    #> (1)_X2          -0.34709866 -0.35175606 -0.34933932 -0.34937267 -0.34993934
    #> (1)_dose         0.45263847  0.46441427  0.44410129  0.44380519  0.44801466
    #> (2)_(Intercept)  0.04416972  0.03946571  0.05916562  0.05869433  0.05266985
    #> (2)_X1          -0.51311892 -0.51003494 -0.51290413 -0.51271583 -0.51036888
    #> (2)_X2           0.70793136  0.70434790  0.70604308  0.70611143  0.70370686
    #> (2)_dose         0.43049029  0.44543353  0.41753398  0.41768779  0.42386430

    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"))
    #> Note: BMA may require extensive computation time in the order of multiple hours
    #> Fitting RC
    #> Fitting ERC
    #> Fitting MCML
    #> Fitting FMA
    #> Fitting BMA
    #> Defining model
    #> Building model
    #> Setting data and initial values
    #> Running calculate on model
    #>   [Note] Any error reports that follow may simply reflect missing values in model variables.
    #> Checking model sizes and dimensions
    #>   [Note] This model is not fully initialized. This is not an error.
    #>          To see which variables are not initialized, use model$initializeInfo().
    #>          For more information on model initialization, see help(modelInitialization).
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> running chain 1...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|
    #> running chain 2...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|

    summary(fit.ameras.clogit)
    #> Call:
    #> ameras(data = data, family = "clogit", Y = "Y.clogit", dosevars = dosevars,
    #>     X = c("X1", "X2"), setnr = "setnr", methods = c("RC", "ERC",
    #>         "MCML", "FMA", "BMA"), deg = 2, doseRRmod = "EXP", CI = c("wald.orig",
    #>         "percentile"), nburnin.BMA = 1000, niter.BMA = 5000)
    #>
    #> Total run time: 399.2 seconds
    #>
    #> Runtime in seconds by method:
    #>
    #>  Method Runtime
    #>      RC     0.5
    #>     ERC   343.6
    #>    MCML     4.3
    #>     FMA     6.2
    #>     BMA    44.6
    #>
    #> Summary of coefficients by method:
    #>
    #>  Method         Term Estimate      SE CI.lowerbound CI.upperbound Rhat  n.eff
    #>      RC           X1  0.54553 0.08896       0.37116      0.719888   NA     NA
    #>      RC           X2 -0.53392 0.11711      -0.76345     -0.304387   NA     NA
    #>      RC         dose  0.68029 0.10131       0.48173      0.878856   NA     NA
    #>      RC dose_squared -0.05146 0.02242      -0.09540     -0.007522   NA     NA
    #>     ERC           X1  0.61917 0.09205       0.43876      0.799584   NA     NA
    #>     ERC           X2 -0.51784 0.11993      -0.75291     -0.282781   NA     NA
    #>     ERC         dose  0.35155 0.08030       0.19416      0.508940   NA     NA
    #>     ERC dose_squared  0.03687 0.01013       0.01701      0.056735   NA     NA
    #>    MCML           X1  0.55083 0.08923       0.37594      0.725708   NA     NA
    #>    MCML           X2 -0.53547 0.11712      -0.76502     -0.305906   NA     NA
    #>    MCML         dose  0.69334 0.09315       0.51077      0.875908   NA     NA
    #>    MCML dose_squared -0.05581 0.01950      -0.09404     -0.017584   NA     NA
    #>     FMA           X1  0.55060 0.08916       0.37658      0.724789   NA     NA
    #>     FMA           X2 -0.53508 0.11710      -0.76360     -0.306706   NA     NA
    #>     FMA         dose  0.69731 0.09408       0.51428      0.883860   NA     NA
    #>     FMA dose_squared -0.05670 0.01995      -0.09734     -0.018434   NA     NA
    #>     BMA           X1  0.54979 0.08586       0.39285      0.724006 1.00 800.00
    #>     BMA           X2 -0.53448 0.11963      -0.75983     -0.299263 1.01 800.00
    #>     BMA         dose  0.68990 0.09957       0.49151      0.876201 1.00 183.00
    #>     BMA dose_squared -0.05469 0.02140      -0.09535     -0.013915 1.00 160.00

    coef(fit.ameras.clogit)
    #>                       RC         ERC        MCML         FMA         BMA
    #> X1            0.54552627  0.61917374  0.55082666  0.55059936  0.54979498
    #> X2           -0.53391889 -0.51784345 -0.53546549 -0.53508245 -0.53448252
    #> dose          0.68029437  0.35155174  0.69333718  0.69730863  0.68989926
    #> dose_squared -0.05145934  0.03687427 -0.05581249 -0.05670202 -0.05468785

    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"))
    #> Note: BMA may require extensive computation time in the order of multiple hours
    #> Fitting RC
    #> Fitting ERC
    #> Warning in ameras.rc(family = family, dosevars = dosevars, data = data, :
    #> WARNING: Hessian was not invertible or inverse was not positive definite,
    #> variance matrix could not be obtained
    #> Fitting MCML
    #> Fitting FMA
    #> Fitting BMA
    #> Defining model
    #> Building model
    #> Setting data and initial values
    #> Running calculate on model
    #>   [Note] Any error reports that follow may simply reflect missing values in model variables.
    #> Checking model sizes and dimensions
    #>   [Note] This model is not fully initialized. This is not an error.
    #>          To see which variables are not initialized, use model$initializeInfo().
    #>          For more information on model initialization, see help(modelInitialization).
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> Compiling
    #>   [Note] This may take a minute.
    #>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
    #> running chain 1...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|
    #> running chain 2...
    #> |-------------|-------------|-------------|-------------|
    #> |-------------------------------------------------------|

    summary(fit.ameras.clogit.lin)
    #> Call:
    #> ameras(data = data, family = "clogit", Y = "Y.clogit", dosevars = dosevars,
    #>     X = c("X1", "X2"), setnr = "setnr", methods = c("RC", "ERC",
    #>         "MCML", "FMA", "BMA"), deg = 1, doseRRmod = "EXP", CI = c("wald.orig",
    #>         "percentile"), nburnin.BMA = 1000, niter.BMA = 5000)
    #>
    #> Total run time: 199.1 seconds
    #>
    #> Runtime in seconds by method:
    #>
    #>  Method Runtime
    #>      RC     0.3
    #>     ERC   154.2
    #>    MCML     1.6
    #>     FMA     3.0
    #>     BMA    40.0
    #>
    #> Summary of coefficients by method:
    #>
    #>  Method Term Estimate      SE CI.lowerbound CI.upperbound Rhat  n.eff
    #>      RC   X1   0.5463 0.08890        0.3721        0.7206   NA     NA
    #>      RC   X2  -0.5231 0.11692       -0.7522       -0.2939   NA     NA
    #>      RC dose   0.4725 0.04268        0.3888        0.5562   NA     NA
    #>     ERC   X1   0.6405      NA            NA            NA   NA     NA
    #>     ERC   X2  -0.4907      NA            NA            NA   NA     NA
    #>     ERC dose   0.3430      NA            NA            NA   NA     NA
    #>    MCML   X1   0.5514 0.08905        0.3769        0.7259   NA     NA
    #>    MCML   X2  -0.5191 0.11710       -0.7486       -0.2896   NA     NA
    #>    MCML dose   0.4575 0.04293        0.3734        0.5417   NA     NA
    #>     FMA   X1   0.5515 0.08918        0.3764        0.7257   NA     NA
    #>     FMA   X2  -0.5186 0.11723       -0.7487       -0.2885   NA     NA
    #>     FMA dose   0.4578 0.04287        0.3745        0.5421   NA     NA
    #>     BMA   X1   0.5468 0.08828        0.3823        0.7188 1.00 806.00
    #>     BMA   X2  -0.5265 0.11725       -0.7787       -0.3116 1.00 907.00
    #>     BMA dose   0.4578 0.04323        0.3715        0.5430 1.00 549.00

    coef(fit.ameras.clogit.lin)
    #>              RC        ERC       MCML        FMA        BMA
    #> X1    0.5463374  0.6405011  0.5514108  0.5515131  0.5467938
    #> X2   -0.5230641 -0.4906780 -0.5190952 -0.5186443 -0.5265150
    #> dose  0.4724948  0.3429971  0.4575164  0.4578075  0.4577984

    traceplot(fit.ameras.clogit.lin)