"No class variable" in a multilevel imputation model

[MuzzJ]

Dear Prof.Van Buuren, thank you very much for your book.
I am glad I have a chance to ask you a question here about it as well!

I have the following problem. It happened that I have only one level-2 variable and it has several missing values(it was an online survey from two organizations and there was a trouble with the data collection, so I am am sure the missing data is MCAR). I try to use 2lonly.pmm to fill in the gaps. I specify this variable as a -2 in pred:
suppressWarnings(suppressMessages( ini <- mice(df, maxit = 0) ))
pred <- ini$pred
pred [,"school"] <- -2

When running the mice function it gives me the following error:
Error in check.cluster(data, predictorMatrix) : Convert cluster variable school to integer by as.integer()

I switched the organization variable to an integer and the error message then became as the following:
Error in .imputation.level2(y = y, ry = ry, x = x, type = type, wy = wy, : No class variable

I added additional variable "city", where both organizations are at, and specified it as a cluster variable in the pred. It didn't help, same error.

Then I tried to specify the full pred matrix in excel and read it from a csv file as a matrix, the error changed to:
Error in predictorMatrix[h, ] : subgroup out of bounds

I consulted your book again and saw the following passage:
"It is conceptually straightforward to extend imputations to higher levels (Yucel 2008). If there are two levels, combine all level-2 predictors with an aggregate (e.g., the cluster means) of the level-1 predictors and the level-1 outcomes. Once we have this, we may choose suitable methods from Chapter 3 to impute the missing level-2 variables in the usual way. No new issues arise."

Is there a way to fix my problem? Or, given the above, is there maybe a way to not treat this variable as a cluster variable and just impute it with logreg in mice since it's the only one variable? Do you mean by aggregate to compute mean values or cluster variables?

Thank you in advance for reading the post and have a great day ahead!

Best Regards,
Yulia

[irenerlover]

I have the same issue! Suspect it's a bug!

[thomvolker]

Multilevel imputation can indeed be tricky, but there is no way of helping you (or determining whether it is a bug) without a reprex. If you are unsure how to write a reprex, please see here: https://reprex.tidyverse.org/articles/reprex-dos-and-donts.html

[irenerlover]

Hi Thom,

I have created some synthetic data for the error that I had:

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(mice)
#> 
#> Attaching package: 'mice'
#> The following object is masked from 'package:stats':
#> 
#>     filter
#> The following objects are masked from 'package:base':
#> 
#>     cbind, rbind

set.seed(20260106)

#----------------------------
# 1) DGP: 2-level long data
#----------------------------
n_id    <- 180
times   <- c(-1, 0, 3, 6, 12, 13)   
n_t     <- length(times)

# subject-level
id_tbl <- tibble(
  subj        = 1:n_id,
  site        = factor(sample(c("A","B","C","D"), n_id, replace = TRUE,
                              prob = c(0.45, 0.25, 0.20, 0.10))),
  stage_grp   = factor(sample(1:4, n_id, replace = TRUE, prob = c(0.35,0.30,0.20,0.15))),
  base_flag1  = factor(rbinom(n_id, 1, 0.35), levels = 0:1, labels = c("No","Yes")),
  base_flag2  = factor(rbinom(n_id, 1, 0.25), levels = 0:1, labels = c("No","Yes")),
  base_age    = round(rnorm(n_id, 52, 14)),
  u0          = rnorm(n_id, 0, 1)  # random intercept
)

# long format
dat <- tidyr::expand_grid(subj = 1:n_id, visit_time = times) %>%
  left_join(id_tbl, by = "subj") %>%
  mutate(
    t01 = (visit_time - min(times)) / (max(times) - min(times)),  # scaled time
    
    # a time-varying latent disease activity affected by stage + subject RE
    act = 0.6*u0 + 0.4*(as.numeric(stage_grp) - 1) + 0.8*t01 + rnorm(n(), 0, 0.7)
  ) %>%
  mutate(
    # Longitudinal continuous vars (to be imputed with 2l.pmm)
    lab_alpha = pmax(0.1, 2.5 + 0.9*act + rnorm(n(), 0, 0.8)),      # positive
    lab_beta  = pmax(1,   55  - 8.0*act + rnorm(n(), 0, 10)),       # like eGFR-ish
    lab_gamma = pmax(5,   35  - 1.8*act + rnorm(n(), 0, 4)),        # like albumin-ish
    lab_delta = pmax(0.05, exp(-0.4 + 0.35*act + rnorm(n(), 0, 0.5))),# like protein-ish
    
    # Baseline-only variable (constant within subj) but stored in long:
    base_pct = pmin(100, pmax(0, 100*plogis(-0.5 + 0.5*u0 + 0.4*(as.numeric(stage_grp)-2) +
                                              rnorm(n(), 0, 0.35))))
  ) %>%
  group_by(subj) %>%
  mutate(
    # enforce baseline-only as constant within subject (take first value)
    base_pct = base_pct[visit_time == min(visit_time)][1]
  ) %>%
  ungroup() %>%
  mutate(
    # Outcomes / auxiliaries (not imputed; just predictors or ignored)
    event_y = rbinom(n(), 1, plogis(-2 + 0.6*act)),
    dur_y   = pmax(0, 12 + 5*act + rnorm(n(), 0, 4))
  )

# introduce missingness (MAR-ish: more missing when act high, and later times)
set.seed(20260106 + 1)
dat <- dat %>%
  mutate(
    pr_miss = plogis(-2 + 0.8*t01 + 0.6*scale(act)[,1])
  ) %>%
  mutate(
    lab_alpha = ifelse(runif(n()) < pr_miss*0.35, NA, lab_alpha),
    lab_beta  = ifelse(runif(n()) < pr_miss*0.30, NA, lab_beta),
    lab_gamma = ifelse(runif(n()) < pr_miss*0.25, NA, lab_gamma),
    lab_delta = ifelse(runif(n()) < pr_miss*0.30, NA, lab_delta),
    base_pct  = ifelse(runif(n()) < 0.18, NA, base_pct),
    base_flag1 = ifelse(runif(n()) < 0.10, NA, as.character(base_flag1)) |> factor(levels = c("No","Yes"))
  )

colnames(dat)[12] = 'P-ET_type'
#----------------------------
# 2) MICE: 2-level setup
#----------------------------

# Predictor matrix
pred <- make.predictorMatrix(dat)

# 2-level structure:
# - cluster variable (subj) must be -2
pred[, "subj"] <- -2

# - do NOT use site as predictor (optional)
pred[, "site"] <- 0

# - do NOT use outcomes as predictors (optional)
pred[, c("event_y", "dur_y")] <- 0

# - baseline-only vars should not be predicted by time directly (optional)
pred[c("base_flag1","base_flag2","base_age","base_pct","stage_grp"), "visit_time"] <- 0

# Methods
meth <- make.method(dat)

# Longitudinal continuous: 2l.pmm
meth[c("lab_alpha","P-ET_type","lab_gamma","lab_delta")] <- "2l.pmm"

# Baseline-only:
meth["base_pct"]   <- "2lonly.norm"   # could be 2lonly.pmm if you prefer
meth["base_flag1"] <- "2lonly.pmm"    # binary factor via pmm works OK often
# (Leave base_flag2/base_age/stage_grp/site/event_y/dur_y as "" = not imputed)
meth[c("base_flag2","base_age","stage_grp","site","event_y","dur_y","visit_time","act","t01","u0")] <- ""

# Post-processing: mimic your group-based constraints + range squeezes
post <- make.post(dat)

# base_pct constraint depending on stage_grp:
# - stage 1 => set to 0
# - stage 2 => squeeze to (0.01, 25)
# - stage 3 => squeeze to (25.01, 75)
# - stage 4 => squeeze to (75.01, 100)
post["base_pct"] <- "
  g <- data[!r[, j], 'stage_grp']
  y <- imp[[j]][, i]

  y[g == '1'] <- 0
  y[g == '2'] <- squeeze(y[g == '2'], c(0.01, 25))
  y[g == '3'] <- squeeze(y[g == '3'], c(25.01, 75))
  y[g == '4'] <- squeeze(y[g == '4'], c(75.01, 100))

  imp[[j]][, i] <- y
"

post["P-ET_type"] <- "
  imp[[j]][, i] <- squeeze(imp[[j]][, i], c(5, 140))
"
post["lab_delta"] <- "
  imp[[j]][, i] <- squeeze(imp[[j]][, i], c(0.05, 12))
"
post["lab_gamma"] <- "
  imp[[j]][, i] <- squeeze(imp[[j]][, i], c(10, 55))
"

# Run imputations
set.seed(3528)
imp <- mice(
  data = dat,
  m = 5,
  maxit = 5,
  predictorMatrix = pred,
  method = meth,
  post = post,
  printFlag = TRUE
)
#> 
#>  iter imp variable
#>   1   1  base_flag1
#> Error in if (sum(type == -2L) != 1L) stop("No class variable"): missing value where TRUE/FALSE needed

imp <- mice(
  data = dat,
  m = 5,
  maxit = 5,
  printFlag = TRUE
)
#> 
#>  iter imp variable
#>   1   1  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   1   2  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   1   3  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   1   4  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   1   5  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   2   1  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   2   2  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   2   3  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   2   4  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   2   5  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   3   1  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   3   2  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   3   3  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   3   4  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   3   5  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   4   1  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   4   2  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   4   3  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   4   4  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   4   5  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   5   1  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   5   2  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   5   3  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   5   4  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#>   5   5  base_flag1  lab_alpha  P-ET_type  lab_gamma  lab_delta  base_pct
#> Warning: Number of logged events: 1

If I used the multilevel imputation (probably related to 2lonly.pmm from traceback), then an error occurred; however, if I imputed the data ignoring the clustering, the mice could run without any error.
Maybe it is also related to the variable name that I provided ('P-ET_type')? But why did it pop up this error instead of Error in terms.formula(formula, data = data) : invalid model formula in ExtractVars?

Thank you.

Best,
Irene

[irenerlover]

One more issue I faced:

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(tidyr)
library(magrittr)
#> 
#> Attaching package: 'magrittr'
#> The following object is masked from 'package:tidyr':
#> 
#>     extract
library(mice)
#> 
#> Attaching package: 'mice'
#> The following object is masked from 'package:stats':
#> 
#>     filter
#> The following objects are masked from 'package:base':
#> 
#>     cbind, rbind
library(miceadds)
#> * miceadds 3.18-36 (2025-09-12 09:54:54)
library(lme4)
#> Loading required package: Matrix
#> 
#> Attaching package: 'Matrix'
#> The following objects are masked from 'package:tidyr':
#> 
#>     expand, pack, unpack
#> 
#> Attaching package: 'lme4'
#> The following object is masked from 'package:mice':
#> 
#>     toenail
#Generate fake data
n_id <- 129
time_pts <- 1:6
clevels <- c("a a", "b b", "c c", "d d")

set.seed(20260105)
# subject-level type (constant across time within subject)
subj_type <- sample(
  clevels,
  size = n_id,
  replace = TRUE,
  prob = c(0.01, 0.15, 0.25, 0.59)   # adjust as you like
)

dat_fake <- tidyr::expand_grid(
  id = 1:n_id,
  time = time_pts
) %>%
  mutate(
    Type = factor(subj_type[id], levels = clevels),
    
    # subject random intercept
    u = rep(rnorm(n_id, mean = 0, sd = 1.0), each = length(time_pts)),
    
    # mean structure over time + type effect
    base = 3.5 + 0.15 * time,
    type_eff = dplyr::case_when(
      Type == "a a" ~  0.2,
      Type == "b b" ~  0.0,
      Type == "c c" ~ -0.1,
      Type == "d d" ~  0.3
    ),
    
    # positive/right-skewed outcome
    lin = base + type_eff + u + rnorm(n(), 0, 0.8),
    measure_a = round(exp(lin / 2), 2),
    
    # missingness (more missing later)
    miss_p = 0.02 + 0.02 * (time - 1),
    measure_a = ifelse(runif(n()) < miss_p, NA_real_, measure_a)
  ) %>%
  select(id, Type, time, measure_a)

str(dat_fake)
#> tibble [774 × 4] (S3: tbl_df/tbl/data.frame)
#>  $ id       : int [1:774] 1 1 1 1 1 1 2 2 2 2 ...
#>  $ Type     : Factor w/ 4 levels "a a","b b","c c",..: 3 3 3 3 3 3 3 3 3 3 ...
#>  $ time     : int [1:774] 1 2 3 4 5 6 1 2 3 4 ...
#>  $ measure_a: num [1:774] 4.13 4.1 6.25 NA 7.77 ...
head(dat_fake, 12)
#> # A tibble: 12 × 4
#>       id Type   time measure_a
#>    <int> <fct> <int>     <dbl>
#>  1     1 c c       1      4.13
#>  2     1 c c       2      4.1 
#>  3     1 c c       3      6.25
#>  4     1 c c       4     NA   
#>  5     1 c c       5      7.77
#>  6     1 c c       6      9.56
#>  7     2 c c       1     22.0 
#>  8     2 c c       2     27.0 
#>  9     2 c c       3     NA   
#> 10     2 c c       4     28.2 
#> 11     2 c c       5     64.4 
#> 12     2 c c       6     37.9

#Set predictor matrix
pred3 <- make.predictorMatrix(dat_fake)
pred3[, c("id")] <- -2
dat_fake$id %<>% as.integer()

#method 
meth3 <- make.method(dat_fake)
meth3["measure_a"] <- "2l.pmm"

#impute with clustering
set.seed(3528)
imp5 <- mice(data = dat_fake, m = 5, predictorMatrix = pred3, method = meth3,
             maxit = 5)
#> 
#>  iter imp variable
#>   1   1  measure_a
#> Error in str2lang(x): <text>:1:20: unexpected symbol
#> 1: dv._lmer ~ 1+Typeb b
#>                        ^

#impute without clustering
set.seed(3528)
imp5 <- mice(data = dat_fake, m = 5, #predictorMatrix = pred3, method = meth3,
             maxit = 5)
#> 
#>  iter imp variable
#>   1   1  measure_a
#>   1   2  measure_a
#>   1   3  measure_a
#>   1   4  measure_a
#>   1   5  measure_a
#>   2   1  measure_a
#>   2   2  measure_a
#>   2   3  measure_a
#>   2   4  measure_a
#>   2   5  measure_a
#>   3   1  measure_a
#>   3   2  measure_a
#>   3   3  measure_a
#>   3   4  measure_a
#>   3   5  measure_a
#>   4   1  measure_a
#>   4   2  measure_a
#>   4   3  measure_a
#>   4   4  measure_a
#>   4   5  measure_a
#>   5   1  measure_a
#>   5   2  measure_a
#>   5   3  measure_a
#>   5   4  measure_a
#>   5   5  measure_a
#No error
lmer(measure_a~Type + (1|id), data = dat_fake)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: measure_a ~ Type + (1 | id)
#>    Data: dat_fake
#> REML criterion at convergence: 4680.551
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  id       (Intercept) 6.049   
#>  Residual             5.236   
#> Number of obs: 719, groups:  id, 129
#> Fixed Effects:
#> (Intercept)      Typeb b      Typec c      Typed d  
#>    10.33600      0.28181     -0.36017      0.02311

#Convert the factor level to number
dat_fake$Type = as.factor(as.numeric(dat_fake$Type))

#run the imputation with clustering again
set.seed(3528)
imp5 <- mice(data = dat_fake, m = 5, predictorMatrix = pred3, method = meth3,
             maxit = 5)
#> 
#>  iter imp variable
#>   1   1  measure_a
#>   1   2  measure_a
#>   1   3  measure_a
#>   1   4  measure_a
#>   1   5  measure_a
#>   2   1  measure_a
#>   2   2  measure_a
#>   2   3  measure_a
#>   2   4  measure_a
#>   2   5  measure_a
#>   3   1  measure_a
#>   3   2  measure_a
#>   3   3  measure_a
#>   3   4  measure_a
#>   3   5  measure_a
#>   4   1  measure_a
#>   4   2  measure_a
#>   4   3  measure_a
#>   4   4  measure_a
#>   4   5  measure_a
#>   5   1  measure_a
#>   5   2  measure_a
#>   5   3  measure_a
#>   5   4  measure_a
#>   5   5  measure_a

So this error is related to the illegal factor levels names ('a a', 'b b', 'c c')?
But when I ran the model with lmer, it didn't show any error.

Thank you.

Irene