Dear Stackoverflow friends, I have a challenging task. I am trying to purposely (for research/teaching) estimate a misspecified multilevel model and retrieve its cluster-specific estimates and CIs without using Bayesian statistics.

Here is the model:

mod <- lme4::lmer(ment_health~MT+(1+MT|clinic), data=dat)

The reason the model is misspecified is that the random effect of treatment (MT) doesn’t truly exist in the population. The data is generated by:

intercept <- 10 # Grand mean
MTeffect <- 0 # Overall, fixed effect
sd_MT <-  0 # (No) variability in cluster-specific coefficients
sigma <-  3  # within-cluster variability 
cor <- 0 # No association between slopes and intercepts
 sd_int <-  5 #  variability in cluster-specific intercept
n <-  10 # sample size in each cluster
n_clusters <-  5 # number of clusters
N <- n*n_clusters # Total sample size
clinic <- rep(1:n_clusters, each = n) # Defining cluster ID
MT <- rep(c(0, 1), length.out = N) 
      # Setting dichotomous predictor: control vs. treat
varmat <- matrix(c(sd_int^2, cor, cor, sd_MT^2), 2, 2) 
           # Variance-covariance matrix
 re <- mvtnorm::rmvnorm(n_clusters, sigma = varmat) 
# Generate cluster-level population information from var-covar
colnames(re) <- c('Intercept', 'MTeffect')
ment_health <- (intercept + re[clinic, 'Intercept'])   +
                (MTeffect  + re[clinic, 'MTeffect'])*MT +  
# basically, this line is redundant bc it's all 0 
# in the null condition
                 rnorm(N, sd = sigma)
  # Putting everything together
  dat <- data.frame(ment_health, MT = factor(MT), clinic, ppt_ID)

Sometimes the model is estimated relatively well with fixed MT effect, sd, and corr close to 0. But, more frequently I get the following warning:

boundary (singular) fit: see help('isSingular')

where the Corr is estimated at 1 or -1.

I want to use bootstrapping or MCMC to solve this. I also want to get the clinic (cluster) specific MT effects and their CIs. I usually do this with the package mixedup, but I am not sure how to retrieve these particular details using methods other than lme4 for estimation.

I know it's a lot, but any help would be tremendously appreciated!

@BenBolker,

Thank you so much for your quick and helpful response! I wanted to follow up with some more detail that didn’t fit in the comments.

I’m hoping to obtain stable estimates of cluster-specific effects, even when the model is knowingly misspecified (in this case, due to the absence of a true random effect for MT).

I understand that the singular fit warning indicates issues with the model specification, particularly when the random effect variance is estimated at or near zero, but I wonder if I can still trust the cluster-specific coefficients (I am not directly interested in the parameter of random effects correlation).

To be even more specific, aside from the point estimates of each clinic (i.e., cluster), I am trying to run an informal significance test on each cluster to test whether its coefficient is statistically different from zero (evaluating Type I errors). This is why I was trying to retrieve the confidence intervals. Thank you for pointing out that these are actually the quantiles of the conditional distributions.

Do you think it is reasonable to use the quantiles of the conditional distributions as indicators of “significance” (i.e., whether they include 0)?

Thank you again for your insights, I truly appreciate your advice!

Udi