Package {TestREnlme}

Pooled within-cluster variance for MM methods

Description

Computes the pooled within-cluster residual variance estimate \hat\sigma^2 = \text{SSE} / \sum_i (n_i - k), where k is the number of parameters in start.

Usage

.Sigma2hat_MM(data, COEFs, Expr, group)

Arguments

Value

A scalar estimate of \sigma^2.

Compute starting values automatically via multi-start NLS

Description

Internal helper used when start is not supplied by the user. Identifies model parameters in Expr (symbols not present in data), and searches for starting values using nls.multstart::nls_multstart(), which tries multiple starting points within \pmrange for each parameter and refits an nls() model on the pooled data (ignoring random effects and grouping structure) from each, keeping the best-converging fit.

Usage

.auto_start(data, Expr, range = 10, seed = 123)

Arguments

Value

A named numeric vector of starting values, suitable for use as the start argument of Dmethod and related functions.

Estimate fixed effects by nonlinear least squares

Description

Estimates the fixed-effects parameter vector by unweighted NLS (first attempt via nls), or by weighted NLS via nlminb when weights is supplied or when the unweighted NLS fails. Used both for initial unweighted estimation and for GLS re-estimation at the end of VLS using the estimated variance components.

Usage

Beta_hat(data, Expr, start, weights = NULL, group = NULL, verbose = 1)

Arguments

Value

A named numeric vector of estimated fixed effects.

Examples

Expr   <- conc ~ Dose * exp(ai2 + ai3 - ai1) *
            (exp(-Time * exp(ai3)) - exp(-Time * exp(ai2))) /
            (exp(ai2) - exp(ai3))
start  <- c(ai1 = -3.22, ai2 = 0.47, ai3 = -2.45)
Beta_hat(as.data.frame(Theoph), Expr, start)

Permutation test for random effects in linear and nonlinear mixed-effects models

Description

Performs a nonparametric permutation test for all random effects or any user-specified subset, using the test statistic of Drikvandi et al. (2013) adapted for nonlinear mixed-effects models in Uwimpuhwe (2026).

Usage

Dhypothesis_test(
  data,
  Expr,
  group,
  random,
  start = NULL,
  bi_out = NULL,
  method = c("VLS", "MM", "MMF"),
  Dhatt = NULL,
  Thatt = NULL,
  MM_base_obj = NULL,
  nperm = 200,
  seed = NULL,
  sig_alpha = 0.05,
  kappa_max = 10000,
  RR_catof = "kappa",
  verbose = 1,
  perm_freq = 10
)

Arguments

Value

An object of class "Dtest", a list with components:

Decision

Character, "Reject H0" or "Do not reject H0".

pvalue

The empirical permutation p-value.

Tobs

The observed test statistic T_\text{obs}.

Tperm

Numeric vector of length nperm containing the permutation statistics T^{(1)}, \ldots, T^{(B)}.

Dhatt

The Dmethod object used.

bi_out

The random effects tested.

plot

A ggplot2 histogram of the permutation null distribution with T_\text{obs} annotated.

References

Uwimpuhwe, G., Drikvandi, R., and Blozis, S.A. (2026). Testing random effects in nonlinear mixed-effects models. Statistics in Medicine. doi: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. Journal of Statistical Software.

Demidenko, E. (2013). Mixed Models: Theory and Applications with R (2nd ed.). 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.

Examples

d      <- as.data.frame(Theoph)
Expr   <- conc ~ Dose * exp(ai2 + ai3 - ai1) *
            (exp(-Time * exp(ai3)) - exp(-Time * exp(ai2))) /
            (exp(ai2) - exp(ai3))
start  <- c(ai1 = -3.22, ai2 = 0.47, ai3 = -2.45)
random <- c("ai1 ~ B1 + bi1", "ai2 ~ B2 + bi2", "ai3 ~ B3 + bi3")
DVLS  <- Dmethod(d, Expr, group = "Subject",
                   random = random, start = start)


