Score-based confidence intervals for random effect covariance parameters: FEV1 example

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:

library(reconf)
library(lme4)
#> Loading required package: Matrix

data(fev1)
dim(fev1)
#> [1] 1994    6
head(fev1)
#>   id   ht     age baseht baseage logfev1
#> 1  1 1.20  9.3415    1.2  9.3415 0.21511
#> 2  1 1.28 10.3929    1.2  9.3415 0.37156
#> 3  1 1.33 11.4524    1.2  9.3415 0.48858
#> 4  1 1.42 12.4600    1.2  9.3415 0.75142
#> 5  1 1.48 13.4182    1.2  9.3415 0.83291
#> 6  1 1.50 15.4743    1.2  9.3415 0.89200

Each subject has between 1 and 12 observations over time.

range(table(fev1$id))
#> [1]  1 12

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.

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

Maximum likelihood estimates

round(fixef(fit), 3)
#> (Intercept)         age          ht     baseage      baseht 
#>      -2.049       0.019       1.636      -0.022       0.273

as.data.frame(VarCorr(fit), order = "lower.tri")[
  , c("grp", "var1", "var2", "vcov")]
#>        grp        var1 var2          vcov
#> 1       id (Intercept) <NA>  1.453194e-02
#> 2       id (Intercept)  age -4.725486e-04
#> 3       id         age <NA>  5.347985e-05
#> 4 Residual        <NA> <NA>  3.496540e-03

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}\).

ci <- vc_ci(fit)
ci
#> 95% confidence intervals for variance-covariance parameters (score, ML)
#> 
#>                           estimate       lower       upper
#> var_(Intercept)|id      0.01453194  0.01061200  0.02034969
#> cov_age.(Intercept)|id -0.00047255 -0.00082375 -0.00024082
#> var_age|id              0.00005348  0.00003678  0.00007862
#> var_Residual            0.00349654  0.00325851  0.00377196

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

vc_ci(fit, parm = "var_age|id")   # random slope variance psi_3
#> 95% confidence intervals for variance-covariance parameters (score, ML)
#> 
#>              estimate      lower      upper
#> var_age|id 0.00005348 0.00003678 0.00007862

Comparison with profile likelihood

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

confint(fit, method = "profile", parm = "theta_")

#>               2.5 %       97.5 %
#> .sig01  0.000000000  0.139517390
#> .sig02  0.006037058  0.007818737
#> .sig03 -0.671833222 -0.353296130
#> .sigma  0.057054602  0.061757510

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}\):

fit_reml <- lmer(
  logfev1 ~ age + ht + baseage + baseht + (age | id),
  data = fev1, REML = TRUE,
  control = lmerControl(optimizer = "bobyqa")
)
vc_test(fit_reml, parm = 3L)
#> 
#>  Score test (REML)
#> 
#> data:  fit_reml
#> S = 117.2, df = 1, p-value < 2.2e-16
#> alternative hypothesis: true var_age|id is not equal to 0

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.

vc_test(fit_reml, parm = c(2L, 3L))
#> 
#>  Score test (REML)
#> 
#> data:  fit_reml
#> S = 127.41, df = 2, p-value < 2.2e-16
#> alternative hypothesis: two.sided
#> null values:
#> cov_age.(Intercept)|id             var_age|id 
#>                      0                      0

References

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