Estimating the incidence and diagnosed proportion of HIV infections in Japan: a statistical modeling study

Surveillance data of HIV and AIDS in Japan

The present study investigated the epidemiological surveillance data of HIV and AIDS in Japan, which is publicly reported by the Committee of AIDS Trends, Ministry of Health, Labor and Welfare, Japan (2018), belonging to the Ministry of Health, Labor and Welfare, Japan. Of the reported datasets, our analyses are focused on Japanese nationals because estimation of infection among foreigners requires accounting for human migration, and the decision of migratory behavior (e.g., leaving Japan) is highly dependent on the diagnosis of HIV infection and AIDS. As of the end of 2017, there were 16,663 HIV infections and 7,587 AIDS cases among Japanese nationals (Committee of AIDS Trends, Ministry of Health, Labor and Welfare, Japan, 2018). As mentioned, HIV diagnoses reflect HIV-infected individuals who undertook voluntary diagnostic testing before the onset of AIDS. An AIDS case indicates a patient who has never been diagnosed with HIV infection prior to an AIDS diagnosis and who meets the clinical diagnostic criteria: (i) confirmed HIV infection and (ii) the presence of one of 23 indicator diseases representing opportunistic infections or tumors. According to the Infectious Disease Law, HIV and AIDS are classified as a category V notifiable disease, and once diagnosed, physicians must notify the case within 7 days of diagnosis. In the present study, the yearly incidence of HIV infections and AIDS diagnoses from 1985 to 2017 was used. The data are structured by sex and also by the most likely route of transmission (e.g., heterosexual, homosexual or intravenous drug use, based on a physician’s interview of patients). The latter information, i.e., the mode of transmission, is discarded because it is believed that a substantial proportion of men having sex with men do not disclose the actual contact and inform physicians that they acquired infection through heterosexual contact (Inoue et al., 2015). Thus, a stratified estimation by sex was conducted. Although the magnitude of the epidemic in Japan is relatively small compared with that in Western industrialized countries, the incidence of HIV infection in Japan is believed to have steadily increased over time, especially among men who have sex with men and young adults (Kihara et al., 2003; Nemoto, 2004).

Derivation of likelihood using a mathematical model

The proposed statistical model is derived from the following partial differential equation (PDE) model, which is referred to as the McKendrick equation (Nishiura & Inaba, 2011; Ejima, Aihara & Nishiura, 2014). Figure 1 shows a compartmental diagram of the data-generating process. Once infected with HIV, individuals who are undiagnosed and in the incubation period experience two different hazards, i.e., the force of HIV diagnosis α(t) that depends on calendar time t and the hazard of illness onset ρ(τ) that depends on the time elapsed since infection τ. Let h (t, τ) be undiagnosed incubating HIV infections at calendar time t and the time since infection τ (i.e., undiagnosed HIV infection without AIDS), the dynamics of HIV diagnosis and illness onset are described by (1)${\displaystyle \left(\frac{\partial }{\partial t}+\frac{\partial }{\partial s}\right)h\left(t,s\right)=-\left(\alpha \left(t\right)+\rho \left(s\right)\right)h\left(t,s\right),}$ with a boundary condition (2)${\displaystyle \lambda \left(t\right):=h\left(t,0\right),}$ where λ(t) represents the HIV incidence (i.e., the number of new HIV infections) at calendar time t. It should be noted that ρ(τ) yields f (τ), the probability density function of the incubation period as follows: (3)${\displaystyle f\left(s\right)=\rho \left(s\right)\mathit{exp}\left(-{\int }_0^s\rho \left(y\right)dy\right),}$

