easystats / easystats/performance

`check_singularity.lme()`: errors on models with more than one level of nesting; misplaced parenthesis

Open
#937 0 comments 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

Dominant language
R
Stars
1.2k
Forks
109
Avg merge
6h 34m
Merged PRs (30d)
8

Description

Two small, independent things in check_singularity.lme(). performance 0.17.1
(= CRAN, = main @ 50552a0), nlme 3.1-170, R 4.6.1.

1. Errors instead of answering, for ~ 1 | A/B

nlme::getVarCov.lme() refuses more than one grouping level
(if (length(obj$groups) > 1) stop("not implemented for multiple levels of nesting")),
so check_singularity() errors out rather than returning TRUE/FALSE:

data(Pixel, package = "nlme")
fit <- nlme::lme(pixel ~ day, random = ~ 1 | Dog / Side, data = Pixel)

nlme::VarCorr(fit)          # works fine
#> Dog =       pdLogChol(1)
#> (Intercept) 647.2931   25.44196
#> Side =      pdLogChol(1)
#> (Intercept) 218.3420   14.77640
#> Residual    233.3463   15.27568

performance::check_singularity(fit)
#> Error in getVarCov.lme(x) : not implemented for multiple levels of nesting

Single-level lme models are fine, and the lmer equivalent
(pixel ~ day + (1 | Dog/Side)) returns FALSE normally.

fit$modelStruct$reStruct looks like the natural replacement: it returns one
block per nesting level, and the blocks are already relative to sigma^2,
so they are dimensionless:

b <- as.matrix(fit$modelStruct$reStruct)
b
#> $Dog
#>             (Intercept)
#> (Intercept)     2.77396
#>
#> $Side
#>             (Intercept)
#> (Intercept)   0.9356995

lapply(b, function(v) as.matrix(v) * fit$sigma^2)   # reproduces VarCorr exactly
#> $Dog  647.2931
#> $Side 218.342

Since this errors rather than misreports, any caller inherits the error unless
it shields itself: insight:::.compute_variances() does, wrapping the call in
.safe(); a direct caller does not.

2. abs() wraps the comparison rather than the variances

L228 of R/check_singularity.R, identical in CRAN 0.17.1 and at HEAD:

any(abs(stats::na.omit(as.numeric(diag(nlme::getVarCov(x)))) < tolerance))

The closing parenthesis of abs() sits after < tolerance, so abs() receives
the logical vector, not the variances:

v <- c(5.4150876, 0.0512695)
v < 1e-5              #> FALSE FALSE
abs(v < 1e-5)         #> 0 0            <- abs() of a logical
any(abs(v < 1e-5))    #> FALSE
any(abs(v) < 1e-5)    #> FALSE          <- presumably intended

It is essentially harmless — abs() of a logical is a 0/1 numeric and any()
coerces back — and the two forms can only differ for a negative variance, which
getVarCov() does not return. But the intent was clearly abs(diag(...)) < tolerance,
as in the sibling methods on L206, L214 and L220. (bbolker made a related remark
about a redundant abs() in the old .merMod line in easystats/insight#878.)

The whole fix is moving that parenthesis: any(abs(...) < tolerance).

Contributor guide

Open the contributing guide

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

Research direction

Start in R/check_singularity.R at the lme handling around L228, and compare the variance extraction with the sibling methods around L206, L214, and L220. Reproduce the nested Dog/Side example with nlme, then verify that check_singularity() returns TRUE or FALSE instead of erroring and that the tolerance comparison applies abs() to the variances.

Written by the indexing model from the issue text.

Assessment

Tech stack
r
Domain
tooling
Issue type
Bug
Difficulty
3/5
Estimated time
1-2 days
Activity status
Active
Clarity
Clearly specified
Newbie friendliness
72/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.