This is an open access article distributed under the terms of the Single-case experimental designs (SCEDs) provide a valuable framework for analysing intervention effects on individuals through repeated measures (Bono & Arnau, 2014; Kazdin, 1982; Shadish & Sullivan, 2011). Their objective is to establish a functional relationship between an intervention and changes in a behavioural outcome, as applied contexts often hinder meeting traditional causal criteria (logical connection, covariation, temporal precedence, and full control of confounding variables; Kazdin, 1977; Manolov et al., 2014; Virués-Ortega & Haynes, 2005). While traditional qualitative methods, like visual inspection of time-series data, can be somewhat subjective and less effective for detecting subtle changes (Busk & Serlin, 1992; Kazdin, 1982; Kratochwill et al., 2013; Manolov & Moeyaert, 2017; Parsonson & Baer, 1986; Van den Noortgate & Onghena, 2003a), quantitative methods, particularly multilevel models, have increasingly addressed these considerations.
Multilevel linear models (MLMs) address the nested structure of the data, improving the precision of effect estimates and enabling the analysis of contextual variables at multiple levels (Hoffman, 2014; Manolov & Moeyaert, 2017; Moeyaert et al., 2020; Van den Noortgate & Onghena, 2003b). They effectively address statistical issues like autocorrelation, which are inherent to SCED data and can undermine the validity of parametric tests (Bono & Arnau, 2014; Gentile et al., 1972; Keselman & Leventhal, 1974). MLMs also complement non-parametric methods in SCED analysis, such as non-overlap indices (PND, PEM, NAP; Parker et al., 2011), which measure improvement across phases but have certain limitations (Rodríguez-Prada & Olmos, 2019). MLMs provide a comprehensive framework for statistical inference in SCEDs (Botella & Caperos, 2019).
MLMs can overcome some SCED challenges, such as model specification issues, limited information for estimating parameters and incorrect standard error estimations from techniques that overlook the data’s hierarchical structure (Rodabaugh & Moeyaert, 2017). By modelling hierarchical structures, MLMs efficiently capture trends and allow flexible error covariance specifications, enhancing robustness and reliability in statistical decisions (Hoffman, 2014). MLMs also analyse fixed effects, shared among participants, and random effects, which account for individual variability in treatment outcomes (Manolov & Moeyaert, 2017). This dual-level approach is particularly well-suited for SCEDs, where individual responses are central to understanding intervention effects, and the hierarchical nature of the data must be accounted for to ensure valid conclusions (Moeyaert et al., 2014; Van den Noortgate & Onghena, 2003a; 2003b). Researchers can assess individual factors influencing clinical outcomes by measuring between-individual variability by explicitly modelling moderators to address why some benefit more from interventions (Moeyaert et al., 2024; Moeyaert & Yang, 2021). MLMs handle complex data structures, specify covariance structures, and model variability, making them suitable for small-sample designs including SCEDs. However, challenges remain in estimating higher-level variance components with few individuals.
Two main approaches to estimate MLMs are frequentist and Bayesian frameworks, differing in parameter estimation and uncertainty quantification. Frequentist methods, typically using maximum likelihood estimation (ML), view parameters as fixed but unknown. Restricted Maximum Likelihood (REML) is preferred for random effects models, as it adjusts for the loss of degrees of freedom in variance component estimation. However, these approaches depend on asymptotic assumptions that require large samples for accuracy (Hoffman, 2014). In small-sample contexts like SCEDs, these assumptions often fail, resulting in biased covariance estimates despite unbiased fixed effect estimates (Baek et al., 2020; Moeyaert et al., 2017; Van den Noortgate & Onghena, 2003a).
Bayesian estimation has emerged as a promising alternative for small-sample designs in SCED studies (Baek et al., 2020; McNeish, 2016). Unlike frequentist methods, it incorporates prior information about parameters, helping to mitigate limitations of datasets with few Level-2 units (Gelman et al., 2013; van de Schoot et al., 2015). Selecting appropriate prior distributions is crucial in estimating Level-2 variance components, as weak or overly restrictive priors can impact results (Baek & Ferron, 2020; Moeyaert et al., 2017). Balancing informative and weakly informative priors is essential to prevent bias and ensure robustness. Additionally, Bayesian methods use advanced algorithms like Markov Chain Monte Carlo (MCMC; Brooks, 1998) to estimate parameters in complex models, enhancing precision and stability (McNeish, 2016; Rindskopf, 2014).
Whether using frequentist or Bayesian frameworks, a major challenge in applying MLMs to SCED data is selecting the best model and its random effect structure. This choice affects the balance between Type I error rates and statistical power, influencing hypothesis testing conclusions (Matuschek et al., 2017). Underparameterized models may oversimplify data and miss key patterns, while overparameterized models can inflate standard errors and reduce statistical power (Hoffman, 2014; Martínez-Huertas et al., 2022; Martínez-Huertas & Olmos, 2022). Thus, model selection is essential for valid and meaningful statistical decisions in SCED data interpretation, in combination with expert knowledge of the clinicians and applied researchers.
A practical approach uses information criteria to compare models, rewarding good fits and penalising overfitting. These criteria arise from Kullback-Leibler (KL) divergence and entropy, estimating model generalizability by prioritising predictive accuracy over simple goodness of fit. Frequentist methods employ the Akaike Information Criterion (AIC; Akaike, 1998) and Bayesian Information Criterion (BIC; Schwarz, 1978), with AIC favouring models that minimise information loss and BIC selecting the most probable model under a Bayesian framework with stricter complexity penalties (Raftery, 1995; Weakliem, 1999). Bayesian approaches use indices like the Watanabe-Akaike Information Criterion (WAIC; Watanabe, 2010) and Leave-One-Out Cross-Validation (LOO; Vehtari et al., 2017) to enhance model evaluation by incorporating posterior distributions. In contrast to frequentist criteria, WAIC prioritises predictive utility over determining the true population model (Nicenboim & Vasishth, 2016), while LOO assesses model performance using cross-validation across data points, which makes it robust to outliers. Both indices select the model with the lowest value, and WAIC and LOO often outperform AIC and Deviance Information Criterion (DIC; Spiegelhalter et al., 2002), offering greater reliability in comparisons (Gelman et al., 2013). In the present simulation study, the performance of all indices is assessed according to their ability to identify the data-generating model.
Aims of the Present Study
This study compares Bayesian and frequentist methods in MLMs for model selection in SCED data via simulation. It focuses on the information criteria AIC, BIC, WAIC, and LOO fit indices to select the true population model and analyze their impact on statistical power and Type I errors. Three objectives are established. The first objective is to evaluate if various Bayesian priors present differences regarding statistical power and Type I error rates. Building on Moeyaert et al. (2017) and Baek et al. (2020), we explore how prior specifications affect robustness and accuracy. The second objective assesses selection accuracy using frequentist (AIC, BIC) and Bayesian (WAIC, LOO) information criteria, aiming for consistency in the Bayesian framework. The third examines how model selection affects statistical power and Type I error rates in intervention effects, identifying optimal estimation methods under unknown population models. In this context, we evaluate the operating characteristics (Type I error and power) of different model selection procedures (AIC, BIC, WAIC, LOO) when the data-generating process is unknown. Thus, this study addresses the gap of analysing model selection's impact from an ecological viewpoint, reflecting researchers’ experiences in applied settings. In practical terms, this paper examines how AIC, BIC, WAIC, and LOO behave in terms of statistical power and Type I error rates when the true data-generating multilevel model is unknown — a common situation in applied SCED studies. Based on these results, the study aims to provide an evidence-based recommendation on which criterion offers the best trade-off between power and Type I error control under model uncertainty.
Method
Simulation Study Design
A Monte Carlo simulation of an SCED AB-design with baseline (A) and intervention (B) phases was conducted. The dependent variables included statistical power, Type I error rate, and proportion of true model selection per information criteria. Statistical power indicates the rate of correctly detecting the intervention effect, while the Type I error rate measures the incorrect identification of an effect when none exists. Correct model selection quantifies how often the information criteria identified the true population model during comparison.
Data Generation and Population Values
SCEDs data were generated based on four models varying in the number of random effects (covariance parameters). The models were:
1. Minimal model: Assumes a constant for the initial value and the intervention effect across individuals. This means that the dependent variable is initially absent or at a floor level for all individuals, and that the intervention effect does not present variation between individuals:
,
where is the response variable for individual j at observation i, represents baseline mean (intercept), is the average intervention effect (slope), is the Level-1 residual error for individual j at observation i and is a dummy-coded predictor indicating whether observation i for individual j belongs to the baseline phase (0) or the intervention phase (1). This model does not include random effects.
2. Partial intercepts model: Includes random intercepts for baseline differences among individuals. This model is relevant in SCED contexts where individuals have different baseline levels, but the intervention is expected to produce a uniform effect:
,
where is the random intercept in the baseline for individual j.
3. Partial slopes model: Incorporates random slopes to represent different intervention effects among individuals, essential in SCED contexts where responses are initially low or absent, or where individual characteristics may influence intervention outcomes:
,
where is the random slope for individual j.
4. Maximal model: Includes both random intercepts and slopes, accounting for variability in both baseline scores and intervention effects. This is the most ecological model because it considers the variability that an individual's learning history and personal factors can hold:
,
where and are random intercepts and slopes for individual j, respectively.
In all the population models, fixed effects were set at , random effects followed normal distributions , , and residual variance was set to . Random effects were uncorrelated (r = 0) to minimise model complexity, as suggested by Van den Noortgate and Onghena (2003a) and Moeyaert et al. (2017). In addition, residuals were assumed to be independent over time, and no autocorrelation or secular nor temporal trends were simulated.
Simulation Conditions
The simulation comprised 144 scenarios from various factor combinations:
1️. Number of individuals (Nj = 3, 5, 7) (Moeyaert et al., 2017; Shadish & Sullivan, 2011).
2️. Number of repeated measurements (RRMMj = 10, 20, 30, 40) (Moeyaert et al., 2017, Shadish & Sullivan, 2011). We added two additional values (10 and 30) to increase design resolution across short-to-moderate series lengths and to assess whether performance changes monotonically with series length.
3️. Effect size for intervention effect : null (d = 0), medium (d = 1.15), and large (d = 2.70), following Ferguson’s (2009) guidelines.
4️. Population model structure, ranging from a Minimal model (no random effects) to increasingly complex specifications (Partial Intercepts, Partial Slopes, and Maximal models). These random structures are in line with previous simulation studies (Martínez-Huertas et al., 2022; Martínez-Huertas & Olmos, 2022; Matuschek et al., 2017).
In each scenario, four model information criteria were evaluated: frequentist criteria (AIC, BIC) and Bayesian criteria (WAIC, LOO), allowing a comparative analysis of selection approaches. Seven different priors were tested to assess their effects on statistical power, Type I error rates for the intervention effect, and the accuracy of model selection using WAIC and LOO. In line with the default parameterisation in brms, weakly informative priors (see Gelman, 2006; McElreath, 2020) were specified on the standard deviations of the random effects (σᵤ0 and σᵤ1), not on the variance components (σ2ᵤ0 and σ2ᵤ1). Because standard deviations are constrained to be positive, Half-Cauchy and Half-normal priors with scale parameters of 10, 20, and 50 were used, alongside a weakly informative Uniform (0, 100) prior. Following the approach outlined by Moeyaert et al. (2017), who derived these values from reanalyses of empirical SCED studies with normally distributed continuous outcomes, our aim was to reflect plausible ranges for variance components in this context. This approach allows the priors to regularize estimation in small-sample settings while still letting the normally distributed outcome data contribute substantially to the results. Increasing the scale parameter in the Half-Cauchy and Half-Normal distributions makes the prior less informative by allowing for greater variability in the variance components. A smaller scale concentrates the prior mass closer to zero, while a larger scale spreads the distribution and assigns more weight to larger variance values, thus reducing the prior's influence. We adopted this approach because one of the main challenges in this context is selecting the appropriate prior distribution. Although future studies should explore this issue in greater depth, our current aim is to increase variability so that the prior parameterization remains sensible while still allowing the data to provide valuable insight.
Conducting the Simulation Study
Five hundred replications were conducted for each of the 144 simulation conditions in R. Frequentist models used the lme4 package (Bates et al., 2015) for models with random effects and nlme package (Pinheiro et al., 2021) for the minimal model, always using REML as the estimation method and Satterthwaite’s approximation to the denominator degrees of freedom (via lmerTest). Bayesian estimation utilised the brms package (Bürkner, 2017), employing the Hamiltonian Monte Carlo (HMC) algorithm with its NUTS extension for efficient sampling of complex models (Hoffman & Gelman, 2014). Bayesian models were fitted using two chains of 1,000 iterations each, with 400 warm-up iterations per chain, yielding a total of 1,200 post–warm-up draws per model. The adapt_delta parameter was set to 0.95 to reduce divergent transitions and enhance sampling reliability. Extracted data included point estimates, standard deviations, p values (frequentist), and posterior means with 95% credible intervals (Bayesian). Information Criteria fit indices (AIC, BIC, WAIC, LOO) and convergence diagnostics were also extracted , with indicating satisfactory convergence). According to this criterion, approximately 99% of the replicas showed appropriate convergence. The convergence dataset is available on our OSF project (see Rodríguez-Prada et al., 2026a).
Data Analysis
Data processing and analysis were performed in R using RStudio. Power and Type I error rates were calculated based on intervention effects. Frequentist decisions used a significance threshold (p < 0.05), while Bayesian ones relied on 95% credible intervals (CIs). Model selection was assessed using AIC, BIC, WAIC, and LOO by comparing each candidate model (minimal, partial intercepts, partial slopes, maximal) to the population model. In addition, analysis of variance (ANOVA) was conducted to examine the influence of simulation parameters on performance indices, and partial eta-squared () was reported as a measure of effect size. Partial eta-squared was used as a relative indicator of the impact of simulation conditions. Effect size magnitudes were not interpreted using conventional cut-off values, as the aim was not to compare absolute effect sizes with previous studies but to examine differences across simulation conditions. Partial eta-squared was selected over eta-squared because it provides an unbiased estimate of the unique variance explained by each factor in multifactorial designs (Richardson, 2011). Consistent with prior methodological simulation studies (e.g., Moeyaert et al., 2017), cut-off points of .01, .06, and .14 were adopted to classify small, medium, and large effects, respectively (Cohen, 1988). This approach allowed us to identify the most influential simulation parameters beyond statistical significance. The dataset, prepared for the dissemination of scientific data, is available on the Open Science Framework at Rodríguez-Prada et al. (2026a). This repository contains the scripts for data processing and analysis, which are also available at Rodríguez-Prada (2025).
Results
Analysis of Different Priors on Power, Type I Error Rate, and Model Selection in Bayesian Methods
Seven Bayesian priors were evaluated for estimating variance components across four population multilevel models under the simulated conditions (Table S1, Rodríguez-Prada et al., 2026b). Minimal differences were observed in statistical power among the priors (, Type I error rates (F(6, 107568) = 0.34, p = 0.916, = < .001), and model selection accuracy (F(6, 428638) = 1.50, p = 0.174; = < .001). Statistical power remained moderate (≈ 0.62), Type I error rates were conservative (≈ 0.033), and WAIC and LOO correctly identified the population model in 82% of cases. The truncated Cauchy prior was selected for further analyses due to slightly higher power and proximity to the nominal Type I error rate of 0.05, reflecting its suitability for this simulation context.
Model Selection Using Information Criteria
Correct model selection rates for frequentist (AIC, BIC) and Bayesian (WAIC, LOO) indices showed minimal overall differences (F(1.64, 117212) = 1557.99, p < .001; = 0.007). However, there was a significant interaction between the index and population model (F(4.93, 117212.53) = 2883.85, p < .001; = .11). Marginal means (Table S2, Rodríguez-Prada et al., 2026b) indicate that WAIC and LOO excel at identifying the maximal model but perform worse with simpler models, selecting them 76%–82% of the time. In contrast, frequentist indices favour simpler models: BIC overwhelmingly selects the minimal model (98%), while AIC demonstrates a more balanced across all population models (Figure 1). Although the interaction effects of the information criteria with both the number of individuals (F(3.28, 117212.53) = 213.47; p < 0.001; = 0.002) and repeated measures (F(4.93, 117212.53) = 88.85; p < 0.001; = 0.001) were small, they reveal noteworthy patterns. More individuals and repeated measures improve correct model selection. AIC and LOO are less sensitive to these factors, while frequentist indices, especially BIC, show greater variability. Overall, Figure 1 shows that BIC and AIC consistently outperform WAIC and LOO in simpler structures (minimal and partial–intercepts models), whereas their performance decreases when slope variability is included in the model. In contrast, WAIC and LOO provide relatively stable, though lower, accuracy across models. Thus, the results suggest that the choice of information criterion has important implications: BIC and AIC tend to favor parsimony, whereas WAIC and LOO yield more balanced but less accurate selections.
Figure 1
Proportion of Correct Selections (Accuracy Rate) for Each Relative Fit Index Based on Population Models
A detailed analysis of errors in BIC, AIC, WAIC, and LOO (i.e., the model selected when the fit index fails to identify the correct one) reveals distinct error patterns. BIC tends to overly penalise complex models and, when incorrect, selects one of the simpler models (minimal or partial-intercepts) in 77.8% of the occasions. In contrast, WAIC and LOO often err by selecting the maximal model nearly half the time. AIC shows the most balanced error distribution, although it slightly prefers the random-intercepts model (Figure 2). When AIC and BIC misidentify the true model, the most frequent error is selecting simpler models (e.g., BIC favours partial-intercepts at 47.8%). In contrast, WAIC and LOO tend to misclassify more complex models, with nearly half of errors involving the maximal specification (46.3% for LOO and 47.9% for WAIC). These results suggest that AIC and BIC are biased toward parsimony, whereas WAIC and LOO are more likely to overfit, highlighting systematic tendencies in model misclassification.
Figure 2
Distribution of Errors: Model Selection When the Information Criteria Fails to Identify the True Model
Power and Type I Error Rate Conditioned on Model Selection
Power and Type I error rates conditioned on model selection were analyzed, simulating scenarios in which the true population model is unknown, as happens in ecological contexts.
On power, ANOVA results showed a significant main effect of the fit index (F(1.27, 60132.31) = 2702.46; p < 0.001; = 0.05), with BIC yielding the highest power, followed by AIC, WAIC and LOO (Table S3, Rodríguez-Prada et al., 2026b). However, for BIC, this increase in power was accompanied by substantially inflated Type I error rates in more complex models. Significant effects were also observed for population model (F(3, 47482) = 11393.98; p < 0.001; = 0.42), effect size of the intervention (F(1, 47482) = 9713.90; p < 0.001; = 0.17), and number of individuals (F(2, 47482) = 2782.18; p < 0.001; = 0.10). Larger effects and sample sizes enhance power, whereas higher error variance in models, as in maximal models, diminishes it (Table S4, Rodríguez-Prada et al., 2026b).
The interaction between the population model and effect size (F(3, 47482) = 2309.91; p < 0.001; = 0.108) showed stable statistical power with large effect sizes across models. However, there were significant declines for moderate effects, particularly in Bayesian methods, reflecting their conservative nature (Figure 3; Table S4, Rodríguez-Prada et al., 2026b). Overall, across all criteria, statistical power is high and stable for simpler models (minimal and partial–intercepts) but declines sharply when slope variability and maximal structures are introduced. The drop is most pronounced for smaller effects (ES = 1.15), where power often falls below .40, particularly for LOO and WAIC. Larger effects (ES = 2.70) mitigate but do not eliminate this decline. These results indicate that model complexity disproportionately reduces power for detecting smaller effects, with Bayesian indices (LOO, WAIC) being more sensitive to this loss than frequentist ones (AIC, BIC).
Figure 3
Effects of the Interaction Between Population Model, Effect Size and Fit Index on Conditioned Power
The model fit index, conditioned on the selected model, significantly impacts the Type I error rate (F(1.49, 35636.95) = 511.52, p < .001, = .02). Although the effect size is small, meaningful differences emerge among indices. Frequentist indices tend to exceed the nominal Type I error rate, while the Bayesian indices demonstrate a more conservative behaviour, staying closer to the nominal value of 0.05 (Table S5, Rodríguez-Prada et al., 2026b). The interaction between population model and model fit indices is statistically significant but small (F(4.47, 35636.95) = 165.28, p < .001, = .02). Examining marginal means (Table S5, Rodríguez-Prada et al., 2026b), the BIC index shows an unacceptable Type I error rate for complex models, reaching 13% for the maximal model. The AIC index is less sensitive than BIC, with a Type I error rate of 8.9% for the maximal model. WAIC and LOO remain stable across all population models, aligning closely with the nominal level in complex models (Figure 4). Thus, error rates remain close to the nominal .05 level for minimal and partial–intercepts models across all criteria. However, as complexity increases, AIC and BIC show inflated error rates, particularly under maximal specifications (AIC ≈ .13, BIC ≈ .09). By contrast, Bayesian fit indices (WAIC and LOO) maintain more stable control of Type I error, even in complex models. These findings suggest that AIC and BIC may increase the risk of false positives when applied to highly parameterized structures, whereas WAIC and LOO provide more conservative performance.
Figure 4
Type I Error Rate Conditioned on Model Selection as a Function of the Population Model
In summary, BIC exhibited the highest values for both statistical power and the Type I error rate. WAIC and LOO yielded the lowest values, remaining close to but not exceeding the nominal 0.05 level for false positives. AIC demonstrated a balanced performance between statistical power and the Type I error rate. Interaction effects highlighted the influence of between-subject factors on both power and the Type I error rate, particularly the intervention effect size, sample size, and population model for power; and only the interaction with the population model for the Type I error rate.
Discussion
Statistical modelling has historically been underutilised in SCEDs within clinical and psychological contexts due to inherent challenges including small sample sizes, limited observations, and reduced data availability. However, advances in methods such as MLMs and Bayesian approaches offer promising solutions (Baek & Ferron, 2013; Moeyaert et al., 2017). These methods are particularly valuable for estimating covariance parameters, enabling researchers to quantify between-individual subject variability and explore its causes (e.g., why some individuals benefit more from an intervention than others). At a descriptive level, multilevel models enable researchers to quantify intervention effects, estimate the magnitude and direction of changes in both levels and slopes, and assess their practical significance (Ferron et al., 2009). Several studies have focused on exploring estimation in terms of bias and standard errors (Baek et al., 2020; Ferron et al., 2010; Moeyaert et al., 2017). However, there has been limited focus on other statistical properties of MLMs. A key objective of this study was to compare the performance of model fit information criteria from two different frameworks (frequentist via REML estimation and Bayesian with varying levels of weakly-informativeness) in a context where the true population model is unknown.1 Bayesian methods showed no significant differences in statistical power, Type I error rate, and model selection among the chosen priors, confirming previous studies (Baek et al., 2020; Moeyaert et al., 2017). The half-Cauchy and half-normal priors, which favor smaller variances, performed slightly better, as suggested in earlier research. However, prior selection is context-dependent, with effectiveness varying by intervention type, outcome characteristics, and research design. Commonly used priors, like half-normal or half-Cauchy distributions, are often based on simulations and may not suit empirical cases where the dependent variable has greater dispersion. Priors should be tailored to the outcome's scale and variability for robust estimation. A stronger integration between simulation and applied studies is needed to improve these decisions. As highlighted by Moeyaert et al. (2017), prior calibration can be informed by reanalyses of empirical SCEDs which helps to ensure that the scale of the priors reflects realistic ranges for variance components. For applied researchers, however, translating prior knowledge into formal prior distributions remains challenging. More empirical work and meta-analytical evidence are needed to provide practical guidance on how to derive informative priors in SCED contexts.
The information criteria (AIC, BIC, WAIC, LOO) showed minor differences in their ability to identify the true data-generating model. This finding partially aligns with those of Moeyaert et al. (2017), which reported equivalence between REML and Bayesian methods for estimating fixed effects, as both frameworks produced similar results. However, interaction effects revealed nuanced behaviours: Bayesian indices tended to favour complex models, while frequentist indices were more accurate with simpler models. This pattern suggests that Bayesian methods may lean toward overparameterisation, whereas frequentist methods might prefer underparameterized models. More recent contributions have extended the scope of MLM research in SCEDs by exploring novel approaches, including the application of generalized linear models (GLMs) and refinements in variance component estimation. Li et al. (2024) recommends using AIC and BIC to select an optimal model in case of count data with overdispersion; but if there is overdispersion and zero-inflation, the recommendation is to use methods with lower penalty for these complex models. However, no Bayesian option was explored in Li et al. (2024). Recent simulation work has further emphasized the importance of model specification for variance components in SCEDs. Li et al. (2022) demonstrated that biased estimates of between-case variance are particularly problematic when the true variance is small, and that unconstrained optimization methods combined with post hoc model selection procedures (e.g., bootstrap-based RLRT) can improve estimation accuracy and inference. These findings highlight that not only the choice between frequentist and Bayesian approaches but also the technical details of variance component estimation and covariance structure specification critically affect the robustness of MLMs applied to SCEDs.
It is worth mentioning that there are some differences between fixed and random effects regarding ecological validity and generalization of models in SCEDs depending on their complexity. On the one hand, simpler models may have lower ecological validity despite their greater generalizability. At the same time, overparameterised models may lead to problems of model fit and reduced predictive validity, particularly when unnecessary fixed effects are added. In the field of SCEDs, including additional fixed effects can enrich the description of treatment effects, and adding more random effects is often crucial for adequately modelling change. But this increases the estimation challenges. Our findings align with previous literature suggesting that Bayesian methods may offer a valuable alternative to handle this kind of complex model (Baek et al., 2020; Moeyaert et al., 2017; Van den Noortgate & Onghena, 2003b). This may be particularly true considering that real SCEDs are generally more complex than those simulated in methodological studies, making the Bayesian framework potentially more suitable for real-world data. In contrast, frequentist approaches may be more appropriate for simpler scenarios as, for example, frequentist methods struggle to estimate covariance parameters in complex models (Moeyaert et al., 2017). AIC performs well across models, making it suitable for straightforward scenarios. In contrast, WAIC and LOO excel in complex models, being less affected by sample size and repeated measures, limitations often faced in SCEDs. BIC’s performance varies and is less reliable with complex models. Each information criterion shows distinct error patterns: WAIC and LOO favor complex models, whereas AIC and BIC lean towards simpler ones, consistent with previous findings (Martínez-Huertas et al., 2022; Martínez-Huertas & Olmos 2022). This divergence may reflect frequentist methods’ limitations in estimating random effects in complex settings, highlighting the importance of understanding the unique strengths of each index in different modelling contexts.
The analysis of power and Type I error rates was focused on emulating real-world conditions researchers face, where the true data-generation model is unknown. When the intervention effect size was large, no significant differences in statistical power were found between Bayesian and frequentist methods, indicating that either can offer reasonable conclusions. However, frequentist methods, especially BIC, exhibited higher Type I error rates in complex scenarios. Higher power driven by inflated Type I error is not desirable. AIC showed better performance with a moderate excess of Type I error (8.9%) and a good balance with statistical power. Bayesian methods, while slightly conservative (Error rate ≈ 0.04), offered more stable and consistent results across conditions. Therefore, despite limitations and considering all these findings, Bayesian methods may be more robust in complex SCED settings than frequentist methods.
Model specification significantly influences statistical performance as it affects the estimation of confidence intervals, standard errors, and effect sizes. For instance, when standardising mean differences, using biased standard errors can result in biased effect size estimates. Model detectability varies based on complexity, information criteria, and estimation methods. This study analysed four models of increasing complexity, ranging from one with no random effects to one with random intercepts and random slopes. In this simulation, simpler models, which assumed no individual variability in intervention effects (an unrealistic assumption), yield lower standard errors from residual variance, enhancing detectability but overlook stochastic dependencies in nested data, assuming independence (Moerbeek, 2004). While fitting well, they may not represent real scenarios. Additional variability from random slopes heightens standard errors of the intervention effect, complicating detection. Increased error variance in complex models might explain lower performance, with detectability improving only with larger intervention effect sizes. Bayesian methods excel in these scenarios, at least using weakly informative priors and fit indices like WAIC or LOO. In fact, models can become more complex with relevant effects like temporal trends (Baek et al., 2020; Moeyaert et al., 2017), making them more intricate than those examined here. The MLMs simulated in this study were simpler, enabling controlled exploration of fit indices’ performance. Unreported simulations revealed estimating the correlation between slopes and intercepts often yielded impossible values or failed. This simplicity may clarify the stronger performance of frequentist methods over Bayesian approaches in some conditions. However, as complexity increases, Bayesian methods are expected to outperform frequentist methods, especially with small sample sizes and few repeated measures.
Bayesian methods have distinct advantages for SCEDs, especially with complex random effects. They have been successfully applied in situations with heterogeneous Level-1 variances (Baek & Ferron, 2013) and autocorrelated error terms with trends (Natesan, 2019), which are challenging for frequentist methods. Bayesian approaches’ flexibility makes them ideal for these complexities. Software like the brms package (Bürkner, 2017) enhances accessibility, enabling researchers to apply Bayesian methods without deep knowledge of computation. Other tools like JAGS, WinBUGS, STAN, or rstan further support users without advanced statistics expertise. The findings emphasise matching estimation methods to model complexity. AIC is suitable for simpler models due to its stability, while Bayesian methods are better for complex models or small samples. In any case, we recommend researchers to perform sensitivity analyses with multiple fit indices to ensure consistent model selection, enhancing their findings' robustness.
Limitations
The present simulation study evaluated the operating characteristics of Bayesian and frequentist methods in MLMs applied to SCEDs when the true population model is unknown. Findings are specific to simulated conditions and need empirical validation. Our findings should be interpreted conditional on the simulated data-generating mechanisms: an AB design with normally distributed outcomes and no explicit autocorrelation or phase-specific trends, nor cross-level interactions. In applied SCEDs, however, outcomes are often counts and time-related structure is common. These features can change the effective model complexity and the penalty-fit trade-off, potentially altering the relative behavior of AIC/BIC versus WAIC/LOO and frequentist versus Bayesian estimations. This is consistent with recent work extending multilevel SCED models to GLMM settings and more realistic covariance structures (e.g., Li et al., 2024, 2025a, 2025b), and motivates future simulations that jointly manipulate outcome distribution, trends, and autocorrelation to evaluate whether the performance patterns reported here generalize to these more ecological scenarios. These decisions limit the scope of our findings, since model complexity of MLM applications often arises precisely from elements such as trends (general or phase-specific) or autocorrelation (Natesan Batley & Hedges, 2021). We deliberately opted for a simple design to keep the conditions more controlled and to examine the performance of model selection criteria in a focused way, particularly regarding random effects and the detection of random slopes. Including both baseline and intervention trends would have substantially increased the complexity of the simulation, making it more difficult to isolate the behaviour of the information. Future research should explore the behavior of information criteria on model selection in more complex designs (e.g., multiple baseline or reversal designs) and advanced covariance structures, such as autoregressive terms, to better reflect real-world data. Using dependent variables that follow normal distributions, uncommon in SCEDs where ceiling and floor effects and heteroscedasticity tend to feature, also limits generalizability. Future studies should investigate the behavior of both frequentist and Bayesian approaches within generalized linear multilevel models (GLMMs). By incorporating variables with asymmetric distributions, such as Poisson distributions, researchers can address this gap and explore small effect sizes that are frequently overlooked in qualitative analyses. Recent research has demonstrated that linear mixed models can be adapted to handle count outcomes in SCEDs, leading to more accurate estimates, particularly in cases of overdispersion or small sample sizes (Declercq et al., 2019). More recently, Li and colleagues have shown how GLMMs can be effectively applied to various count data scenarios, including zero-inflated and overdispersed distributions (Li, 2024; Li et al., 2025b; Li et al., 2024). They have also provided step-by-step tutorials designed for applied researchers (see Li et al., 2025a). These contributions emphasize the importance of moving beyond normally distributed outcomes to enhance the ecological validity and robustness of statistical inferences in SCED data.
Regarding the Bayesian fit indices, WAIC and LOO, both derived from the conditional likelihood, were chosen for their accessibility and compatibility with the evaluated models. However, future work could benefit from a broader exploration of Bayesian model comparison strategies. Greater emphasis should be placed on marginal rather than conditional likelihood-based approaches, which some authors argue are more suitable for model comparison in Bayesian frameworks (Ariyo et al., 2022; Merkle et al., 2019). Additionally, recent methodological advances demonstrate the growing applicability of Bayes factors within different SCEDs, such as ABAB, alternating treatments and changing criterion designs (Yamada & Okada, 2024, 2025). Incorporating Bayes factors in future analyses could offer complementary insights into model selection decisions. The best results of these information criteria were observed under optimal conditions: more individuals, repeated measures, and larger effect sizes, consistent with previous research (Baek et al., 2020; Moeyaert et al., 2017). While MLMs effectively capture complex realities, their mathematical and statistical demands underscore the importance of clear model specification. Bickel (2007) observes that these models perform best when grounded in robust theories and literature. Thus, advancing MLMs should also promote theoretical development, particularly within psychology. It is important to remember that statistical decisions should never be made without considering the wider context and the insights provided by clinical expertise. An analysed effect might not reach the traditional thresholds of statistical significance but can still be socially relevant, aligning with the concept of social validity (Kazdin, 1977; Snodgrass et al., 2023), and ways of control the antecession of the intervention before the change is fundamental to reassure the efficacy of an intervention (Perone, 1999). Future research should investigate the contingency relations between the decisions proposed by these models and those made by applied researchers, to understand how they are related.
Conclusions and Recommendations
MLM is a particularly versatile option for analysing SCED data, offering general and individual effect quantification, statistical testing, and the ability to model complex effects. This study highlights the strengths of Bayesian methods for highly complex models, where they maintain stable power and Type I error rates across varying effect sizes, sample sizes, and repeated measures. For simpler models, the frequentist REML approach performs equally well or better, particularly when AIC is used for model selection. When protection against Type I errors is paramount, as in assessing therapeutic change, Bayesian methods provide a reliable and robust framework, making them a valuable tool for applied researchers in SCEDs. Given that the data-generating model is typically unknown in applied SCEDs, the practical contribution of this study is to characterise the pros and cons of common model selection fit indices. Under our simulation conditions, REML estimations with AIC offered the most favourable balance between power and Type I error rates in simpler population models. But Bayesian estimations provided a more stable positive performance when model complexity increased. Thus, our recommendation is to treat model selection as risk management under uncertainty, triangulating across plausible random-effects structures, and using sensitivity analyses across AIC/BIC and WAIC/LOO to ensure that conclusions are not criterion-dependent.