In this section, we describe the mathematical framework we use to estimate the prevalence of CHB in populations of interest. We developed a state-transition model of CHB progression based on the natural history of the disease, building on earlier work in [14]. The model includes several late-stage clinical health states, such as DC and HCC, which are more readily diagnosed than the asymptomatic stages of CHB. Indeed, we assume that the symptoms of these late-stage clinical health states are severe enough that they will be diagnosed within 1 year. Using a Bayesian calibration method based on the Metropolis-Hastings algorithm, we used publicly available diagnosis data for these late-stage complications, in combination with diagnosis data for CHB itself, to back-calculate the prevalence and undiagnosed fraction of CHB.

We developed a state transition model to describe the infection, progression, and treatment process, following the natural history of CHB. A schematic of this model is presented in Figure 1. We assume that people infected with CHB will remain in exactly one state for the duration of a given calendar year. Once per year, individuals may transition to a different state based on probabilities derived from the parameters in Table 1 and Multimedia Appendix 1. Some of these state transitions, such as the transition from an undiagnosed state to a diagnosed state, correspond to observable events recorded by public health officials. Our goal is to infer the number of individuals in each state by calibrating the model to reproduce the time evolution of these observables.

The structure of our state-transition model is depicted in Figure 1. For clarity, the early stages in the top portion of the diagram are organized along x-, y-, and z-axes based on their nature. The states are organized in rows and denoted by U, D, and T, corresponding to the cascade of care, specifically, the undiagnosed and untreated states, for those who have CHB but are unaware of their infection; the diagnosed and untreated states, for those with a positive CHB diagnosis who have not initiated treatment; and diagnosed and treated states, for those receiving treatment for their CHB infection. Each of these basic states is further subdivided into sub-stages, organized into columns, that represent the natural history of CHB. These sub-stages are denoted using a subscript. The subscript 1 denotes chronic hepatitis B e-antigen (HBeAg) + hepatitis B (hepatitis B surface antigen [HBsAg]: positive, HBeAg: positive, viral load: high, ALT level: high). The subscript 2 denotes Inactive Hepatitis B (HBsAg: positive, HBeAg: negative, viral load: low, ALT level: normal). The subscript 3 denotes Chronic HBeAg-Hepatitis B (HBsAg: positive, HBeAg: negative, viral load: high, ALT level: high).

Finally, each of these states can exist with or without compensated cirrhosis (CC), which, when present, is represented with a superscript C. In order to cover the whole natural history of CHB, we also include 2 states, U0 and D0, corresponding to the immune tolerant state (HBsAg: positive, HBeAg: positive, viral load: high, ALT level: normal).

People infected with CHB may progress into advanced stages of liver disease, which are depicted in the bottom half of Figure 1. These advanced stages are categorized into the following states: HCC, DC, liver transplant (LT), and post–liver transplant (PLT).

We define the LT state to last at most 1 year. Transplant recipients who have survived more than 1 year after their liver transplantation enter the PLT state. Additionally, we assume that individuals in these later states have a probability of mortality greater than the background all-cause mortality. This is represented by arrows to a liver death state, LD. Background mortality is represented by an arrow to a death state, d. A full list of model states is included in Multimedia Appendix 2.

Finally, we assume that in the year t, there are u(t) new CHB infections. These new infections include observed infections that are diagnosed within the first year, in which case they will be directed to the state D0, as well as unobserved cases that remain undiagnosed and go to U0.

The transition probabilities of the model depend on a mix of known and unknown parameters. The known parameters are obtained from the literature and are presented in Multimedia Appendix 1. The unknown parameters are presented in (Table 1). We will estimate these unknown parameters using our calibration procedure, which is described in later sections of this paper.

Because CHB is initially asymptomatic, we assume that the probability of receiving a diagnosis is constant across all the undiagnosed states U_i. In principle, these diagnosis probabilities can vary significantly from year to year. As a simplifying approximation, we assume that this trend is linear over the period being modeled, with a diagnosis probability of dl in the first year being modeled, tinit, and a final diagnosis probability of dr in the last year modeled, tfin. For notational convenience, we denote this time-dependent quantity dl/r(t). For undiagnosed states with compensated cirrhosis, we add an annual probability of CC-specific diagnosis, dcc, to allow for the fact that compensated cirrhosis may likely produce symptoms, which increases the probability of diagnosis. Thus, the diagnosis probability for an individual in one of the UiC states would be dl/r(t)+dcc. Finally, we apply a different diagnosis probability for the first year of infection, daa. In any given year, daau(t) new cases will be added to the D0 state, with the remaining (1−daa)u(t) cases being added to the U0 state. This first-year diagnosis probability allows us to account for CHB-infected individuals who seek screening following high-risk behavior or immigrants who disclose their health status on entering the country.

