MMLN: An R Package for Mixed-Effects Multinomial Logistic-Normal Regression and Model Diagnostics

MMLN: An R Package for Mixed-Effects Multinomial Logistic-Normal Regression and Model Diagnostics Eric A. E. Gerber1 1Northeastern University, 360 Huntington Ave, Boston, MA 02115 Abstract Multinomial outcomes arise in numerous fields---from sports and species counts to genomics---yet existing software often focuses on simpler fixed-effects or purely multinomial (logistic) frameworks. The MMLN package introduces a suite of functions to fit more complex multinomial logistic-normal regression models, including incorporation of random effects, and evaluate the fit of all multinomial regression models using the squared Mahalanobis distance residuals [6]. The MMLN() function fits mixed-effects multinomial logistic-normal models via MCMC sampling, while the MDres() function computes the squared Mahalanobis distance residuals to comprehensively evaluate model adequacy. Users can visualize or formally test these using quantilequantile plots and Kolmogorov-Smirnov tests. We describe the design and usage of the functions provided in the MMLN package and demonstrate the package's capabilities for modeling multinomial data by integrating flexible modeling tools, summaries, visualization, and robust diagnostics. Key Words: multinomial regression, mixed effects models, residuals, diagnostics 1. Introduction Multinomial logistic-normal (MLN) models extend the classical multinomial framework by embedding the probability vectors of the multinomial data in a latent Gaussian space, more robustly accounting for overdispersion among categories than traditional multinomial logit or hierarchical multinomial-Dirichlet models. Bayesian hierarchical models are powerful tools for fitting both fixed-effect and mixed-effect MLN models used for capturing overall and group-level covariance structures among the multinomial categories. Posterior predictive checks provide essential diagnostics for assessing model adequacy but have historically been difficult to derive or unintuitive for multinomial models of any type. Squared Mahalanobis residuals calculated using samples from predictive distributions have been developed to address this gap [6]. The MMLN package implements both fixed-effects and mixed-effects MLN models using Gibbs sampling with flexible Metropolis-Hastings updates, accompanied by tools for residual analysis in R. The specific MLN functions add to the lexicon of established multinomial regression models, while the MDres() function is simple to use for assessing model fits under any multinomial regression framework, from models as complex as MLN to as simple as basic multinomial logit. 2. Multinomial Logistic-Normal Models and Diagnostics Define the general form for a nominal outcome regression model: 𝑦𝑦𝑖𝑖∼ ℳ 𝐽𝐽(𝑛𝑛𝑖𝑖,𝜋𝜋𝑖𝑖) where ℳ 𝐽𝐽 represents the multinomial distribution with 𝐽𝐽 categories, exposure 𝑛𝑛𝑖𝑖> 1, and probability vector 𝜋𝜋𝑖𝑖. The commonly fit model to these data, the multinomial logistic, struggles to account for overdispersion and cannot account for positive correlations between the outcome counts of the different categories. The multinomial Dirichlet model has been shown as one option for accounting for overdispersion [7]. However, the multinomial logistic-normal is a more flexible alternative, which does a better job of accounting for positive correlations between categories [1]. Multinomial logistic-normal models arise from modeling the probability vector using the inverse additive logistic ratio transformation, where a multivariate normal noise term is added to the log odds, in other words (and equivalently): (𝑙𝑙𝑙𝑙𝑙𝑙(𝜋𝜋𝑖𝑖1 𝜋𝜋𝑖𝑖𝐽𝐽 ), ⋯,𝑙𝑙𝑙𝑙𝑙𝑙(𝜋𝜋𝑖𝑖(𝐽𝐽−1) 𝜋𝜋𝑖𝑖𝐽𝐽 )) ∼ 𝑁𝑁𝐽𝐽−1 𝜋𝜋𝑖𝑖=𝑎𝑎𝑙𝑙𝑟𝑟−1(𝑍𝑍𝑖𝑖+𝜀𝜀𝑖𝑖) 𝜀𝜀𝑖𝑖∼ 𝑁𝑁𝐽𝐽−1(0, Σ) 2.1 Fixed-Effects MLN The FMLN() function fits the purely fixed effects version of multinomial logistic-normal regression [5], where: 𝜋𝜋𝑖𝑖=𝑎𝑎𝑙𝑙𝑟𝑟−1(𝑋𝑋𝑖𝑖𝛽𝛽+𝜀𝜀𝑖𝑖) The model is fit via a Gibbs sampler with Metropolis or Metropolis-Hastings (depending on proposal distribution) proposals for the latent 𝑊𝑊𝑖𝑖=𝑋𝑋𝑖𝑖𝛽𝛽+𝜀𝜀𝑖𝑖. The algorithm iteratively proceeds: (1) Sample all 𝑊𝑊𝑖𝑖|𝑌𝑌 𝑖𝑖,𝑋𝑋𝑖𝑖,𝛽𝛽,Σ ∈ ℝ𝐽𝐽−1 via Metropolis-Hastings (latent variables; depending on proposal distribution) (2) Sample 𝛽𝛽|𝑊𝑊,𝑋𝑋,Σ ∼ 𝑁𝑁 (fixed effects) (3) Sample Σ|𝑊𝑊,𝑋𝑋,𝛽𝛽 ∼ Inv-Wishart (residual covariance) 2.2 Mixed-Effects MLN To accommodate more complex data structures, including group level variation (for 𝑚𝑚= 1, ⋯,𝑀𝑀 groups), the MMLN() function introduces estimation of random intercepts, where: 𝜋𝜋𝑚𝑚𝑖𝑖 =𝑎𝑎𝑙𝑙𝑟𝑟−1(𝑋𝑋𝑚𝑚𝑖𝑖𝛽𝛽+𝜓𝜓𝑚𝑚+𝜀𝜀𝑚𝑚𝑖𝑖) and 𝜓𝜓𝑚𝑚∼ 𝑁𝑁𝐽𝐽−1(0, Φ) while the rest of the model parameterization remains the same. This mixed effects multinomial logistic-normal model [5] is fit via a Metropolis-within-Gibbs sampler extended from the one used to fit the fixed effects version: (1) Sample all 𝑊𝑊𝑚𝑚𝑖𝑖|𝑌𝑌 𝑚𝑚𝑖𝑖,𝑋𝑋𝑚𝑚𝑖𝑖,𝛽𝛽,𝜓𝜓𝑚𝑚,Σ ∈ ℝ𝐽𝐽−1 via Metropolis-Hastings (latent variables; depending on proposal distribution) (2) Sample 𝜓𝜓𝑚𝑚|𝑊𝑊,𝑋𝑋,𝛽𝛽,Σ,Φ ∼ 𝑁𝑁 (random effects) (3) Sample Φ|Ψ ∼ Inv-Wishart (random effects covariance; Ψ is the full matrix of random effects) (4) Sample 𝛽𝛽|𝑊𝑊,𝑋𝑋,Ψ,Σ ∼ 𝑁𝑁 (fixed effects) (5) Sample Σ|𝑊𝑊,𝑋𝑋,𝛽𝛽 ∼ Inv-Wishart (residual covariance) 2.3 Mahalanobis Residuals Recently, randomized quantile residuals for binary outcomes [3] were extended for use in nominal multinomial modeling frameworks [6]. These residuals take the form of squared Mahalanobis distances in the transformed additive log-ratio space from the observed data to the samples from a predictive distribution under any model fit. These residuals, implemented in the MDres() function, are not specific to any multinomial model. Generally, for each observation 𝑖𝑖, 𝐾𝐾 samples of predicted counts 𝑦𝑦𝑖𝑖𝑖𝑖 are generated from a fitted multinomial model to form the sampling distribution of the additive log-ratio transformed vectors: 𝑤𝑤𝑖𝑖𝑖𝑖 =𝑎𝑎𝑙𝑙𝑟𝑟(𝑦𝑦𝑖𝑖𝑖𝑖 ∗), 𝑘𝑘= 1,2, … , 𝐾𝐾 Squared Mahalanobis distances of these model-generated log-odds, {𝑤𝑤𝑖𝑖}, as well as for the observed log-odds, 𝑤𝑤𝑖𝑖 𝑜𝑜𝑜𝑜𝑜𝑜 are calculated: 𝑀𝑀𝐷𝐷𝑖𝑖𝑖𝑖 2= ((𝑤𝑤𝑖𝑖𝑖𝑖 − 𝑤𝑤𝑖𝑖)𝑇𝑇Σ ^ 𝑤𝑤𝑖𝑖 −1((𝑤𝑤𝑖𝑖𝑖𝑖 − 𝑤𝑤𝑖𝑖)) and 𝑀𝑀𝐷𝐷𝑖𝑖 2(𝑜𝑜𝑜𝑜𝑜𝑜)= ((𝑤𝑤𝑖𝑖 𝑜𝑜𝑜𝑜𝑜𝑜 − 𝑤𝑤𝑖𝑖)𝑇𝑇Σ ^ 𝑤𝑤𝑖𝑖 −1((𝑤𝑤𝑖𝑖 𝑜𝑜𝑜𝑜𝑜𝑜 − 𝑤𝑤𝑖𝑖)) where 𝑤𝑤𝑖𝑖 is the sample mean of model-generated log-odds and Σ ^ 𝑤𝑤𝑖𝑖 is their sample covariance matrix. A percentile for the observed distance, 𝑀𝑀𝐷𝐷𝑖𝑖 2(𝑜𝑜𝑜𝑜𝑜𝑜) relative to the empirical cdf of the model-generated distances, 𝐹𝐹 ^ 𝐾𝐾(𝑀𝑀𝐷𝐷𝑖𝑖 2) is calculated. A uniform random variable is then generated where the minimum, 𝑎𝑎𝑖𝑖, and maximum, 𝑏𝑏𝑖𝑖, depend on the observed value's ordered location among the modelbased 𝑀𝑀𝐷𝐷𝑖𝑖 2: (1) If 𝑀𝑀𝐷𝐷𝑖𝑖 2(𝑜𝑜𝑜𝑜𝑜𝑜)≤ 𝑚𝑚𝑖𝑖𝑛𝑛(𝑀𝑀𝐷𝐷𝑖𝑖 2) 𝑢𝑢𝑖𝑖∼ 𝒰𝒰(𝑎𝑎𝑖𝑖= 0, 𝑏𝑏𝑖𝑖=𝐹𝐹 ^ 𝐾𝐾(𝑚𝑚𝑖𝑖𝑛𝑛(𝑀𝑀𝐷𝐷𝑖𝑖 2))) (2) If 𝑀𝑀𝐷𝐷𝑖𝑖 2(𝑜𝑜𝑜𝑜𝑜𝑜)>𝑚𝑚𝑎𝑎𝑚𝑚(𝑀𝑀𝐷𝐷𝑖𝑖 2) 𝑢𝑢𝑖𝑖∼ 𝒰𝒰(𝑎𝑎𝑖𝑖=𝐹𝐹 ^ 𝐾𝐾(𝑚𝑚𝑎𝑎𝑚𝑚(𝑀𝑀𝐷𝐷𝑖𝑖 2)), 𝑏𝑏𝑖𝑖= 1) (3) Otherwise, 𝑢𝑢𝑖𝑖∼ 𝒰𝒰(𝑎𝑎𝑖𝑖=𝐹𝐹 ^ 𝐾𝐾(𝑀𝑀𝐷𝐷 ~ 𝑖𝑖 2), 𝑏𝑏𝑖𝑖=𝐹𝐹 ^ 𝐾𝐾(𝑀𝑀𝐷𝐷𝑖𝑖 2(𝑜𝑜𝑜𝑜𝑜𝑜))) where 𝑀𝑀𝐷𝐷 ~ 𝑖𝑖 2=𝑚𝑚𝑎𝑎𝑚𝑚(𝑀𝑀𝐷𝐷𝑖𝑖 2<𝑀𝑀𝐷𝐷𝑖𝑖 2(𝑜𝑜𝑜𝑜𝑜𝑜)) If the data fit the model, these percentiles are distributed𝒰𝒰(0,1). These percentiles are backtransformed to standard normal values to serve as residuals, 𝑟𝑟𝑖𝑖=Φ−1(𝑢𝑢𝑖𝑖). Normal quantile-quantile and residual plots are then used to assess fit. 3. The MMLN R Package 3.1 Package Structure The package is organized into four main R scripts: • mln_helpers.R: utility functions • mln_functions.R: core MCMC samples FMLN() and MMLN() • multi_res: Mahalanobis residuals MDres(), summary and plotting helpers • real_data_examples.R: vignettes for applying the functions to real data sets 3.2 Function Overview Table 1 gives the description for the three primary functions of the package, as well as the three most important helper functions. There are several other functions which are discussed as needed. 3.3 Using FMLN The FMLN() function takes as its primary arguments the count matrix 𝑌𝑌 and input matrix 𝑋𝑋 of fixed-effects covariates. Parameters also include the total number of MCMC iterations, burn-in (number of initial iterations to discard), thinning interval, scaling factor for Metropolis-Hastings proposal covariance, settings for the prior distributions on the fixed effects and residual covariance matrix, and choice of proposal distribution for the Metropolis-Hastings. The verbose argument allows for printing of progress updates. _______________________________________________________________________________ res_f <- FMLN( Y = sim$Y, X = sim$X, n_iter = 2000, burn_in = 500, thin = 2, proposal = "normbeta", verbose = TRUE ) _______________________________________________________________________________ 3.4 Using MMLN The MMLN() function behaves similarly to the FMLN(), though now requires the definition of the 𝑍𝑍 random effects design matrix, currently only supporting random intercepts for group-level observations. All other arguments remain the same, though the prior settings also now account for the inclusion of the random effects. _______________________________________________________________________________ res_m <- MMLN( Y = sim$Y, X = sim$X, Z = sim$Z, n_iter = 2000, burn_in = 500, thin = 2, proposal = "normbeta", verbose = TRUE ) _______________________________________________________________________________ Table 1: Primary MMLN Package Functions Function Description Fit fixed-effects MLN model via MH-Gibbs sampling Fit mixed-effects MLN model with group-level random intercepts Generate traceplots and posterior summary tables Calculate DIC from log-likelihood samples Simulate posterior predictive counts for model checking Compute Mahalanobis residuals FMLN MMLN plot_trace_and_summary compute_dic sample_posterior_predictive MDres 3.5 Diagnostic Tools Trace plots and posterior Markov chain summaries can be displayed with the plot_trace_and_summary() function after passing one of the posterior chain objects returned by either FMLN() or MMLN() through the simplify2array() function. By default, the trace plots are displayed in groups of four. _______________________________________________________________________________ beta_chain_array <- simplify2array(res_m$beta_chain) trace_stats <- plot_trace_and_summary(beta_chain_array, "beta") trace_stats _______________________________________________________________________________ The Deviance Information Criterion (DIC) [2] for comparing model fits is computed via the compute_dic() function after using the returned posterior chains and the true data counts to estimate the log likelihood functions using the dmnl_loglik() function. _______________________________________________________________________________ ll_chain <- sapply(res_m$w_chain, function(W) dmnl_loglik(W, sim$Y)) W_hat <- alr(compress_counts(sim$Y) / rowSums(sim$Y)) ll_hat <- dmnl_loglik(W_hat, sim$Y) dic_res <- compute_dic(ll_chain, ll_hat) _______________________________________________________________________________ Finally, the squared Mahalanobis residuals can be computed for any set of predictive distribution samples for each observation using the MDres() function. The function has a summary() class method which prints out the results of the Kolmogorov-Smirnov test for normality as a formal test of model fit and displays the normal quantile-quantile plot of the residuals for a convenient graphical assessment. _______________________________________________________________________________ Y_pred_list <- lapply(seq_along(res_m$w_chain), function(i) { sample_posterior_predictive(X = sim$X, beta = res_m$beta_chain[[i]], Sigma = res_m$sigma_chain[[i]], n = sim$n, Z = sim$Z, psi = res_m$psi_chain[[i]], mixed = TRUE ) }) resids <- MDres(sim$Y, Y_pred_list) summary(resids) _______________________________________________________________________________ 3.5.1 Example Output The MMLN package also includes several vignettes for demonstrating the implementation and utility of the models and diagnostic output on both simulated and real data. One vignette involves helper function, run_pollen_models(), that shows the residuals ability to capture the well-established existence of overdispersion [7] in pollen count data. There is also a simulate_mixed_mln_data() function which will simulate data from the MMLN model. As an example, we simulate data under the MMLN, then fit those data with both the MMLN() and FMLN() functions. The example Figure 1, and Kolmogorov-Smirnov test results presented demonstrate one use case of the Mahalanobis residuals and its summary class method. Figure 1: Example QQ-plots of summary(MDres) output for FMLN (left) and MMLN (right) models fit to MMLN data. _______________________________________________________________________________ > resids <- MDres(observed_counts, fitted_counts_list) > summary(resids) Kolmogorov-Smirnov test for normality of Mahalanobis residuals: Asymptotic one-sample Kolmogorov-Smirnov test D = 0.13429, p-value = 0.05429 alternative hypothesis: two-sided _______________________________________________________________________________ 4. Discussion and Future Work The MMLN package equips users with flexible tools for modeling multinomial outcomes in the presence of overdispersion, together with comprehensive diagnostics via squared Mahalanobis residuals. The modular design and simple interfaces facilitate usage and application to a wide range of data. Future extensions are planned to include handling of more robust random effects, as the infrastructure of the mixed effects model should be easily extended: 𝑊𝑊𝑚𝑚𝑖𝑖 =𝑋𝑋𝑚𝑚𝑖𝑖𝛽𝛽+𝑍𝑍𝑚𝑚𝑖𝑖𝜓𝜓𝑚𝑚+𝜀𝜀𝑚𝑚𝑖𝑖 𝑣𝑣𝑣𝑣𝑣𝑣(𝜓𝜓𝑚𝑚)∼ 𝑁𝑁(𝐽𝐽−1)𝑞𝑞(0, Φ) where, given 𝑞𝑞 group-level random covariates, the log-odds latent variables have the multivariate normal distribution: 𝑊𝑊𝑚𝑚𝑖𝑖|𝑋𝑋𝑚𝑚𝑖𝑖,𝛽𝛽,𝜓𝜓𝑚𝑚,Σ ∼ 𝑁𝑁(𝐽𝐽−1)(𝑋𝑋𝑚𝑚𝑖𝑖𝛽𝛽+𝑍𝑍𝑚𝑚𝑖𝑖𝜓𝜓𝑚𝑚,Σ) and, unconditionally: 𝑣𝑣𝑣𝑣𝑣𝑣(𝑊𝑊 𝑚𝑚)|𝑋𝑋𝑚𝑚𝑖𝑖,𝛽𝛽,Σ ∼ 𝑁𝑁(𝐽𝐽−1)𝑛𝑛𝑖𝑖(𝑣𝑣𝑣𝑣𝑣𝑣(𝑋𝑋𝑚𝑚𝛽𝛽), 𝑄𝑄𝑚𝑚 −1) where 𝑄𝑄𝑚𝑚 −1 = (𝑍𝑍𝑚𝑚⊗ 𝐼𝐼(𝐽𝐽−1))Φ(𝑍𝑍𝑚𝑚⊗ 𝐼𝐼(𝐽𝐽−1))𝑇𝑇+ (𝐼𝐼𝑛𝑛𝑚𝑚⊗ Σ) However, the addition of additional random effects drastically increases computation cost, and will thus require more robust implementation, perhaps by leveraging the Rcpp package [4] for integrating R and C++. In the future, support for alternative priors for the parameters of the Bayesian model may also be included. References [1] J. Aitchison. The Statistical Analysis of Compositional Data. Journal of the Royal Statistical Society: Series B (Methodological), 44(2): 139-177, 1982. [2] D. J. Spiegelhalter, N. G. Best, B. P. Carlin, and A. Van Der Linde. Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society Series B: Statistical Methodology, 64(34):583-639, 2002. [3] K. P. Dunn and G. K. Smyth. Randomized quantile residuals. Journal of Computational and Graphical Statistics, 5:1-10, 1996. [4] D. Eddelbuettel and R. Francois. Rcpp: Seamless R and C++ integration. Journal of Statistical Software, 40(8):1-18, 2011. [5] E. A. E. Gerber and B. A. Craig. A mixed effects multinomial logistic-normal model for forecasting baseball performance. Journal of Quantitative Analysis in Sports, 17(3):221239, 2021. [6] E. A. E. Gerber and B. A. Craig. Residuals and diagnostics for multinomial regression models. Statistical Analysis and Data Mining: An ASA Data Science Journal, 17(1):e11645, 2024. [7] J.E. Mosimann. On the Compound Multinomial Distribution, the Multivariate 𝛽𝛽Distribution, and Correlations among Proportions. Biometrika, 49(1-2):65-82, 1962.