Multilevel Modeling for Causal Inference in Oral Health Research

Abstract

Multilevel modeling (MLM) is increasingly applied in oral health research to reflect the hierarchical structure of dental data, in which teeth are nested within individuals and individuals are embedded within social and environmental contexts. By partitioning both within- and between-individual variability, MLM improves estimate precision, reduces bias, and strengthens causal interpretation, thereby enhancing the evidence base for preventive strategies, clinical practice, and public health policy. This review synthesizes recent methodological developments and ongoing challenges in the application of MLM to oral health research, with particular emphasis on causal inference. Key issues include reciprocal relationships between individual- and contextual-level factors, establishing temporality, collinearity among contextual variables, population mobility, and the selection of appropriate geographic or institutional levels. We also highlight methodological advances, such as integrating propensity score weighting within multilevel frameworks, which offer promising approaches to improving causal assessment in observational studies. Accumulating evidence demonstrates the importance of contextual determinants, such as school and neighborhood environments, in shaping oral health outcomes across the life course. However, causal inference regarding these distal exposures remains underdeveloped. Persistent threats to validity include measurement error in contextual indicators, unmeasured multilevel confounding, and limited consideration of cross-level interactions. Using a school-based example, we illustrate how food insecurity mediates the association between household income and active dental caries, highlighting the value of analytical approaches that jointly address individual- and contextual-level determinants. Overall, MLM provides a robust framework for clarifying causal pathways, improving analytical rigor, and guiding interventions that target both behavioral and structural drivers of oral health inequalities.

Keywords

epidemiology statistics dental public health social determinants inverse probability weighting mixed-effects modeling

Introduction

Multilevel modeling (MLM), also known as hierarchical linear or mixed-effects modeling, has transformed the analysis of complex, nested data. Its origins trace back to the mid-20th century, when researchers in education and sociology recognized a fundamental challenge: observations were rarely independent. Students were nested within classrooms, patients within hospitals, repeated measures within individuals. Ignoring these hierarchies often led to underestimated variability and biased inferences using traditional regression methods (Snijders and Bosker 2012; Hox et al 2017). Groundbreaking work by Goldstein and colleagues in the 1970s and 1980s formalized variance components models, enabling the partitioning of variability across multiple levels (Bryk and Raudenbush 1992; Goldstein 1995). Raudenbush and Bryk (2002) further advanced the field by developing accessible methods and software that captured group-level differences through random intercepts and slopes. By the 1990s, MLM had evolved into a flexible framework capable of accommodating longitudinal data, repeated measures, generalized outcomes, and complex covariance structures, enabling the simultaneous estimation of fixed and random effects while distinguishing within- and between-cluster variability (Verbeke and Molenberghs 2000). The development of user-friendly software—including HLM, MLwiN, and R packages such as lme4 and nlme—broadened its adoption across disciplines, cementing MLM as a cornerstone of modern statistical analysis (Bates et al 2015; Charlton et al 2025).

This review synthesizes recent methodological developments and ongoing challenges in the application of MLM to oral health research, with particular emphasis on causal inference. Using a school-based example, we illustrate the use of MLM for causal inference by integrating multilevel propensity score (PS) methods into hierarchical frameworks.

MLM in Epidemiology and Oral Health

In epidemiology, MLM has become indispensable for analyzing health outcomes influenced by factors at multiple levels—from individual behaviors to household, community, and policy contexts. By accounting for clustering (eg, patients within clinics or individuals within neighborhoods), MLM disentangles contextual from compositional effects, providing critical insight into social determinants of health and inequalities (Diez-Roux 2000). Context refers to the natural or built environmental characteristics surrounding individuals, encompassing physical, political, social, cultural, and economic features, as well as their interactions with individuals. A contextual effect is defined as a change in health status caused by a factor operating at a higher, contextual level (Diez Roux 2002).

In oral health research, MLM has been particularly powerful. Outcomes such as dental caries, periodontitis, and tooth loss result from a complex interplay of individual behaviors (eg, oral hygiene, diet), biological factors (eg, salivary composition, genetics), and environmental determinants (eg, fluoridation, socioeconomic context) (Sabbah et al 2007; Listl et al 2015; Peres et al 2019). MLM has quantified the independent impact of neighborhood deprivation on caries prevalence (Antunes et al 2006) and evaluated school- or clinic-level interventions while accounting for intracluster correlation (ICC) (Tadakamadla et al 2022). Longitudinal MLM further enables examination of oral health trajectories across the life course, capturing within-person changes over time while acknowledging higher-level influences (Celeste et al 2020).

Three main considerations justify MLM in oral health research. First, oral health outcomes are inherently clustered (eg, teeth within individuals), violating the independence assumption of standard single-level analyses. Clustering is quantified via ICC for continuous variables or variance partitioning coefficients for categorical outcomes, representing the percentage of total variability at the cluster level (Goldstein et al 2002). Even small ICCs (1%–5%) can bias associations due to ecological fallacy (Diez-Roux 1998). Ignoring such dependencies leads to underestimated standard errors, overly narrow confidence intervals (CIs), and unreliable P values. Second, MLM allows the explicit inclusion of meaningful higher-level predictors (eg, school-level factors) rather than merely adjusting for clustering, offering deeper insights into contextual effects (Snijders and Bosker 2012). Third, MLM naturally handles unequal cluster sizes—common in oral health data (eg, varying numbers of teeth per child)—by providing appropriate standard errors and minimizing risks of under- or overestimation that can occur in simpler adjustments (Maas and Hox 2005).

