Extending Bayesian back-calculation to estimate age and time specific HIV incidence

CD4-based multi-state back-calculation methods are key for monitoring the HIV epidemic, providing estimates of HIV incidence and diagnosis rates by disentangling their inter-related contribution to the observed surveillance data. This paper, extends existing approaches to age-specific settings, permitting the joint estimation of age- and time-specific incidence and diagnosis rates and the derivation of other epidemiological quantities of interest. This allows the identification of specific age-groups at higher risk of infection, which is crucial in directing public health interventions. We investigate, through simulation studies, the suitability of various bivariate splines for the non-parametric modelling of the latent age- and time-specific incidence and illustrate our method on routinely collected data from the HIV epidemic among gay and bisexual men in England and Wales. Electronic supplementary material The online version of this article (10.1007/s10985-019-09465-1) contains supplementary material, which is available to authorized users.


Introduction
Quantification of HIV incidence and prevalence is key to HIV surveillance and the design and evaluation of targeted interventions. Direct measurement of these quantities is, however, infeasible: infection times are unobserved and, due to the long asymptomatic incubation period, a large proportion of infections remain undiagnosed. Therefore, a number of statistical approaches have been developed to estimate HIV burden from routinely collected surveillance data. The back-calculation method, initially proposed by Gail (1987, 1988) still plays a key role in the monitoring of HIV and other long incubation diseases (Deuffic-Burban et al. 2007;Sweeting et al. 2007;van Sighem et al. 2015). The idea underlying this approach is that the infection process can be reconstructed from time series data on disease endpoint events and knowledge of the distribution of the time between infection and the endpoints of interest. For HIV, in a discrete time formulation, this is formally expressed by the convolution equation (Becker et al. 1991): where a i is the expected number of new AIDS diagnoses in the ith interval, (t i−1 , t i ], h i 0 is the expected number of new infections in (t i 0 −1 , t i 0 ], and f i−i 0 is the probability of an AIDS diagnosis in the (i − i 0 )th interval after infection.
A further extension, again aimed at providing a better insight into the epidemic, has been to estimate time-and age-specific infection rates, making use of the information that age at infection is the strongest predictor of HIV progression. In principle, this would entail the specification of a latent time-and age-(denoted i 0 and j 0 ) specific bivariate surface h i 0 , j 0 , which is particularly challenging to estimate. To avoid modelling a two-dimensional infection surface, Verdecchia and Mariotto (1995), Greenland (1996) and Wand et al. (2009) applied age-independent back-calculation models to diagnosis data stratified by birth-cohorts, deriving age-dependent incidence estimates through the combination of the estimates resulting from each cohort. Also to simplify the problem, some authors (Becker and Marschner 1993;Becker et al. 2003) used a multiplicative model, h i 0 , j 0 = h i 0 h j 0 which however cannot capture different time-trends across age-groups. The non-parametric bivariate step-function in Rosenberg (1995) added flexibility to the infection surface model and Marschner and Bosch (1998) improved parameter identifiability by imposing (thin plate spline) smoothing at the corner-points. The level of smoothing in Marschner and Bosch (1998) was however isotropic (i.e. equal in both the time and age dimensions), which may not always be appropriate, for instance when time and age are measured on different scales. All these approaches solely considered counts of AIDS diagnoses as endpoint data, with the exception of Becker et al. (2003) who additionally incorporated HIV diagnoses.
In this paper we reconsider the problem of age specific estimation but in the context of a CD4 count based multi-state model, extending the work of Aalen et al. (1997), Sweeting et al. (2005) and Birrell et al. (2012) to estimate age-time dependent HIV incidence, CD4 state specific diagnosis probabilities and number of undiagnosed infections. We also extend the work of Marschner and Bosch (1998) by adopting bivariate splines (Wood 2006a) to model the incidence surface as a continuous function of age and time and we further investigate tensor product splines (Eilers and Marx 2003;Wood 2006b), allowing for differential smoothing in the time and age dimensions. In contrast to earlier age-dependent back-calculation approaches (e.g. Marschner and Bosch 1998), the use of a multiplicity of data introduces the complication that the back-calculation model cannot be expressed as a generalised linear model (GLM), so that the inferential problem becomes non-standard. We propose a Bayesian approach to estimation, which can more easily tackle the non-standard nature of the problem and automatically allows propagation of uncertainty to all the derived quantities of interest.
Section 2 describes the motivating application and describes the data available from England and Wales, with relevant notation. Section 3 introduces the CD4-staged back-calculation model and links this model to the data. Section 4 looks at a range of spline models suitable for smoothing the incidence surface and the merits of each spline are examined in the simulation study of Sect. 5. In Sect. 6, appropriate model parameterisations are then used to estimate age-stratified HIV incidence in England and Wales over the last 20 years. We conclude with a discussion in Sect. 7.