for s>0. It is well known that the McKendrick equation can be solved along the characteristic line, i.e., (4)${\displaystyle h\left(t,s\right)=\lambda \left(t-s\right)\mathit{exp}\left(-{\int }_{t-s}^t\alpha \left(x\right)dx-{\int }_0^s\rho \left(y\right)dy\right),}$ for t −s>0. Eqs. (1) and (4) indicate that the incidence of HIV diagnosis at calendar time t, u (t), is written as (5)${\displaystyle u\left(t\right)={\int }_0^t\alpha \left(t\right)h\left(t,s\right)ds={\int }_0^t\alpha \left(t\right)\lambda \left(t-s\right)\mathit{exp}\left(-{\int }_{t-s}^t\alpha \left(x\right)dx-{\int }_0^s\rho \left(y\right)dy\right)ds,}$ and similarly, the incidence of AIDS cases at time t, a (t), is (6)${\displaystyle a\left(t\right)={\int }_0^t\rho \left(s\right)h\left(t,s\right)ds={\int }_0^t\rho \left(s\right)\lambda \left(t-s\right)\mathit{exp}\left(-{\int }_{t-s}^t\alpha \left(x\right)dx-{\int }_0^s\rho \left(y\right)dy\right)ds.}$ Equations (5) and (6) read similarly to the so-called extended backcalculation (Hall et al., 2008), which is derived from a competing risk model (Marschner, 1994; Cui & Becker, 2000). The abovementioned process can be used as the generalization.

Statistical model and estimation

The datasets are reported in a discrete time interval (i.e., year); thus, here I discretized models Eqs. (5) and (6) as (7)${\displaystyle u_t=\sum _{s=1}^t{\lambda }_{t-s}{\alpha }_t\prod _{x=t-s+1}^{t-1}\left(1-{\alpha }_x\right)\prod _{y=1}^{s-1}\left(1-{\rho }_y\right),}$

and (8)${\displaystyle a_t=\sum _{s=1}^t{\lambda }_{t-s}{\rho }_s\prod _{x=t-s+1}^{t-1}\left(1-{\alpha }_x\right)\prod _{y=1}^{s-1}\left(1-{\rho }_y\right).}$

