Coronavirus disease (COVID-19), declared a global pandemic on March 11, 2020, by the World Health Organization (WHO), disproportionately impacted aging adults. In Canada, senior adults aged 65 and older experienced significantly higher death rates and severe complications compared to younger populations [1], [2], [3]. These outcomes are driven by age-related vulnerabilities, underlying comorbidities, and structural determinants of health, such as socioeconomic status and racial or ethnic disparities [4], [5], [6].

In epidemiological studies, spatial modelling is used to detect geographical variation in disease burden and identify influencing factors [7], [8]. Multiple health outcomes measured at the same spatial locations are often correlated, and their spatial structures may show similar patterns. Traditional approaches often model single health outcomes independently, missing the potential interdependence among related events. Joint spatial modelling addresses this by analyzing multiple outcomes simultaneously, improving the precision of risk estimates while accounting for common underlying factors [9], [10]. This method is notably valuable when outcomes are spatially correlated or when underreporting might obscure incidence patterns [11]. Understanding variations and trends for these health outcomes at the local level, where health services are planned, organized, and delivered, is important in addressing health inequalities during a pandemic [12]. It is, therefore, necessary to provide knowledge about the spatial distribution of risk at a local level with a precise risk estimate to be able to recognize with confidence the areas with high incidence and understand risk factors to improve upon the existing intervention and prevention processes [7].

Although joint spatial models have gained attention, only a few studies have applied them to understand COVID-19 outcomes at a local level, particularly among older adults. Prior research in Sweden used Bayesian joint modelling to assess COVID-19 incidence, intensive care unit (ICU) admission, and mortality for hotspot identification [13], while a study in Ontario employed cointegration analysis to jointly examine hospitalization and deaths over time [14]. However, these studies either focused on broader populations, lacked spatial specificity or were limited to temporal trends. There remains a gap in understanding how severe COVID-19 outcomes cluster geographically among seniors and what shared risks contribute to these patterns.

This study addresses these gaps by applying a Bayesian shared component model to jointly examine three COVID-19 severity-weighted outcomes, deaths, multiple hospitalizations, and single hospitalizations, among Ontario residents aged 65 and older. Using data from January 2020 to March 2022, we identify spatial variations in relative risk across forward sortation areas (FSAs) and assess shared and outcome-specific socioeconomic and demographic risk factors. Additionally, we explore whether rural seniors experienced elevated risks during the pandemic. Unlike previous studies, our approach provides a spatially explicit joint model focused on an older adult population.