## Permutation test (slow -- use nperm = 200 minimum for real use)
H      <- Dhypothesis_test(d, Expr, group = "Subject",
                           random = random, start = start,
                           Dhatt = DVLS, nperm = 20, seed = 1)
H$pvalue
H$plot

Non-parametric approach to estimate variance components in a linear and nonlinear mixed-effects model

Description

Computes the scaled variance covariance matrix \hat D_* and the error variance \hat\sigma^2 using one of three nonparametric estimators: Variance Least Squares ("VLS"), Method of Moments ("MM"), or Method of Moments with First-Order Approximation ("MMF"). The method also estimate weighted fixed effects

Usage

Dmethod(
  data,
  Expr,
  group,
  random,
  start = NULL,
  method = c("VLS", "MM", "MMF"),
  MM_base_obj = NULL,
  kappa_max = 10000,
  RR_catof = "kappa",
  Beta_nls = NULL,
  verbose = 1,
  is_permuting = FALSE
)

Arguments

Value

An object of class "Dmethod", a list with components:

Dhat

The estimated scaled covariance matrix of random effects \hat D_*.

Sigma2

The estimated error variance \hat\sigma^2.

Beta

The estimated fixed-effects vector \hat\beta.

method

The estimation method used.

R

Per-subject Jacobian matrices R_i(\hat\theta_0).

Ypred

Per-subject fitted values under fixed effects only.

residuals

Per-subject residuals \hat e_i.

IDconvegence

Subject inclusion/exclusion summary (MM/MMF only; NULL otherwise).

Additional internal objects are stored as an "internal" attribute (attr(object, "internal")) and used by package methods. They are not part of the user interface.

References

Uwimpuhwe, G., Drikvandi, R., and Blozis, S.A. (2026). Testing random effects in nonlinear mixed-effects models. Statistics in Medicine. doi: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. Journal of Statistical Software.

Demidenko, E. (2013). Mixed Models: Theory and Applications with R (2nd ed.). John Wiley & Sons

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.

Examples

d      <- as.data.frame(Theoph)
Expr   <- conc ~ Dose * exp(ai2 + ai3 - ai1) *
            (exp(-Time * exp(ai3)) - exp(-Time * exp(ai2))) /
            (exp(ai2) - exp(ai3))
start  <- c(ai1 = -3.22, ai2 = 0.47, ai3 = -2.45)
random <- c("ai1 ~ B1 + bi1", "ai2 ~ B2 + bi2", "ai3 ~ B3 + bi3")
DVLS   <- Dmethod(d, Expr, group = "Subject",
                  random = random, start = start)
## method defaults to VLS; explicit: method = "VLS"
DVLS[c("Dhat", "Sigma2", "Beta")]

Build model expressions for MM/MMF estimation

Description

Constructs the set of model expressions needed by MM_base: the full expression with random effects substituted, the second-stage expression incorporating subject-level covariates from Q, and the nlsList-style formula for per-subject NLS fitting.

Usage

Expressions(data, group, random, Expr, start, seed = 123)

Arguments

Value

A list with components Expr, Expr1, Expr2, Expr_MM_all0, random0, start, and start_MM_all.

Compute first-stage quantities for MM and MMF estimators

Description

Performs the first-stage nonlinear least squares (NLS) fit for each subject individually, and computes the per-subject Jacobian matrices R_i^T R_i and condition numbers \kappa_i. The returned object can be passed to Dmethod via the MM_base_obj argument to avoid recomputing the first stage when calling both method = "MM" and method = "MMF" on the same data.

Usage

MM_base(
  data,
  Expr,
  group,
  random,
  start = NULL,
  kappa_max = 10000,
  RR_catof = "kappa",
  verbose = 1
)

Arguments

Value

A list of class "MM_base" with components:

data

The (re-ordered, group-renumbered) data.frame used internally.

group

Name of the grouping variable.

start