Motivating application
The methodology developed in this paper is motivated by the surveillance data routinely collected by Public Health England (PHE) to monitor the HIV epidemic among men-who-have-sex-with-men (MSM) in England and Wales (see Fig. 1). Available data include diagnoses of HIV over time classified in two groups, according to the presence or absence of AIDS related symptoms within 3 months of the initial HIV diagnosis. These will be loosely expressed as diagnoses of AIDS and HIV respectively. Information on the CD4 cell counts around diagnosis (i.e. taken within 3 months of HIV diagnosis) is also available for a large, and increasing, proportion of the new HIV diagnoses (see Fig. 1a, b).
A Bayesian back-calculation analysis of this type of data (Birrell et al. 2012) collected over the whole epidemic history (i.e. 1978-2015), resulted in the estimated yearly number of new HIV infections levelling off at approximately 3000 (see Fig. 2), following a steady increase over the period 2007-2013. However, stratification of new diagnoses by age (Fig. 1c) reveals heterogeneous trends, questioning whether the apparent plateau in incidence might mask contrasting trends in different age-groups, suggesting the need for age specific incidence estimates.
Formally, assume the epidemic period (t 0 , t T ] is split into T disjoint, consecutive intervals (t i−1 , t i ], i = 1, . . . , T . Similarly the age-range (a 0 , a A ] is subdivided into Year of Diagnosis   Birrell et al. (2013) A disjoint, consecutive groups (a j−1 , a j ], j = 1, . . . , A. We will, in places, refer to

