| Title: | Confidence Intervals for Random Effect Parameters in Linear Mixed Models |
|---|---|
| Description: | Computes score-based confidence intervals and hypothesis tests for variance and covariance parameters of random effects in linear mixed models fitted with 'lme4'. The intervals invert a profile score statistic on the variance-covariance scale, with nuisance parameters estimated on an extended parameter set. The intervals have near-nominal coverage uniformly in the parameter, including near the boundary of the parameter space, where Wald and likelihood ratio intervals can be unreliable. The methods are described in Shedden and Ekvall (2026) <doi:10.48550/arXiv.2610.04181>, building on Ekvall and Bottai (2026) <doi:10.1093/biomet/asaf079>, Zhang, Ekvall, and Molstad (2025) <doi:10.1093/biomet/asaf010>, and Ekvall and Bottai (2022) <doi:10.1214/22-AOS2177>. |
| Authors: | Karl Oskar Ekvall [aut, cre, cph], Yiqiao Zhang [ctb], Matias Shedden [ctb] |
| Maintainer: | Karl Oskar Ekvall <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.2.0 |
| Built: | 2026-10-23 18:46:48 UTC |
| Source: | https://github.com/koekvall/reconf |
Repeated measurements of forced expiratory volume in one second (FEV1) on a random sample of 300 girls from Topeka, Kansas, in the Six Cities Study of Air Pollution and Health, with one to twelve measurements per girl.
fev1fev1
A data frame with 1994 rows and 6 variables:
subject identifier.
height (meters).
age (years).
height at the first measurement (meters).
age at the first measurement (years).
natural logarithm of FEV1 (liters).
Fitzmaurice, G. M., Laird, N. M., and Ware, J. H. (2011). Applied Longitudinal Analysis, 2nd edition. Wiley. https://content.sph.harvard.edu/fitzmaur/ala2e/. Data courtesy of D. W. Dockery.
Dockery, D. W., Berkey, C. S., Ware, J. H., Speizer, F. E., and Ferris, B. G. (1983). Distribution of FVC and FEV1 in children 6 to 11 years old. American Review of Respiratory Disease 128, 405–412.
Print method for confidence intervals
## S3 method for class 'vc_ci' print(x, digits = 4, ...)## S3 method for class 'vc_ci' print(x, digits = 4, ...)
x |
A |
digits |
Number of significant digits. Default 4. |
... |
Unused. |
x, invisibly.
Returns the intervals as a data frame with the column names of the broom package.
## S3 method for class 'vc_ci' tidy(x, ...)## S3 method for class 'vc_ci' tidy(x, ...)
x |
A |
... |
Unused. |
A data frame with columns term, estimate,
conf.low, and conf.high.
Computes confidence intervals for covariance parameters of a linear mixed model by inverting a signed profile score or likelihood ratio statistic. See the references for the theory.
vc_ci( object, parm = NULL, level = 0.95, statistic = c("score", "rlrt"), expected = TRUE, known = NULL, nonneg = TRUE, onestep = FALSE, step_size = NULL, num_points = 500L, method = c("auto", "q_side", "n_side", "spectral"), ... )vc_ci( object, parm = NULL, level = 0.95, statistic = c("score", "rlrt"), expected = TRUE, known = NULL, nonneg = TRUE, onestep = FALSE, step_size = NULL, num_points = 500L, method = c("auto", "q_side", "n_side", "spectral"), ... )
object |
a fitted linear mixed model of class |
parm |
parameters, as names (see Value) or as indices in the order
of |
level |
confidence level. |
statistic |
|
expected |
if |
known |
parameters, as for |
nonneg |
if |
onestep |
if |
step_size |
initial step length of the search (see Details), on the
scale of the parameter. If |
num_points |
maximum number of steps in each direction. The search
also stops at distance |
method |
computational method for the likelihood: |
... |
arguments passed to the optimizer |
The interval at level for a parameter is
where is the signed statistic with the other parameters profiled
out and is the standard normal quantile. The
parameter set is the set of covariance parameters for which the covariance
matrix of the response is positive definite. It includes negative
variances, and the nuisance parameters are maximized over it.
The bounds are found by a search outward from the estimate in each
direction, with the nuisance parameters at each step computed from those
at the previous step. Step lengths are chosen by secant extrapolation of
the statistic, each at most twice the previous one. If the covariance
matrix of the response is not positive definite at the start of a step,
the step is halved, up to 20 times, after which a starting point is
constructed as in vc_test. A bound is where the
statistic crosses , located to within
step_size. If there is no crossing within the search distance, the
bound is -Inf or Inf and a warning is issued, reporting
the statistic at the end of the search. As a variance grows, the
statistic tends to a finite limit, which for a grouping factor with few
levels can be above ; the upper bound is then
infinite for any search distance. If the statistic at the estimate
exceeds the critical value, the search starts instead from the maximizer
over the parameter set, with a warning.
With statistic = "rlrt", the search starts from the maximizer over
the parameter set. A warning is issued if that maximizer has a negative
variance; the interval for that variance can then lie entirely below zero.
The likelihood is computed by one of three methods. "q_side" uses
sparse matrices, where is the number of random
effects. "n_side" uses dense matrices. For models
with one random-effect variance, "spectral" uses one
eigendecomposition and operations per evaluation.
"auto" uses "spectral" or "n_side" if
and more than 10 percent of the entries of the
random-effects design matrix are nonzero, and "q_side" otherwise.
A matrix of class vc_ci with one row per parameter and
columns estimate, lower, and upper. For an lme4
fit the estimates are those of lme4, and rows are named
"var_x|g" for the variance of the random effect x with
grouping factor g, "cov_y.x|g" for the covariance of
x and y, and "var_Residual" for the error
variance, following the names lme4 uses for standard deviations
and correlations.
Shedden, M. and Ekvall, K. O. (2026). Score-based confidence intervals for variance-covariance parameters in linear mixed models. arXiv:2610.04181.
library(lme4) fit <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy) # 95% intervals for all covariance parameters vc_ci(fit) # The random slope variance, by name or by index vc_ci(fit, parm = "var_Days|Subject") vc_ci(fit, parm = 3)library(lme4) fit <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy) # 95% intervals for all covariance parameters vc_ci(fit) # The random slope variance, by name or by index vc_ci(fit, parm = "var_Days|Subject") vc_ci(fit, parm = 3)
Returns the log-likelihood of a linear mixed model, or its restricted log-likelihood for a restricted likelihood fit, as a function of the covariance parameters.
vc_loglik(object, method = c("auto", "q_side", "n_side", "spectral"))vc_loglik(object, method = c("auto", "q_side", "n_side", "spectral"))
object |
a fitted linear mixed model of class |
method |
computational method for the likelihood: |
The likelihood is maximized over the fixed effects. Outside the
parameter set, the set of covariance parameters for which the
covariance matrix of the response is positive definite, the value is
-Inf and the score and information are zero. With prior weights
in an lme4 fit, the value differs from logLik() by
.
A function of psi, the covariance parameters in the order
of vc_model, with logical arguments get_val,
get_score, and get_inf selecting what to compute and
expected selecting the expected or the observed information, all
TRUE by default. It returns a list with components value,
score, and information, the last two named by the
parameters; components not requested are NULL.
library(lme4) fit <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy) ll <- vc_loglik(fit) ll(c(600, 10, 35, 650))library(lme4) fit <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy) ll <- vc_loglik(fit) ll(c(600, 10, 35, 650))
Builds the object that vc_ci, vc_test, and
vc_loglik work with: the likelihood of a linear mixed model
as a function of its covariance parameters, with quantities precomputed
once, and the parameter names and estimates. Those functions also accept
a fitted model and build this object themselves; building it once avoids
repeating the precomputation across calls.
## S3 method for class 'lmerMod' vc_model(x, method = c("auto", "q_side", "n_side", "spectral"), ...) vc_model(x, ...) ## S3 method for class 'vc_model' print(x, digits = 4, ...)## S3 method for class 'lmerMod' vc_model(x, method = c("auto", "q_side", "n_side", "spectral"), ...) vc_model(x, ...) ## S3 method for class 'vc_model' print(x, digits = 4, ...)
x |
a fitted linear mixed model of class |
method |
computational method for the likelihood; see
|
... |
unused. |
digits |
number of significant digits printed. |
For an lme4 fit the likelihood is the restricted likelihood if
the model was fitted by REML, and the estimates are those of lme4. The
parameters are ordered as in
as.data.frame(VarCorr(x), order = "lower.tri"), with the error
variance last. Prior weights and offsets are supported; as in
lme4, observation then has error variance
, with the error variance parameter.
An object of class "vc_model". Printing it shows the
covariance parameters, in the order that indices in vc_ci and
the other functions refer to, and their estimates.
library(lme4) fit <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy) m <- vc_model(fit) m vc_ci(m, parm = "var_Days|Subject") vc_test(m, parm = "var_Days|Subject")library(lme4) fit <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy) m <- vc_model(fit) m vc_ci(m, parm = "var_Days|Subject") vc_test(m, parm = "var_Days|Subject")
Computes the score test statistic for a hypothesis about the covariance parameters of a linear mixed model.
vc_test( object, parm = NULL, null_value = 0, expected = TRUE, known = NULL, method = c("auto", "q_side", "n_side", "spectral"), ... )vc_test( object, parm = NULL, null_value = 0, expected = TRUE, known = NULL, method = c("auto", "q_side", "n_side", "spectral"), ... )
object |
a fitted linear mixed model of class |
parm |
the tested parameters, as names or indices as in
|
null_value |
values of the tested parameters under the null
hypothesis, recycled to the length of |
expected |
if |
known |
parameters, as for |
method |
computational method for the likelihood: |
... |
arguments passed to the optimizer |
The statistic is the profile score statistic: the score for the
tested parameters, evaluated at the maximizer of the likelihood with the
tested parameters at their null values and the parameters in
known at their estimates, and standardized by the efficient
information. The maximization is over the other covariance parameters
and, for a maximum likelihood fit, the fixed effects. It starts from the
estimates in the model, with the covariances of a random effect
whose variance is tested at zero set to zero. If the covariance matrix
of the response is not positive definite at this start, free entries
of are changed to make its blocks weakly diagonally
dominant, and the error variance is doubled until it is.
An object of class "htest" with components
statistic; parameter, the degrees of freedom, which is the
number of tested parameters; p.value, from the chi-squared
distribution; null.value, the null values named as in
vc_ci; alternative; method; and
data.name.
Shedden, M. and Ekvall, K. O. (2026). Score-based confidence intervals for variance-covariance parameters in linear mixed models. arXiv:2610.04181.
library(lme4) fit <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy) # All random-effect parameters zero vc_test(fit) # Random slope variance zero vc_test(fit, parm = "var_Days|Subject")library(lme4) fit <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy) # All random-effect parameters zero vc_test(fit) # Random slope variance zero vc_test(fit, parm = "var_Days|Subject")