This multilevel perspective is further supported by ecological models that recognize the broader determinants of social, behavioral, clinical, and immune phenotypes. As illustrated in Figure 1, oral health outcomes are shaped not only by individual behaviors and biological factors but also by the layered influences of family, community, societal structures, and the environment. For instance, individual choices such as smoking, sugar intake, and use of fluoridated toothpaste interact with social network influences—such as family norms and access to health education—and broader societal elements such as public health campaigns, health care policies, and economic conditions. In addition, the human genome and microbiome contribute to the transmission and acquisition of clinical and immune phenotypes, reinforcing the need for a systems-level approach with multiple levels to understanding oral health inequalities and outcomes across the life course. MLM, by incorporating random intercepts and slopes, cross-level interactions, and time-varying covariates, provides a robust framework to examine these interactions, identify inequalities, and guide interventions.

Figure 1.

Multilevel factors shaping clinical and immune phenotypes.

Critical Methodological Considerations in MLM in Oral Health Research

After defining the hierarchy of multilevel data to disentangle within-cluster and between-cluster effects, sample size planning is critical. Limited higher-level units reduce the power and reliability of variance components and random effects estimates, whereas moderate-to-large cluster sizes aid within-cluster estimation (Hox et al 2017). Intracluster correlation (ICC), average cluster size at each level, and cluster size variation must be considered alongside the usual parameters in sample size determination. ICC quantifies the clustering effect by measuring the proportion of total outcome variance between clusters (Hox et al 2017). A high ICC for “teeth within patients” implies that knowing the outcome of one tooth from the same patient strongly predicts the outcomes for other teeth. High ICC requires more (or smaller) clusters to maintain statistical power. The average cluster size is the typical number of observations per cluster and inflates the total required sample. For 2-level models, the compute design effect is $1 + (m - 1) I C C$ (where m is the average cluster size), which inflates the single-level sample size (Snijders and Bosker 2012). For 3- or higher-level models, simulation- or Monte Carlo–based power analysis is often needed (Snijders and Bosker 2012).

Missing data (eg, missing measurements or patient dropouts) can bias results. Multilevel models commonly assume missing at random, but the pattern (within vs between clusters) matters, and missing not at random is plausible if linked to unobserved disease progression. Appropriate handling requires multilevel multiple imputation, maximum likelihood under the assumed missing-data mechanism, or joint modeling, which is complex in multilevel contexts (Enders et al 2016; Grund et al 2016).

Model building should start simple, adding complexity only when theory and data justify it (Cheung et al 2019). Variance components models (the simplest form) include only intercepts and random terms to partition total variance across levels (Snijders and Bosker 2012). They assess variance at each level and decide which levels should be retained for analysis. Nonsignificant levels can be dropped, although retention may be warranted for data structure or comparability (Wan et al 2009). Explanatory variables (predictors, confounders) that may associate with or influence the outcome can be added next. Random intercept models allow cluster-varying baselines while keeping the effects of explanatory variables (fixed effects) constant across clusters (Albandar and Goldstein 1992). Random slope models extend this by allowing effects of independent variables to vary between clusters, yielding random slopes (Snijders and Bosker 2012) but demand larger samples and more computation.

Linear mixed models are widely used to model continuous, approximately normally distributed outcomes (Bates et al 2015). Generalized linear mixed models handle nonnormal outcomes via suitable link functions and distributions (Hox et al 2017). Logistic or probit links are appropriate for binary outcomes such as the presence/absence of arrested caries lesions (Mabangkhru et al 2020). Multilevel Poisson or negative binomial models, with or without zero-inflation to handle excess zeros, suit count outcomes such as the number of caries lesions (Li et al 2021; Piovesan and Bernabe 2025), and multilevel survival models address time-to-event outcomes, such as the longevity of restorations (Wierichs et al 2025). When studying the progression of oral diseases, trajectories may be nonlinear over time; growth-curve models or nonlinear mixed models can be used (Snijders and Bosker 2012; Hox et al 2017).

Depending on model complexity, either frequentist or Bayesian approaches for parameter estimation can be used (Snijders and Bosker 2012; Hox et al 2017). Frequentist approaches (eg, lme4, glmer in R; PROC GLIMMIX in SAS) are efficient but falter with complex random-effects structures (Bates et al 2015). Bayesian approaches (eg, brms, rstanarm in R; MCMC in MLwiN) manage small samples and complexity via priors, requiring careful specification and more computation (Bürkner 2017; Muth et al 2018). Many estimation procedures are iterative, such as iterative generalized least squares, which alternates coefficient and covariance estimation until convergence (Goldstein 1986). Nonconvergence or unstable estimates may signal an overly complex or poorly identified model, suggesting the need for simplification (Hedeker and Gibbons 2006). Convergence diagnostics and checks for identifiability problems are therefore essential.

When models include continuous predictors, centering (subtracting either the grand mean or the group mean) is an important consideration. Group-mean (within-group) centering of level 1 predictors isolates within-group effects and separates them from between-group effects, facilitating interpretation of within-cluster associations. When between-group effects are of primary interest, grand mean centering is often preferable (Enders and Tofighi 2007). Explicit reporting of centering choices is needed to avoid misinterpretation.

Postfitting diagnostics include residual plots at each level, assessment of the distribution of random effects, and influence diagnostics to identify outlying or influential clusters or units (Loy and Hofmann 2014; Cho et al 2022). When multiple models are considered, fit can be compared using Akaike information criterion, Bayesian information criterion, or likelihood ratio tests for nested models (Burnham and Anderson 2004). More advanced issues, including cross-classified and multiple membership structures, cross-level interactions, longitudinal data with time-varying effects, and other extensions, go beyond the scope of this overview but can be found elsewhere (Snijders and Bosker 2012; Hox et al 2017).

