Confidence intervals

library(ameras)
#> Loading required package: nimble
#> nimble version 1.4.2 is loaded.
#> For more information on NIMBLE and a User Manual,
#> please visit https://R-nimble.org.
#> 
#> Attaching package: 'nimble'
#> The following object is masked from 'package:stats':
#> 
#>     simulate
#> The following object is masked from 'package:base':
#> 
#>     declare

Introduction

To compute confidence intervals, first fit the model using ameras, then use the confint method to attach confidence intervals. Several types of confidence intervals are supported, which should be supplied to the type argument of confint, see below. When confint is called with methods containing at least one of RC, ERC, and MCML and at least one of FMA and BMA, type should be a vector of length 2 with one method for RC, ERC and MCML and one for FMA and BMA.

Regression calibration, extended regression calibration, and Monte Carlo maximum likelihood

For (extended) regression calibration and Monte Carlo maximum likelihood, Wald and profile likelihood intervals can be obtained. When a parameter transformation \(\mathbf\theta = h(\mathbf\eta)\) is used, type="wald.transformed" yields the CI at significance level \(\alpha\) of \(h(\mathbf\eta \pm z_{1-\alpha/2} \mathbf V)\) where \(z_{1-\alpha/2}\) is the \(1-\alpha/2\)-quantile of the standard normal distribution and \(\mathbf V\) is the vector of standard deviations estimated using the inverse Hessian matrix (returned by optim), and type="wald.orig" uses the delta method to obtain the CI \(h(\mathbf\eta)\pm z_{1-\alpha/2} \mathbf V_*\) where \(\mathbf V_*\) is the vector of standard deviations estimated using \(J H^{-1} J^T\) with \(J\) the Jacobian of the transformation (obtained with transform.jacobian) and \(H\) is the Hessian. When no transformation is used, type="wald.orig" should be used. The third option is type="proflik", which uses the profile likelihood to compute confidence bounds. For this, the profile log (partial) likelihood for parameter \(\theta_p\) is defined as \[ PL_p (\theta_p^*) = \max_{\mathbf\theta: \theta_p = \theta_p^*} \ell (\mathbf \theta), \] where \(\ell\) is the log (partial) likelihood. Next, profile confidence intervals \((\theta_p^l, \theta_p^h)\) are obtained for parameter \(\theta_p\) at significance level \(\alpha\) by solving \(-2 \{PL_p(\theta_p^*) - \ell(\hat{\mathbf{\theta}})\}=\chi^2_{1,1-\alpha}\) using the bisection method, with \(\hat{\mathbf{\theta}}\) the maximum likelihood estimate. Note that profile likelihoods are more computationally intensive to obtain. For this reason, confint offers the option to only determine them for the exposure-related parameters, which is the default setting. To obtain profile likelihood intervals for all parameters, use parm = "all".

To illustrate, we determine the three types of confidence intervals for a regression calibration analysis using the example data, using a linear excess relative risk model with the default exponential transformation (see Parameter transformations).

data(data, package="ameras")
dosevars <- paste0("V", 1:10)
fit.ameras <- ameras(Y.binomial~dose(V1:V10, model="ERR")+X1+X2, data=data, 
                            family="binomial", methods=c("RC"))
#> Fitting RC

fit.ameras.waldorig <- confint(fit.ameras, type="wald.orig")
fit.ameras.waldtransformed <- confint(fit.ameras, type="wald.transformed")
fit.ameras.proflik <- confint(fit.ameras, type="proflik", parm="all")
#> Obtaining profile likelihood CI for (Intercept)
#> Obtaining profile likelihood CI for X1
#> Obtaining profile likelihood CI for X2
#> Obtaining profile likelihood CI for dose