Model specification
The data previously described arise as a result of three distinct, interlinked processes: infection, disease progression and diagnosis. Figure 3 shows the structure of a discretetime non-homogeneous population-level CD4 count multi-state model that explicitly specifies the contribution of these three processes to the dynamics of an infected population.
The infection process is modelled by a two dimensional non-homogeneous Poisson process (e.g. Rosenberg 1995), with time (u) and age (v) dependent infection rate λ (u, v). Then the expected number of new infections in (t i−1 , After infection, individuals are subject to competing disease progression and diagnosis pressures, represented by movements to undiagnosed states (1, . . . , K ) with lower CD4 counts and to the absorbing diagnosis states (K + 1, . . . , 2K + 1), respectively.
Given the discrete time framework, the progression and diagnosis processes are expressed in terms of probabilities. Diagnosis probabilities are state, time-and agespecific, to reflect the weakening of the immune system and to allow for the impact of testing campaigns over time, possibly targeted at specific age groups. Denote d k,i, j the probability of being diagnosed from the kth undiagnosed state in the ith time and jth age group. For those infected in the j 0 th age group and remaining undiagnosed, let q j 0 k denote the probability of progressing from the kth to the (k + 1)th state in the same interval.

Model dynamics
Previous work (Aalen et al. 1997;Sweeting et al. 2005;Birrell et al. 2012) has characterised the number of infected individuals in the disease states through a Markov chain.
However, in the model of Fig. 3, progression depends on the age at infection. The number of individuals in a state, at a given time and age group no longer depends solely on the number of individuals at the previous time and age, unless the infected population is also stratified by age at infection. This substantially complicates the model dynamics, now described by progression and transition matrices, Q j 0 i, j and D j 0 i, j respectively, both depending on the age group at infection j 0 . The K × K matrix Q j 0 i, j specifies the probabilities of moving between the undiagnosed states of the model in the ith time interval and jth age group, for individuals infected in the j 0 th age group. Its (k, l)th entry is expressed as: i, j has the (k, l)th entry giving the corresponding probability of moving from the undiagnosed state k to the diagnosed states K + l: Note that the dynamics are slightly different for individuals in state K as progression to AIDS is assumed to always result in a diagnosis. The above matrices reflect the assumption that the time intervals are sufficiently small so that at most one transition event can happen and that diagnosis events occur before progression. We further assume that the time intervals and age groups are of equal width. Let e i, j,K ) T denote the K × 1 vector of the expected number of individuals in the undiagnosed states in the ith time interval and jth age group who were infected in the j 0 th age group. Similarly, μ with entries giving the corresponding expected numbers of new diagnoses in the absorbing states K + 1, . . . , 2K + 1. These are the result of the recursive equations: . 1 j=A is an indicator function, equal to one if j = A and zero otherwise. The starting values of the recursion, when j = j 0 are defined so that: for j 0 = 1, . . . , A, i = 2, . . . , T and we further let e j 0 1, j 0 = (h 1, j 0 , 0, . . . , 0) T for j 0 = 1, . . . , A. Note that although the Ath age-group, in accordance with the data is cumulative, it is assumed that infections do not occur at ages greater than a A . From (4) to (7), the expected number of undiagnosed individuals e i, j and the expected number of new diagnoses μ i, j in the ith time interval and jth age group are obtained by summing over the infection age-groups j 0 : Note that the dynamic equations discussed here can be appropriately modified when, as it may happen in practice, data are available on a coarser time scale, or uneven time and age scales or they might not be collected from the beginning of the epidemic (see Sect. 1 of the online resource).

Likelihood
The aim is to estimate the expected number of new time-and age-specific infections H = {h 1,1 , . . . , h T ,A }, to which we refer as the incidence surface (or simply incidence), and the diagnosis probabilities D = {d 1,1 , . . . , d T ,A }, when the progression probabilities Q = {q 1 , . . . , q A } are assumed to be known from external cohort studies. The components of H and D could be treated as free parameters, however a more parsimonious parameterisation can be achieved by introducing parameters θ and δ respectively, so that H ≡ H(θ ) and D ≡ D(δ). Note that all the quantities defined in Sect. 3.2 become dependent on these parameters. For notational convenience, this dependency will be suppressed, e.g.
By the properties of the non-homogeneous Poisson process (Cox and Isham 1980) characterising the infection process, the number of arrivals into the diagnosis state k in (t i−1 , t i ] and (a j−1 , a j ] results in a set of independent Poisson random variables with means μ i, j,k [Eq. (9)]. Hence the likelihood of HIV and AIDS diagnoses is given by independent Poisson random variables, Y H i, j and Y A i, j : for i = 1, . . . , T and j = 1, . . . , A, where the means are The contribution of the subsample of HIV diagnoses with a linked CD4 count is included based on the assumption that the distribution of the available CD4 counts is representative of the CD4 count distribution for all individuals. As the number of new HIV diagnoses in the ith time interval and jth age group is the sum of K independent Poisson random variables with means μ i, j,1 , . . . , μ i, j,K , the distribution of the number of diagnoses in the states {K + 1, . . . , 2K } conditional on their sum is multinomial: where The likelihood, expressed in terms of θ and δ, is proportional to: 4 Bivariate smoothing methods