The Gap in Causal Inference Using MLM in Oral Health Research

Over the past 2 decades, the usage of MLM has surged, with more than 100 studies in 2024 leveraging this approach to analyze clustered data in oral health (Fig 2). However, the application of MLM for causal inference—estimating the true causal impact of exposures on oral health outcomes—remains limited. In 2024, only approximately 20 studies used MLM for causal analyses despite the prevalence of hierarchical data in this field (Fig 2). Notably, while the overall use of MLM in oral health research has increased markedly, its application for causal inference has grown much more slowly. The proportion of MLM studies explicitly addressing causal relationships has declined by an average of 4% annually (average annual percentage change = −4% [95% CI, −7.6% to −0.3%]), highlighting a critical gap in methodological advancement for robust causal inference in the field.

Figure 2.

Temporal trends in multilevel modeling (MLM) and causal inference applications in oral health research, 1990 to 2024. Left panel: Annual number of publications using MLM (red line) compared with those applying MLM for causal inference (blue line). Right panel: Proportion of MLM publications explicitly addressing causal inference, with the fitted trend line showing the average annual percentage change of −4% (95% confidence interval, −7.6% to −0.3%).

Causal inference addresses the fundamental challenge of missing counterfactual outcomes—what would have happened to an individual under a different exposure (Rosenbaum and Rubin 1983). In public health, where randomized controlled trials are often infeasible due to ethical or practical constraints, PS methods have become increasingly vital for causal inference. To emulate the balance of randomization, traditional regression methods are based on assessing the association of a potential confounding factor with the outcome. Another class of adjustment is based on assessing the association of potential confounding with the main exposure, rebalancing the exposure groups based on weights or scoring methods (Haneuse 2021). The PS, defined as the conditional probability of receiving a treatment or exposure given baseline characteristics, balances covariate distributions across exposure groups, with the assumption of no unmeasured confounders and positivity (Rosenbaum and Rubin 1983; Austin 2011). Initially developed for baseline treatments, PS methods have evolved to accommodate time-varying exposures and confounders, enhancing their applicability in longitudinal settings (Austin 2020).

Common PS-based approaches include PS matching, PS stratification, and inverse probability weighting (IPW). These methods are typically implemented after PS estimation and balance assessment (eg, standardized mean differences <0.1), and they have been increasingly applied in dental research for estimating the causal effects of interventions, risk factors, and care utilization patterns (Burgette et al 2016; Shetty et al 2016). Among these, IPW stands out for preserving the full sample while directly targeting the population average treatment effect. In contrast, PS matching often leads to substantial sample loss due to incomplete or poor-quality matches, potentially reducing precision and generalizability. PS stratification retains most of the sample by dividing data into strata based on PS quantiles but generally achieves weaker covariate balance than IPW or PS matching, especially with few strata. Sensitivity analyses remain essential to evaluating robustness against violations of key assumptions including no unmeasured confounders and positivity.

Causal inference also encompasses mediation analysis to decompose mechanistic pathways. In multilevel settings, mediation analysis—often via traditional MLM or multilevel structural equation modeling—allows separation of within- and between-cluster indirect effects but faces challenges such as level conflation, confounding across levels, and power limitations with small higher-level samples (Zhang et al 2009; Preacher et al 2010; Preacher 2011). While promising for elucidating hierarchical pathways in oral health, multilevel mediation remains underused in the field.

While PS techniques have advanced considerably, their extension to multilevel structures is still underdeveloped in oral health research (Chang and Stuart 2022), even though applications in clustered data have been explored in other fields (Su and Cortina 2009; Lee et al 2021). Multilevel models incorporating fixed or random effects at the cluster level reduce bias in PS estimation and treatment effect inference compared with single-level approaches (Arpino and Mealli 2011; Thoemmes and West 2011; Fuentes et al 2022). Ignoring clustering can lead to underestimated variance and unreliable causal inferences (Cafri et al 2019).

Challenges of Making Causal Inferences in Multilevel Studies on Contextual Determinants of Oral Health

Contextual variables are traditionally classified according to their measurement (Morgenstern 1982, 1995). Aggregate variables are means or proportions derived from individual-level measures, such as neighborhood mean income or the percentage of residents below a defined income threshold. Global variables, in contrast, are measured directly at the group level that cannot be reduced to individual units—such as community water fluoridation in a city. Health policies are a typical example, as they are defined within specific political or administrative boundaries (eg, cities, states, or countries).

Even with careful classification, causal inference about contextual determinants in multilevel oral health studies remains challenging. A key difficulty is distinguishing true contextual effects from compositional confounding—where group-level associations primarily reflect aggregated individual characteristics—and from unmeasured confounding across levels, including residential self-selection and time-varying exposures (Diez Roux 2004). Endogeneity frequently occurs because individuals sort nonrandomly into contexts based on socioeconomic position, behaviors, or health status, violating core exchangeability assumptions required for valid causal claims. Aggregate variables (eg, neighborhood mean income or school-level oral health knowledge) are particularly vulnerable to misinterpretation as independent contextual influences. Global variables offer stronger potential for causal interpretation but still face issues such as spillover effects, boundary interference, limited between-unit variation, and small numbers of higher-level clusters, resulting in reduced statistical power and potential bias in variance component estimates.

