Healthcare access has been central to health policy and decision-making (Anderson, 2014;Arifkhanova, 2017;Blount & Miller, 2009) as it plays a pivotal role in shaping overall health outcomes across the lifespan (Jerant et al., 2012).

As early as 1978, it has been recognized that healthcare access is hindered by many barriers, including financial, geographic, cultural, and community-driven barriers (WHO, 2021). These access barriers exhibit varying levels of relevance, depending on the type of healthcare system (e.g., universal versus market-based), population at-risk (e.g., paediatric versus adult), services provided (e.g., outpatient versus inpatient), and conditions addressed (e.g., mental health, asthma, diabetes) (Serban, 2019). Consequently, modelling the healthcare access must account for its complexities across all levels: people, processes, organizations, and ecosystems (Rouse & Serban, 2014).

Specifically, healthcare access modelling accounts for the interdependence between the supply (i.e., caseload) of healthcare services available from multiple organizations and individual providers, and the demand for healthcare services of the population at risk seeking access to the services. Thus, healthcare access captures the dynamics between demand and supply of specific services within a healthcare system with its characteristics and constraints. Such dynamics reflect the opportunity or potential to receive care, which can vary across sub-populations and communities (Serban, 2019).

Common approaches to match provider-level caseload and individual or community-level demand for services employ optimization modelling, wherein these models account for system constraints, access barriers, and population preferences (Gentili et al., 2016;Z. Li et al., 2015;Nobles et al., 2014;Patrick et al., 2008;F. Wang, 2012). Optimization models for access measurement commonly take the perspective of a centralized planner determining all the actions within the system and seeking to attain a minimal overall opportunity cost across all communities to achieve a system-wide objective such as reduced wait time or cost-effectiveness. Such an approach would provide a "best-case scenario" of potential access as it optimizes the allocation of resources across the entire system. However, an optimization approach typically fails to account for service inefficiencies, for example, due to constraints and dynamics in the appointment scheduling, variations or seasonality in the demand across the period of interest, or uncertainties in noshow-ups. Such inefficiencies reduce access to care; hence, optimization models tend to overestimate access with implications on decision-making and health policy.

Matching demand and caseload of healthcare services has also been considered within appointment scheduling systems (Ala et al., 2022, Ala et al. 2023;Cayirli & Gunes, 2014;Cayirli & Veral, 2003;Gupta & Denton, 2008;Homem de Mello et al., 2022). Commonly applied to management of patient arrivals in a healthcare system, the primary goal of appointment scheduling is to efficiently allocate resources within healthcare organizations, such as clinics or hospitals, aiming to reduce wait times and improve service availability. Appointment scheduling is adopted in strategic, tactical, and operational decisionmaking (Ahmadi-Javid et al., 2017) towards system improvement performance, mainly for primary and secondary levels of healthcare systems. However, appointment scheduling has not been applied to model patient arrivals at the public health level, in which regional or national health policies are planned (S. N. Roy et al., 2020). This particularly extends to measuring healthcare access, which can be represented as the arrival of patients choosing among multiple potential providers or organizations within a healthcare network, along with demand constraints and preferences that may not align with the available caseload. Providers in the network compete for the demand for care while subject to constraints on providing care and under various sources of uncertainties, including no-show-ups, cancellations, and provider's limit in availability (Agrawal et al., 2023).

Discrete event simulation has demonstrated efficacy in capturing the complex and stochastic system dynamics inherent in many systems and is subsequently used to inform decision-making, particularly within the domain of healthcare service delivery (Günal & Pidd, 2010;Kadri et al., 2014;Philip et al., 2023;S. Roy et al., 2021; S. N. Roy et al., 2020;Ruiz et al., 2024). Such systems are subject to stochastic inputs, including variations in provider availability, patient demand, and patient behaviours (Goodwin et al., 2024;Harper & Mustafee, 2019;Liu et al., 2010;Patrick et al., 2008). Given the flexibility of discrete event simulation for capturing such complex dynamics, simulation-optimization frameworks are frequently utilized to facilitate decision-making under uncertainty in these settings. Recent literature in this area has explored the co-integration of digital twins with discrete event simulation to generate improved representations of real-world stochastic systems (Anyene et al., 2025;Santos et al., 2022). These advanced approaches can integrate real-time data streams into simulation models (e.g., real-time simulation), enabling them to better reflect and control system disturbances (Harper et al., 2023;Marmor et al., 2009;Mousavi et al., 2011;Munavalli et al., 2021;Z. Wang & Hu, 2024). Consequently, when coupled with optimization procedures, this capability supports informed shortterm, or even real-time, operational decisionmaking particularly valuable for applications such as emergency service planning and staffing. However, for public health, real-time data for healthcare service delivery is often not achievable, and longer-term strategic planning, not short-term decisions, is the focus. Decision-making under uncertainty through simulation in the public health domain is the problem we tackle in this paper.

We introduce an access modelling approach using a discrete-event simulation (DES) assigning the demand for healthcare services to the caseload provided by all available providers in a healthcare network to measure or estimate access at the community level. The simulation model captures the dynamics between caseload and demand by weekly appointment scheduling, accounting for system and provider-level constraints and individual-level preferences and access barriers under uncertainties in both caseload and demand. We develop simulation input distributions to account for such uncertainties in caseload and demand inputs informed by data, which, in turn, are used to quantify how such uncertainties impact decision-making. A hierarchical Bayesian input modelling is adopted, where the uncertainty about the distributions of the input parameters is modelled using posterior distributions derived from prior information and data acquired from multiple sources.