The starting values used.

Beta_nls

Named numeric vector of pooled (first-pass) fixed-effects estimates from Beta_hat.

random

The random formulas, aligned with start.

re_names

Named character vector mapping random-effect names to subject-specific parameter names.

Tmatrix

Matrix of per-subject scaled covariance terms (one row per subject, k^2 columns).

Ti

List of per-subject k \times k covariance matrices derived from Tmatrix.

ddi

List of per-subject k \times k second-stage deviation outer-product matrices (eq.~8.20).

Lamda

Minimum eigenvalue of the standardised scatter matrix, used to bias-correct the MM/MMF variance estimates.

sigma2

Pooled residual variance across retained subjects.

Beta_GB

Named numeric vector of non-random-effect fixed-effects estimates.

m

Total number of second-stage fixed-effects parameters (k * ncol(Q)).

k

Number of random-effects parameters.

id, uid, nud

Subject index vector, unique subject IDs, and number of unique subjects.

IDworks

Integer vector of subject IDs retained after filtering.

Q

Subject-level design matrix.

Aai

List with elements Ai (per-subject design matrices) and a0i (per-subject parameter estimates).

kappa

Named numeric vector of per-subject condition numbers.

IDconvegence

List with elements ID_all, ID_used, and ID_supressed (itself a list with highRR and Unconverged).

Examples

d      <- as.data.frame(Theoph)
Expr   <- conc ~ Dose * exp(ai2 + ai3 - ai1) *
            (exp(-Time * exp(ai3)) - exp(-Time * exp(ai2))) /
            (exp(ai2) - exp(ai3))
start  <- c(ai1 = -3.22, ai2 = 0.47, ai3 = -2.45)
random <- c("ai1 ~ B1 + bi1", "ai2 ~ B2 + bi2", "ai3 ~ B3 + bi3")
mb <- MM_base(d, Expr, group = "Subject", random = random, start = start)

Skill acquisition data

Description

Response latency data from a study of quantitative skill acquisition on a learning task (Blozis 2004). Log-transformed response latencies are recorded for N = 204 subjects across J = 11 trial blocks, stored in wide format. Used to illustrate fitting a nonlinear mixed-effects model with a subject-level covariate (working memory) incorporated through the second-stage model, and to demonstrate reshaping wide-format longitudinal data into the long format required by Dmethod and Dhypothesis_test.

Usage

SkillAcq

Format

A data frame with 204 rows (one per subject) and 13 columns:

id

Subject identifier.

ly1

Log-transformed response latency at trial block 1.

ly2

Log-transformed response latency at trial block 2.

ly3

Log-transformed response latency at trial block 3.

ly4

Log-transformed response latency at trial block 4.

ly5

Log-transformed response latency at trial block 5.

ly6

Log-transformed response latency at trial block 6.

ly7

Log-transformed response latency at trial block 7.

ly8

Log-transformed response latency at trial block 8.

ly9

Log-transformed response latency at trial block 9.

ly10

Log-transformed response latency at trial block 10.

ly11

Log-transformed response latency at trial block 11.

wm2

Subject-level working-memory covariate.

Details

SkillAcq is stored in wide format, as commonly encountered in practice. Before use with Dmethod it must be reshaped to long format, with one row per subject-trial observation; see Examples below.

The nonlinear mixed-effects model used for this dataset is

Y_{ij} = a_{i1} - (a_{i1} + a_{i0}) \exp(a_{i2} T_{ij}) + \varepsilon_{ij},

with

a_{ik} = \beta_{0k} + \beta_{1k} wm2_i + b_{ik}, \quad k \in \{0, 1, 2\},

where a_{i0}, a_{i1}, and a_{i2} represent, respectively, the subject-specific initial performance offset, lower asymptote, and learning rate, each with a regression component on wm2 and a subject-specific random effect. The questions of interest are whether all three random effects are necessary, and whether one or more can be removed to obtain a more parsimonious model; see Dhypothesis_test.