These challenges are especially pronounced in oral health research due to the hierarchical and clustered nature of the data (eg, teeth within individuals), unequal cluster sizes, and complex mediating pathways involving biological, behavioral, and psychosocial factors (Newton and Bower 2005). Robust causal inference therefore requires explicit frameworks (eg, potential outcomes approach), rigorous sensitivity analyses for unmeasured confounding, and advanced methods such as marginal structural models, g-computation, instrumental variables, or natural experiments when observational data permit. These approaches are essential to move beyond associational evidence toward policy-relevant causal conclusions in dental public health.

Empirical Application: Using MLM for Causal Inference

This illustration aimed to promote causal inference in oral health research by integrating multilevel PS methods (eg, IPW with random intercept models) into hierarchical modeling frameworks. This is particularly important given that conventional MLM alone does not establish causal relationships between exposures and outcomes. Mediation analysis is not the primary methodological focus of this study. Although a 2-level model was demonstrated here (as individual-level outcomes are more common than tooth-level analyses in most epidemiological and intervention studies in this field), the proposed approach is readily generalizable to 3- or higher-level hierarchical structures (eg, teeth within individuals within communities or schools). This allows researchers to address more complex hierarchical structures common in oral health data while preserving robust causal estimation.

We illustrated the application of PS weighting in multilevel analysis using a school-based cross-sectional study conducted in Monteiro, a small town in northeast Brazil in 2023. This study recruited 463 schoolchildren aged 12 years from 14 schools to evaluate whether food insecurity (moderate and severe insecurity vs security and mild insecurity) (Segall-Corrêa et al 2014; Reichenheim et al 2016) mediates the effect of equivalized household income (dichotomized at median income) on the number of surfaces with active caries (a count outcome), as detailed in Appendix Figure S1. In this study, data were clustered at the school level with sex and total number of tooth surfaces as level 1 confounders and school type as a level 2 confounder. Sensitivity analyses were performed using (1) continuous log10-transformed equivalized family income and (2) an alternative dichotomization of food insecurity (any food insecurity vs food security) to assess the robustness of findings to information loss from dichotomization.

Model Selection and Analysis

Given evidence of overdispersion in the count of surfaces with active caries (dispersion parameter = 1.3; deviance-based likelihood ratio test of negative binomial vs Poisson: P < 0.001), we fitted a negative binomial regression model (Appendix Table S1). To account for clustering, variance components analysis of this 2-level random-intercept negative binomial model (students within schools) indicated that approximately 19% of the variability in active caries surfaces was attributable to school-level differences, supporting the multilevel approach (model 2 in Appendix Table S1).

PS Weighting

PSs were generated to estimate stabilized IPWs for exposure (dichotomized equivalized household income) and mediator (dichotomized food insecurity) for student i from school j:

Exposure of student i from school j:

I P W_{i j}^{E} = \frac{P (E = e_{i j})}{P (E = e_{i j} | X = x_{i j})}

Mediator of student i from school j:

I P W_{i j}^{M} = \frac{P (M = m_{i j} | E = e_{i j})}{P (M = m_{i j} | E = e_{i j}, X = x_{i j})}

The final $I P W_{i j} = I P W_{i j}^{E} \times I P W_{i j}^{M}$ .

Here, $e_{i j}$ and $m_{i j}$ were the dichotomized equivalized household income and dichotomized food insecurity status for student i from school j, respectively. $x_{i j}$ were confounders depending on the different logistic regression models (Arpino and Mealli 2011).

M0: Single-level logistic regression with level 1 (sex, total number of tooth surfaces) and level 2 (school type) confounders.

M1: Single-level logistic regression with level 1 confounders only (sex, total number of tooth surfaces).

M2: Two-level random-intercept logistic regression with level 1 confounders and school-specific intercepts.

M3: Two-level fixed-effects logistic regression with level 1 confounders and school dummy variables for constant school-specific effects.

Here, the total number of tooth surfaces refers to the total count of erupted tooth surfaces in the permanent dentition, which varies across individuals due to differences in tooth eruption and accounts for differential opportunity for caries development. This was included to adjust for the varying at-risk surface area, separated from the outcome of active caries on those surfaces (Appendix Methods S1).

Students were weighted by IPW to generate a pseudo-population with balanced confounders. Figure 3 shows that level 1 confounders (sex, total number of tooth surfaces) were balanced (standardized mean difference < 0.1) in the IPW standardized pseudo-populations, but level 2 confounder (school type) remained unbalanced. This imbalance stemmed from limited within-school variation, as private schools were more accessible to higher-income families, a contextual effect that IPW could not fully adjust. Therefore, school type’s influence on active caries was included in the outcome regression as a specific contextual effect when assessing the effect of household income.

Figure 3.

Love plot: covariate balance before and after standardization using inverse probability weighting.

Causal Effect Estimates

Compared to the analysis without specific contextual effect (Table 1), Table 2 shows stable estimates of the total effect (TE), controlled direct effect (CDE) of income on active carries, proportion mediated by food insecurity using a single-level outcome model, and random-intercept outcome model with specific contextual effect of school type, across different IPW-balanced pseudo-populations (generated using M0–M4). Compared with the single-level outcome model, the random-intercept model provided more robust estimates of TE, CDE, and specific contextual effect with a wider 95% CI regardless of the exposure model for IPW generation. Schoolchildren from lower income families had ~40% higher prevalence of tooth surfaces with active caries, with at least 10% mediated by food insecurity. The prevalence of surfaces with active caries among schoolchildren from public schools was almost triple compared with those from private schools. When the exposure (family income) and mediator (moderate to severe food insecurity) interaction was considered, the random intercept model showed that minimal differences in TE and natural direct effect regardless of methods for IPW generation. However, the proportion mediated was much lower when both cluster and contextual factors were not considered for IPW generation (M1). Conclusions remained robust by sensitivity analysis using food insecurity at different cutoff (any food insecurity vs food security) as a mediator (Appendix Table S2) or log10 transformed equivalized family income as exposure (Appendix Table S3).