We motivate the modelling approach for estimating and making inferences on access to psychosocial services for the population of Medicaid-insured children in Georgia, USA. Psychosocial services encompass a wide range of treatments addressing both the psychological and social factors affecting children's mental health, such as assessments, psychotherapy, psycho-education, and psychosocial rehabilitation (Bitsko et al., 2022;NAMI, 2025, Xie et al., 2025). Specifically, we focus on the impact of the workforce availability on meeting the demand for psychosocial services with comparison across communities, identifying communities to be targeted for intervention to expand the workforce. The workforce shortage for psychosocial services has been identified as a crisis in healthcare delivery (for National Council for Behavioral Health, 2017;Hawkins, 2018aHawkins, , 2018b)). The Health Resources and Services Administration has projected a shortage of at least 50,000 behavioural health providers by 2030 (for U.S. Department of Health and Human Services, Health Resources and Services Administration, National Center for Health Workforce Analysis, 2018). However, the lack of access to psychosocial services does not affect the population evenly. Fewer than 15% of children in poverty needing psychosocial care receive services, and even fewer complete treatment (Hodgkinson et al., 2017;Santiago et al., 2013). Many of these children are insured through public insurance, including Medicaid or Children's Insurance Program (CHIP); however, such programs face additional access barriers due to lack of availability of providers, transportation challenges, and lack of parents' availability to take time off work, among others. Thus, in this paper, we will apply the simulation access model to make inference on access to psychosocial services.

The contributions of this paper are multi-fold: 1) To establish a DES framework for modelling the dynamics between demand and caseload, differentiating caseload by care settings and accounting for varying population stratification in estimating demand; 2) To implement the access simulation model accounting for uncertainty in caseload and demand using a Bayesian hierarchical model; and 3) To demonstrate the value of the simulation approach by evaluating interventions to improve access to psychosocial services for Medicaid-insured children.

The remainder of this paper is organized as follows. Section 2 introduces the simulation model, describing the flow process of the weekly arrivals of new and continuing patients seeking care from available providers. Section 3 introduces the details of the hierarchical Bayesian model and the distributional assumptions within the context of the motivating application. This section also introduces the data sources informing the hyperparameters of the prior distributions for both the demand and caseload. Section 4 demonstrates the numerical results from the simulation model, with the focus on policies to improve access. Section 5 concludes this paper.

In this section, we describe the implementation details of the DES model within a motivating application: access to paediatric psychosocial services. While we demonstrate the applicability of the DES model within a specific application, the framework is general and applicable to other service settings or study populations.

For the ease of presentation, we present relevant notation in Table 1. We differentiate between model parameters (Table 1(a)), schedule, and prioritization, and define the process logic governing the interactions among patients and providers. The model inputs, including patient arrival, patient demand, and provider caseload, are simulated samples of real-world data, modelled by an hierarchical Bayesian approach described in Section 3.

The DES model tracks multiple patient states (Table 1(d)) and summary statistics (Table 1(c)), where its output is further used in deriving several outcome measures (Table 1(e)). We evaluate system effectiveness and efficiency by patient accessibility and provider utilization, and system equity by identifying service shortage areas.

We aim to model access to psychosocial services at the community level within a region of interest over the study period (i.e., year 2018). We denote by J the number of provider locations or practices delivering psychosocial services, and denote by C the number of communities for which demand for services is generated. Each community serves as a demand location, approximated by its geographical centre. We model access through making appointment scheduling decisions on a weekly basis. For each week t with t ¼ 1; . . . ; T, the j-th provider with j ¼ 1:; ; :J has the capacity to provide A t:j appointments (herein caseload). In each week t with t ¼ 1; . . . ; T, I t;c number of patients from the c-th community with c ¼ 1; . . . ; C begin requesting appointments. Each patient i from community c requests a total of D c;i appointments across the treatment period.

In the context of the motivating application, patients consist of Medicaid-insured children from Georgia, with census tracts serving as proxies to communities. Census tracts are small, homogeneous, and contiguous geographic subdivisions defined by the US Census Bureau (US Census Bureau, 2022). The DES model matches demand estimated at the census tract level to available caseload from each practice location every week.

At the core of the DES model stands the appointment scheduling logic, emulating the complex matching dynamics between demand and caseload, derived from deterministic or stochastic processes of patients (e.g., treatment state, no-show, arrival order) and of providers (e.g., available caseload), along with constraints on travel distance, yearly caseload, and prioritization. Below we highlight these processes supporting the appointment scheduling decisions at the beginning of each week t at both individual and system levels (Figure 1).

We present the processes for assigning care to individual patient i from community c starting at the beginning of week t. We assume the patient either starts requesting his/her first appointment in week t stochastically (i.e., generated by I t;c ) defined as system entry, or may continue receiving services that week from previously assigned provider j defined as system re-entry. We refer to patients labelled system entry as new patients and to the patients labelled system reentry as assigned patients; together they make up the patients in the appointment scheduling system. A patient may exit the system in week t after the system entry or re-entry if that patient's demand cannot be met, or may exit the system at the end of the requested treatment if that patient's demand has been fully met. We differentiate the system exit as without or with re-entry, based on whether the same patient enters the system again at the beginning of week t þ 1. Patients interact with the appointment scheduling system, transitioning among eight different patient states in Table 1(d) as depicted in Figure A1 in Appendix A. We also describe each state in detail below.

At the beginning of week t, new patient i from community c, with UNASSIGNED status, enters the system and starts requesting appointments (Figure 1(a)). The DES seeks to match the requested appointments with a provider who has caseload available to accommodate the requests.

