CODEWITHSUNDEEP
rmixed-modelssingular
kris
kris

Reputation: 153

Contrasting results about singularity in random intercept logistic model

I estimated a random intercept logistic regression model in R using the lme4::glmer function. The model converges without issues. However, when producing post-estimation goodness of fit statistics with performance package, this is what I obtain:

> library(performance)
> check_singularity(model_base)
[1] FALSE
> r2(model_base)
boundary (singular) fit: see help('isSingular')
# R2 for Mixed Models

  Conditional R2: 0.777
  Marginal R2: 0.687

So while check_singularity does not detect issues, the r2 command produces the boundary (singular) fit: see help('isSingular') warning. Can you help me understand why this might be the case?

Upvotes: 1

Views: 43

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226712

You haven't given us a reproducible example, but I can confirm that r2() (specifically, r2_nakagawa(), which calls insight::get_variance(), which calls insight:::get_variance.merMod with a null model) does try to fit the null model (only an intercept in the fixed effects). If the full model wasn't singular and the null model is, then I'd expect to see the message you got.

I wrote some simulation code to generate cases where the full model is non-singular and the null model is singular, and found some, but can't replicate your message ...

One thing that's a bit confusing is that check_singularity(), r2_nakagawa(), and lme4::isSingular() (which does lme4's internal checks) may all use different tolerances for checking ... (all three of these methods detect singularity within a tolerance, not exact [floating-point 0] singularity)

library(lme4)
library(performance)

sfun0 <- function(x, theta = 0.5) {
    cat(".")
    dd <- data.frame(x = rnorm(100), g = factor(rep(1:10, each = 10)))
    dd$y <- simulate( ~ x + (1|g),
                     newdata = dd,
                     family = gaussian,
                     newparams = list(beta = c(0, 1), theta = theta, sigma = 1))[[1]]
    m1 <- lmer(y ~ x + (1|g), data = dd)
    return(m1)
}

check_sing <- function(m1) {
    m0 <- update(m1, . ~ . - x)
    sapply(list(m0, m1), check_singularity, tolerance = 1e-4)
}

nsim <- 100
res <- matrix(NA, ncol = 2, nrow = nsim,
              dimnames = list(NULL, c("m0_sing", "m1_sing")))
for (i in 1:nsim) {
    set.seed(100+i)
    res[i, ] <- check_sing(sfun0())
}

with(as.data.frame(res), table(m0_sing, m1_sing))
which(res[, "m0_sing"] & !res[, "m1_sing"])
which(!res[, "m0_sing"] & res[, "m1_sing"])

set.seed(183)
## match default check_singularity() tolerance
performance::r2_nakagawa(sfun0(), tolerance = 1e-5)

set.seed(124)
performance::r2(sfun0())

Upvotes: 4

Related Questions