Table 1.

Average Treatment Effect of Low Equivalized Household Income via IPW-Weighted Regression without Specific Contextual Effect (Dichotomized Food Insecurity as Mediator: Moderate to Severe Food Insecurity vs Security or Mild Insecurity).

IPW	TE (95% CI)	CDE/NDE (95% CI)	Proportion Mediated
Single-level outcome model without EM interaction
Traditional	1.82 (1.38, 2.40)	1.72 (1.30, 2.29)	12.2%
M0	1.64 (1.20; 2.24)	1.56 (1.12; 2.16)	12.5%
M1	1.86 (1.37; 2.52)	1.76 (1.28; 2.43)	11.6%
M2	1.84 (1.35; 2.50)	1.74 (1.26; 2.40)	11.9%
M3	1.57 (1.11; 2.22)	1.48 (1.05; 2.09)	15.8%
Random intercept outcome model without EM interaction
Traditional	1.51 (1.13, 2.03)	1.46 (1.08, 1.96)	9.8%
M0	1.52 (1.14; 2.02)	1.45 (1.08; 1.96)	13.5%
M1	1.55 (1.16; 2.08)	1.49 (1.11; 2.01)	10.9%
M2	1.54 (1.15; 2.06)	1.49 (1.10; 2.00)	9.3%
M3	1.50 (1.13; 2.00)	1.42 (1.06; 1.91)	16.0%

CDE, controlled direct effect; CI, confidence interval; EM, exposure and mediator; IPW, inverse probability weight; M0, single-level logistic regression with level 1 (sex, total number of surfaces) and level 2 (school type) confounders; M1, single-level logistic regression with level 1 confounders only; M2, 2-level random-intercept logistic regression with level 1 confounders and school-specific intercept; M3, 2-level fixed-effects logistic regression with level 1 confounders and school dummy variables for constant school-specific effects; TE, total effect; NDE, natural direct effect; NIE, natural indirect effect.

CDE, NDE when there is no EM interaction, only NDE is reported when there is an EM interaction. All effect estimates are expressed in terms of prevalence ratios. Proportion mediated = (NDE × (NIE – 1))/(NDE × NIE – 1).

Table 2.

Average Treatment Effect of Low Equivalized Household Income via IPW-Weighted Regression with Specific Contextual Effect (Dichotomized Food Insecurity as Mediator: Moderate to Severe Food Insecurity vs Security or Mild Insecurity).

IPW	TE (95% CI)	CDE/NDE (95% CI)	Specific Contextual Effect (95% CI)	Proportion Mediated
Single-level outcome model without EM interaction
Traditional	1.41 (1.05, 1.89)	1.36 (1.01, 1.84)	2.84 (1.81, 4.47)	12.2%
M0	1.47 (1.07, 2.02)	1.41 (1.01, 1.95)	2.79 (1.72, 4.55)	12.8%
M1	1.46 (1.06, 1.99)	1.41 (1.02, 1.94)	2.70 (1.67, 4.37)	10.9%
M2	1.44 (1.05, 1.98)	1.39 (1.00, 1.92)	2.74 (1.69, 4.43)	11.4%
M3	1.40 (0.98, 1.98)	1.33 (0.94, 1.88)	2.93 (1.75, 4.92)	17.5%
Random intercept outcome model without EM interaction
Traditional	1.38 (1.03, 1.84)	1.33 (0.99, 1.79)	2.97 (1.55, 5.69)	13.2%
M0	1.41 (1.05, 1.88)	1.36 (1.01, 1.83)	2.75 (1.33, 5.66)	12.2%
M1	1.40 (1.05, 1.88)	1.36 (1.01, 1.84)	2.82 (1.50, 5.32)	10.0%
M2	1.40 (1.04, 1.87)	1.35 (1.00, 1.82)	2.85 (1.52, 5.35)	12.5%
M3	1.40 (1.05, 1.86)	1.33 (0.99, 1.79)	2.84 (1.37, 5.89)	17.5%
Random intercept outcome model with EM interaction
Traditional	1.36 (0.97, 1.70)	1.32 (0.93, 1.66)	3.80 (2.28, 5.31)	11.1%
M0	1.41 (1.00, 1.76)	1.36 (0.94, 1.72)	3.01 (1.83, 4.27)	12.2%
M1	1.40 (1.02, 1.75)	1.37 (0.97, 1.72)	3.12 (1.95, 4.29)	7.5%
M2	1.39 (1.00, 1.75)	1.35 (0.97, 1.71)	3.15 (1.87, 4.31)	10.3%
M3	1.40 (0.99, 1.74)	1.34 (0.95, 1.71)	3.04 (1.75, 4.24)	15.0%

CDE, controlled direct effect; CI, confidence interval; EM, exposure and mediator; IPW, inverse probability weight; M0, single-level logistic regression with level 1 (sex, total number of surfaces) and level 2 (school type) confounders; M1, single-level logistic regression with level 1 confounders only; M2, 2-level random-intercept logistic regression with level 1 confounders and school-specific intercepts; M3, 2-level fixed-effects logistic regression with level 1 confounders and school dummy variables for constant school-specific effects; NDE, natural direct effect; NIE, natural indirect effect; TE, total effect.