Due to travel barriers, patient i can only travel δ miles for an appointment, restricting the available providers. Among this set of nearby providers, the DES implementation checks whether some provider j DES Model Input and Relevant Sets Set J Number of providers: provider j for j ¼ 1; . . . ; J C Number of communities: community c for c ¼ 1; . . . ; C I Number of patients: patient i for i ¼ 1; . . . ; I Model Input D c;i Number of demanded appointments from patient i in community c I t;c Number of patients in community c that start requesting appointments in week t A t;j Number of available appointments by provider j in week t (c) DES Model Statistics t Index of the current simulated week t � Index of a future simulated week T i;t appt List of appointments scheduled for patient i at week t M c;i Number of missed appointments of patient i in community c until week t R c;i Number of remaining demand of patient i in community c until week t, after previous served appointments (d) DES Patient States UNASSIGNED Patient without provider assigned, enters to request first appointments SCHEDULED

(1) UNASSIGNED patient got assigned to provider who can serve immediately (2) Patient is scheduled to re-enter for subsequent future appointments PENDING

(1) UNASSIGNED patient got assigned to provider who can serve in near future (2) Patient is scheduled to re-enter for the first appointment in near future UNSERVED UNASSIGNED patient exits the system, unable to be assigned to any provider MISSED Patient exits the system upon missing too many scheduled appointments RESCHEDULED

(1) Patient re-enters per schedule to be served by the assigned provider (2) Patient is scheduled to re-enter for subsequent future appointments UNMET

(1) Patient exits the system because assigned provider can no longer accommodate patient's remaining demand (2) Patient may re-enter until the previously scheduled appointments are exhausted COMPLETED Patient exits the system after all demand are fulfilled.

Patient Accessibility %unservedP Patients who are unable to receive any service %unmetP Patients who are unable to receive complete service %unmetD Counterpart to %unmetP, quantify incomplete demand Provider Utilization %utilization Provider utilization Shortage Area %shortage State-wide shortage area proportion short Whether a community is a shortage area can accommodate T N appt appointments (Figure 1(d)) in the current week t or some future week t � . However, due to limiting provider caseload A t � ;j , patient i may not be assigned to a provider, causing him/her to exit the system as UNSERVED. Otherwise, provider j is identified to accommodate patient i's requested appointments with the service starting at week t (i.e., fulfil one demand immediately; SCHEDULED) or at some future week t � (PENDING). SCHEDULED or PENDING patients are then assigned to provider j, reentering the system based on his/her scheduled appointments in future weeks T i;tþ1 appt .

At the beginning of week t, we denote by T i;t appt the set of weeks wherein an assigned patient i has been assigned to provider j for appointments (Figure 1(b)). Regardless whether the patient shows up or not, week t is removed from the patient's scheduled appointments.

If patient i shows up to the appointment (with probability 1 À p), then the remaining demand R c;i is decreased by one. Patient i may exit the system without re-entry if there is no remaining demand (COMPLETED) or re-enter the system for future appointments. The simulator may schedule more appointments for patient i (herein reschedule), after checking whether the assigned provider j can accommodate these further requests (Figure 1(d)). Failure to do so indicates that the system cannot fully fulfil patient i's demand (UNMET). Note that an UNMET patient i will not be rescheduled, but can re-enter the system until all his/her scheduled appointments in T i;tþ1 appt are fully met.

On the other hand, if patient i does not show up to his/her appointment (with probability p), M c;i is incremented by 1. If patient i misses M appointments, patient i will exit the system without re-entry as MISSED. Otherwise, patient i's state remains unchanged while the simulator may reschedule for future appointments (Figure 1(d)).

The DES model stipulates provider's service processes by specifying (1) the order to process patients entering in the same week (prioritization); (2) whether provider j can accommodate patient i (allocation); and (3) how to adjust the provider caseload when patients exit the system (deallocation).

At the system level, provider prioritization takes place when multiple providers nearby can accommodate a new patient. Additionally, when patients arrive without a prespecified order, the simulator needs to determine which patients should be served first (patient prioritization). For provider prioritization, the simulator matches each patient with the provider that has the highest yearly caseload. Patient prioritization is determined by the policies defined in Table 1(a) dictating which subgroup of patients to be processed first; Policy 1, 2, and 3 respectively prioritize the assigned patients, new patients, and patients with fewer demand. As depicted in Figure 1(c), the simulator first groups and orders the patients' entries and re-entries each week based on the prioritization policy. Among each group, the arrivals are randomized and served by the provider according to the first-in-first-out principle.

At an individual level, during week t, the allocation process of provider j checks whether patient i can be accommodated as depicted in Figure A2 in Appendix A. Given patient i's demand and state, provider j evaluates how many appointments are needed to reserve for the patient (denoted by K), how many gap weeks are allowed when accommodating patient i (denoted by T gap ), and which week to begin service (denoted by t � ). Based on these decisions, provider j checks the availability in week t � ; A t � ;j , using a recursive procedure illustrated in Figure 1(d) to determine whether the provider can accommodate patient i, which will then be input back into the patient's processes (Figure 1(a-b)). If provider j can accommodate patient i, the simulator further retrieves the list of future appointments for patient i, T i;t appt , and assigns patient i to provider j by decreasing provider's caseload A t � ;j .

show up for his/her scheduled appointment in week t, or when patient i exits the system as MISSED, the simulator frees up his/her assigned provider j's caseload accordingly (Figure 1(e)). In case of a no-show, A t;j increases by one since provider j will not serve patient i in week t. In case that patient i exits as MISSED, A t � ;j increases by one for all future weeks t � in patient i's scheduled appointment list T i;tþ1 appt , reflecting that provider j no longer serves patient i.