Bivariate splines
To parameterise H(θ) we employ bivariate splines. In general terms, given a vector of n observations y = (y 1 , . . . , y n ) T with associated two dimensional covariates : R 2 → R used to smoothly model the (x, y) relationship. Splines are constructed from a set of basis functions {b 1 (x i ), · · · , b p (x i )} and a related p × 1 vector of parameters θ . For any x, the spline takes values g(x) = j θ j b j (x) and can be expressed as a generalised linear model, where the data arise from a distribution of the exponential family, and with n × p design matrix X, with (i, j)th entry X i, j = b j (x i ). Estimation of θ is typically carried out through minimisation of a penalised loglikelihood criteria: where l(θ | y) is the log-likelihood of the data and N s is the number of p × p matrices S s chosen to penalise the roughness of the resulting spline curve. A large number of parameters can be specified to guarantee flexibility, with any induced overfitting effect counteracted by the scaling, through the smoothing parameters λ s , of the penalty term. Large λ s values favour smooth curves over more volatile ones. Closed form, and numerical, solutions are available for obtaining optimalθ andλ s if y arise from a distribution from the exponential family and can be expressed as a GLM (Wood 2006a). Note that (13) can be re-interpreted from a Bayesian perspective, as a sum of the loglikelihood and a log-prior giving a log-posterior distribution. Specifically, the penalty term is equivalent to a zero-mean multivariate Normal prior for θ , with p × p precision matrix N s s=1 λ s S s . Flat priors are implicitly assigned to the λ s , though alternatives could be chosen. Table 1 summarizes all the splines considered in what follows. Two main types of bivariate splines exist: thin plate splines and tensor product splines. Thin plate splines (Green and Silverman 1994) are defined by a bivariate spline basis obtained by introducing a set of knot points κ = {κ 1 , · · · , κ p } (see tps in Table 1 as well as Sect. 2.2.2 of the online resource). Roughness is quantified by the Laplacian integral: that imposes isotropic smoothing (i.e. equal smoothing in the x 1 and x 2 dimensions) and can be conveniently expressed in a quadratic form θ T S 1 θ (in this case N p = 1, see Sect. 2.1.1 of the online resource for details). Thin plate splines may be sensitive to the choice of κ, hence Wood (2003) proposed thin plate regression splines that avoid specifying the location of the knots (see Sect. 2.1.2, and 2.2.3 of the online resource for further details). Here a slight modification, due to Marra and Wood (2011), is implemented (see tprs in Table 1 as well as Sect. 2.1.3 and 2.2.4 of the online resource). Tensor product splines are constructed by defining two univariate splines, with design matrices X (1) and X (2) , (of dimension n 1 × p 1 , and n 2 × p 2 ) and roughness matrices S (1) and S (2) (of dimension p 1 × p 1 and p 2 × p 2 ). The bases for the joint bivariate spline are then obtained by multiplying the basis functions of the marginal splines. Tensor product splines allow for differential smoothing in the two dimensions (N p = 2), by applying the univariate penalty matrices marginally (see Sect. 2.2.5 of the online resource). Eilers and Marx (2003) constructed tensor product splines from marginal cubic B-spline and measuring marginal roughness via a first order difference penalty squared (see ptensbs in Table 1). Wood (2006b) extended their approach to handle any type of marginal spline, such as thin plate regression splines, and any type of marginal penalty, such as the integrated second derivative squared (see ptenstprs in Table 1). The different penalty measures of ptensbs and ptenstprs imply that, in absence of information, the marginal splines revert towards a flat and linear trend respectively.