Much like the diagnosis probability, we also assume that CHB incidence, μ, is a linear function of time. Following the same notational convention as earlier, the number of new cases in a given year is given by u(t)=μl/r(t)⋅Pop(t), where Pop(t) is the population of the cohort being studied. We note that several of the parameter extrema, such as a CHB incidence of 0.1, are set implausibly high to grant the model the flexibility to explore all plausible parameter values.

The final calibration parameter, rhcc0, is an adjustment parameter that aims to account for uncertainty in the probability of developing HCC. The annual probabilities of developing HCC, as reported in the literature, are compiled in Multimedia Appendix 1.

Finally, we note that these model parameters, both known and unknown, represent the underlying probabilities associated with any basic state transition. In principle, simulating the full state transition model represented by Figure 1 would also require knowledge of joint probabilities. In order to obtain a tractable model, we make the simplifying assumption that each of these basic events is independent.

We highlight some key simplifying assumptions made during the development of the model (Textbox 1).

For each of the model states (U0, D0, U1, ..., HCC, DC, LT, and PLT), the state-transition model generates a time series for the expected number of individuals in each state. We denote these by U0(t), D0(t), ..., and PLT(t), for t=tinit,...,tfin. Additionally, we record the number of individuals transitioning between states. These simulated time series—the number of individuals in each state and the number of individuals transitioning between states—are used to calculate 4 observables, which are used to calibrate the model. The first of these is the number of people receiving treatment, T1(t)+T3(t)+T1C(t)+T3C(t). Based on clinical guidelines [15], treatment is mainly offered to infected individuals in HBeAg+ and HBeAg- states, and thus our model does not include an inactive treatment state, T2. The remaining observables are the number of CHB diagnoses, the number of HCC diagnoses, and the number of late HCC diagnoses, which we define as any transition from an undiagnosed state, Ui-/C, to the HCC state. These diagnosis events correspond to transitions between undiagnosed states to diagnosed states. The details of these computations can be found in (Multimedia Appendix 2).

The observables described earlier are compared to 3 data sources: the annual number of HCC diagnoses, obtained from Statistics Canada [16], which we will denote YHCCTotal(t); the annual number of CHB diagnoses, obtained from the Public Health Agency of Canada’s (PHAC) Notifiable Diseases Online [17], denoted YCHB(t); and CHB treatment data, which is based on prescription data obtained from IQVIA via collaboration with PHAC, denoted YTreat(t). The CHB and treatment data can be directly compared to their corresponding observables, but the other 2 observables require additional assumptions on the data before they can be used.

The HCC figures from Statistics Canada include all instances of HCC, regardless of cause. Following [14], the estimated number of CHB-related HCC cases is 9% of total HCC cases, YHCC(t) = 0.09⋅YHCCTotal(t).

In [18], an estimated 1.07% of CHB diagnoses (198 of 18,434) occur within ±6 months of an HCC diagnosis. Thus, we assume the number of late HCC diagnoses per year is given by

Reporting HBV cases is mandatory. Therefore, we can be confident that the Notifiable Diseases Online database accurately reflects the number of HBV diagnoses. However, although HBV diagnoses must be reported, the chronic versus acute distinction is not reportable in all provinces and territories. Unspecified cases were apportioned based on the chronic-to-acute ratio observed in the specified cases. We note that although the rate of HBV diagnosis per se may be biased (not all people who are infected will seek medical attention, and therefore the Notifiable Diseases database might systematically underrepresent populations less likely to interact with the health care system), the structure of our mathematical model accounts for this possibility. Bias or underreporting in the diagnoses of late-stage complications represents a greater source of uncertainty in the model. As such, the HCC data from Statistics Canada represents the most significant source of data-related uncertainty in our analysis. The data for model calibration are available in Multimedia Appendix 3.

The model contains several unknown quantities that need to be inferred via model calibration. These come in 2 broad categories: the unknown parameters from Table 1, and the initial population sizes for each of the model states.

If the initial population distributions in the year tinit deviate too much from the true values, the calibration converges to a spurious local optimum. In order to avoid this, we use a 2-stage calibration process. We begin by assuming a prevalence of 0.4%, in line with earlier seroprevalence results [3], and distribute these among the states in the model based on percentages used in previous modeling studies [19]. We then calibrated the model against data from a subset of provinces—Alberta, British Columbia, Ontario, and Saskatchewan—whose diagnosis, treatment, and HCC data were the most complete and selected the population distribution that yielded the closest match to their data. These population distributions then served as our initial distribution in all subsequent calibrations.

In order to obtain an initial estimate of incidence, μl/r, we assume the system is in approximate equilibrium and apply the same 0.4% seroprevalence estimate. This yields a preliminary estimate of μl/r≈0.00004.