Source

Blozis, S. A. (2004). Structured latent curve models for the study of change in multivariate repeated measures. Psychological Methods, 9(3), 334–353. https://doi.org/10.1037/1082-989X.9.3.334

Used as a worked example in Uwimpuhwe, G., Drikvandi, R. and Blozis, S. A. (in preparation). TestREnlme: An R Package for Testing Random Effects in Nonlinear Mixed-Effects Models. Journal of Statistical Software.

Examples

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

## Model: Y_ij = ai1 - (ai1 + ai0) * exp(ai2 * t_ij) + e_ij
## with aik = B0k + B1k * wm2_i + bik, k in {0, 1, 2}
Expr_learn <- Y ~ ai1 - (ai1 + ai0) * exp(ai2 * t)
random_learn <- c("ai0 ~ B00 + B10 * wm2 + bi0",
                  "ai1 ~ B01 + B11 * wm2 + bi1",
                  "ai2 ~ B02 + B12 * wm2 + bi2")


## Estimate variance components (VLS) and test all three random effects
DVLS_learn <- Dmethod(qrt1, Expr_learn, group = "id",
                      random = random_learn, start = NULL)
DVLS_learn[c("Dhat", "Sigma2", "Beta")]

H_learn <- Dhypothesis_test(qrt1, Expr_learn, group = "id",
                            random = random_learn, start = NULL,
                            Dhatt = DVLS_learn, nperm = 200, seed = 1)
H_learn$pvalue

Compute the observed test statistic

Description

Computes T_\text{obs} = N^{-1} \sum_i \text{tr}(R_i \hat D_* R_i^T) from a fitted Dmethod object.

Usage

Tstat(Dobj, bi_out = NULL)

Arguments

Value

A scalar, the value of T_\text{obs}.

Compute Jacobian matrix and fixed-effects predicted values

Description

Evaluates the Jacobian R_i(\hat\theta_0) of the nonlinear model function with respect to the random-effects parameters, and computes fitted values under the fixed-effects-only model.

Usage

ZandYPred(data, Expr, group, random, Beta_nls = NULL, start = NULL)

Arguments

Value

A list with components:

R

Named list of per-subject Jacobian matrices R_i(\hat\theta_0).

Ypred

Named list of per-subject predicted vectors f_i(A_i\hat\beta, \hat\gamma).

residuals

Named list of per-subject residual vectors \hat e_i = Y_i - f_i.

Bootstrap standard errors for variance component estimates

Description

Computes bootstrap standard errors for the fixed effects \hat\beta, random effect variance components \hat D_*, and error variance \hat\sigma^2 using one of two bootstrap strategies documented in the literature for mixed-effects models.

Usage

bootstrap_se(
  Dobj,
  nboot = 200,
  type = c("case", "residual"),
  seed = NULL,
  verbose = 1
)

Arguments

Value

A list of class "Dboot" with components:

Beta

Matrix with columns Estimate and SE for the fixed effects.

Dhat

Matrix with columns Estimate and SE for each variance/covariance component (flattened via row:column naming).

Sigma2

Matrix with columns Estimate and SE for the error variance.

Boot_Beta

Matrix of bootstrap fixed-effects estimates (nboot x length(Beta)).

Boot_Dhat

Array of bootstrap variance component matrices (k x k x nboot).

Boot_Sigma2

Numeric vector of bootstrap error variance estimates (length nboot).

type

The bootstrap strategy used.

nboot

Number of bootstrap samples requested.

nfail

Number of bootstrap samples that failed.

References

Carpenter, J.R., Goldstein, H. and Rasbash, J. (2003). A novel bootstrap procedure for assessing the relationship between class size and achievement. Journal of the Royal Statistical Society C, 52, 431–443.

Thai, H.T., Mentre, F., Holford, N.H.G., Veyrat-Follet, C. and Comets, E. (2013). A comparison of bootstrap approaches for estimating uncertainty of parameters in linear mixed-effects models. Pharmaceutical Statistics, 12, 129–140.