Splines within back-calculation
The bivariate splines are used to model the log-incidence surface γ = (γ 1,1 , . . . , γ T ,A ) T , where γ i, j = log(h i, j ), by letting γ = Xθ . X denotes the design matrix corresponding to the chosen type of spline. Irrespectively of the parameterisation of D(δ), back-calculation cannot be expressed as a GLM: the likelihood (Eq. 13) includes Poisson and Multinomial terms so that a single link function cannot be specified, and the expected number of diagnoses μ i, j is a non-linear function of θ and δ. In this case, standard algorithms to estimate the spline parameters θ and λ s within a GLM penalised likelihood context cannot be implemented. Although the penalised likelihood can be For tensor product splines, the spline basis is a product of two univariate splines in the time and age dimensions denoted by subscripts (1) and (2) numerically maximized, estimation of λ s and quantification of uncertainty become computationally prohibitive (Brizzi 2018, Sect. 6.4.4). Even Expectation Maximization (EM) based algorithms, often used for back-calculation (Becker et al. 1991;Becker and Marschner 1993;Marschner and Bosch 1998) cannot be efficiently employed, as the derivatives of the likelihood are not analytically tractable. An alternative Bayesian approach (Wood 2016) offers a number of advantages allowing direct estimation of both model and smoothing parameters and automatic quantification of uncertainty. Moreover external sources of information (e.g. underreporting rates, see Birrell et al. 2012) can easily be incorporated and implementation can be achieved using standard software for Bayesian analysis.

Simulation study
Here we investigate the most appropriate type of spline model. To do this, we carried out a simulation study starting from the age-dependent back-calculation model described and p i, j respectively, are then obtained through a generalisation of the dynamical equations described in (5) and (9) (see Sect. 3.1 in the online resource). These are used to simulate data according to (10-12), taking n i, j to be equal to the n i, j , i.e. the number of samples observed in the last 20 years of this study (Fig. 1a). Three data-generating bivariate incidence surfaces are derived by assuming: and v i, j is the proportion of h i occurring among age groups (a j−1 , a j ], with 52 j=1 v i, j = 1, for all i. Three plausible time profiles h i are considered. These are identical until the most recent 3 years when they differ to allow an increasing, a constant and a decreasing trend in incidence (see Fig. 4a). The v i, j are constructed such that, in all the three time profiles, the mean age at infection shifts linearly from age 43, in (t 0 , t 1 ], to 33, in (t 19 , t 20 ]. The resulting age-specific time profiles of the incidence surfaces are shown in Fig. 5. To limit the computational burden of the simulation study, the diagnosis probabilities used to generate the data (Fig. 4b) are assumed to be independent of age, i.e. d k,i, j ≡ d k,i for all j. The values specified are available in the online resource, Sect. 3.

Study design
For each of the three incidence surfaces considered, 50 sets of simulated data were generated. Estimation of the incidence surface was then carried out for each simulated dataset, with the incidence surface modelled using each of the four splines discussed in Table 1. All splines considered have 80 parameters. For thin plate splines, knots are located at intervals of 2 years in the time dimension, and every 6.5 years in the age dimension (i.e. for a total of 10 and 8 knots in the time and age dimension respectively). For each of the two marginal splines of a tensor product we specified 10 and 8 parameters in the time and age dimension respectively (for a total of 80 parameters), using equidistant knots. The weakly informative priors imposed on the reparameterised coefficients are available in Sect. 3.5 of the online resource. The smoothing parameters λ s have a crucial role as they determine the roughness of the estimated incidence curve. To reflect a lack of prior knowledge and a weak preference towards smooth curves, diffuse half-t prior distributions with 2 degrees of freedom and scale parameter 200 are chosen so that 95% of the prior density lies in the [0, 400] region (Gelman 2006).
Alongside the smoothed infection process, a model for the diagnosis process D(δ) is also specified. Diagnosis probabilities are expressed on a logistic scale, i.e. δ k,i = log d k,i 1−d k,i , using a first order random walk: A total of 600 scenarios are considered, where each scenario refers to a combination of the data-generating incidence surface, a spline model for γ (e.g. tprs) and a simulated dataset. Inference is carried out using Stan (version 2.14), which employs Hamiltonian Monte Carlo methods (Hoffman and Gelman 2014;Carpenter et al. 2017). Each posterior estimate is obtained using three chains of 2000 iterations with burn-in of 1000. Splines are implemented via the R package mgcv (Wood 2017), and the reparameterisations discussed in Wood (2016) are implemented for computational efficiency (see online resource, Sect. 2.3). The weakly informative priors imposed on the reparameterised spline coefficients, δ k,1 and σ 2 k parameters are available in Sect. 3.5 of the online resource. The approximate running time per scenario is 10 h. Codes are available at https://github.com/frbrz25/Thesis_Codes.