CDE = NDE when there is no EM interaction, only NDE is reported when there is an EM interaction. All effect estimates are expressed in terms of prevalence ratios. Proportion mediated = (NDE × (NIE – 1))/(NDE × NIE – 1).

Comparison with Traditional Methods

When compared with the traditional confounders adjustment without the exposure–mediator interaction, the IPW weighted random-intercept outcome models yielded the most robust results, with a mild difference between traditional regression and IPW weighted regression (Table 2). However, this difference could be more pronounced in cohort studies with multiwave data and feedback effect from outcome to exposure as IPW weighted regression will not block this path. In addition, IPW weighted regression also facilitates the analysis of time-varying exposures and confounders and enables the evaluation of counterfactual outcomes, which is not possible using traditional regression.

Limitations

The empirical example used cross-sectional data with known temporal ordering (family income precedes caries development), supporting causal inference with assumptions of no unmeasured confounding. While longitudinal designs are stronger, cross-sectional analyses can provide valid estimates when temporal precedence is justified (Savitz and Wellenius 2023). Results should be interpreted cautiously, as the primary goal is to illustrate multilevel PS methods.

MLM accounts for clustering (eg, children within schools), but it has limitations with few higher-level units (here, only 14 schools). Simulation studies indicate that <20 to 30 clusters can bias variance components and standard errors, particularly in nonlinear models such as logistic regression, although fixed-effect estimates remain relatively unbiased when level 1 sizes are adequate (Maas and Hox 2005; Snijders and Bosker 2012). In such cases, single-level models with cluster-robust standard errors often provide more reliable inference for fixed effects. MLM may be less suitable when the ICC is negligible (<1% to 2%), clusters are few and substantively important, or dependencies are nonhierarchical. Common misuses include skipping diagnostics (convergence, singular fits, overdispersion), overinterpreting higher-level variance as causal cluster effects, or relying on MLM alone for causal inference without adequately addressing confounding (eg, via PS methods). To assess robustness, we compared single-level and 2-level models with PS weighting: M0 (single-level with level 1 and level 2 confounders), M1 (single-level with level 1 confounders only), M2 (2-level random-intercept with level 1 confounders), and M3 (2-level fixed effects with level 1 confounders and school dummies). Across these approaches, treatment effect estimates were qualitatively consistent, although MLM models yielded slightly wider CIs due to explicit modeling of clustering. These findings support the robustness of our results despite the limited number of clusters. Nevertheless, variance components should be interpreted cautiously, and future work with more clusters is recommended.

Conclusions and Recommendations

MLM represents an essential analytical approach for oral health research, offering a powerful means to rigorously analyze hierarchical data structures that conventional single-level regression cannot adequately capture. By explicitly accounting for clustering and variability across multiple levels, MLM not only improves the precision of estimates but also deepens insight into the complex interplay of individual, contextual, and societal determinants of oral health. Despite its potential, the use of MLM for causal inference remains limited; embedding causal frameworks such as PS weighting within multilevel models can yield more robust and credible estimates of exposure effects.

Key recommendations:

Strengthen methodological integration: Leverage MLM to surpass the limitations of single-level approaches and apply IPW within multilevel structures to enhance causal inference.

Mitigate contextual confounding: Exploit MLM’s capacity to model contextual determinants, enabling more accurate estimation of causal effects.

Advance longitudinal causal analyses: Use MLM to examine time-varying exposures, confounders, and mediators, supporting stronger evidence for causal pathways in oral health research.

Authors’ Contributions

H. Li, contributed to conception and design, data acquisition, drafted and critically revised the manuscript; M.C.M. Wong, R.K. Celeste, contributed to conception and design, data analysis and interpretation, drafted and critically revised the manuscript; M.A. Peres, contributed to conception and design, data acquisition, analysis, and interpretation, drafted and critically revised the manuscript. All authors gave final approval and agreed to be accountable for all aspects of the work.

Supplemental Material

sj-docx-1-jdr-10.1177_00220345261452997 – Supplemental material for Multilevel Modeling for Causal Inference in Oral Health Research

Supplemental material, sj-docx-1-jdr-10.1177_00220345261452997 for Multilevel Modeling for Causal Inference in Oral Health Research by H. Li, M.C.M. Wong, R.K. Celeste and M.A. Peres in Journal of Dental Research

Footnotes

Acknowledgements

Generative AI was used in the production of Figure 1 to enhance clarity and quality of the figure.

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

ORCID iDs

H. Li

M.C.M. Wong

R.K. Celeste

M.A. Peres

Data Statement

Data sharing does not apply to this article.

A supplemental appendix to this article is available online.

References

Albandar

Goldstein

1992. Multi-level statistical models in studies of periodontal diseases. J Periodontol. 63(8):690–695. https://doi.org/10.1902/jop.1992.63.8.690

Antunes

Peres

de Campos Mello

Waldman

EA.

2006. Multilevel assessment of determinants of dental caries experience in Brazil. Community Dent Oral Epidemiol. 34(2):146–152. https://doi.org/10.1111/j.1600-0528.2006.00274.x

Arpino

Mealli

2011. The specification of the propensity score in multilevel observational studies. Comput Stat Data Anal. 55(4):1770–1780. https://doi.org/10.1016/j.csda.2010.11.008

Austin

. 2011. An introduction to propensity score methods for reducing the effects of confounding in observational studies. Multivariate Behav Res. 46(3):399–424. https://doi.org/10.1080/00273171.2011.568786

Austin

. 2020. Advances in propensity score analysis. Stat Methods Med Res. 29(3):641–643. https://doi.org/10.1177/0962280219899248

Bates