There is no prior notion as to the shape of the epidemic curve (i.e., the frequency of transmission) over time. Thus, the incidence of HIV infection in year t, λ_t, is modeled as a step function: (9)${\displaystyle {\lambda }_t=\left\{\begin{array}{c}{\lambda }_1\mathrm{for}t<1989,\hfill \\ {\lambda }_2\mathrm{for}1989\le t<1993,\hfill \\ ⋮\hfill \\ {\lambda }_9\mathrm{for}2013\le t,\hfill \end{array}\right.}$ such that the yearly incidence can be directly dealt with as the parameter. The yearly probability of diagnosis in year t, α_t, is similarly modeled as (10)${\displaystyle {\alpha }_t=\left\{\begin{array}{c}{\alpha }_1\mathrm{for}t<1989,\hfill \\ {\alpha }_2\mathrm{for}1989\le t<1993,\hfill \\ ⋮\hfill \\ {\alpha }_9\mathrm{for}2013\le t.\hfill \end{array}\right.}$

The probability mass function of the incubation period is assumed as known, and in discrete time, this is written as ${\rho }_s{\prod }_{y=1}^{s-1}\left(1-{\rho }_y\right)$. As is widely assumed for HIV infection, the incubation period is modeled using the Weibull distribution. Using the property of Weibull distribution with the scale parameter η and shape parameter k, the discrete Weibull model is connected to the continuous version as (11)${\displaystyle {\rho }_s=1-\frac{exp\left(-{\left(\frac{t+1}{\eta }\right)}^k\right)}{exp\left(-{\left(\frac{t}{\eta }\right)}^k\right)},}$

and (12)${\displaystyle \prod _{y=1}^{t-1}\left(1-{\rho }_y\right)=exp\left(-{\left(\frac{t}{\eta }\right)}^k\right).}$

Using the abovementioned model, undiagnosed HIV infections at the end of year t are computed as (13)${\displaystyle x_t=\sum _{s=1}^t{\lambda }_{t-s}\prod _{x=t-s+1}^{t-1}\left(1-{\alpha }_x\right)\prod _{y=1}^{s-1}\left(1-{\rho }_y\right).}$

The diagnosed proportion of HIV infections is calculated either as ∑(a + u)∕∑(x + a + u) or ∑u∕∑(x + u), taking the summations over time. The former calculates the proportion of diagnosed HIV-positive individuals out of the cumulative number of HIV-positive individuals. This calculation has the drawback of including patients with AIDS who have already died by the year of calculation. As of 2017, it has been reported that a total of 2,321 cases resulted in death (Iwamoto et al., 2017). Alternatively, the latter calculates the fraction of individuals who are HIV positive but have not yet developed AIDS out of the cumulative number of HIV-positive individuals but including undiagnosed individuals, considering that the incubation period in most cases of HIV infection is now considerably extended by ART. The drawback of the latter calculation is that patients with AIDS who have survived and have received ART are excluded; thus, the calculated proportion may not be strictly in line with the target figure in the first goal of the 90-90-90 initiative. Therefore, when estimating the undiagnosed number of HIV infections and the diagnosed proportion at the end of 2017, both calculations are made, and the former is adjusted by subtracting 2,321 AIDS deaths from the cumulative count of AIDS cases.

To quantify the proposed system of equations, we estimate parameters λ_t and α_t by means of the maximum likelihood method. Considering that HIV infections are generated as the nonhomogenous Poisson process, the resulting HIV diagnoses and AIDS cases would also follow Poisson distributions. The likelihood function of HIV diagnoses is (14)${\displaystyle L_1=\mathit{constant}\times \prod _{t=1985}^{2017}E{\left(u_t\right)}^{r_t}\mathit{exp}\left(-E\left(u_t\right)\right),}$ where r_t denotes the reported (observed) number of HIV diagnoses in year t in the surveillance record. Similarly, the likelihood of new AIDS diagnoses is (15)${\displaystyle L_2=\mathit{constant}\times \prod _{t=1985}^{2017}E{\left(a_t\right)}^{w_t}\mathit{exp}\left(-E\left(a_t\right)\right),\left(15\right)}$ where w_t denotes the reported number of new AIDS diagnoses in year t. Consequently, the total likelihood L is given by (16)${\displaystyle L=L_1{L}_2.}$

Maximum likelihood estimates of parameters are obtained by minimizing the negative logarithm of Eq. (16). As mentioned above, the incubation period distribution is assumed as known, and to address the uncertainty, three different estimates are derived from published studies (Boldson et al., 1988; Brookmeyer & Goedert, 1989; Munoz & Xu, 1996). A widely cited estimate by Brookmeyer & Goedert (1989) was derived from the study of patients with hemophilia over 20 years of age with η = 11.6 and k =2.5, resulting in a median incubation period of 10.0 years. Boldson et al. (1988) investigated a cohort of AIDS cases in San Francisco with η = 14.3 and k =2.5, yielding a median incubation period of 12.3 years. The estimate by Munoz & Xu (1996) was obtained from the Multicenter AIDS Cohort Study with η = 10.0 and k =1.3, and the median incubation period is 7.5 years. All three estimates have been used in the present study to address uncertainty with respect to the incubation period. In addition to Eqs. (14) and (15), we have also explored the over-dispersed likelihood function, employing the negative binomial distribution with time-independent dispersion parameter for HIV and AIDS counts, respectively, (Althaus, 2015) and compared the Akaike Information Criterion (AIC) against Poisson distributed likelihood, as part of sensitivity analysis.

The 95% confidence interval (CI) of parameters was derived from the profile likelihood. The 95% CI of model estimates (e.g., the number of undiagnosed HIV infections and the proportion diagnosed) was derived using a parametric bootstrap method. In the bootstrapping exercise, model parameters were resampled from a multivariate normal distribution with vectors of mean θ and standard deviation σ. The latter vector was derived from the covariance matrix, taking diagonal elements of the inverse Hessian matrix (${\sigma }^2=\text{diag}\left(H^{-1}\left(\theta \right)\right)$). For each set of parameters, the model solution is obtained, and 1,000 times of parameter resampling results in a simulated distribution of model solutions. By taking the 2.5th and 97.5th percentile points of the simulated distribution, the 95% CI is obtained. All statistical data were analyzed using R version 3.1 (Comprehensive R Archive Network) (R Core Team, 2016) and JMP version 12.0.1 statistical software (SAS Institute Inc., Cary, NC, USA).

Ethical considerations

In the present study, the analyzed data are publicly available (Committee of AIDS Trends, Ministry of Health, Labor and Welfare, Japan, 2018). As such, the datasets used in our study are deidentified and fully anonymized in advance, and the analysis of publicly available data with no identifying information does not require ethical approval.