Assessment
For the mth (m = 1, . . . , 600)  The Predictive Mean Squared Error (PMSE) is the mean of squared errors between the data-generating and the estimated incidence curves. For the mth scenario this has the expression: The distribution of PMSE( H m ) can be evaluated for different splines, with lower PMSE( H m ) values indicating that the data-generating incidence curve is more accurately estimated. PMSE( D m k ) for the diagnosis probabilities, from the kth state can analogously be defined.
Convergence of HMC chains of 1000 (post burn-in) iterations is assessed using thê R statistics of Gelman and Rubin (1992). Figure 6 shows that all the spline models considered (tps, tprs, ptenstprs, ptensbs) reasonably reproduce the time profile of the flat incidence surface, except in the first and last 3 years of the epidemic, where estimates diverge. Estimates in the initial period are sensitive to the choice of the initial expected number of infected individuals π . A further sensitivity analysis (see Brizzi 2018) showed that these estimates are only affected by π for, at most, a period of 7 years.

Results
In the last 3 years, the time profiles of the incidence surface are overestimated in the majority of the scenarios under each spline model (especially ptenstprs and ptensbs). This is induced by an incorrect attribution of recent diagnoses to an increase in incidence (Fig. 6) resulting in a consistent under-estimation of the diagnosis probabilities from state 1 in most recent years. There is also increased variability across the estimates at this time, a common feature of back-calculation.
The age-specific time profiles of the incidence are adequately estimated for all age ranges and spline models. In Fig. 7 aggregated incidence over the 15-24 and 25-34 age ranges are accurate, even in the later years. Estimates in the 35-44 and 45+ age-ranges are more volatile, due to fewer diagnoses occurring in these age-groups. Figure 8 shows PMSE( H m ) for each of the splines. Among thin plate splines, tprs outperform tps (similar findings were obtained in Wood 2003). Among tensor product splines, ptensbs outperforms ptenstprs with the PMSE( H m ) distributions of the tprs and ptensbs being similar. The incidence time profiles, estimated by tprs and ptensbs (Fig. 7), only differ in the latest years, with the estimates from tprs visibly less biased, but more volatile (especially in most recent years). This different performance is attributable to the assumptions on the behaviour of the splines in most recent years, for which data are only weakly informative. The linear trend of the tprs splines, occasionally results in extreme estimates, whereas the ptensbs flattens out. In this simulation study, tprs and ptensbs splines perform similarly. Note that the time intervals and age groups are both measured on a yearly scale and hence the isotropy assumption appears to hold. This assumption is hardly testable in practice  1 (b, d, f, h) for the different splines: tps (a, b); tprs (c, d); ptenstprs (e, f); ptensbs (g, h). Black dashed lines represent the values used in data generation and does not apply in situations where data are collected on an uneven time and age scale; as ptensbs splines do not rely on isotropy, they may be preferred to tprs splines. All of the above conclusions consistently apply when also considering the increasing and decreasing incidence profiles for recent infections. Additionally there was no detectable difference in the goodness-of-fit achieved by the different spline incidence models (see Brizzi 2018, Appendix G.2.1.).
A further sensitivity analysis (see Brizzi 2018, Sect. 4.6.4) revealed that incidence estimates are robust to the specified weakly informative prior for the smoothing parameters λ s .

Application to the MSM-HIV epidemic in England and Wales
As an illustration, we apply the back-calculation model described in Sect. 5, to the data introduced in Sect. 2. Specifically, we focus on reconstructing incidence from the mid 1990s when CD4 data started to become more reliable, including a total of 45,972 diagnoses from 1995 to 2015. Individuals are assumed to have seroconverted between 15 and 66 years of age. The expected number of undiagnosed individuals at the beginning of 1995 (i.e. π ) and progression probabilities Q are set as in Sects. 3.2 and 3.3 of the online resource.