In principle, the early diagnosis parameter, daa, can take any value between 0 and 1, though in practice we know that values close to 1 are unrealistic. We restrict our search range to 0≤daa≤0.3 [18].

Following [10], we use a 2-stage calibration algorithm, which is designed to account for the uncertainty in the literature-derived parameters. Let vlit= [q01,...,dld4] be a set of the literature-derived model parameters, ie, a set of parameters sampled from the ranges in Multimedia Appendix 1. Similarly, let vcal= [dl,...,rhcc0,U0(tinit),...,PLT(tinit)] denote a set of the calibration parameters, ie, a set of the unknown model parameters from Table 1 and a set of initial population estimates. Taken together, a single pair (vlit,vcal) determines the initial conditions and transition probabilities, and thus fully specifies the behavior of the state transition model. The goal of our calibration is to generate pairs that yield a good fit with the available data.

Our 2-stage calibration algorithm is described as follows:

Stage 1: We begin by generating k sets of the literature-derived parameters,

[v1lit,…,vklit], by sampling uniformly from the ranges in Multimedia Appendix 1. For each sample vilit, we calibrate the model using the Metropolis-Hastings MCMC algorithm [20], using the likelihood function described in Multimedia Appendix 2, until the Markov chain converges to a stationary distribution, P(vcal|Y,vilit). Discarding the burn-in, we obtain a set of n samples for the unknown parameters, [v1cal,…,vncal], sampled from this distribution.

Stage 2: The samples generated in the first stage are conditioned on a choice of vilit. Thus, in order to obtain a full set of model parameters, we perform a second sampling step. We first sample uniformly from the literature-derived parameters obtained in step 1, [v1lit,…,vklit], and then sample uniformly from the associated distribution obtained via the MCMC algorithm. The result is a single combined set of parameters, (vlit,vcal), sampled from the distribution P(vlit,vcal|Y).

For each pair (vlit,vcal) derived in this way, the model produces a time series for each of the model states, U0(t), D0(t), ..., PLT(t). In order to describe the typical behavior of our model, we sample multiple pairs (vlit,vcal) and compute the means of the associated time series, U0¯(t), D0¯(t), ..., PLT¯(t), as well as the corresponding standard deviations.

We obtain a prevalence estimate by summing over all model states and dividing by the total population size.

The undiagnosed fraction is computed as the ratio of all undiagnosed states and the total CHB.

U

Using this 2-stage calibration algorithm, we first estimate the national prevalence. We then conduct a stratified analysis for male and female populations. Finally, we conduct a one-way sensitivity analysis to explore the uncertainty associated with these estimates.

