The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

TestREnlme: Testing Random Effects in Linear and Nonlinear Mixed-Effects Models

Germaine Uwimpuhwe, Reza Drikvandi, Shelley A. Blozis

Introduction

TestREnlme provides a nonparametric framework for testing random effects in linear mixed-effects models (LMEMs) and nonlinear mixed-effects models (NLMEMs) with additive random error. The package implements the permutation tests of Uwimpuhwe, Drikvandi and Blozis (2026, Statistics inMedicine, doi:[10.1002/sim.70605](https://doi.org/10.1002/sim.70605)), allowing users to test all random effects jointly or selected subsets of random effects without relying on distributional assumptions. The package provides three distribution-free variance component estimators:

Key features of TestREnlme:

The package can be installed from CRAN using:

install.packages("TestREnlme")
library(TestREnlme)

Model Formulation

The NLMEM is \[ Y_i = f_i(a_i, \gamma) + \varepsilon_i, \qquad a_i = A_i\beta + b_i, \] where \(f_i(\cdot)\) is a known nonlinear function, \(\gamma\) is a vector of common fixed parameters, and \(\beta\) contains fixed effects paired with random effects structure. The random effects \(b_i\) have mean zero and covariance \(D_* = \sigma^2 D\), and the random errors \(\varepsilon_i\) have mean zero and variance \(\sigma^2\). if \(f_i(\cdot)=X_i\beta+Z_ib_i\) the NLMEM reduces to the LMEM.

The model is defined through the following arguments:

Theophylline Pharmacokinetic Data

Data and Model

The built-in Theoph dataset contains serum concentration measurements for 12 patients. The one-compartment pharmacokinetic model is: \[ Y_{ij} = \frac{D_i \exp(a_{i2} + a_{i3} - a_{i1}) \big[\exp(-e^{a_{i3}} T_{ij}) - \exp(-e^{a_{i2}} T_{ij})\big]} {\exp(a_{i2}) - \exp(a_{i3})} + \varepsilon_{ij} \] where \(a_{i1} = \beta_1 + b_{i1}\), \(a_{i2} = \beta_2 + b_{i2}\), and \(a_{i3} = \beta_3 + b_{i3}\) are the subject-specific log clearance, log absorption rate, and log elimination rate, respectively.

Exploratory Plot

d      <- as.data.frame(Theoph)
p <- plot_profiles(d, group = "Subject", time = "Time",
                   response = "conc",
                   title = "Theophylline: individual profiles")
## Colour by subject using ggplot2 customisation
p + ggplot2::aes(colour = .data$`.group`) +
    ggplot2::scale_colour_manual(values = rainbow(12), guide = "none")
Individual concentration profiles with mean (black and bold).
Individual concentration profiles with mean (black and bold).

Testing All Random Effects

Step 1: Model preparation

Expr   <- conc ~ Dose * exp(ai2 + ai3 - ai1) *
            (exp(-Time * exp(ai3)) - exp(-Time * exp(ai2))) /
            (exp(ai2) - exp(ai3))

random <- c("ai1 ~ B1 + bi1",
            "ai2 ~ B2 + bi2", 
            "ai3 ~ B3 + bi3")

start  <- c(ai1 = -3.22, ai2 = 0.47, ai3 = -2.45)

Step 2: Variance components estimation

## VLS (default: no subject-specific estimation)
DVLS <- Dmethod(d, Expr, group = "Subject",
                random = random, start = start,
                verbose = 0)
summary(DVLS)
#> 
#> TestREnlme: Variance component estimates
#> Method : VLS 
#> 
#> Fixed effects (Beta):
#>        B1        B2        B3 
#> -3.265220  0.454294 -2.520843 
#> 
#> Error variance (Sigma2):
#>   0.530631 
#> 
#> Variance-covariance matrix (Dhat):
#>          bi1      bi2      bi3
#> bi1 0.043393 0.024829 0.011664
#> bi2 0.024829 0.289258 0.015266
#> bi3 0.011664 0.015266 0.003404

When both MM and MMF are needed on the same data, the expensive first-stage NLS fits can be computed once via MM_base() and shared between both methods, substantially reducing computation time:

mb   <- MM_base(d, Expr, group = "Subject",
                random = random, start = start,
                kappa_max = 1e4, verbose = 1)
DMM  <- Dmethod(d, Expr, group = "Subject",
                random = random, start = start,
                method = "MM", MM_base_obj = mb, verbose = 0)


DMMF <- Dmethod(d, Expr, group = "Subject",
                random = random, start = start,
                method = "MMF", MM_base_obj = mb, verbose = 0)
summary(DMM)
#> 
#> TestREnlme: Variance component estimates
#> Method : MM 
#> 
#> Fixed effects (Beta):
#>        B1        B2        B3 
#> -3.266823  0.447620 -2.529429 
#> 
#> Error variance (Sigma2):
#>   0.490269 
#> 
#> Variance-covariance matrix (Dhat):
#>           bi1       bi2       bi3
#> bi1  0.067649 -0.016485  0.036102
#> bi2 -0.016485  0.478569 -0.001676
#> bi3  0.036102 -0.001676  0.019373
#> 
#> Subject summary:
#>   Total           : 12 
#>   Used            : 12 
#>   High-kappa excl.: 0 
#>   Non-converged   : 0
summary(DMMF)
#> 
#> TestREnlme: Variance component estimates
#> Method : MMF 
#> 
#> Fixed effects (Beta):
#>        B1        B2        B3 
#> -3.267043  0.441197 -2.530169 
#> 
#> Error variance (Sigma2):
#>   0.490269 
#> 
#> Variance-covariance matrix (Dhat):
#>           bi1       bi2       bi3
#> bi1  0.090767 -0.022465  0.048952
#> bi2 -0.022465  0.640604 -0.003296
#> bi3  0.048952 -0.003296  0.027208
#> 
#> Subject summary:
#>   Total           : 12 
#>   Used            : 12 
#>   High-kappa excl.: 0 
#>   Non-converged   : 0
## Compare variance component estimates across methods
plot_variance(list(VLS = DVLS, MM = DMM, MMF = DMMF), component = "diagonal")

Step 3: Bootstrap standard errors.

Once variance components have been estimated, uncertainty in all model parameters can be assessed using bootstrap_se():

BVLS <- bootstrap_se(DVLS, nboot = 200, seed = 42)
print(BVLS)

Step 4: Permutation tests

We consider three permutation tests: (i) testing all random effects jointly, (ii) testing \(b_{i1}\) and \(b_{i3}\) while retaining \(b_{i2}\), and (iii) testing \(b_{i3}\) alone while retaining \(b_{i1}\) and \(b_{i2}\). Since all tests use the same variance estimate, Dhatt is computed once via Dmethod() and reused across all calls.

## Test 1: H0: D = 0 (all RE absent)
HVLS <- Dhypothesis_test(d, Expr, group = "Subject",
                          random = random, start = start,
                          Dhatt = DVLS, nperm = 200, seed = 1,
                          verbose = 0)
cat("Test 1 - p-value:", HVLS$pvalue, "| Decision:", HVLS$Decision, "\n")
#> Test 1 - p-value: < 0.001 | Decision: Reject H0

## Test 2: H0: d11 = d13 = d33 = 0 given  bi2 retained in the model
HVLS13 <- Dhypothesis_test(d, Expr, group = "Subject",
                            random = random, start = start,
                            bi_out = c("bi1", "bi3"),
                            Dhatt = DVLS, nperm = 200, seed = 1,
                            verbose = 0)
cat("Test 2 - p-value:", HVLS13$pvalue, "| Decision:", HVLS13$Decision, "\n")
#> Test 2 - p-value: < 0.001 | Decision: Reject H0

## Test 3: H0: d33 = 0 given bi1 and bi2 retained in the model
HVLS3 <- Dhypothesis_test(d, Expr, group = "Subject",
                           random = random, start = start,
                           bi_out = "bi3",
                           Dhatt = DVLS, nperm = 200, seed = 1,
                           verbose = 0)
cat("Test 3 - p-value:", HVLS3$pvalue, "| Decision:", HVLS3$Decision, "\n")
#> Test 3 - p-value: 0.91 | Decision: Do not reject H0

Permutation histograms

library(ggpubr)
ggarrange(
  plot_perm_hist(HVLS,   title = "H0: D = 0"),
  plot_perm_hist(HVLS13, title = "H0: d11=d13=d33=0"),
  plot_perm_hist(HVLS3,  title = "H0: d33=0"),
  nrow = 1
)
Permutation null distributions for the three tests.
Permutation null distributions for the three tests.

Test 1 and Test 2 strongly reject \(H_0\) (\(p < 0.001\)). Test 3 fails to reject (\(p =\) 0.91), suggesting the elimination rate \(K_{ei}\) can be treated as a common fixed parameter across subjects.

Reduced Random Effect Strureture.

Since \(b_{i3}\) is not required, we refit the model retaining only \(b_{i1}\) and \(b_{i2}\). The parameter \(B_3\) now appears directly in Expr2 as a common fixed parameter.

Step 1: Model preparation

Expr2   <- conc ~ Dose * exp(ai2 + B3 - ai1) *
             (exp(-Time * exp(B3)) - exp(-Time * exp(ai2))) /
             (exp(ai2) - exp(B3))
random2 <- c("ai1 ~ B1 + bi1",
             "ai2 ~ B2 + bi2")
start2 <- c(ai1 = -3.22, ai2 = 0.47, B3 = -2.45)

Step 2: Variance components estimation

DVLS_C <- Dmethod(d, Expr2, group = "Subject",
                  random = random2, start = start2,
                  verbose = 0)
summary(DVLS_C)
#> 
#> TestREnlme: Variance component estimates
#> Method : VLS 
#> 
#> Fixed effects (Beta):
#>        B1        B2        B3 
#> -3.263142  0.454111 -2.517851 
#> 
#> Error variance (Sigma2):
#>   0.567943 
#> 
#> Variance-covariance matrix (Dhat):
#>          bi1      bi2
#> bi1 0.022460 0.027285
#> bi2 0.027285 0.307774

Diagnostic plots

Observed versus fitted values can be visualised using plot_fitted(): In the example below, we set subjects = NULL (the default), meaning that all subjects are plotted. Alternatively, a subset of subjects can be selected for visualisation.

plot_fitted(DVLS_C, time = "Time", subjects = NULL, overlay = TRUE)
Observed vs fitted values for all 12 subjects (reduced model).
Observed vs fitted values for all 12 subjects (reduced model).

Calling plot() on a Dmethod object produces all relevant diagnostic plots (profiles, fitted values, raw and standardised residuals; and condition-number plot for VLS based objects):

## Shows: profiles, fitted, residuals (raw + standardised)
## For MM/MMF also shows condition-number plot
plot(x=DVLS_C, time = "Time")

Condition-number diagnostics

For MM or MMF fits, plot_condition() shows the per-subject condition numbers with the exclusion threshold. In the reduced model, all 12 subjects are well below the threshold kappa_max = 1e4:

DMM_C <- Dmethod(d, Expr2, group = "Subject",
                 random = random2, start = start2,
                 method = "MM", verbose = 0)
plot_condition(DMM_C, kappa_max = 1e4)
Condition numbers for the reduced MM model. All subjects are below the exclusion threshold.
Condition numbers for the reduced MM model. All subjects are below the exclusion threshold.

Step 3: Bootstrap standard errors.

Bootstrap standard errors follow the same codes as in the full model, replacing DVLS with DVLS_C in the bootstrap_se() call.

Step 4: Permutation

Both retained random effects are highly significant:

## Test for the need of all random effects in the reduced model
HC <- Dhypothesis_test(d, Expr2, group = "Subject",
                        random = random2, start = start2,
                        Dhatt = DVLS_C, nperm = 200, seed = 1,
                        verbose = 0)
cat("Testing of all RE (bi1, bi2) - p-value:", HC$pvalue, "| Decision:", HC$Decision, "\n")
#> Testing of all RE (bi1, bi2) - p-value: < 0.001 | Decision: Reject H0

# ## Test for the need of bi1 given the presence of bi2 in the reduced model
HC1 <- Dhypothesis_test(d, Expr2, group = "Subject",
                        random = random2, start = start2,bi_out = "bi1",
                        Dhatt = DVLS_C, nperm = 200, seed = 1,
                        verbose = 0)
cat("Testing bi1 give bi2 in the model- p-value:", HC1$pvalue, "| Decision:", HC1$Decision, "\n")
#> Testing bi1 give bi2 in the model- p-value: < 0.001 | Decision: Reject H0

Skill Acquisition Data

The SkillAcq dataset, available in the TestREnlme Package, contains log-transformed response latencies for \(N = 204\) subjects across \(J = 11\) trial blocks, with a subject-level working memory (WM) covariate wm2. The nonlinear learning model is: \[ Y_{ij} = a_{i1} - (a_{i1} + a_{i0})\exp(a_{i2}\,T_{ij}) + \varepsilon_{ij} \] where \(a_{i0}\) is the initial offset, \(a_{i1}\) the asymptote, and \(a_{i2}\) the learning rate. Each parameter has a WM regression component: \(a_{ik} = \beta_{0k} + \beta_{1k} WM2_i + b_{ik}\).

## Reshape from wide to long format
library(reshape2)
data("SkillAcq")
qrt1 <- melt(SkillAcq, id.vars = c("id", "wm2"),
             variable.name = "ly", value.name = "Y")
qrt1$t <- as.numeric(sub("ly", "", qrt1$ly))


#model formulation
Expr_learn  <- Y ~ ai1 - (ai1 + ai0) * exp(ai2 * t)
random_learn <- c("ai0 ~ B0 + B01 * wm2 + bi0",
                  "ai1 ~ B1 + B11 * wm2 + bi1",
                  "ai2 ~ B2 + B21 * wm2 + bi2")

We first fit the full with MM to assess subject-specific convergence, and the starting values are automatically compute.

## MM estimator 
DMM_all <- Dmethod(data=qrt1, Expr=Expr_learn, group = "id",
                   random = random_learn, method = "MM")
summary(DMM_all)
#> 
#> TestREnlme: Variance component estimates
#> Method : MM 
#> 
#> Fixed effects (Beta):
#>        B0        B1        B2       B11       B01       B21 
#> -2.567417  2.054550 -0.255008 -0.246641  0.378805  0.104123 
#> 
#> Error variance (Sigma2):
#>   0.005694 
#> 
#> Variance-covariance matrix (Dhat):
#>           bi0       bi1       bi2
#> bi0  0.104305 -0.051163 -0.007414
#> bi1 -0.051163  0.038786 -0.006110
#> bi2 -0.007414 -0.006110  0.007466
#> 
#> Subject summary:
#>   Total           : 204 
#>   Used            : 143 
#>   High-kappa excl.: 24 
#>   Non-converged   : 37

Full Model (Three Random Effects)

The full model excludes 61 subjects (24 high-kappa + 37 non-converged); this is also the case for the MMF approach. The figure below, produced using plot_condition(), displays the condition-number diagnostic, which illustrates this practical situation where the default method VLS should be used (as it uses all 204 subjects) and not MM nor MMF. The subsequent analysis should therefore proceed using VLS following Steps 1 to 4 as shown in the theophylline example above.


## Condition number plot -- subjects above threshold excluded
set.seed(123)
NS <- sort(sample(1:204, 20))
plot_condition(DMM_all, kappa_max = 1e4)
Skill acquisition data: the per-subject condition numbers from MM for the full three-random-effects model, sorted in decreasing order for converged subjects. The subjects above the dashed line (red dots) were excluded from MM second-stage estimation, and the subjects below (blue dots) were retained and used in the analysis.
Skill acquisition data: the per-subject condition numbers from MM for the full three-random-effects model, sorted in decreasing order for converged subjects. The subjects above the dashed line (red dots) were excluded from MM second-stage estimation, and the subjects below (blue dots) were retained and used in the analysis.

Summary of Key Functions

Function Purpose
Dmethod() Estimate \(\hat{D}_*\) and \(\hat{\sigma}^2\) (VLS, MM, MMF)
Tstat() Compute observed test statistic \(T_{\rm obs}\)
Dhypothesis_test() Run permutation test; returns \(p\)-value and histogram
MM_base() Pre-compute shared first-stage for MM and MMF
bootstrap_se() Bootstrap standard errors for all estimates
plot_profiles() Individual profile plot with mean
plot_fitted() Fitted curves overlaid on observed data
plot_residuals() Residual diagnostic plots
plot_condition() Condition number diagnostic (MM/MMF)
plot.Dmethod() S3 plotting method for objects: produces profiles, fitted values, residuals, and condition numbers
plot_variance() Compare variance estimates across methods
plot_perm_hist() Permutation null distribution histogram

References

Uwimpuhwe, G., Drikvandi, R., and Blozis, S.A. (2026). Testing random effects in nonlinear mixed-effects models. Statistics in Medicine. https://doi.org/10.1002/sim.70605

Uwimpuhwe, G., Drikvandi, R. and Blozis, S. A. (in preparation). TestREnlme: An R Package for Testing Random Effects in Nonlinear Mixed-Effects Models.

Demidenko, E. (2013). Mixed Models: Theory and Applications with R. Wiley.

Drikvandi, R., Verbeke, G., Khodadadi, A. and Nia, V. P. (2013). Testing multiple variance components in linear mixed-effects models. Biostatistics, 14(1), 144–159.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.