Initial investigations
As in Sect. 5 a yearly scale for both time and age (i.e. T = 21, A = 52) is assumed. Incidence is modelled using a ptensbs spline and diagnosis probabilities through a random walk on a logistic scale, independent of current age. As before, models are implemented in Stan, using four chains of 2000 iterations, the first 1000 of which were burn-in. The resulting posterior sample of 4000 iterations was obtained in approximately 8 h. Figure 9a is a plot of the estimated incidence surfaces obtained by sequentially including an additional year of data from 2010 to 2015. Letĥ model to make the diagnosis probabilities dependent on current age, may allow the model to better adapt to recent changes in the data. We consider three models for the age-dependence of diagnosis probabilities and two alternative time scales, using six models in total.
1. YAID: yearly model, with age-independent diagnosis probabilities (as discussed in the previous section). Let i = 2, . . . , T , j = 1, . . . , A, k = 1, . . . , 4 and: with initial condition: where m k and σ 0,k are known fixed constants, whereas σ k are estimated. 2. YADD0: As YAID, except for an additive term α j in (17) and (18). This term has the interpretation of an age-specific linear time trend in the logistic diagnosis probabilities and α j is estimated, after imposing a N(0, 1) prior on it. 3. YADD1: As YADD0, except that we use an age and state specific α j,k time trend in the logistic diagnosis probabilities. 4. QAID: As YAID, but using a quarterly time scale. 5. QADD0: As YADD0, but using a quarterly time scale. 6. QADD1: As YADD1, but using a quarterly time scale. Figure 9 displays the sensitivity of the estimated time profile of the incidence surface to the sequential addition of further years of data and demonstrates that the QAAD1 model is the most robust among the models considered. Estimates of incidence, both at population and at age-specific level in the most recent years are only slightly revised when further years of data are added, suggesting that the estimated trends in incidence are not artificial.
All the models are consistent with the HIV diagnosis data, as judged by the goodness of fit to the observed data (Brizzi 2018, Appendix H.4.). However the YAAD1 and QAAD1 seem to better fit the AIDS and CD4 count data, especially in the 15-24 and 45+ age ranges. The number of diagnoses with CD4 count in the (200, 500] range, between 2005 and 2015, only increase in the 15-24 age range. A model with a stateand-age dependent starting value for the diagnosis probabilities, allows this feature of the data to be captured. For all age ranges, the posterior-predictive distribution of CD4 count data include all data points, but credible intervals are wide. Although overfitting may be an issue, as suggested by the noisy fit to CD4 data, QAAD1 successfully achieves robust incidence estimates. Figure 10 shows the results obtained from the QAAD1 model. Figure 10a plots the expected number of infections over time; incidence has steadily increased from 2007 onwards, even though a plateau is reached in the latest years. However, age-specific back-calculation reveals that this plateau hides a sharp increase in the expected number of infections for the 25-34 age-group from 2007 (Fig. 10b). Incidence has remained approximately constant in this period for the other age ranges. As a result the distribution of age at infection has shifted towards younger ages: in 2000, 17%, 42%, 30%, 11% of individuals were respectively newly infected in age ranges 15-24, 25-34, 35-44 and 45+, compared to 19%, 45%, 24%, 12% of individuals in 2015 (Fig. 10c).

Illustration results
Similarly, Fig. 10d shows that the diagnosis probabilities from state 1 vary with age, and are estimated to be higher for the 25-34 and the 35-44 age ranges.
The age-dependent back-calculation model further reveals that underlying a constant trend in the expected number of undiagnosed infections (Fig. 10e) Fig. 10 Posterior mean (and related 95 % credible intervals) of a number of relevant quantities for the chosen final model, QAAD1: a incidence surface time profile; b incidence surface time profile, stratified by age range; c proportion of incidence in each age range over time; d diagnosis probabilities from state 1; e expected number of undiagnosed individuals; f expected number of undiagnosed individuals, by age range years, there is a sharp increase in the expected number of undiagnosed individuals living with HIV in the 25-34 age range, and a sharp decrease in the 35-44 age range (Fig. 10f).
It is further interesting to note that age-dependent incidence estimates are reassuringly in agreement with results obtained from the simpler age-independent model (Fig. 2), apart from the 1995-1998 period. Over these years the incidence estimates are highly sensitive to the choice of π , the age specific distribution of the infected undiagnosed population in 1995, which is instead estimated from historical data in the age-independent model. Results from the two models become consistent after 1998, when results are no longer influenced by the specification of π .

