MultiLevelOptimalBayes (MLOB) is designed for estimating two-level latent variable models, particularly in small sample settings. This is especially useful in psychology, education, and other fields with hierarchical or nested data structures. We present the R package MultiLevelOptimalBayes (MLOB) for estimating between-group effects in multilevel latent variable models. MLOB employs a regularised Bayesian estimator devised by Dashuk, Hecht, Luedtke, Robitzsch, and Zitzmann (2025a), which was subsequently enhanced for additional covariates by the same authors. This estimator chooses prior parameters to minimise the mean squared error (MSE) of the between-group effect by effectively balancing bias and variance. The regularised Bayesian estimator provides MSE-optimal estimations due to the mean-variance tradeoff, especially in scenarios of small sample sizes and poor intraclass correlation (ICC). The MLOB software supports imbalanced group sizes through integrated data-balancing methods and offers comprehensive inference, including point estimates, standard errors, p-values, and confidence intervals for both primary regressors and covariates. To gain comprehensive understanding, we initially examine the theoretical underpinnings of the regularised Bayesian estimator (Dashuk et al. 2025a, 2025b), followed by a discussion of its implementation in MLOB, namely the core function mlob(). We illustrate the application of mlob() using real datasets. Consequently, we provide researchers in psychology, education, and related disciplines a robust, user-friendly instrument for dependable multilevel latent variable estimation, particularly in contexts characterised by small sample sizes and low ICCs.
The core function mlob() estimates the between-group
coefficient (beta_b) using a regularized Bayesian approach,
and applies a data balancing procedure if the groups are unbalanced.
Below is the signature for the mlob() function. This
shows the arguments an theit default values you can pass, but
note that this chunk is not meant to be executed.
mlob(
formula,
data,
group,
balancing.limit = 0.2,
conf.level = 0.95,
jackknife = FALSE,
punish.coeff = 2,
...
)Arguments:
formula: A formula (e.g., Y ~ X + C) where
Y is the outcome, X is the context variable of interest, and C
represents covariates.data: A data.frame containing all variables in the
formula and the grouping variable.group: A string naming the grouping variable.balancing.limit: Proportion (0-1) of the dataset that
can be removed to balance group sizes. Default is 0.2.conf.level: Confidence level for confidence intervals.
Default is 0.95 (95% CI).jackknife: Logical. If TRUE, standard errors and CIs
are computed via jackknife resampling. Default is FALSE.punish.coeff: Multiplier for penalizing removal of
entire groups during balancing. Higher values discourage full-group
deletions.Balancing Procedure: The mlob() function also verifies whether the data is balanced, that is each group consist of exactly the same number of individuals. If the data is unbalanced, the balancing procedure comes into effect and identifies the optimal number of individuals and groups to delete based on the punishment coefficient. If the amount of data to be deleted is more than the threshold (balancing.limit), the regularized Bayesian estimator is not calculated and mlob() produces an error. This forces the user to increase the balancing limit manually and warns that the results should be interpreted with caution. # Examples
result_iris <- mlob(
Sepal.Length ~ Sepal.Width + Petal.Length,
data = iris,
group = "Species",
conf.level = 0.99,
jackknife = FALSE
)
summary(result_iris)
#> Call:
#> mlob(Sepal.Length ~ Sepal.Width + Petal.Length, data = iris, group = Species, conf.level = 0.99, jackknife = FALSE)
#>
#> Summary of Coefficients:
#> Estimate Std. Error Lower CI (99%) Upper CI (99%) Z value
#> beta_b 0.2957937 0.30457475 -0.4887389 1.080326 0.9711694
#> gamma_Petal.Length 0.4679522 0.05038678 0.3381645 0.597740 9.2872029
#> Pr(>|z|) Significance
#> beta_b 0.331464
#> gamma_Petal.Length 0.000000 ***
#>
#>
#> For comparison, summary of coefficients from unoptimized analysis (ML):
#> Estimate Std. Error Lower CI (99%) Upper CI (99%) Z value
#> beta_b 0.6027440 0.87389866 -1.6482698 2.853758 0.6897184
#> gamma_Petal.Length 0.4679522 0.05038678 0.3381645 0.597740 9.2872029
#> Pr(>|z|) Significance
#> beta_b 0.4903713
#> gamma_Petal.Length 0.0000000 ***
#>
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Note:
#> The standard error from unoptimized ML estimation is about 186.9% larger than the standard error obtained through our optimization procedure,
#> meaning that the optimized estimates are more accurate.
#> Concerning the estimates themselves, the unoptimized ML estimates may
#> differ greatly from the optimized estimates and should not be reported.
#> As the optimized estimates are always at least as accurate as the
#> unoptimized ML estimates,
#> please use them and their corresponding standard errors (first table of
#> output) for interpretation and reporting.
#> For more information, see Dashuk et al. (2025a).result_chick <- mlob(
weight ~ Time,
data = ChickWeight,
group = "Diet",
punish.coeff = 1.5,
jackknife = FALSE
)
print(result_chick)
#> Call:
#> mlob(weight ~ Time, data = ChickWeight, group = Diet, conf.level = 0.95, jackknife = FALSE)
#>
#> Coefficients
#> beta_b
#> -0.0171691
#>
#> Standard_Error
#> beta_b
#> 0.05542791
#>
#> Confidence_Interval (95%)
#> Lower Upper
#> beta_b -0.1258058 0.0914676
#>
#> Z value
#> beta_b
#> -0.3097556
#>
#> p value
#> beta_b
#> 0.7567468
summary(result_chick)
#> Call:
#> mlob(weight ~ Time, data = ChickWeight, group = Diet, conf.level = 0.95, jackknife = FALSE)
#>
#> Summary of Coefficients:
#> Estimate Std. Error Lower CI (95%) Upper CI (95%) Z value Pr(>|z|)
#> beta_b -0.0171691 0.05542791 -0.1258058 0.0914676 -0.3097556 0.7567468
#> Significance
#> beta_b
#>
#>
#> For comparison, summary of coefficients from unoptimized analysis (ML):
#> Estimate Std. Error Lower CI (95%) Upper CI (95%) Z value Pr(>|z|)
#> beta_b 2.209348 7.137118 -11.77915 16.19784 0.3095574 0.7568975
#> Significance
#> beta_b
#>
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Note:
#> The standard error from unoptimized ML estimation is about 12776% larger than the standard error obtained through our optimization procedure,
#> meaning that the optimized estimates are more accurate.
#> Concerning the estimates themselves, the unoptimized ML estimates may
#> differ greatly from the optimized estimates and should not be reported.
#> As the optimized estimates are always at least as accurate as the
#> unoptimized ML estimates,
#> please use them and their corresponding standard errors (first table of
#> output) for interpretation and reporting.
#> For more information, see Dashuk et al. (2025a).Interpretation of the results for the ChickWeight dataset
All chicks are weighed at the same time points, so Time
hardly varies between the diet groups: its estimated between-group
variance is not positive, and mlob() warns that the
estimates of beta_b and their standard errors are not
meaningful for these data. The example illustrates the balancing
procedure and this warning; for a substantive analysis, the predictor
needs to vary between groups.
result_mtcars <- mlob(
mpg ~ hp + wt + am + hp:wt + hp:am,
data = mtcars,
group = "cyl",
balancing.limit = 0.35
)
summary(result_mtcars)
#> Call:
#> mlob(mpg ~ hp + wt + am + hp:wt + hp:am, data = mtcars, group = cyl, balancing.limit = 0.35, conf.level = 0.95)
#>
#> Summary of Coefficients:
#> Estimate Std. Error Lower CI (95%) Upper CI (95%) Z value
#> beta_b -0.022150377 0.02200285 -0.06527516 0.02097441 -1.0067051
#> gamma_wt -5.109432015 7.60951564 -20.02380861 9.80494458 -0.6714530
#> gamma_am 5.194266715 12.06095092 -18.44476271 28.83329614 0.4306681
#> gamma_hp:wt 0.007257284 0.04403082 -0.07904153 0.09355610 0.1648228
#> gamma_hp:am -0.047559777 0.10255003 -0.24855413 0.15343458 -0.4637715
#> Pr(>|z|) Significance
#> beta_b 0.3140765
#> gamma_wt 0.5019320
#> gamma_am 0.6667097
#> gamma_hp:wt 0.8690834
#> gamma_hp:am 0.6428115
#>
#>
#> For comparison, summary of coefficients from unoptimized analysis (ML):
#> Estimate Std. Error Lower CI (95%) Upper CI (95%) Z value
#> beta_b -0.044388444 0.06405843 -0.16994066 0.08116377 -0.6929368
#> gamma_wt -5.109432015 7.60951564 -20.02380861 9.80494458 -0.6714530
#> gamma_am 5.194266715 12.06095092 -18.44476271 28.83329614 0.4306681
#> gamma_hp:wt 0.007257284 0.04403082 -0.07904153 0.09355610 0.1648228
#> gamma_hp:am -0.047559777 0.10255003 -0.24855413 0.15343458 -0.4637715
#> Pr(>|z|) Significance
#> beta_b 0.4883492
#> gamma_wt 0.5019320
#> gamma_am 0.6667097
#> gamma_hp:wt 0.8690834
#> gamma_hp:am 0.6428115
#>
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Note:
#> The standard error from unoptimized ML estimation is about 191.1% larger than the standard error obtained through our optimization procedure,
#> meaning that the optimized estimates are more accurate.
#> Concerning the estimates themselves, the unoptimized ML estimates may
#> differ greatly from the optimized estimates and should not be reported.
#> As the optimized estimates are always at least as accurate as the
#> unoptimized ML estimates,
#> please use them and their corresponding standard errors (first table of
#> output) for interpretation and reporting.
#> For more information, see Dashuk et al. (2025a).The output is an object of class mlob_result, which
contains:
beta_b and gamma values)The mlob_result object supports a comprehensive set of
S3 methods that follow standard R conventions, making it easy to work
with results in familiar ways. Here are all available methods:
# Get a basic result for demonstration
result <- mlob(weight ~ Time, data = ChickWeight, group = 'Diet', jackknife = FALSE)
# Print method - displays coefficients, standard errors, confidence intervals, Z-values, and p-values
print(result)
#> Call:
#> mlob(weight ~ Time, data = ChickWeight, group = Diet, conf.level = 0.95, jackknife = FALSE)
#>
#> Coefficients
#> beta_b
#> -0.01861521
#>
#> Standard_Error
#> beta_b
#> 0.06168138
#>
#> Confidence_Interval (95%)
#> Lower Upper
#> beta_b -0.1395085 0.1022781
#>
#> Z value
#> beta_b
#> -0.3017963
#>
#> p value
#> beta_b
#> 0.7628074# Summary method - comprehensive summary with significance stars and comparison to unoptimized ML
summary(result)
#> Call:
#> mlob(weight ~ Time, data = ChickWeight, group = Diet, conf.level = 0.95, jackknife = FALSE)
#>
#> Summary of Coefficients:
#> Estimate Std. Error Lower CI (95%) Upper CI (95%) Z value
#> beta_b -0.01861521 0.06168138 -0.1395085 0.1022781 -0.3017963
#> Pr(>|z|) Significance
#> beta_b 0.7628074
#>
#>
#> For comparison, summary of coefficients from unoptimized analysis (ML):
#> Estimate Std. Error Lower CI (95%) Upper CI (95%) Z value Pr(>|z|)
#> beta_b 2.520974 8.362585 -13.86939 18.91134 0.3014587 0.7630648
#> Significance
#> beta_b
#>
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Note:
#> The standard error from unoptimized ML estimation is about 13458% larger than the standard error obtained through our optimization procedure,
#> meaning that the optimized estimates are more accurate.
#> Concerning the estimates themselves, the unoptimized ML estimates may
#> differ greatly from the optimized estimates and should not be reported.
#> As the optimized estimates are always at least as accurate as the
#> unoptimized ML estimates,
#> please use them and their corresponding standard errors (first table of
#> output) for interpretation and reporting.
#> For more information, see Dashuk et al. (2025a).# Extract coefficients as a data frame
coef(result)
#> beta_b
#> 1 -0.01861521
# Extract standard errors
se(result)
#> beta_b
#> 0.06168138
# Extract variance-covariance matrix (diagonal only)
vcov(result)
#> beta_b
#> 0.003804592
# Extract confidence intervals
confint(result)
#> 2.5% 97.5%
#> beta_b -0.1395085 0.1022781
# Extract confidence intervals for specific parameters
confint(result, "beta_b")
#> 2.5% 97.5%
#> beta_b -0.1395085 0.1022781
# Extract confidence intervals with different confidence level
confint(result, level = 0.99)
#> 0.5% 99.5%
#> beta_b -0.1774959 0.1402655# Convert results to a data frame format
as.data.frame(result)
#> Estimate Std. Error Lower CI (95%) Upper CI (95%) Z value
#> beta_b -0.01861521 0.06168138 -0.1395085 0.1022781 -0.3017963
#> Pr(>|z|)
#> beta_b 0.7628074
# Get dimensions (number of parameters)
dim(result)
#> [1] 1 1
# Get number of parameters
length(result)
#> [1] 1
# Get parameter names
names(result)
#> [1] "beta_b"# Update model with new parameters (e.g., different confidence level)
updated_result <- update(result, conf.level = 0.99)
summary(updated_result)
#> Call:
#> mlob(weight ~ Time, data = data, group = Diet, conf.level = 0.99, jackknife = FALSE)
#>
#> Summary of Coefficients:
#> Estimate Std. Error Lower CI (99%) Upper CI (99%) Z value
#> beta_b -0.01861521 0.06168138 -0.1774959 0.1402655 -0.3017963
#> Pr(>|z|) Significance
#> beta_b 0.7628074
#>
#>
#> For comparison, summary of coefficients from unoptimized analysis (ML):
#> Estimate Std. Error Lower CI (99%) Upper CI (99%) Z value Pr(>|z|)
#> beta_b 2.520974 8.362585 -19.01962 24.06157 0.3014587 0.7630648
#> Significance
#> beta_b
#>
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Note:
#> The standard error from unoptimized ML estimation is about 13458% larger than the standard error obtained through our optimization procedure,
#> meaning that the optimized estimates are more accurate.
#> Concerning the estimates themselves, the unoptimized ML estimates may
#> differ greatly from the optimized estimates and should not be reported.
#> As the optimized estimates are always at least as accurate as the
#> unoptimized ML estimates,
#> please use them and their corresponding standard errors (first table of
#> output) for interpretation and reporting.
#> For more information, see Dashuk et al. (2025a).You can discover all available methods for mlob_result
objects using:
Here’s a practical example showing how to use multiple methods together:
# Run analysis
result <- mlob(weight ~ Time, data = ChickWeight, group = 'Diet', jackknife = FALSE)
# Get basic information
cat("Number of parameters:", length(result), "\n")
#> Number of parameters: 1
cat("Parameter names:", paste(names(result), collapse = ", "), "\n")
#> Parameter names: beta_b
# Extract key statistics
coefficients <- coef(result)
standard_errors <- se(result)
confidence_intervals <- confint(result, level = 0.99)
# Create a custom summary table
custom_summary <- data.frame(
Parameter = names(result),
Estimate = as.numeric(coefficients),
SE = as.numeric(standard_errors),
CI_Lower = confidence_intervals[, 1],
CI_Upper = confidence_intervals[, 2]
)
print(custom_summary)
#> Parameter Estimate SE CI_Lower CI_Upper
#> 1 beta_b 0.05792323 0.1159205 -0.2406683 0.3565148All these methods follow standard R conventions, making your
mlob_result objects compatible with existing R workflows
and familiar to users of other statistical packages.
While MultiLevelOptimalBayes provides a robust solution
for regularized estimation in two-level models, users should be aware of
the following limitations:
Local Grid Search:
The prior variance is searched within five standard deviations of the ML
estimate of the between-group variance of the predictor (bounded below
by 1% of its total variance). The minimum can lie on the boundary of
this region, especially for few groups and a small ICC of the predictor;
the estimate is then the best one within the searched region rather than
the global minimum.
Assumption of Equal Group Sizes:
The estimator assumes equal group sizes to simplify the model. While
averaging group sizes is a proposed solution, the method does not yet
handle unbalanced group sizes natively, but find the optimal reduction
to the balanced size.
Jackknife Computational Cost:
Jackknife resampling improves interval coverage in small samples. Note
that it may be computationally demanding for larger samples.
Limited Model Scope:
Currently, MLOB only handles two-level models with continuous outcomes.
Extensions to support generalized linear mixed models, three- and
more-level structures, or multivariate outcomes are not yet
available.
Dashuk, V., Hecht, M., Luedtke, O., Robitzsch, A., & Zitzmann, S. (2025a). An Optimally Regularized Estimator of Multilevel Latent Variable Models, with Improved MSE Performance. https://doi.org/10.1017/psy.2025.10045
Dashuk, V., Hecht, M., Lüdtke, O., Robitzsch, A., & Zitzmann, S. (2025b). Estimating context effects in small samples while controlling for covariates: an optimally regularized Bayesian estimator for multilevel latent variable models. https://doi.org/10.1007/s41237-025-00264-7
Luedtke, O., Marsh, H. W., Robitzsch, A., et al. (2008).
The multilevel latent covariate model: A new, more reliable approach to
group-level effects in contextual studies.
https://doi.org/10.1037/a0012869