Mächler

Bolker

Walker

2015. Fitting linear mixed-effects models using lme4. J Stat Softw. 67(1):1–48. https://doi.org/10.18637/jss.v067.i01

Bryk

Raudenbush

. 1992. Hierarchical linear models for social and behavioral research: applications and data analysis methods. Sage. p 16–61.

Burgette

Preisser

Rozier

. 2016. Propensity score weighting: an application to an Early Head Start Dental Study. J Public Health Dent. 76(1):17–29. https://doi.org/10.1111/jphd.12106

Bürkner

P-C

. 2017. Brms: an r package for bayesian multilevel models using stan. J Stat Software. 80(1):1–28.

10.

Burnham

Anderson

. 2004. Multimodel inference: understanding AIC and BIC in model selection. Sociol Methods Res. 33(2):261–304. https://doi.org/10.1177/0049124104268644

11.

Cafri

Wang

Chan

Austin

. 2019. A review and empirical comparison of causal inference methods for clustered observational data with application to the evaluation of the effectiveness of medical devices. Stat Methods Med Res. 28(10–11):3142–3162. https://doi.org/10.1177/0962280218799540

12.

Celeste

Eyjólfsdóttir

Lennartsson

Fritzell

2020. Socioeconomic life course models and oral health: a longitudinal analysis. J Dent Res. 99(3):257–263. https://doi.org/10.1177/0022034520901709

13.

Chang

Stuart

. 2022. Propensity score methods for observational studies with clustered data: a review. Stat Med. 41(18):3612–3626. https://doi.org/10.1002/sim.9437

14.

Charlton

Rasbash

Browne

Healy

Cameron

2025. MLwiN version 3.14. Centre for Multilevel Modelling, University of Bristol.

15.

Cheung

Wong

MCM

Lam

Corbet

Leung

. 2019. Multilevel modelling for the analysis of longitudinal periodontal data. In: Leung

Goh

, editors. Periodontitis: symptoms, prevention and treatment options. Nova Science. p 231–278.

16.

Cho

De Boeck

Naveiras

Ervin

2022. Level-specific residuals and diagnostic measures, plots, and tests for random effects selection in multilevel and mixed models. Behav Res Methods. 54(5):2178–2220. https://doi.org/10.3758/s13428-021-01709-z

17.

Diez-Roux

. 1998. Bringing context back into epidemiology: variables and fallacies in multilevel analysis. Am J Public Health. 88(2):216–222. https://doi.org/10.2105/ajph.88.2.216

18.

Diez-Roux

. 2000. Multilevel analysis in public health research. Annu Rev Public Health. 21:171–192. https://doi.org/10.1146/annurev.publhealth.21.1.171

19.

Diez Roux

. 2002. A glossary for multilevel analysis. J Epidemiol Community Health. 56(8):588–594. https://doi.org/10.1136/jech.56.8.588

20.

Diez Roux

. 2004. Estimating neighborhood health effects: the challenges of causal inference in a complex world. Soc Sci Med. 58(10):1953–1960. https://doi.org/10.1016/s0277-9536(03)00414-3

21.

Enders

Mistler

Keller

. 2016. Multilevel multiple imputation: a review and evaluation of joint modeling and chained equations imputation. Psychol Methods. 21(2):222–240. https://doi.org/10.1037/met0000063

22.

Enders

Tofighi

2007. Centering predictor variables in cross-sectional multilevel models: a new look at an old issue. Psychol Methods. 12(2):121–138. https://doi.org/10.1037/1082-989x.12.2.121

23.

Fuentes

Lüdtke

Robitzsch

2022. Causal inference with multilevel data: a comparison of different propensity score weighting approaches. Multivariate Behav Res. 57(6):916–939. https://doi.org/10.1080/00273171.2021.1925521

24.

Goldstein

. 1986. Multilevel mixed linear-model analysis using iterative generalized least-squares. Biometrika 73(1):43–56. https://doi.org/10.1093/biomet/73.1.43

25.

Goldstein

. 1995. Multilevel statistical model. Edward Arnold. p 19–80.

26.

Goldstein

Browne

Rasbash

2002. Partitioning variation in multilevel models. Underst Stat. 1(4):223–231. https://doi.org/10.1207/S15328031US0104_02

27.

Grund

Lüdtke

Robitzsch

2016. Multiple imputation of multilevel missing data: an introduction to the R package pan. Sage Open. 6(4). https://doi.org/10.1177/2158244016668220

28.

Haneuse

. 2021. Regression analysis part II: model fitting and assessment. In: Lash

VanderWeele

Haneuse

Rothman

, editors. Modern epidemiology. 4th ed. Philadelphia (PA): Wolters Kluwer. p. 467–496.

29.

Hedeker

Gibbons

2006. Longitudinal data analysis. Wiley-Interscience.

30.

Hox

Moerbeek

van de Schoot

2017. Multilevel analysis: techniques and applications. 3rd ed. Routledge.

31.

Lee

Nguyen

Stuart

. 2021. Partially pooled propensity score models for average treatment effect estimation with multilevel data. J R Stat Soc Series A 184(4):1578–1598. https://doi.org/10.1111/rssa.12741

32.

, et al. 2021. Trend on dental caries status and its risk indicators in children aged 12 years in China: a multilevel analysis based on the repeated national cross-sectional surveys in 2005 and 2015. BMC Public Health. 21(1):2285. https://doi.org/10.1186/s12889-021-12262-x

33.

Listl

Galloway

Mossey

Marcenes

2015. Global economic impact of dental diseases. J Dent Res. 94(10):1355–1361. https://doi.org/10.1177/0022034515602879