The simulator monitors the processes on a weekly basis, outputting multiple weekly statistics such as (1) number of patients in each patient state; and (2) caseload for providers after scheduling decisions; see Table 1(e). These statistics can then be aggregated to derive the outcome measures at different resolutions (i.e., temporal, geographic). In the context of the motivating application, we focus on deriving such outcomes towards generating decision-making insights for population healthcare access to psychosocial services considering the trade-offs between system effectiveness (patient accessibility), efficiency (provider utilization), and equity (shortage areas). Common access measurements in the existing literature include provider shortage, healthcare service utilization and travel distance (Gentili et al., 2016;Joseph & Phillips, 1984;Nobles et al., 2014;Serban, 2019), which assist public service agencies such as the Georgia Department of Community Health and the Georgia Budget and Policy Institute to evaluate healthcare access and strategize the resources (GBPI, 2016;Georgia DCH, 2025). Compared to classic access models, using a DES approach can capture novel and granular access measures, derived from individual-level statistics. Formulations for these access measures together with aggregation across different resolutions are provided in Table A1 in Appendix A.

Both new and assigned patients alike experience access barriers when requesting services from providers with limited caseload quantified by the measures described below.

New patients might exit the system if no nearby provider can accommodate their initial appointment requests. The metric %unservedP evaluates the patient-level inaccessibility to the system, computed as the proportion of UNSERVED patients among new patients.

Assigned patients might exit the system if their assigned providers cannot accommodate all of their demand. The metric %unmetP evaluates the patient-level inadequacy of receiving complete psychosocial treatments, computed as the proportion of UNMET patients among assigned patients. Similar to %unmetP, %unmetD is the appointment-level measure taking the amount of unfulfilled demand into consideration.

Due to the inevitable inefficiency in the appointment scheduling as well as the patients' no-show behaviours, not all providers' caseload is utilized. We thus quantify the utilization of provider caseload as an indicator of system efficiency.

Initially, each provider j has A t;j caseload available to serve patients each week t. Towards the end of each week t, after allocating appointments per patients' requests as well as deallocating appointments due to no-shows, we record the remaining caseload, denoted as L t;j . The metric %utilization is computed by dividing the difference between A t;j and L t;j by A t;j .

Shortage areas arise due to the significant imbalance between patient demand and provider caseload in specific areas or communities. The more such service areas exist in a region of interest, the higher the systematic disparities in access to care.

For community c, short θ ðcÞ indicates whether community c experiences service shortage. Specifically, we adopt a hypothesis testing framework described in Section 4.4.1 to identify short θ ðcÞ by assessing %unmetDðcÞ against an acceptable rate θ for all communities c ¼ 1; . . . ; C simultaneously. The metric, %shortage θ , is computed as the proportion of communities experiencing service shortage.