summary(fit.ameras.waldorig)
#> Call:
#> ameras(formula = Y.binomial ~ dose(V1:V10, model = "ERR") + X1 + 
#>     X2, data = data, family = "binomial", methods = c("RC"))
#> 
#> Total run time: 0.2 seconds
#> 
#> Runtime in seconds by method:
#> 
#>  Method Runtime
#>      RC     0.2
#> 
#> Summary of coefficients by method:
#> 
#>  Method        Term Estimate      SE CI.lowerbound CI.upperbound
#>      RC (Intercept)  -1.0641 0.08788       -1.2363       -0.8918
#>      RC          X1   0.4409 0.07628        0.2914        0.5904
#>      RC          X2  -0.3360 0.09544       -0.5230       -0.1489
#>      RC        dose   0.8508 0.14517        0.5663        1.1353
summary(fit.ameras.waldtransformed)
#> Call:
#> ameras(formula = Y.binomial ~ dose(V1:V10, model = "ERR") + X1 + 
#>     X2, data = data, family = "binomial", methods = c("RC"))
#> 
#> Total run time: 0.2 seconds
#> 
#> Runtime in seconds by method:
#> 
#>  Method Runtime
#>      RC     0.2
#> 
#> Summary of coefficients by method:
#> 
#>  Method        Term Estimate      SE CI.lowerbound CI.upperbound
#>      RC (Intercept)  -1.0641 0.08788       -1.2363       -0.8918
#>      RC          X1   0.4409 0.07628        0.2914        0.5904
#>      RC          X2  -0.3360 0.09544       -0.5230       -0.1489
#>      RC        dose   0.8508 0.14517        0.6050        1.1827
summary(fit.ameras.proflik)
#> Call:
#> ameras(formula = Y.binomial ~ dose(V1:V10, model = "ERR") + X1 + 
#>     X2, data = data, family = "binomial", methods = c("RC"))
#> 
#> Total run time: 0.2 seconds
#> 
#> Runtime in seconds by method:
#> 
#>  Method Runtime
#>      RC     0.2
#> 
#> Summary of coefficients by method:
#> 
#>  Method        Term Estimate      SE CI.lowerbound CI.upperbound pval.lower
#>      RC (Intercept)  -1.0641 0.08788       -1.2391       -0.8971    0.05248
#>      RC          X1   0.4409 0.07628        0.2908        0.5917    0.04878
#>      RC          X2  -0.3360 0.09544       -0.5256       -0.1490    0.04822
#>      RC        dose   0.8508 0.14517        0.6009        1.1784    0.05040
#>  pval.upper
#>     0.05153
#>     0.04853
#>     0.04902
#>     0.04952

Frequentist and Bayesian model averaging

For frequentist and Bayesian model averaging methods, the options are percentile which uses equal-tailed quantiles, and hpd which computes highest posterior density intervals using HPDinterval from the coda package, using either the FMA samples or Bayesian posterior samples.

Again, we use the example data to illustrate.

set.seed(123)
fit.ameras2 <- ameras(Y.binomial~dose(V1:V10, model="ERR")+X1+X2, data=data, 
                            family="binomial", methods=c("FMA"))
#> Fitting FMA

fit.ameras.hpd <- confint(fit.ameras2, type="hpd")
fit.ameras.percentile <- confint(fit.ameras2, type="percentile")

summary(fit.ameras.hpd)
#> Call:
#> ameras(formula = Y.binomial ~ dose(V1:V10, model = "ERR") + X1 + 
#>     X2, data = data, family = "binomial", methods = c("FMA"))
#> 
#> Total run time: 0.8 seconds
#> 
#> Runtime in seconds by method:
#> 
#>  Method Runtime
#>     FMA     0.8
#> 
#> Summary of coefficients by method:
#> 
#>  Method        Term Estimate      SE CI.lowerbound CI.upperbound
#>     FMA (Intercept)  -1.0576 0.08770       -1.2284       -0.8856
#>     FMA          X1   0.4429 0.07652        0.2939        0.5937
#>     FMA          X2  -0.3378 0.09572       -0.5260       -0.1530
#>     FMA        dose   0.8447 0.14494        0.5640        1.1317
summary(fit.ameras.percentile)
#> Call:
#> ameras(formula = Y.binomial ~ dose(V1:V10, model = "ERR") + X1 + 
#>     X2, data = data, family = "binomial", methods = c("FMA"))
#> 
#> Total run time: 0.8 seconds
#> 
#> Runtime in seconds by method:
#> 
#>  Method Runtime
#>     FMA     0.8
#> 
#> Summary of coefficients by method:
#> 
#>  Method        Term Estimate      SE CI.lowerbound CI.upperbound
#>     FMA (Intercept)  -1.0576 0.08770       -1.2290       -0.8860
#>     FMA          X1   0.4429 0.07652        0.2938        0.5935
#>     FMA          X2  -0.3378 0.09572       -0.5250       -0.1518
#>     FMA        dose   0.8447 0.14494        0.5616        1.1298