Relating Recent Infection Prevalence to Incidence with a Sub-population of Non-progressors

We present a new analysis of relationships between disease incidence and the prevalence of an experimentally defined state of `recent infection'. This leads to a clean separation between biological parameters (properties of disease progression as reflected in a test for recent infection), which need to be calibrated, and epidemiological state variables, which are estimated in a cross-sectional survey. The framework takes into account the possibility that details of the assay and host/pathogen chemistry leave a (knowable) fraction of the population in the recent category for all times. This systematically addresses an issue which is the source of some controversy about the appropriate use of the BED assay for defining recent HIV infection. Analysis of relative contributions of error arising variously from statistical considerations and simplifications of general expressions indicate that statistical error dominates heavily over all sources of bias for realistic epidemiological and biological scenarios. Numerical calculations validate the approximations made in analytical relations.


Introduction
Reliable estimation of disease incidence (rate of occurrence of new infections) and prevalence (the fraction of a population in an inflected state) are central to the determination of epidemiological trends, especially for the allocation of resources and evaluation of interventions. Prevalence estimation is relatively straightforward, for example by cross-sectional survey. Incidence estimates are notoriously problematic, though potentially of crucial importance. An approximate measure of incidence in a population is required for the proper planning of sample sizes and costing for clinical trials and other population based studies. Repeated follow-up of a representative cohort is often touted as the 'gold standard' for estimating incidence, but is costly, time intensive and still prone to some intrinsic problems. For example, there may be bias in the factors determining which subjects are lost to follow-up. Furthermore, ethical considerations demand that a cohort study involve substantial support for subjects to avoid becoming infected, which may make the cohort unrepresentative of the population of interest.
Numerous methods have been proposed for inferring incidence from single or multiple cross-sectional surveys rather than following up a cohort [1,2,3,4,5,6,7,8,9,10]. A central idea in most of these [1,3,4,5,6,8,10] is to count the prevalence of a state of 'recent infection', which naturally depends on the recent incidence. The relationship between the two is in general not simple and depends in detail on the recent population dynamics as well as distributions which capture the inter-subject variability of progression through stages of infection, as they are observed by the specific laboratory assays used in the test for recent infection (TRI). For this approach to be sensible, a working definition of 'recent infection' must be calibrated, for example by repeatedly following up subjects over a period during which they become infected. This is effectively as much effort as one measurement of incidence by follow-up. The calibration is then used to infer incidence from each of many subsequent independent cross-sectional surveys.
Owing to the devastating impact of the HIV epidemic, and the many challenges of research and intervention design, the problem of estimating HIV incidence has attracted considerable interest in recent years. The prospect of using a TRI is in principle very attractive. Given the range of values of incidence likely to be observed in populations with a major epidemic (say 1-10% per annum) a mean definition of 'recent' of approximately half a year is desirable to yield reasonable statistical confidence for sample sizes of a few thousand. The BED assay is currently the leading candidate for such a test, but controversy has arisen about the possibility of conducting a reliable calibration. This stems from the fact that a subset of individuals (approximately 5% [8,10], potentially variable between viral and host populations) fail to progress above any statistically useful threshold set on the assay in the definition of 'recent' infection. This subset of individuals, who consequently remain classified as 'recently infected' for all times, poses a problem to which there is currently no consensus remedy.
We present a new analysis of the interaction between epidemiological trends and a model of inter-subject variability of progression through an experimental category of 'recent infection'. Our model yields simple formulae for inference even when a fraction of the population fails to progress out of the recent category. The only physiological assumption required to deal with the non-progressors is that their survival after infection is the same as the progressors. This assumption is also implicit in previous work on using the BED assay to estimate incidence.
A key conceptual point about our analysis, which distinguishes it from all others of which we are aware, is that we confront the fundamental limitation of what can be inferred from a cross-sectional survey. In particular, even perfect knowledge of an instantaneous population state does not uniquely determine the instantaneous incidence. At best, a weighted average of recent incidence can be inferred. Although the discrepancies between this weighted average and instantaneous incidence can be shown to be small compared to statistical errors (for our application), it can in principle be systematically incorporated into estimation of trends from multiple cross-sectional surveys.
The article is organized as follows. In Section 2 we first develop a basic continuous time model defined by a time dependent incidence and susceptible population, a distribution of times after infection spent under the threshold on a TRI and a distribution of post-infection survival times. We note that there is in principle no specific relationship between the instantaneous incidence and the prevalence of individuals who are infected and under the threshold. At best, one obtains a relationship between the prevalence of under-threshold individuals and a convolution of the recent incidence with a specific weighting function which is implied by the use of a TRI. This relationship in principle includes all moments of the distribution of the waiting times that individuals spend under the threshold. We show that, for realistic rates of variation in the susceptible population, only the mean of the waiting time distribution is needed, and a simple expression for a weighted average of the incidence is obtained. The basic model is extended to allow some fraction of individuals (specified by a new parameter) to be assigned infinite waiting times under the threshold of the TRI. This leads to only very minor modifications of the previous expression for weighted incidence, namely a systematic 'subtraction' of over-counted 'not recently infected' individuals which are included in the experimental category 'under threshold'. This subtraction is similar, but not identical, to that proposed in [10].
Section 3 explores the consequences of designing a cross-sectional survey with a sample size N based on the relations derived in Section 2. Using a systematic expansion of the incidence estimator in powers of 1/ √ N (derived in the Appendix), we note consistency of the estimator (no bias in the limit of large N ) and derive an approximate expression for its relative variance. These expressions facilitate error estimation both from a study design and data analysis point of view. On calibration, we note that trends, as opposed to absolute values, for incidence can be obtained without any information about the distribution of finite waiting times under the threshold. However, an estimate of the fraction of non-progressors is essential. A key observation is that, for realistic population dynamics and sample sizes, statistical error is much larger than bias.
Numerical simulations are presented in Section 4. These demonstrate that the approximate statistical analysis of Section 2 is essentially correct, with discrepancies of the size expected from the 1/ √ N expansion. In the conclusion, we note that the framework presented here is quite general and is applicable to any TRI, as long as any finite probability for non-progression can be calibrated, survival is the same for progressors and non-progressors, and there is no 'relapse' from over to under the recency threshold. It may be possible to modify the analysis to relax these requirements. We point out that the relationship between the present analysis and other proposals for using the BED TRI for measuring HIV incidence should be more systematically investigated. Some preliminary work in this direction has already been produced [11].
2 Relating the Prevalence of 'Recent Infection' to Incidence We now outline a quite general approach to relating the key demographic, epidemiological and biological processes which are relevant to the estimation of incidence from cross-sectional surveys of the prevalence of 'recent infection'. This refines the naive intuition that a high prevalence of 'recently infected' individuals means a high incidence.

The Basic Model
A test for recent infection, such as the CDC STARHS algorithm, is typically obtained through the administration of two assays of different sensitivity. The more sensitive test distinguishes infected from healthy individuals and the less sensitive test, applied to the infected individuals, distinguishes 'recent' from 'long' established infection.
Consider an assay which yields a quantitative result, the value of which typically increases with time from infection. The BED assay is of this type. The quantitative result is a normalized optical density (ODn), which is an increasing function of the proportion of HIV-1 specific IgG. The hypothetical ongoing observation of an individual might lead to a curve similar to that displayed in Figure 1A. Such an assay becomes the less sensitive component of a test for recent infection when we declare a threshold value and define 'observed to be recently infected' to be a test value under the threshold.
In practice, there is inevitably inter-individual variation in these progression curves. Plotting the curves for multiple individuals on a single graph would lead to something like Figure 1B. Clearly, the category 'observed to be recently infected' is not sharply defined by a time boundary, and we now adopt the more precise labels under threshold (U) and over threshold (O). The variability of times spent in the under-threshold category, conditional on being alive long enough to reach the threshold, is captured by a distribution of waiting times f U|A which may look something like that shown in Figure 1C.
It is now possible to construct the basic epidemiological model shown in Figure 2A. Since our analysis will focus on a variety of survival functions S(t), we shall refer to the susceptible population as the healthy population H(t). Upon infection, individuals move from the healthy population to the under-threshold infected population U (t). Those that live long enough, reach the threshold after a waiting time distributed according to f U|A , and enter the over-threshold population O(t). We denote by τ U|A a waiting time generated by the density f U|A . The corresponding cumulative probability function is given by while the probability of 'survival' (persistence) in the population U, conditional on being alive, is and the mean waiting time is Analogously, we define f A , τ A , F A , S A and E [τ A ] in order to capture survival times (how long individuals remain alive after the moment of infection). We assume that survival time and waiting time to threshold are independent in this model. Hence, the probability, at a time delay ∆t after infection, of being simultaneously alive and under the threshold on the assay is Similarly, the probability of being simultaneously alive and over the threshold is Hence, the mean time spent in the category U, accounting for both progression and mortality, is E [τ U ]. New infections are generated by a non-homogeneous Poisson process with an intensity (probability per unit time of new arrivals) λ(t). Let the instantaneous incidence be given by I(t). Then, in a period dt around time t, the expected number of new cases dC is given by We can now write down numerous expressions resulting from the model. For example, the expected number of historically accumulated cases up until time t is given by The expected populations of infected persons under and over the threshold at time t are and Our goal is to relate I for recent times to instantaneous values of H, U and O. We wish to emphasize the cautionary note that there is fundamentally a loss of information when one tries to characterize the history of a population based on observations made at a single time point. The recent historical course of a population, and even instantaneous values of state variables which are rates, like incidence, are in general not inferable from counting data obtained in a single survey. This is due to the fact that counts are, unavoidably, convolutions of historical epidemiological variables, as in (8) and (9) above. Any attempt to derive incidence estimates from the counting of infections accumulated in the recent past faces this problem, and at best some sort of weighted average of the recent values of incidence can be inferred without additional assumptions.
In general, a well defined construction of an estimate for incidence, based on data obtained in a survey conducted at time t, will be some sort of weighted average of past values where W (s, t) is a statistical weight arising from a convolution of population history and biology. Since our goal is to estimate incidence from a count of recently infected individuals, a natural weighting function is one that reflects the relative contributions to this count made by infections from different times in the recent past. Hence, we consider since W (s, t) is proportional to the probability that individuals are 1. available for being infected at time s < t, and 2. still alive and classified as under the threshold at time t, if infected at time s.
Using (11) as the weighting function leads to an expression for the incidence given by The numerator in this expression is an instantaneous state variable, while the denominator in principle involves data from the entire history of the system as well as full knowledge of the survival function S U . We will shortly show how to obtain a systematic approximation of the denominator, but a few remarks are in order about the practical meaning of this weighted average. In the case of constant incidence, the weighted average is the instantaneous value. In the case of a narrowly peaked distribution f U|A , a constant rate of change of I and a constant healthy population, the weighted average is approximately equal to the instantaneous incidence at a time E [τ U ] /2 prior to the cross-sectional survey. If trends are fitted to the results of multiple cross-sectional surveys, this time lag could be more systematically accounted for.

A Simple Expression for Incidence
A simplified expression for weighted incidence in terms of sample and calibration data is now derived. We express the healthy population using the expansion and use the identity which follows directly from integration by parts. It then follows that the weighted incidence (12) can be expressed as If the healthy population is approximately constant for the times where the weight W is non-vanishing, we obtain the simple relation which gives a weighted recent incidence in terms of instantaneous state variables (H and U ) and the expected waiting time in the under-threshold category. Expectation values of the form E [τ n U ] are not state variables and should in principle be measured independently of a particular cross-sectional survey. Usually this would be accomplished in a calibration cohort follow-up study. Thus, after calibrating some of these expectation values, we can deal with a truncated expansion for H without further assumptions about the behavior of I.
Some cautionary comments on calibration are, however, necessary. It seems unlikely that accurate measurements for terms of higher order than just E [τ U ] are practically possible for the case of the BED assay. Finding the non-leading terms H i (for i ≥ 1) in the expansion of H will also require considerable demographic research.
These considerations mean that it is most likely that we will be constrained to use the simple expression (16) even if the healthy population is not approximately constant over the times where W is non-vanishing. The key question then is: how severe is the bias introduced by using the simple formula under realistic non-constancy of the healthy population?
In order to explore this issue systematically we consider a non-constant healthy population given by H(t) = H 0 e αt . This has a conveniently tunable degree of failure to conform with the constancy assumption required for (16). When α = 0 we have a constant number of healthy individuals, while a value of α = ln(x) means the population grows by a factor of x in one year. We can provide a survival function S U (t) for time measured in years, roughly inspired by the ODn progression on the BED assay, by specifying f U to be a Weibull distribution with scale parameter l = 0.57 and shape parameter k = 8. We now numerically evaluate the denominator of (15) and compare it to the denominator of the simple formula (16). Note that this bias calculation is independent of the actual time dependance in I.
In Figure 3, the bias in the naive denominator (reported as a fraction of the unbiased denominator) is shown as a function of α, reported as the annual percentage growth. Note that for a population growing at 4% per year, the bias is about 1.1%. As we shall see later when analyzing a slightly more complex model of a TRI, this is small compared with the statistical error that arises as a result of using realistic sample sizes. Thus, bias arising from the non-constancy of the healthy population is not a key concern unless there is very dramatic variation in H. The bias calculation demonstrated here is also applicable to the more complex model that now follows.

Modeling Non-progressing Individuals
A complicating factor for the BED assay is the fact that a small number of individuals utterly fail to progress beyond any practical ODn threshold used to define recency. This is due to individual variation in biochemical details such as immune response, for example. The non-progression phenomenon leads to a long term accumulation of apparently recently infected individuals, as classified by the TRI, even though many of them have been infected for a long time. There is currently no consensus on how to handle this complication. We now generalize the previous analysis to the situation where some individuals fail to progress to the over-threshold category. Consider the simple model captured in Figure 2B. At the moment of infection, individuals transition from the healthy population to either a non-progressing population (NP) or to a progressing under-threshold population (P U ). The probability of non-progression is P NP , and hence the probability of progression is 1 − P NP . Those individuals in P U wait for a stochastic delay after which they move into the progressing over-threshold category P O . In the previous model, f U|A was the distribution of waiting times governing the transition, but since the waiting times for non-progressing individuals are infinite, f U|A cannot be normalized. Therefore, in order to specify the transition times from P U to P O in terms of a normalized density, we introduce the density of waiting times in the state of being a progressor and under the threshold, conditional on being a alive, denoted by f PU|A . Then S U|A (t), S PU|A (t) and P NP are related by The difficulty is that the TRI will classify as 'recently infected' all the individuals in the NP and P U categories even though some potentially large number in NP are long infected. This can systematically be addressed by the following two key steps.
Firstly, we assume that the same survival function S A is applicable to both progressing and non-progressing individuals. This is true if the differences between individuals which account for progression versus non-progression do not translate into significant differences in post-infection survival. This assumption has also been made, at least implicitly, in previous work on use of the BED assay for estimating HIV incidence (for example, see [11] for analysis of [8]). Its applicability should in principle be tested, but we are unaware of any evidence suggesting that it is false.
Secondly, we introduce two artificial categories by separating the non-progressing population into 'recently infected' (NP R ) and 'long infected' (NP L ) sub-populations. Individuals entering the NP R sub-population are assigned a waiting time drawn from f PU|A after which they transition to the NP L category. Note that this is a book-keeping device used for convenience and, unlike the assumption about survival, does not rely on any property of disease progression. It is now possible to provide a sensible definition for the class of 'recently infected individuals' (R) which has a population given by Note that, since both P U and NP R now have the same exit waiting times, the distribution of waiting times for R is given by f R|A = f PU|A , with corresponding survival function S R|A . These two steps lead to the model in Figure 2C. It is now possible to recycle our preceding analysis and write down expected counts in these new classes. Survival in the state of being simultaneously alive and recently infected, is given by S R (t) = S A (t)S R|A (t), and hence for the progressing populations we obtain and Note the similarity with expressions for U (t) and O(t) in the basic model. For the non-progressing populations we obtain and For convenience we define and note that and, more importantly, These equations express the symmetry between the progressing and non-progressing sub-populations of Figure 2C. Substituting (19) and (21) into (18), we can write It is appropriate to use a weighting scheme analogous to the one used in the basic model since W (s, t) is now proportional to the probability that individuals are alive and classified as recently infected at time t if they become infected at time s, regardless of whether they are progressors or non-progressors. Then the weighted incidence, denoted I W , is given by The populations of under -threshold (U) and over -threshold (O) individuals are related to the populations defined in Figure 2C by and Using the above two equations and (18) and (25) this means that the population of recent infections is related to the under-threshold and over-threshold populations by Performing the same expansion technique as before and assuming a slowly varying healthy population gives the simple expression This expresses the incidence in terms of the calibration parameters E [τ R ] and ǫ (equivalently P NP ), and the state variables H, U and O.
All that has changed, as a result of allowing non-progressors into the model, is the shift in the numerator from U to R = U − ǫO. The same bias calculations as before apply immediately, but there is an increase in statistical sensitivity. This can be understood by noting that the gross error in U becomes the dominant part of a fractional error in R and that R is smaller than U .

Statistics and Calibration
The population models of the preceding section are expected to be in ever closer correspondence to a real population as the population size increases.
It should be noted that we do not address issues relating to experiment design and selection bias.

Statistical Fluctuations
As noted in the preceding section, in a relatively established epidemic where the smallest class is U, the major source of fluctuations in I est is N U . Crudely speaking then, we can estimate the reproducibility by blaming all the statistical uncertainty on the measurement of N U , which has a standard deviation σ(N U ) = N P U (1 − P U ) ≈ √ N P U . This leads to the 'back of the envelope' formula for the relative standard deviation given by σ(I est ) However, the trinomial counts in the estimator (33) fluctuate and are correlated since they are constrained to add up to the sample size N . In the Appendix we demonstrate how these counts can be modeled by two independent draws (α 1 and α 2 ) from a standard normal distribution. We obtain a particular incidence estimator by inserting the counts, as functions of α 1 and α 2 into (33). Organizing the resulting expression into a natural expansion in powers of 1/ √ N gives where γ 1 = P H (1 − P H ) and The leading term, is just the estimator evaluated at the expectation values of the counts. The O(1/N ) term, which we have omitted, contains a term proportional to α 2 1 . This means that there is in principle a finite sample size bias in the estimator, which is however suppressed by O(1/N ) relative to the dominant term, as is borne out by numerical calculations in Section 4. The retained sub-leading term (of order O(1/ √ N)) is distributed according to a bivariate normal distribution. Thus, to this order, there is no bias and we can read off the likelihood of observing a value of I est , where f (·) is the standard normal density and Comparison with numerical simulation suggests that this approximate result is essentially indistinguishable from the exact answer in the regime of realistic values for N and the proportions (P H , P U and P O ), given that in practice one uses the sample proportions as estimates of the population proportions. Figure 4 plots the relations (34) and (39) for different values of P NP and fixed values of P H , P O and P U . Comparison with Figure 3 confirms that the truncation of the expansion of the healthy population to a constant term (the crucial step in obtaining the simple incidence relation on which the estimator is based) leads to a bias that is small compared to realistic statistical errors. The close correspondence between the two curves suggests that the simple formula (34) should be sufficiently accurate to choose a sample size for an intended study.

Calibration
Aside from the sample counts N H , N U and N O , all other quantities in relations of the kind derived in the previous section, such as E [τ R ] in the estimator, should be regarded as parameters that need to be estimated independently of a cross-sectional survey used to infer incidence. Even in the more general case, where the healthy population is allowed to vary considerably over the time when the weighting function is non-vanishing, calibration consists only of estimating P NP and expressions of the form E [τ n R ]. We have already remarked that for the BED assay it will probably not be possible to obtain reasonable estimates of E [τ n R ] for values of n other than 1 and that it appears that only this term is really needed for practical purposes.
Note that if one wishes merely to estimate trends in incidence, as opposed to absolute incidence values, then it is not necessary to have an estimate for E [τ R ] at all, since it is just an overall factor. If the overall scale of incidence estimates is to be known, considerable effort will need to be invested in the estimation of E [τ R ]. Note that this is the mean time progressing individuals spend under the threshold, with mortality accounted for. However, surveys conducted at different times will not yield comparable values of I est ∝ (U − ǫO)/H unless ǫ (equivalently P NP ) is known with some accuracy, since it appears in one of two terms in the numerator. Consider two surveys which use the same point estimate where ǫ 0 is the real value and ǫ 1 is the error due to methodological and statistical factors. The first survey obtains values of N O . This leads to an estimate of the difference between the two incidences of Knowledge of the exact value ǫ 0 leads to from which we see that the error in ∆I est , due to the error in ǫ, is The direction and magnitude of error depend in detail on many factors, such as population renewal and long term post-infection survival. While it is not possible to summarize all the effects that may be produced by imperfect estimation of P NP , in Section 4 we conduct a number of numerical simulations which demonstrate the kind of bias that may arise.

Numerical Simulations
In this section we briefly outline two numerical simulations. The first serves to test the accuracy of the bias and standard deviation estimates derived from the truncation of the systematic expansion of the stochastic estimator I est (α 1 , α 2 ). For each of 10,000,000 iterations, two standard normal variables were drawn. Counts for healthy, under-threshold and over-threshold sub-populations, within a sample of size N = 5000, were generated according to the procedure provided in the Appendix. This algorithm incorporates the assumption that the trinomial distribution of the sample proportions can be approximated as normal, but involves no truncation of the estimator in powers of 1/ √ N . The resulting ensemble of point estimates I est produced suitably converged estimates for the mean and standard deviation.
The entire procedure of the preceding paragraph was repeated for values of P O ∈ [0.1, 0.5]. For each value of P O , the value of P U was varied to produce an incidence in the range [0.01, 0.2]. The average fractional discrepancy between the mean incidence and I est (0, 0) was 0.0003, with values ranging from 0.00001 to 0.0008, confirming that the intrinsic sampling bias is of order 1/N . The average fractional discrepancy between the observed standard deviation of the incidence and the approximate expression (39) was 0.0006, with values ranging from 0.00004 to 0.001, which is also consistent with the truncation at O(1/ √ N ). A second simulation demonstrates the use of the simple expression (33) for incidence estimation. Arrival times of new infections were generated according to a non-homogeneous Poisson process with intensity given by λ(t) = H(t)I(t) as described in more detail below. Newly infected individuals were initially classified as under the recency threshold of a TRI. A fraction 1 − P NP progressed to the over-threshold category according to waiting times generated by f R|A . Weibull functions were used for the waiting time distributions f R|A and f A . Unique individuals were drawn from the population at intervals, to produce counts N H , N U and N O , and hence estimates for incidence.
To demonstrate the incidence estimation process, a 50 year population scenario was produced. Figure 5 shows the Weibull distributions used for f A and f R|A (top) and the corresponding survival curves S A and S R|A (bottom). The Weibull shape and scale parameters for the distributions were chosen to give approximately realistic values for the mean and standard deviations for the window period and infected life expectancy, as detailed in Table 1. The healthy population was set to H(t) = 100, 000 + 5, 000t, with t measured in years. The incidence was set at 0.01 (hazard per person per year) for the first ten years, climbing linearly to 0.1 over the next ten years, then remaining at this high level for a further ten years, followed by ten years of linear decline to 0.03 and maintained at this level for the last ten years of the simulation.  Figure 6 shows output from this simulation. The input incidence parameter is indicated as the dotted instantaneous incidence curve. A sample of 5, 000 individuals was surveyed every year, and an incidence estimate was produced using the simple estimator (33) with exact values of E [τ R ] and P NP , i.e., assuming perfect calibration. These point estimates are indicated as estimated incidence values, using '+' symbols. The combined effects of the previously noted time convolution, in I W , as well as stochastic departure from means in the simulated population, make the input incidence parameter an unrealistic target for simulated incidence measurements. Thus, the solid weighted incidence line has been displayed, which uses full knowledge of all population members' classification into H, U or O, inserted into (28) with full knowledge of the denominator, (both the non-constant H(t) and the exact S R ). This is essentially all that the incidence estimation algorithm can be asked to reproduce. A two standard deviation envelope around the weighted incidence line, calculated from (39) using knowledge of the full population, is shown as two dashed lines.
In Section 3 it was shown that incidence trends can be extracted without E [τ R ] calibration, while an estimate for P NP is vital. We now explore the extent to which the accuracy of the estimate of P NP affects the ability to determine a trend in incidence.
Population fractions for H, U and O were extracted at six times from the population simulation described above and are shown in Table 2. Four instances of incidence trend estimation were simulated by selecting the time intervals (15, 20), (20, 30), (30, 35) and (40, 50). We considered the trends that would be observed if incidence were measured at the beginning and end of each of these intervals using (37). In order to focus on the bias introduced by imperfect estimation of P NP , rather than sample size effects, we assumed perfect knowledge of P H , P U and P O . For each of these intervals, we calculated an incidence estimate at the beginning and end, as a function of the estimated value of P NP (the true value being 0.05), assuming E [τ R ] is known exactly. We also calculated the estimated fractional change in incidence over each interval. Note that the fractional change does not depend on E [τ R ]. The results are shown in Figure 7, where the four intervals (15, 20), (20, 30), (30, 35) and (40, 50) are referred to as scenarios A, B, C and D, respectively.
In each case, the effect of the error in the estimation of P NP is quite different, as can be understood by considering how (43) is impacted by the system history. Note that case B and case D both simulate intervals over which incidence is approximately constant, but the impact (on the estimated incidence change) of incorrect estimation of P NP does not even agree in sign. At a time of 40 years, the incidence estimate becomes negative when the P NP estimate is higher than 0.0882. This breakdown of the model results in the divergence of the fractional change in estimated incidence over the interval in scenario D. In short, incorrect estimates of P NP can lead to the fundamental breakdown of the inference scheme. This makes sense, as P NP impacts the long term accumulation of individuals in the P O category.

Discussion and Conclusion
We have presented a detailed analysis of relations between recent incidence in a population and counts of 'recently infected' individuals. These are in principle complex convolutions involving the epidemiological history as well as all information about the distribution of waiting times in the recently infected category. When the healthy population undergoes realistically modest variation on the time scale of the definition of recency implied by the TRI, we obtain simplified forms which incur very little error. The simplified relations form the basis of estimators which are shown to have considerably more variance than bias under realistic demographic and epidemiological assumptions.
Noting that the assumptions of our model are the least restrictive of any BED based HIV incidence estimation method of which we are aware, we now consider its limitations. We have only modeled one direction of progression of individuals from an experimentally defined state of 'recent' infection to 'non-recent' infection. The reverse apparently occurs for BED optical density in some terminal stage AIDS patients. This process constitutes a substantial complication, and further work is required to investigate how it may be incorporated into an analytical model of the kind developed here. It may be worth exploring previous suggestions [8] to use additional information, from questionnaires or other assays, to remove end-stage patients from the observed recent count. We have also not considered the possibility that calibration parameters are functions of time, for example as a result of substantial vaccine uptake in a population. Even more subtle is a point raised under calibration, namely that imperfect estimation of the non-progressing fraction of the population leads to a complex bias in incidence estimates. This bias is dependent on factors not even present in the incidence estimator, such as long term survival post infection.
A key observation is that, for the purposes of estimating incidence from a TRI, there is no fundamental obstacle posed by having a known fraction of individuals fail to progress over the recency threshold, as long as their distribution of survival times from infection is the same as in the progressing population. In fact, an accurate estimate of this non-progressing fraction alone, is sufficient (and necessary) to infer trends in incidence. This fraction could possibly be estimated for the BED assay from historical records, since there are many viable samples in storage with supporting clinical information indicating long-infected status. However, as demonstrated in the calculations of Section 4, a suitably large error in the estimate of P NP can render TRI based incidence estimates meaningless. A calibration of the mean finite waiting time is required in order to estimate absolute values of incidence. Prospective follow-up is probably the only practical way to estimate this parameter.
In contrast to our model, which has only two calibration parameters, the model of McDougal et al. [8] appears to have three (sensitivity, short-term specificity and long-term specificity). In a separate short note [11] we demonstrate that, under their own assumptions, these parameters can be reduced to ours. This has two advantagesour parameters are easier to calibrate and assuming independence of their parameters would lead to incorrect error estimates.
Besides the explicit assumptions noted, the analysis presented here is quite general. Tests for recent infection continue to be of interest, and new assays are likely to be developed both for HIV infection and other important diseases. In summary, we have presented a simple incidence estimator and a detailed analysis of its likelihood function, which can inform design of appropriate calibration studies and cross-sectional incidence estimation surveys, and can also form the basis of systematic inference algorithms for processing the data obtained from such surveys.
where the counts N X = P X N fluctuate trinomially around their meansN X = P X N , with standard deviations σ X = N P X (1 − P X ), or, alternatively, the realized sample proportions P X fluctuate multinomially around their means P X , with standard deviations σ P X = P X (1 − P X )/N . In order to account for fluctuations in the sample counts, as well as their correlations, we express the three counts as the result of two independent random draws. We assume the counts are sufficiently large so that binomial distributions can be approximated by normal distributions-which will be the case if the survey is to have any reasonable accuracy.
• Draw a number α 1 from a standard normal distribution and set where Defining γ 1 allows us to keep track of powers of 1/ √ N .
• Draw a number α 2 from a standard normal distribution and set and where are the probabilities of an individual being in U and O given that they are not in H, and is the standard deviation of N U (and N O ) subject to α 1 having been determined, so that The sample proportions can now be written as: Note that this procedure works independently of the order in which the values of N H , N U and N O are assigned from the random α's. We can explicitly insert this decomposition into the expression for I est to obtain an expression for the estimator in terms of two independent draws from a standard normal distribution. Keeping only the O(1) and O(1/ √ N ) terms explicit leads to where Two key observations immediately follow: • The proposed I est is a consistent estimator in the sense that terms with bias are at least O(1/N ) suppressed relative to the dominant term.
• To order O(1/ √ N ), I est is Gaussian with relative variance

Bias in Simple Incidence Formula
Annual percentage growth of susceptible population Relative error in simple incidence Figure 3: Fractional error in the simple incidence relation (16) versus the full relation (15), as a function of growth rate of the healthy population, quoted as a percentage annual growth. The scenario is defined by: S(t) = S(0)e αt and f U is a Weibull distribution with scale parameter l = 0.57 and shape parameter k = 8. The parameter α is varied to produce deviation from a constant healthy population (α = 0) in which limit equation (16) is exact.    Table 1.

Monte Carlo Experiment
Years Incidence Instantaneous incidence Weighted incidence (total population) Estimated incidence (sample) Two standard deviation envelope Figure 6: Full stochastic simulation of population with epidemic, individual survival times, and annual sampling of 5, 000 individuals. The healthy population was set to H(t) = 100, 000 + 5, 000t with time measured in years. The instantaneous incidence parameter is the dotted curve. The target of the estimates is the weighted incidence (solid line), which was calculated explicitly as per (28) from all the known inputs. This is flanked by a two standard deviation envelope (dashed lines). Simulated estimated incidence (+ symbols) were obtained by using sample counts in the simple estimator (33). The calibration parameters E [τ R ] and P NP were assumed to be known exactly.

Incidence estimate
Incidence at 40 years Incidence at 50 years Figure 7: Absolute incidence estimates and estimated fractional incidence changes for four pairs of successive times from Table 2.