This study uses publicly deidentified data collected by Public Health Agency of Canada and Statistics Canada. The analysis procedure has been approved by the Research Board at the University of Waterloo (#43995).

We conducted a second calibration stratified by sex. Again, the algorithm identified a set of model parameters and initial conditions that closely match the observed statistics on CHB diagnosis and HCC diagnosis for the male population (Figure 3) and the female population (Figure 4). In order to produce these fits, the stratified analyses suggest a slightly higher prevalence of 0.663% (95% CI 0.572%‐0.816%) in the male population, while the proportion undiagnosed in 2021 was estimated to be 60.0% (95% CI 54.8%‐67.0%; Figure 3). In the female population, we obtain a slightly lower prevalence of 0.382% (95% CI 0.347%‐0.424%), while the proportion undiagnosed in 2021 was estimated to be 48.8% (95% CI 44.2%‐54.0%; Figure 4). The male and female fits had deviance information criteria of 9.0 and 73.2, respectively. The differences in prevalence obtained from this calibration process are consistent with the known sexual dimorphism in CHB [21].

In addition to stratifying by sex, we also attempted an analysis based on age cohort. However, because late-stage complications are highly skewed toward older individuals, the dataset for the younger cohort was too small and noisy to make reliable inferences. We intend to address this in future work with better health administrative data.

We performed a 1-way sensitivity analysis of the calibrated model, varying the parameters in Table 1 to determine which components of the model caused the greatest variation in the estimates of the prevalence and undiagnosed fraction. Parameters were varied by ±50%, and we measured the resulting change in the prevalence and diagnosed fraction, averaged over the time period studied (2007‐2021). The results are presented in Figures 5 and 6.

As expected, the yearly CHB incidence has the largest impact on prevalence. Other parameter changes impact prevalence indirectly by changing the mortality rate. This effect was seen most strongly in the probability of liver death associated with late-stage complications. Changing the diagnosis probabilities or HCC incidence also increased prevalence, presumably because of their downstream effects on CHB-related mortality. Similarly, increasing the treatment uptake and treatment effectiveness increased prevalence due to the associated drop in mortality. The effects of daa on prevalence were negligible.

Changes in the diagnosis probabilities had the largest effect on the undiagnosed fraction, followed by the probability of liver death, CHB incidence, and daa, which had comparable effects. Changes in events downstream of diagnosis, such as HCC incidence or treatment effectiveness, did not change the overall prevalence enough to have a significant impact on the undiagnosed proportion.

CHB has a long incubation period, often progressing undetected until the onset of late-stage complications. This creates significant uncertainties in the prevalence and undiagnosed proportion and poses challenges to policymakers who wish to measure the efficacy of public health interventions or track progress toward elimination targets. We have proposed a model-based method to estimate these quantities. Using limited CHB and HCC diagnosis data, we estimate prevalence in 2021 to be 0.483% (95% CI 0.443%‐0.535%). These results are in line with previous 0.4% seroprevalence estimates [3]. Additionally, the proposed method provides estimates of the fraction of CHB cases that remain undiagnosed, in this case 54.3% (95% CI 50.9%‐58.9%). Again, this is in line with previous estimates. The same Canadian seroprevalence study indicated that 46% of respondents who tested positive for a present HBV infection knew of their infected status.

These numbers are comparable to estimates obtained in other countries. In 2022, approximately 205,549 individuals in Australia were living with CHB, representing 0.78% of the population. It is estimated that 27.9% of these remain undiagnosed [22]. The prevalence of hepatitis B varies across European Union/European Economic Area countries, with estimates ranging from 0.1% in Ireland to 4.4% in Romania [23]. The proportion of undiagnosed infections also varies widely, ranging between 45% and 85%, highlighting gaps in national testing programs [24]. Although our results are focused on Canada, other high-income, immigrant-receiving countries like the United Kingdom and Australia may also benefit from applying the methods developed in this paper.

A 1-way sensitivity analysis identifies CHB incidence as the greatest factor influencing overall prevalence, followed by the probability of liver death in late-stage clinical health states, followed by treatment efficacy against cirrhosis. The undiagnosed fraction was sensitive to changes in diagnosis probabilities and CHB incidence, but treatment efficacy did not impact prevalence enough to have a significant effect on the undiagnosed proportion.

The estimation method we have presented has several advantages. For one, it is systematic: the analysis we have performed can, in principle, be extended to any cohort for which data are available. Additionally, the 2-step calibration method accounts for uncertainty in the model parameters that are derived from the literature, yielding statistically robust estimates. However, a few important limitations constrain its use. We made the simplifying assumption that the CHB incidence is a linear function of time, but the 1-way sensitivity analysis identified CHB incidence as the parameter with the largest impact on prevalence. Therefore, the assumption of linearity greatly limits the time resolution of the model. As such, the time series presented in this paper should be taken as indicators of overall prevalence trends, rather than a detailed year-by-year accounting. In particular, the model incidence as currently formulated cannot accommodate sudden changes in vaccine uptake or immigration rates. This latter limitation may be addressed by immigrant sub-cohort modeling, though we currently lack the data required to perform this analysis. There is also considerable uncertainty surrounding HCC incidence. For one, the HCC data we use for calibration include HCC from all causes. Following [14], we estimate that hepatitis B accounts for 9% of HCC cases, but this estimate does not account for annual variation. Moreover, as with all cancers, HCC incidence is known to be highly age-dependent, but the available HCC data were not detailed enough to calibrate a state transition model that accounts for age. This latter limitation makes it impossible to reliably stratify our analysis by age cohort. Fortunately, both of these HCC-related limitations can be addressed with more accurate health administration data.

With these limitations in mind, our results point to some notable public health concerns. Despite Canada’s 2016 endorsement of ambitious elimination targets, our 2021 CHB prevalence estimate does not indicate a significant decrease over previous estimates from 2013. Hepatitis B continues to be a substantial health burden for at least a hundred thousand Canadians. Of these, a significant proportion remains undiagnosed and untreated, putting them at increased risk of serious complications. Our analysis suggests that measures to decrease CHB incidence, such as vaccination and screening, have the greatest potential to decrease the overall disease burden. Based on our current data limitations, further stratified analysis is not possible. Our next step will be accessing health administrative data to overcome this barrier, with the aim of producing more reliable and detailed estimates.

We have presented a method to estimate the prevalence of CHB, using a state transition model based on the natural history of the disease. The estimates are obtained via a 2-stage calibration process that incorporates multiple sources of uncertainty, yielding results that are robust to perturbations of the model. This calibration process relies only on diagnosis data, which are already collected by the health care system, making it an attractive alternative to seroprevalence studies and other methods that require the collection of novel data. Conversely, a lack of high-quality health administrative data significantly limits the accuracy of the model. This is especially true among small populations and younger age cohorts, as younger individuals will typically not have developed a large and statistically robust pattern of late-stage complications. With more data, these estimates can further be stratified by region or by age cohort, allowing them to inform the targeting of public health interventions.