In this section, we describe the hierarchical Bayesian input models we adopt to generate the DES model inputs in Table 1(

The primary data sources include ( 1 The study population consists of Medicaid-insured children, aged 3 to 17, with a current mental health (MH) diagnosis. We consider only MH diagnoses that are prevalent in youth population (e.g., anxiety disorder, depressive disorder); disorders with rare representation are prone to misdiagnosis and/or require different levels of care than prevalent disorders. We refer to children in this study population as patients herein.

More details on data sources and derivations are provided in Appendix A

The DES model matches patient demand with provider caseload under granular temporal and geographical resolution (i.e., every week, each community).

We model the number of patients entering for appointments (patient arrival) and the number of appointments needed for patients (appointment requested) separately, incorporating the hyperparameters listed in Table A2(a). The hyperparameters are estimated using the Medicaid claims data following the procedures described in Appendix A. I t,c . We assume a multinomial distribution for the number of patient entries per week from each community following the three-step hierarchical framework below (Figure 2(c)).

(1) Model for Yearly Patient Arrival I: We denote the (stochastic) count of patients entering for appointments between 2018-01-01 and 2018-12-31 by I. Although a Poisson distribution is a sensible choice to model randomness in I, it lacks flexibility in modelling as its mean and variance are equal. The patient count data in the healthcare domain tends to be over-dispersed, with the variance being typically larger than the mean. Therefore, we adopted the negative binomial distribution to correct for such over-dispersion. Hence, I,NBðr i ; pi Þ. (2) Model for Weekly Patient Arrival I t : The Tdimensional random vector, I T ¼ ðI 1 ; . . . I T Þ, represents the number of patients entering for appointments weeks 1 � t � T in 2018. We assume that I T follows a Multinomial distribution with total sample size I and T-dimensional proportion vector D T : I T ,MultiNðI; D T Þ. A Dirichlet distribution is imposed as the conjugate prior for D T to account for uncertainty about I T . Hence, D T ,Dirðr 0 þ ĨT Þ. (3) Model for Weekly Patient Arrival from Communities I t;c : The C-dimensional random vector, I t;C ¼ ðI t;1 ; . . . ; I t;C Þ, denotes the number of patients entering from communities 1 � c � C during week t. We assume that I t;C follows a Multinomial distribution with total sample size I t and C-dimensional proportion vector I t;C ,MultiNðI t ; D t;C Þ. A Dirichlet distribution is imposed as the conjugate prior for D t;C accounting for uncertainty about I t;C . Hence, D t;C ,Dirðl 0 t þ Ĩt;C Þ.

We adopt a truncated geometric distribution to model the number of appointments requested by patient i (Figure 2(b)). That is, D c;i ,Geoð 1 ωÞ, with D c;i bounded between ωlb and ωub . The geometric distribution implies that each treatment service has probability ω of being the last appointment independently of all previous appointments. While in reality, future visits depend on the present outcomes, we believe that at population-level, the geometric distribution with constant success probability can adequately capture the macro-level variation for the access system.

We model the provider caseload distribution by incorporating the hyperparemeters listed in Table A2(b). The hyperparameters are estimated from the Medicaid claims, NPPES repository, and APEX dataset following procedures described in Appendix 9.2.2.

The weekly caseload for all providers is assumed to follow a Multinomial distribution, which can be simulated by the two-step hierarchical framework in Figure 2(d).

(1) Model for Yearly Provider Caseload A j : The Jdimensional random vector, A J ¼ ðA 1 ; . . . ; A J Þ, represents the caseload of all J providers between 2018-01-01 and 2018-12-31. We assume that each A j follows a Poisson distribution: A j ,Poiðλ j Þ with λ j being their expected yearly caseload. A Gamma distribution is assumed as the conjugate prior for λ j ,Gammaðα j ; βj Þ. (2) Model for Weekly Provider Caseload A t;j : For each provider j, the T-dimensional random vector, A T;j ¼ ðA 1;j ; . . . ; A T;j Þ, denotes the weekly caseload of provider j for all T weeks. We assume that A T;j follows a Multinomial distribution with total sample size A j and Tdimensional proportion vector S T;j , i.e., A T;j ,MultiNðA j ; S T;j Þ. The conjugate prior for S T;j is a Dirichlet distribution, i.e., S t;j ,Dir

In this section, we demonstrate the applicability of the DES model by investigating access to psychosocial services for Medicaid-insured children in Georgia, USA. Nonetheless, the DES model is applicable to other geographic divisions or states in the USA because the data sources, in particular the NPPES database, are publicly available. We explore various factors that impact healthcare access through different analyses followed by relevant decision-making insights. In Section 4.2, we examine the outcome measures annually, to analyse the impact of different DES model parameters on state-wise accessibility and efficiency. In Section 4.3, we consider a more granular temporal resolution to provide insights into how the parameters contribute to the weekly outcome fluctuations. In Section 4.4, we refine the geographic resolution to the census-tract level, focusing on the disparities in service shortage and the effectiveness of policy interventions. We establish several contributions as follows. First, we demonstrate how the DES is able to model the interactions between Medicaid-insured children requesting appointments for psychosocial services and MH-specialized providers available to provide such services on a weekly basis. Specifically, the DES captures granular information on the behaviour of the access outcomes of interest at multiple resolutions, providing decision-making insights across different temporal (weekly/monthly/yearly) and geographic (census tract/county/state) scales. Second, we evaluate different scheduling policies specified through the DES parameters, suitable for comparative analyses on system effectiveness, efficiency, and equity. Lastly, we show the importance of accounting for data uncertainty through the hierarchical Bayesian model, facilitating more robust decision making.

Table 2 displays 36 scenarios, labelled 1.1 through 1.24 and 2.1 through 2.12, based on which we analyse healthcare access from different viewpoints. Each scenario represents one set of DES parameters focusing on specific access behaviours while keeping several parameters constant (Table B1 Appendix B). For instance, scenario 1.1 specifies the parameters as

Þ and adopts policy 1 for prioritization. We analyse the outputs from the DES model over a 52-week horizon from Z ¼ 20 independent replications with model inputs sampled from the hierarchical Bayesian models described in Section 3. The simulation begins with providers having no appointments allocated and no patients in the system. Therefore, we run the model with T warm weeks before collecting statistics for the simulation horizon and with T wrap weeks after, mitigating the biases pertaining to simulation initiation and termination. With %unmetP and %utilization as performance measures, the statistics derived from the DES model stabilize

We first present how the DES model facilitates uncovering important access insights involving model parameters used in appointment scheduling. We explore 24 scenarios (Table 2(a)) combining different schedule and priority policies (see Table 1(a)). Access effectiveness and efficiency are assessed by patient accessibility (%unservedP, %unmetP, %unmetD) and provider utilization (%utilization) respectively, with statistics reported in Table 3. We include estimated mean, median, and interquartile range (IQR) for these four measures.

Recall policy 1 prioritizes assigned patients, policy 2 new patients, and policy 3 patients with fewer demand. Also recall that a higher %unservedP measure indicates larger proportion of patients remaining inaccessible to service, a higher %unmetP (or %unmetD) measure implies more patients receiving inadequate service, and a higher %utilization demonstrates more efficient utilization of providers' caseload.

We extract three major insights inferred from statistically analysing the four outcome measures.

Across the 24 scenarios, each with 20 replications, we observe strong correlations among the four outcome measures displayed in Figure B2 in Appendix B. Particularly the negative correlation (r ¼ À0:98) between %unmetP and %unservedP highlights the trade-off between serving new versus assigned patients. Inspecting %unmetP and %unservedP in Table 3 strengthens this conclusion as high %unmetP consistently accompanies low %unservedP in all scenarios. The one-way analysis of variance (ANOVA) test on %unmetP and %unservedP on all scenarios (α = 0.05; p-values � 0) followed by the Tukey-HSD test (α = 0.05) reinforces this trade-off. For example, scenarios 1.4 and 1.8 that adopt policy 1 yield significantly lower %unmetP (52%) with higher %unservedP (18%) among all scenarios, whereas scenario 1.9 or 1.13 that adopt policy 2 exhibit the exact opposite results (83% and 1%, respectively). Insuring accessibility among assigned patients hinders that for new patients, likely due to limited provider availability.

A strong negative correlation exists between %unmetP and %utilization (r ¼ À0:99), suggesting that providers achieve higher utilization by assuring treatment adequacy as they complete appointments requested from assigned patients. The one-way ANOVA test on %unmetP, %unmetD, and %utilization (α = 0.05; p-values � 0) together with Tukey-HSD test (α = 0.05) further underlines this insight as the pairwise significance pattern displayed in Figure 3 illustrates the similarities among these three measures. Additionally, scenarios 1.4/1.8 yield significantly lower %unmetP (52%) or %unmetD (51%) with higher %utilization (84%) whereas scenarios 1.9/1.13 exhibit the opposite pattern (83%, 70%, 49% respectively). Under limited caseload, the appointment scheduling system could attain better efficiency if providers prioritize the completion of existing treatments to the patients they are already assigned.

appt or T E gap . We determine priority policies, T N appt , and T E gap as significant factors that impact patient accessibility and provider utilization according to the factorial ANOVA (α ¼ 0:05). Comparing across Table 3, among the three priority policies, we observe policy 1 consistently outputs low %unmetP and %unmetD regardless of other model parameters (Figure B3 Appendix B). High T N appt (e.g., 3) or T E gap (e.g., 1) leads to the same conclusion despite different priority policies. Furthermore, increasing T N appt tends to yield greater reduction in %unmetP (3% vs. 13%) and %unmetD (5% vs. 22%) compared to increasing T E gap (Figure 4). In order to assure better treatment adequacy among patients (and subsequently higher efficiency), providers should consider incorporating policy 1 or increasing T N appt .

The DES model can assist in the discovery of detailed temporal insights originated from appointment scheduling as model statistics fluctuate with time (Figure 5(a)). Specifically, the time-dependent variability in patient accessibility help identify periods with large uncertainty. Comparing against different appointment scheduling parameters informs the direction to control variation in the pursuit of service consistency.

Reusing the scenarios in Table 2(a) with T E gap ¼ 1, we analyse weekly %unmetD over a year and report the temporal variations in Table B2 in Appendix B for the 12 scenarios. The mean/median/IQR on the row labels are computed across Z ¼ 20 replications, which are not to be confused with the mean/median/IQR of the weekly %unmetD presented in the table.

Across the 12 scenarios, we observe different temporal variations in treatment adequacy uncertainty measured as time-dependent IQR of weekly %unmetD (Figure 5 (b-d)). Attempting to reserve the maximum number of appointments (T E appt ¼ 6, T N appt ¼ 3) in advance leads to greater variability regardless of priority policies, which translates into greater uncertainty for patients' treatment continuity throughout the year. Such phenomenon is likely due to fluctuations in provider availability throughout the year. Providers may consider reducing T E appt or T N appt to ensure consistency in providing adequate services to patients. Together with insight #1.3, we identify a mean-variance trade-off in treatment adequacy introduced by T N appt . That is, increasing T N appt improves treatment adequacy on average (lower mean %unmetD within a year) whereas reducing T N appt ensures consistent estimation on the proportion of patient receiving adequate service each week (lower variance for weekly %unmetD).

We observe the highest variability in treatment continuity under scenario 1.8, particularly around the week of Thanksgiving (t ¼ 48), likely due to limited caseload during holidays while assigned patients continue demanding appointments. In order to control such variability, providers can dynamically adjust the appointment scheduling prior to Thanksgiving, such as switching to policy 3 (Figure 5

The DES model can additionally reveal geographical disparities in healthcare access, guiding interventions to improve system equity.

Common interventions for geographic communities experiencing service shortage (i.e., shortage area) include increasing provider caseload. To analyse shortage area, we examine 12 scenarios (Table 2(b)) across different scheduling parameters and three interventions (i.e., Intervention 1/2/3 increases caseload by 0%/10%/20% for all provi- ders), assuming prioritization of the assigned patients and T E gap ¼ 1. We identify shortage area based on short θ ðcÞ and %shortage θ with θ ¼ 50%. We investigate this rudimentary caseload increase; however, targeted interventions for improving access in specific communities can also be considered.

We follow a systematic approach to identify shortage communities in Georgia characterized by the most pronounced unmet demand for paediatric psychosocial services (%unmetD). By running Z ¼ 20 replications of the DES model we obtain Z realizations of %unmetDðcÞ for C communities from which we apply C hypothesis tests to determine whether yearlyaggregated %unmetDðcÞ for community c with c ¼ 1:; ; ; C exceeds the access standard θ. For each community c ¼ 1; . . . ; C, we perform Mann Whitney U Tests under the following null and alternative hypotheses:

Because C hypothesis tests are evaluated simultaneously, we incorporate the Benjamini-Hochberg correction (Benjamini & Hochberg, 1995) to account for multiplicity by controlling the false discovery rate. Based on the spatial structure (rural vs. urban) arising from dependence between geographic communities, we adopt a structure-adaptive Benjamini-Hochberg algorithm (i. e., SABHA; α ¼ 0:05) (Li & Barber, 2018) to derive short θ ðcÞ, using data from the Health Resources and Services Administration (HRSA, 2024) HRSA_Shortage to inform the proportion of communities experiencing service shortage. Following this approach, short θ ðcÞ takes the value 1 when SABHA rejects H c 0 for community c (shortage area), and the value 0, otherwise.

Recall that the census tracts serve as proxies to communities in the motivating application. Using θ ¼ 50%, we report the count of shortage census tracts together with %shortage across 12 scenarios in Table B3(a) in Appendix B. Additionally, examining all three interventions (Table B3(b) Appendix B), we categorize census tracts (herein shortage category) into tracts (1) experiencing no service shortage with and without the increase in caseload, (2) experiencing service shortage before increased caseload (Intervention 1), (3) experiencing no service shortage after Intervention 2 or 3, and (4) experiencing service shortage despite the interventions. We visualize ðT E appt ; T N appt Þ ¼ ð6; 3Þ in Figure 6(c-g), with other scenarios visualized in Figures B4-B6 in Appendix B. 4.4.2.1. Insight #3.1: High T N appt Greatly Improves Access Equity. With increased caseload, we see consistent reduction in shortage areas across different combinations of scheduling parameters (Figure 6 (a-b)), with ðT E appt ; T N appt Þ ¼ ð6; 3Þ showing the greatest reduction (38.33% to 12.39% to 9.63%). T N appt ¼ 3 mitigates more shortage area ( > 20%) compared to T N appt ¼ 1 ( < 3%), suggesting the use of large T N appt to improve equity. 4.4.2.2. Insight #3.2: Rural Community Faces Greater Challenge in Access Equity. The breakdown of shortage categories uncovers the characteristic of community-wide access equity in response to interventions. Under ðT E appt ; T N appt Þ ¼ ð6; 3Þ, 61.47% (mostly urban) of the census tracts never experience service shortage and 7.49% (mostly rural) remains shortage area despite increased caseload (Figure 6(c)). Although the composition greatly varies under different scheduling parameters (Figure 6(b)), rural communities consistently experience service shortage even after intervention.

In this case study, our objective is to illustrate the wide range of analyses performed at different geotemporal resolutions, and showcase the insights derived accordingly. We use high-level (i.e., yearly, state-wise) resolutions to identify the overall access patterns of the system with respect to different parameters. For example, we observe a trade-off between new and assigned patients' accessibility. The high accessibility for assigned patients resembles continuity of treatments, or treatment adequacy, which we found to be highly correlated with high provider utilization, or system efficiency. This adequacy of treatment can be achieved by prioritizing the assigned patients when scheduling the appointments or increasing the parameters T N appt or T E gap . We then increase the temporal resolution to inspect the detailed weekly outcome fluctuations attributed to different model parameters and identify special time periods within the year (e.g., Thanksgiving) that requires attention to ensure adequate health access. Finally, we increase the geographic resolution to identify disparities in access across communities in Georgia. Despite increases in caseload, rural communities still experience a greater service shortage compared to urban communities.

The methodology developed in this study establishes a reliable framework for modelling the dynamics between demand and supply of healthcare services, specifically psychosocial services, accounting for uncertainties from various sources, such as temporal fluctuations and patient no-show patterns. The DES model logic represents a significant advancement, providing a tool that not only measures access but also accommodates the inherent uncertainties within the healthcare system.

The outcomes of the simulation case study offer valuable insights into the effectiveness, efficiency, and equity within the healthcare network, across different temporal and geographical scales. Both appointment schedule parameters and priority policies contribute to the outcomes of patient accessibility and provider utilization. Incorporating input uncertainties enriches the outcome measurements from single estimation to quantifiable distribution, assisting in statistical inference. For example, when providers prioritize their assigned patients, we consistently observe higher provider utilization together with fewer patients receiving inadequate treatments, compared to prioritizing new patients. Similar observations arise if providers reserve more appointments for new patients or do not require patients arriving for consecutive appointments. These consistent conclusions provide insights for adjusting scheduling strategies to ensure effective and efficient care delivery.

The simulation outcomes generate geographical distribution of shortage area, accompanied with interventions such as increasing provider caseload, outlining the map of service disparities. Input uncertainties again provide confidence in identifying specific census tracts with pronounced unmet demand, through hypothesis testing and access mapping. Compared against a simple deterministic approach to classify census tracts as shortage versus non-shortage area, such statistical framework over multiple hypothesis testing not only supports integrating known geographic prior on structural disparities but also achieves reliable conclusion for decision-making. These findings provide a foundation for targeted interventions and strategic resource allocation to enhance access to paediatric psychosocial services. Furthermore, the observed patterns of geographical shortage area distribution, particularly those census tracts experiencing service shortage even after increased caseload, emphasize the need for a more equitable distribution of healthcare resources. The implications of these disparities impact the overall effectiveness and efficiency of the healthcare delivery system. Addressing these regional imbalances is imperative for ensuring comprehensive and equitable healthcare access for all, irrespective of geographic location.

Overall the insights derived from this research hold relevance for policymakers, healthcare administrators, and practitioners seeking evidence-based approaches to improve healthcare accessibility. By integrating rigorous methodology with practical applications, this study contributes to the ongoing discourse on healthcare delivery optimization and underscores the importance of targeted interventions to address regional disparities. As the healthcare landscape evolves, these findings will provide a foundation for informed decision-making, ultimately working towards the goal of enhancing healthcare access for vulnerable populations.

Nonetheless, we acknowledge the limitations of this study. The proposed DES model is a simplified representation of a complex system, through which we aim to capture essential population-level dynamics for evaluating access to psychosocial services. Several improvements, such as incorporating complicated realistic distributions to model patient visit behaviours or introducing adaptive scheduling strategies in response to temporal demand variations, could be worthwhile for future research. Additionally, the proposed framework can be applied to model psychosocial service access to any geographic division within the USA, however, scaling the model to a larger region requires computationally efficient implementations, which is a future research direction.

(3) A t;j represents the caseload available from provider j based on model input; (4) L t;j represents the left-over caseload from provider j after appointment scheduling in week t;

We use P t2T measureðt; cÞ to compute temporal aggregation; P c2C measureðt; cÞ for geographic aggregation; and P t2T;c2C measureðt; cÞ for both. (A) 2 Hierarchical Bayesian (B) 2.1 Data Source Primary Data Source Detail: Recall the primary data sources include Medicaid claims, NPPES repository, and APEX dataset. The Medicaid dataset provides comprehensive patient-level records for the Medicaid-insured population, including personal information such as age and residential zip codes. Furthermore, it also provides detailed care-level and provider-level information, such as service date, procedure (Current Procedural Terminology or CPT) codes, diagnostic (International Classification of Diseases Tenth Revision or ICD-10) codes, provider identification (i.e., National Provider Identifier or NPI), and care setting (e.g., clinic or home). The NPPES repository comprises information on all healthcare providers in the United States identified by NPIs, including details such as practice location, entity type (i.e., individual or organization), and taxonomy indicative of their specialities. The APEX dataset details psychosocial services provided in public participant schools in Georgia, providing information such as school name, counts of provided services, and service month.

Auxiliary Data Source and Relevant Derivation: Some auxiliary data sources include (1) a set of ICD-10 codes to identify patients with prevalent mental health conditions set of relevant CPT codes to (2) a set of relevant CPT codes to extract information on (outpatient) psychosocial services, including psychotherapy, assessment, psychoeducation, and rehabilitation (APA, 2022); (3) place of service (POS) and type of service (TOS) codes in search for case settings; (4) zip code crosswalk file (PD & R, 2018) to estimate residential census tract from zip code, and (5) geospatial information for geocoding and travel distance calculation. The set of CPT codes are derived by applying a frequency procedure on 2016-2018 Medicaid claims data and complemented with additional well-established national psychotherapy codes (Cummings et al., 2022).

We provide the complete list of hyperparameters used in hierarchical Bayesian input modelling in Table A2, for both demand (Table A2(a)) and caseload (Table A2(b)).

We associate occurrences of demand with treatment periods, consisting of temporally ordered (appointment) sequences, each of which signifies a patient visiting a provider (identified by NPI) on

Table A1. Outcome measure Formulations.. Measure Formula %unservedPðt; cÞ P i2I 1ðt¼Ec;i ¼Xc;i Þ P i2I 1ðt¼Ec;iÞ %unmetPðt; cÞ P i2I 1ðt¼Xc;i;0 < Fc;i < Dc;i Þ P i2I 1ðt¼Xc;i ;0 < Fc;i Þ %unmetDðt; cÞ P i2I 1ðt¼Xc;i;0 < Fc;i < Dc;i Þ�ðDc;iÀFc;i Þ P i2I 1ðt¼Xc;i;0 < Fc;iÞ�Dc;i %utilizationðtÞ P j2J 1ðAt;j ÀLt;j Þ P j2J 1ðAt;j Þ %shortage θ P c2C 1ðshortθðcÞ¼1Þ jCj Table A2. Hierarchical Bayesian model hyperparameters.

Appointment Scheduling Analysis: We visualize the pairwise correlation between %unmetP and the other three (%unmetD, %utilization, and %unservedP) in Figure B2. We present the outcome changes with respect to T N appt and T E gap in 10.2. Temporal Analysis: We report temporal analysis statistics as Table B2.

We report statistics as Table B3. We present the detailed visualizations for scenarios in Figure B4 (T E appt ¼ 4, T N appt ¼ 1), Figure B5 (T E appt ¼ 4, T N appt ¼ 3), and Figure B6 (T E appt ¼ 6, T N appt ¼ 1).

Table B1. Constant DES parameter. Parameter Value Description δ 30 Patients should not travel more than 30 miles to get services p 0.1 10% probability that a patient doesn't show up for appointment T N gap 0 New patient should be assigned consecutive appointments T N buffer 2 New patient may be delayed up to 2 weeks to start services M 3 Patients removed from the system upon 3 missed appointments T warm 52 52 weeks warm-up T 52 52 weeks simulation (i.e., one year) T wrap 52 52 weeks wrap-up Table B2. Temporal Analysis report Table. Across-replication mean, median, and IQR of the weekly mean, median, and IQR of %unmetDðtÞ. High (red) values are highlighted for IQR. T E appt T N appt %unmetDðtÞ Mean %unmetDðtÞ Median %unmetDðtÞ IQR Policy 1 Policy 2 Policy 3 Policy 1 Policy 2 Policy 3 Policy 1 Policy 2 Policy 3 4 1 64.40 69.11 64.79 64.41 69.13 64.80 0.77 0.71 0.94 3 51.05 55.85 53.30 51.03 55.85 53.33 1.38 0.96 1.33 Mean 6 1 64.41 69.12 64.78 64.42 69.13 64.77 0.84 0.70 0.98 3 50.84 55.44 53.50 50.92 55.46 53.50 1.66 0.98 1.26 4 1 64.72 69.09 65.49 64.73 69.11 65.48 0.79 0.71 0.93 3 51.38 55.90 53.67 51.31 55.74 53.66 1.31 0.92 1.22 Median 6 1 64.61 69.14 65.40 64.62 69.17 65.32 0.83 0.66 1.01 3 50.96 55.77 54.02 51.23 55.74 54.13 1.51 0.94 1.09 4 1 2.41 2.01 2.26 2.46 1.98 2.28 0.31 0.26 0.36 3 6.82 4.26 5.30 6.82 4.40 5.30 0.48 0.30 0.66 IQR 6 1 2.26 1.95 2.47 2.23 1.83 2.53 0.32 0.29 0.38 3 7.12 5.06 5.18 7.00 5.09 5.12 0.71 0.35 0.57 Table B3. Shortage areas and intervention Analysis report Table. (a) Number of shortage census tracts with %shortage. (b) Number of census tracts in each shortage category. (a) Shortage Area T E appt T N appt #short c ð%shortageÞ Intervention 1 Intervention 2 Intervention 3 4 1 1627 (82.93) 1597 (81.40) 1577 (80.38) 3 755 (38.48) 429 (21.87) 346 (17.64) 6 1 1626 (82.87) 1601 (81.60) 1570 (80.02) 3 752 (38.33) 243 (12.39) 189 (9.63) (b) Shortage Category T E appt T N appt #short c ð%Þ ; f1g f1; 2=3g f1; 2; 3g 4 1 321(16.36) 44(2.24) 34(1.73) 1563(79.66) 3 1207(61.52) 302(15.39) 131(6.68) 322(16.41) 6 1 322(16.41) 40(2.04) 43(2.19) 1557(79.36) 3 1210(61.67) 467(23.8) 138(7.03) 147(7.49)