Discussion
Back-calculation plays an important role in the monitoring of HIV incidence, based on routinely collected surveillance data. The contributions of Odd Aalen to the backcalculation literature, through pioneering the exploitation of newly available sources of data Farewell et al. 1994), particularly within a multi-state model (Aalen et al. 1997;Sweeting et al. 2005), have been fundamental. These ideas have been central to the development of multi-state back-calculation, where the incorporation of information on CD4 count data around HIV diagnosis, has also enabled estimation of trends in diagnosis probabilities and, consequently, trends in the number of undiagnosed infections (Birrell et al. 2012). In this paper, we have proposed a further extension of this CD4-staged back-calculation model, which allows the joint estimation of age-and time-specific HIV incidence, as well as age-and time-specific diagnosis probabilities. This insight into the HIV epidemic is extremely valuable for targeting and evaluating interventions aimed at reducing HIV prevalence and transmission.
Existing approaches to smoothing incidence over time and age used strong multiplicative assumptions or step functions, which require the arbitrary definition of corner-points. We have thoroughly investigated spline models for smoothing incidence jointly over time and age at infection at a finer level of detail (52 yearly age groups, 80 quarterly time periods). Bivariate splines allow the capture of age-and timeinteractions in a continuous manner, with tensor product splines permitting differential smoothing in the two dimensions. Results from the simulation study show that tensor product splines, constructed from marginal cubic B-splines measuring roughness with first order difference penalty squared (ptensbs), are particularly suitable.
Any back-calculation model provides very uncertain incidence estimates over the most recent period, which are the most crucial to inform public health decision making. This is still true, to some extent, for the model we propose here, motivating a further extension of the backcalculation to incorporate additional data on biomarkers indicative of recent infection (Ndawinz et al. 2011;Yan et al. 2011). Since 2009, PHE has introduced the routine application of Recent Infection Testing Algorithms (RITA) to new HIV diagnoses, allowing the identification of 'recent' infections (Aghaizu et al. 2014). In principle, the proposed multi-state back-calculation framework could be extended through the addition of undiagnosed states for newly infected individuals to include RITA data. In practice, this poses some challenges: many new diagnoses are not RITA tested; and an increase in the number of states will result in both model complexity and computational demand. The approach proposed here already requires long running times, and consideration of a reduced and/or more coarse time scale is often necessary to achieve implementation within an acceptable computational budget. Running times are of the order of 10 and 80 h for a yearly and quarterly time scale, respectively, even when only considering the last 20 years of the epidemic. Despite being faster to implement, yearly models produce estimates that are substantially less stable (e.g. to the addition of further years of data) than the respective quarterly estimates. To successfully incorporate RITA data, future research would benefit from focusing on more computationally efficient inferential approaches than used here.
We have presented a method to estimate age-specific trends in incidence, diagnosis and undiagnosed prevalence of HIV using information from routine surveillance. New diagnosis and early infection biomarker data are becoming increasingly available worldwide, even in less developed countries. Our approach is of value in those countries, as our model can be easily adapted to accommodate limited historical data. By assuming an initial distribution π for the infected individuals across the undiagnosed states at a convenient starting point, age-specific incidence can still be estimated with results sensitive to the choice of π only for the initial years. In countries with established surveillance systems as England and Wales, our approach represents an insightful new tool to guide the targeting of test and treat and pre-exposure prophylaxis strategies (Volz et al. 2018) and to support their evaluation.