---
title: "Confidence intervals"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Confidence intervals}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(ameras)

```

# 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](transformations.html)).

```{r fits1, eval = identical(Sys.getenv("NOT_CRAN"), "true")}
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"))

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")

summary(fit.ameras.waldorig)
summary(fit.ameras.waldtransformed)
summary(fit.ameras.proflik)
```

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

```{r fits2, eval = identical(Sys.getenv("NOT_CRAN"), "true")}
set.seed(123)
fit.ameras2 <- ameras(Y.binomial~dose(V1:V10, model="ERR")+X1+X2, data=data, 
                            family="binomial", methods=c("FMA"))

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

summary(fit.ameras.hpd)
summary(fit.ameras.percentile)
```