Examples

d      <- as.data.frame(Theoph)
Expr   <- conc ~ Dose * exp(ai2 + ai3 - ai1) *
            (exp(-Time * exp(ai3)) - exp(-Time * exp(ai2))) /
            (exp(ai2) - exp(ai3))
start  <- c(ai1 = -3.22, ai2 = 0.47, ai3 = -2.45)
random <- c("ai1 ~ B1 + bi1", "ai2 ~ B2 + bi2", "ai3 ~ B3 + bi3")
DVLS   <- Dmethod(d, Expr, group = "Subject",
                  random = random, start = start, method = "VLS")
DVLS[c("Dhat", "Sigma2", "Beta")]

BootSE <- bootstrap_se(DVLS, nboot = 20, type = "case", seed = NULL, verbose = 1)

Plot method for Dmethod objects

Description

By default all relevant diagnostic plots are produced and arranged in a grid. Use the which argument to select specific plots.

Usage

## S3 method for class 'Dmethod'
plot(x, time, which = NULL, ...)

Arguments

Value

Invisibly returns a list of ggplot2 objects.

Plot method for Dtest objects

Description

Displays the permutation null distribution histogram with the observed test statistic annotated.

Usage

## S3 method for class 'Dtest'
plot(x, ...)

Arguments

Value

Invisibly returns the ggplot2 object.

Plot per-subject condition numbers

Description

Displays the condition numbers \kappa_i from a "MM_base" object or the MM_base_obj slot of a "Dmethod" object, sorted in decreasing order. A horizontal reference line at kappa_max marks the exclusion threshold.

Usage

plot_condition(Dobj, kappa_max = 10000, log_scale = TRUE)

Arguments

Value

A ggplot2 object.

Plot fitted curves overlaid on observed data

Description

Draws the model-predicted trajectories on top of the observed data for each subject. When a "Dmethod" object is supplied, subject-specific EBLUP-adjusted predictions are used; otherwise only the population fixed-effects curve is shown.

Usage

plot_fitted(Dobj, time, subjects = NULL, overlay = TRUE, ncol = 4)

Arguments

Value

A ggplot2 object.

Histogram of the permutation null distribution

Description

Regenerates or customises the permutation histogram stored in a "Dtest" object. The observed statistic T_\text{obs} is shown as a vertical dashed line, and the empirical p-value and rejection decision are annotated on the plot.

Usage

plot_perm_hist(
  Htest,
  bins = 30,
  fill = "steelblue",
  line_col = "red",
  title = NULL
)

Arguments

Value

A ggplot2 object.

Plot raw individual trajectories

Description

Produces a spaghetti plot of the raw observed trajectories for all subjects (or a selected subset), with an optional group-level mean profile and subject highlighting.

Usage

plot_profiles(
  data,
  group,
  time,
  response,
  subjects = NULL,
  mean_profile = TRUE,
  highlight = NULL,
  title = "Individual profiles"
)

Arguments

Value

A ggplot2 object.

Examples

d <- as.data.frame(Theoph)
plot_profiles(d, group = "Subject", time = "Time", response = "conc")

Plot residuals versus fitted values

Description

Produces a residuals-versus-fitted-values scatter plot with a horizontal reference line at zero and a loess smoother to aid detection of systematic misfit or heteroscedasticity.

Usage

plot_residuals(Dobj, time, type = c("response", "standardised", "subject"))

Arguments

Value

A ggplot2 object.

Compare variance component estimates across methods

Description

Produces a side-by-side dot plot (or bar chart) of the estimated variance components \hat d_{*jj} for each random effect, comparing results from multiple estimation methods.

Usage

plot_variance(
  Dhatt_list,
  component = c("diagonal", "full"),
  title = "Variance component estimates"
)

Arguments

Value

A ggplot2 object.