Estimating a (purposely) misspecified multilevel model in R using frequentist statistics with MCMC/BS and getting cluster-specific effects and CIs

Question

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

Answer 1

I did some experimentation. Looks like the bias is small (although you have to decide whether you think it is; we can't use relative bias as a metric since the true cluster effects are zero in this case).

This assumes that by "cluster effect" you mean the deviation of the clusters from the population-level estimate; things are different and slightly harder if you want to estimate the effect for each cluster (i.e. the sum of the population-level, fixed, effect and the cluster's deviation from the population mean). In this case coef() will give you the sum, but the standard error of this value isn't straightforward to get. (You can get the predicted value in the treatment group for each cluster with its standard error, but the standard deviation of the estimated control vs treatment difference in each cluster is harder ...)

To test this, I ran your simulation code many times and computed summary statistics. The results:

  level     bias  rmse coverage
  <chr>    <dbl> <dbl>    <dbl>
1 1     -0.0214  0.948    0.497
2 2      0.00767 0.954    0.504
3 3      0.00614 0.959    0.47 
4 4     -0.0145  1.01     0.468
5 5      0.0221  1.07     0.458

As stated, the bias is fairly small but the coverage is way too low (should be around 0.95 for an interval of ±2 SE). I don't know if I made a mistake or if this has to do with singular fits ...

Set up a function to simulate data, run the model, and store the estimates, true values, and SEs (true values are all zero in this example, but this is useful if we want to generalize):

library(lme4)
set.seed(101)
fun2 <- function() {
    sim <- simfun()    ## see simfun def below: modified version of OP's code
    mod <- suppressMessages(
        lme4::lmer(ment_health~MT+(1+MT|clinic), data=sim$dat))
    ans <- subset(broom.mixed::tidy(mod, effects = "ran_vals"),
                  term=="MT1", select =c(level, estimate, std.error))
    cbind(ans, true = sim$true_vals)
}

I used tidyverse out of laziness, but this could be done perfectly well in base R as well ...

library(tidyverse)
## map_dfr function needs to take an argument (which we ignore ...)
sims <- map_dfr(1:1000, \(i) fun2(), .id = "sim")
(sims
    |> mutate(in_ci = between(estimate, true-2*std.error, true+2*std.error),
              error = estimate - true)
    |> group_by(level)
    |> summarise(bias = mean(error),
                 rmse = sqrt(mean(error^2)),
                 coverage = mean(in_ci))
)

OP's code with simulation parameters as default arguments; returns a list with both the simulated data and the true cluster-level effects.

simfun <- function(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)
    list(dat = dat, true_vals = re[,'MTeffect'])
}