---
title: "Score-based confidence intervals for random effect covariance parameters: FEV1 example"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Score-based confidence intervals for random effect covariance parameters: FEV1 example}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
# lme4's profile confidence intervals take ~25 s, which is too slow for the
# time budget of a routine package check. Compute them only when NOT_CRAN is
# set (locally and in CI); otherwise show the precomputed result below.
run_profile <- !identical(Sys.getenv("NOT_CRAN"), "")
```

## Overview

This vignette computes score-based confidence intervals for the three
covariance parameters of the random effects in a model for longitudinal lung
function data: two variances and one covariance. The profile likelihood
intervals of `lme4` are on the standard deviation and correlation scale;
`reconf` gives intervals on the variance and covariance scale.

## Data

The `fev1` dataset contains 1994 repeated measurements of forced expiratory
volume in one second (FEV1) for 300 girls from Topeka, Kansas, described in
Fitzmaurice, Laird, and Ware (2011). Each row is one observation; the
variables are:

| Variable   | Description |
|------------|-------------|
| `id`       | Subject identifier |
| `age`      | Age at measurement (years) |
| `ht`       | Height at measurement (meters) |
| `baseage`  | Baseline age (first visit) |
| `baseht`   | Baseline height (first visit) |
| `logfev1`  | Log FEV1 (natural log, liters) |

```{r load-data}
library(reconf)
library(lme4)

data(fev1)
dim(fev1)
head(fev1)
```

Each subject has between 1 and 12 observations over time.

```{r visits-per-subject}
range(table(fev1$id))
```

## Model

We fit the linear mixed model

$$
\text{logfev1}_{it} = \beta_1 + \beta_2\,\text{age}_{it} + \beta_3\,\text{ht}_{it}
  + \beta_4\,\text{baseage}_i + \beta_5\,\text{baseht}_i
  + U_{i1} + U_{i2}\,\text{age}_{it} + E_{it},
$$

where $(U_{i1}, U_{i2})^T \sim N(0, \Psi)$ independently across subjects,
$E_{it} \sim N(0, \psi_4)$ independently of the random effects, and

$$
\Psi = \begin{pmatrix} \psi_{1} & \psi_{2} \\ \psi_{2} & \psi_{3} \end{pmatrix}.
$$

The random effect covariance parameters $(\psi_{1}, \psi_{2}, \psi_{3})$
are of primary interest here.

```{r fit-model}
fit <- lmer(
  logfev1 ~ age + ht + baseage + baseht + (age | id),
  data = fev1,
  REML = FALSE,
  control = lmerControl(optimizer = "bobyqa")
)
```

## Maximum likelihood estimates

```{r mle}
round(fixef(fit), 3)

as.data.frame(VarCorr(fit), order = "lower.tri")[
  , c("grp", "var1", "var2", "vcov")]
```

The estimated random slope variance $\hat{\psi}_{3}$ is small.

## Score-based confidence intervals

`vc_ci()` computes score-based 95% confidence intervals for all
covariance parameters, including the error variance.  Each row is named by
the corresponding variance or covariance and shows the estimate alongside
the bounds; the random-effect parameters appear in lower-triangular order,
so indices 1, 2, 3 refer to $\psi_{1}$, $\psi_{2}$, $\psi_{3}$.

```{r ci-all}
ci <- vc_ci(fit)
ci
```

The argument `parm` selects parameters by name or index, and only those
intervals are computed.

```{r ci-single}
vc_ci(fit, parm = "var_age|id")   # random slope variance psi_3
```

## Comparison with profile likelihood

For comparison, `lme4::confint()` with `method = "profile"` gives profile
likelihood intervals on the standard deviation and correlation scale:

```{r lr-compare, eval = run_profile}
confint(fit, method = "profile", parm = "theta_")
```

```{r lr-compare-static, eval = !run_profile, echo = FALSE}
# Precomputed output of the chunk above (shown when the profile computation
# is skipped for speed).
cat("              2.5 %       97.5 %",
    ".sig01  0.000000000  0.139517390",
    ".sig02  0.006037058  0.007818737",
    ".sig03 -0.671833222 -0.353296130",
    ".sigma  0.057054602  0.061757510", sep = "\n")
```

The `.sig01` and `.sig03` rows are standard deviations ($\sqrt{\psi_{1}}$
and $\sqrt{\psi_{3}}$); `.sig02` is the correlation
$\psi_{2}/\sqrt{\psi_{1}\psi_{3}}$.

## Score tests

`vc_test()` with a single index tests that parameter against zero,
with the other parameters as nuisance parameters.  For the random slope
variance $\psi_{3}$:

```{r score-test-one}
fit_reml <- lmer(
  logfev1 ~ age + ht + baseage + baseht + (age | id),
  data = fev1, REML = TRUE,
  control = lmerControl(optimizer = "bobyqa")
)
vc_test(fit_reml, parm = 3L)
```

Several indices give a joint test.  That the random effects for age
contribute nothing is $H_0: \psi_{2} = \psi_{3} = 0$.  When
$\psi_{3} = 0$, the slope $U_{i2}$ is identically zero, which forces
$\psi_{2} = 0$ as well, so both are in the null.

```{r score-test}
vc_test(fit_reml, parm = c(2L, 3L))
```

## References

- Fitzmaurice, G. M., Laird, N. M., and Ware, J. H. (2011).
  *Applied Longitudinal Analysis*, 2nd edition. Wiley.