34.

Loy

Hofmann

2014. HLMdiag: a suite of diagnostics for hierarchical linear models in R. J Stat Softw. 56(5):1–28. https://doi.org/10.18637/jss.v056.i05

35.

Maas

CJM

Hox

. 2005. Sufficient sample sizes for multilevel modeling. Methodology. 1(3):86–92. https://doi.org/10.1027/1614-2241.1.3.86

36.

Mabangkhru

Duangthip

Chu

Phonghanyudh

Jirarattanasopha

2020. A randomized clinical trial to arrest dentin caries in young children using silver diamine fluoride. J Dent. 99:103375. https://doi.org/10.1016/j.jdent.2020.103375

37.

Morgenstern

1982. Uses of ecologic analysis in epidemiologic research. Am J Public Health. 72(12):1336–1344. https://doi.org/10.2105/ajph.72.12.1336

38.

Morgenstern

. 1995. Ecologic studies in epidemiology: concepts, principles, and methods. Annu Rev Public Health. 16:61–81. https://doi.org/10.1146/annurev.pu.16.050195.000425

39.

Muth

Oravecz

Gabry

2018. User-friendly bayesian regression modeling: a tutorial with rstanarm and shinystan. Quant Meth Psychol. 14(2):99–119. https://doi.org/10.20982/tqmp.14.2.p099

40.

Newton

Bower

EJ.

2005. The social determinants of oral health: new approaches to conceptualizing and researching complex causal networks. Community Dent Oral Epidemiol. 33(1):25–34. https://doi.org/10.1111/j.1600-0528.2004.00190.x

41.

Peres

, et al. 2019. Oral diseases: a global public health challenge. Lancet. 394(10194):249–260. https://doi.org/10.1016/s0140-6736(19)31146-8

42.

Piovesan

ETA

Bernabe

2025. Municipal oral healthcare coverage and dental caries among Brazilian children. Caries Res. [published online ahead of print]. https://doi.org/10.1159/000548560

43.

Preacher

. 2011. Multilevel sem strategies for evaluating mediation in three-level data. Multivariate Behav Res. 46(4):691–731. https://doi.org/10.1080/00273171.2011.589280

44.

Preacher

Zyphur

Zhang

2010. A general multilevel SEM framework for assessing multilevel mediation. Psychol Methods. 15(3):209–233. https://doi.org/10.1037/a0020141

45.

Raudenbush

Bryk

. 2002. Hierarchical linear models: applications and data analysis methods. Sage. p 16–100.

46.

Reichenheim

, et al. 2016. A model-based approach to identify classes and respective cutoffs of the brazilian household food insecurity measurement scale. J Nutr. 146(7):1356–1364. https://doi.org/10.3945/jn.116.231845

47.

Rosenbaum

Rubin

. 1983. The central role of the propensity score in observational studies for causal effects. Biometrika. 70(1):41–55. https://doi.org/10.1093/biomet/70.1.41

48.

Sabbah

Tsakos

Chandola

Sheiham

Watt

. 2007. Social gradients in oral and general health. J Dent Res. 86(10):992–996. https://doi.org/10.1177/154405910708601014

49.

Savitz

Wellenius

. 2023. Can cross-sectional studies contribute to causal inference? It depends. Am J Epidemiol. 192(4):514–516. https://doi.org/10.1093/aje/kwac037

50.

Segall-Corrêa

Marín-León

Melgar-Quiñonez

Pérez-Escamilla

2014. Refinement of the Brazilian household food insecurity measurement scale: recommendation for a 14-item EBIA. Rev Nutr. 27(2):241–251. https://doi.org/10.1590/1415-52732014000200010

51.

Shetty

, et al. 2016. Methamphetamine users have increased dental disease: a propensity score analysis. J Dent Res. 95(7):814–821. https://doi.org/10.1177/0022034516640478

52.

Snijders

Bosker

2012. Multilevel analysis: an introduction to basic and advanced multilevel modeling. 2nd ed. Sage.

53.

Cortina

2009. What do we gain? Combining propensity score methods and multilevel modeling. Paper presented at: APSA 2009 Annual Meeting; Toronto. https://ssrn.com/abstract=1450058

54.

Tadakamadla

Rathore

Mitchell

Johnson

Morawska

2022. Protocol of a cluster randomised controlled trial evaluating the effectiveness of an online parenting intervention for promoting oral health of 2-6 years old Australian children. BMJ Open. 12(10):e056269. https://doi.org/10.1136/bmjopen-2021-056269

55.

Thoemmes

West

. 2011. The use of propensity scores for nonrandomized designs with clustered data. Multivariate Behav Res. 46(3):514–543. https://doi.org/10.1080/00273171.2011.569395

56.

Verbeke

Molenberghs

2000. Linear mixed models for longitudinal data. Springer. p 1–92.

57.

Wan

, et al. 2009. Effects of smoking on healing response to non-surgical periodontal therapy: a multilevel modelling analysis. J Clin Periodontol. 36(3):229–239. https://doi.org/10.1111/j.1600-051x.2008.01371.x

58.

Wierichs

, et al. 2025. Practice-based study on CAD/CAM inlays, onlays, and crowns: longevity and risks. J Dent Res. 104(11):1218–1227. https://doi.org/10.1177/00220345251332882

59.

Zhang

Zyphur

Preacher

. 2009. Testing multilevel mediation using hierarchical linear models: problems and solutions. Organizational Research Methods. 12(4):695–719. https://doi.org/10.1177/1094428108327450

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

4.39 MB