INTRODUCTION
General overview
In 1993 the World Health Organization (WHO) declared tuberculosis (TB) a global public health emergency, being the only disease thus far to warrant that designation. Although hospitals have been established and chemotherapy has been developed to combat TB, bringing considerable reduction in incidence to developed nations, historical data calculated by the WHO indicate that there have not been great effects on the global problem since the time of Koch. Currently, TB is responsible for more human deaths than any other single infectious agent, representing 26% of all preventable deaths and 7% of all deaths [Reference Pio and Chaulet1].
TB resumption has been attributed to several factors, such as the increase in drug resistance, the HIV/AIDS pandemic (at the beginning of the 1980s), the increase of injecting drug users, changes in social structure, the increase of immigrants from high prevalence nations to developed ones, the ageing of the world's population, the active transmission in environments of human accumulation (e.g. prisons, hospitals, homeless shelters), and the dismantling of health-care systems [Reference Ducati2]. Although TB became a re-emerging disease in European and North-American nations, TB is not an emergent nor re-emerging public health problem in developing countries such as Brazil, but rather a long lasting one [Reference Ruffino-Netto3].
In order to facilitate the comprehension of the various components involved in the interaction between these factors, Ruffino-Netto [Reference Ruffino-Netto4] proposed an expression which reflects the TB burden, represented by the components social inequality, prevalence of HIV-positive individuals, percentual default of treatment, prevalence of primary resistance plus acquired resistance, migration, age of the population, adequate health services, directly observed treatment short-course, educational level, nutrition level, human resources for TB control and degree of political participation of the population.
From the historical data of the observed numbers of TB cases, especially in developed countries, we observe a trend of declining incidence starting at the beginning of 1960 up to 1980, where there was a change in this trend. During the period between 1980 and 1990, we observe an increase in the incidence rates of TB cases; after 1990, we again observe a trend of declining incidence. That is, we have the presence of two change-points in the rates of TB, especially for developed countries.
The case of New York City
The incidence (notification cases) of TB disease in New York City (NYC) between 1970 and 2000 presents three trends (see Table 1): a first period (1970–1979) where the trend of declining incidence was probably associated with good control programmes; a second period (1979–1992) where there was an increase in incidence rates [Reference Coker5, Reference Wallace6], possibly associated with a systematic dismantling of public-health infrastructure of control programmes, social disruption (including homelessness, drug abuse, poverty and housing overcrowding), and mainly caused by the HIV epidemic; a third period (1992–2000) where again there is a decline in incidence rates. It is important to remember the many factors associated with this third period, i.e. implementation of directly observed therapy, broader chemotherapy regimens for patients with TB or suspected multidrug-resistant TB and improved therapeutics for the care of HIV-infected individuals [Reference Paolo and Nosanchuk7].
Table 1. Number of tuberculosis cases in NYC from 1970 to 2000
Source: New York City Department of Health and Mental Hygiene, Bureau of Tuberculosis Control Information Survey.
Table 1 gives the yearly numbers and the accumulated numbers of TB cases in NYC. From Table 1, we see decreasing numbers of TB cases from 1970 to 1978, where there is a minimum. From 1978 to 1992, we observe increasing numbers of TB cases, where there is a maximum number of cases in 1992. From 1992 to 2000, we observe decreasing numbers of cases. That is, we have two change-points for the numbers of cases (see Fig. 1). It is interesting to note that the use of powerful antiviral drugs against HIV commenced around 1990.
Fig. 1. Number of tuberculosis cases in New York City, 1970–2000.
To model the number of TB cases in NYC during the period 1970–2000, we consider the use of a point process to count the numbers of TB cases in each year starting in 1970. In this way, we considered a stratified sample of size n=6721 representing 10% of the total number of TB cases (259 cases in 1970, 257 cases in 1971, 227 cases in 1972 and so on), where for each year, we used an uniform distribution to have the times (in days) for the occurrence of each case since 1 January 1970 until 30 December 2000, i.e. with a total time of observation equal to T=11 323 days.
For this dataset, we assume a non-homogeneous Poisson process (NHPP) in the presence of two change-points considering a Bayesian approach using Markov chain Monte Carlo (MCMC) methods (see e.g. Gelfand & Smith [Reference Gelfand and Smith8]). The use of Bayesian methods has been considered by many authors for analyses of homogeneous or non-homogeneous Poisson processes in the presence of change-points (see e.g. Raftery & Akman [Reference Raftery and Akman9] considering the presence of a change-point in homogeneous Poisson processes, or Ruggeri & Sivaganesan [Reference Ruggeri and Sivaganesan10] considering any number, random or fixed, of change-points in NHPP assuming power-law intensity functions).
The paper is organized as follows: in the Methods section, we introduce the likelihood function and a Bayesian analysis for the model; in the Results section, we introduce the analysis for the NYC data, and finally, in the Discussion, we present some concluding remarks.
METHODS
The likelihood function
Let N(t) be the cumulative number of TB cases that are observed during the interval (0, t) and assume that N(t) is modelled by a NHPP with intensity function λ(t)=dm(t)/dt=dE[N(t)]/dt, where m(t) is the mean value function (see e.g. Cox & Lewis [Reference Cox and Lewis11]). Different parametrical forms could be assumed for the intensity function λ(t) (increasing, decreasing, bathtub shape, unimodal, among many others, see e.g. Musa & Okumoto [Reference Musa and Okumoto12] or Muldholkar et al. [Reference Mudholkar, Srivastava and Friemer13]). We assume power-law processes (PLP) in the presence of two change-points with intensity function for the overall process given by
where θ=(α1, α2, α3, β1, β2, β3, ζ1, ζ2).
Equivalently, letting m j(t)=m(t|θj), the corresponding mean value function is given by
where m 1(t)=(t/α1)β1, m 2(t)=(t/α2)β2 and m 3(t)=(t/α3)β3.
Observe that the intensity function λj(t) in equation (1) is constant for βj=1, decreases for βj<1 and increases for βj>1, j=1, 2, 3. This process is related to the Weibull probability model [Reference Kuo and Yang14] (αj, βj), j=1, 2, 3.
Assuming that the data are observed up to a total time T, where the epochs of occurrence of cases are denoted by t i, i=1, … , n, 0<t 1 <t 2< … 
<…<t n<T, the likelihood function for θ in the presence of two change-points ζ1 and ζ2 is given by
where λj(t) is given in equation (1) and m j(t) is given in equation (2) for j=1, 2, 3. To justify the likelihood function [equation (3)], observe that N(s+t) – N(s) given θ has a Poisson distribution P(m(s+t|θ) – m(s|θ)) for t>0 and independent increments [Reference Raftery and Akman9]. Thus, the sampling distribution for the between occurrence times, say U i, has density 
, 
, and so on. In this way, we obtain the likelihood of the data D T={n; t 1, …, 
, …, t n, T} in the presence of two change-points. Moreover, observe that homogeneous Poisson processes in the presence of one change-point is a special case of equation (3) [Reference Raftery and Akman9].
A Bayesian analysis
For a Bayesian analysis of the PLP with intensity function given in equation (1) in the presence of two change-points ζ1 and ζ2, we assume uniform prior distributions for αj and βj given by
for j=1, 2, 3, where a j, 
are known hyperparameters, b 11 and b 13 are assumed to be equal to 0 and b 21 and b 23 are assumed to be equal to 1 in order to have decreasing intensity functions in the intervals 0<t<ζ1 and ζ2<t<T; b 12 is assumed to be equal to 1 to have increasing intensity function in the interval ζ1<t<ζ2, and a j and b 22 are assumed to have large values (non-informative prior distributions for αj, j=1, 2, 3). We also assume uniform prior distributions for the change-points ζ1 and ζ2, given by
where c ℓ and d ℓ are known hyperparameters, ℓ=1, 2 is assumed to have ζ1<ζ2. We further consider prior independence among the parameters.
The joint posterior distribution for θ is given by [Reference Box and Tiao15]
where D T={n; t 1, …, t n; T}, 0<αj<a j, 0<βj<b j, c ℓ<ζl<d ℓ, j=1, 2, 3 and ℓ=1, 2.
To simulate samples for the joint posterior distribution [equation (6)], we could consider standard MCMC methods such as the Gibbs sampling algorithm [Reference Gelfand and Smith8] or the Metropolis–Hastings algorithm [Reference Smith and Roberts16]. In this case, we need all full conditional posterior distributions Π(θj|θ(j), D T), j=1, 2, …, K and θ(j)=(θ1, …, θj−1, θj+1, …, θK). A great computational simplification is given by WinBugs software [Reference Spiegelhalter, Thomas and Best17], where we only need to specify the joint distribution for the data and the prior distributions for the parameters.
RESULTS
For a Bayesian analysis of the NYC TB data, we assumed the uniform prior distributions [equation (4)] for αj and βj with a j=100, j=1, 2, 3, b 21=b 23=1, b 11=b 13=0 (related to decreasing functions), b 12=1 and b 22=10 (related to an increasing function between the first and second change-points). We also assumed prior distributions [equation (6)] for the change-points ζ1 and ζ2 with c 1=2558, d 1=4383, c 2=7671 and d 2=9131 (number of days since 1 January 1970). This choice of prior distributions, especially for the change-points ζ1 and ζ2 are based on medical knowledge of the epidemic, or if it is known that the first change-point is between 1977 and 1982 and the second change point is between 1991 and 1995. We are assuming in the prior distribution [equation (5)] that the two intervals do not overlap. Using WinBugs software and considering a burn-in sample of size 40 000, we simulated a Gibbs sample of size 100 000 choosing every 50th sample for each parameter to have approximately uncorrelated samples, i.e. obtaining a final Gibbs sample of size 1200 to get the posterior summaries for each parameter. The WinBugs code is given in the Appendix.
Table 2 gives the posterior summaries for each parameter. Convergence of the Gibbs sampling algorithm was monitored by checking the plots of the simulated samples for each parameter to verify if a stationary distribution was obtained by the 1200 simulated Gibbs samples and also using other existing methods to check the convergence of the Gibbs sampling algorithm [Reference Gelman and Rubin18].
Table 2. Posterior summaries for the parameters
Figure 2 shows the plots of the marginal posterior distributions for the change-points ζ1 and ζ2 approximated by the simulated Gibbs samples. Inferences for the change-points are of great interest to epidemiologists.
Fig. 2. Marginal posterior distributions for the change-points. (a) ζ1; (b) ζ2.
Considering the Monte Carlo estimators for the posterior means of α1, α2, α3, β1, β2, β3, ζ1 and ζ2 given in Table 2, we obtain Bayesian estimators for the mean value function m(t) given by equation (2), i.e.
Table 3 gives Monte Carlo Bayesian estimators for m(t) based on the 1200 simulated Gibbs samples and the observed accumulated numbers of TB cases for each year. Figure 3 shows the plot of the estimated mean vale function and the observed accumulated number of TB cases against the years (in days). We observe a good fit for the PLP in the presence of two change-points.
Fig. 3. Mean value function. —, Estimated; ○, observed.
Table 3. Estimators for the mean value function and observed accumulated numbers
DISCUSSION
There was a great increase in TB prevalence in NYC during the 1980s and at the start of the 1990s, with a peak of 3811 cases in 1992. In 1978 a very low number of TB cases in NYC (1307 cases) was observed, following a long period of decreasing numbers in the prevalence of the disease. We observe that the proposed model was well fitted to the data of TB cases in NYC, and it is straightforward to implement this model in WinBugs.
The presence of more than one change-point is common in many applications of medical counting data. Considering the NYC TB data of Table 1, the results obtained in the present study could be easily extended for other epidemiological datasets, where we could have the presence of a finite number of change-points. In this way, the likelihood function [equation (6)] could be easily generalized to accommodate more than two change-points. We usually have great difficulty in obtaining classical inference results for the parameters of NHPP in the presence of change-points and the use of MCMC methods is a suitable way of obtaining Bayesian inferences for this family of models. Using WinBugs software greatly simplifies obtaining the posterior summaries of interest.
Other parametrical forms for the intensity functions [equation (1)] could be considered in place of PLP. In this case we could consider other usual intensity functions commonly used in software reliability studies, e.g. Gompertz growth, logistic growth, etc. [Reference Musa, Iannino and Okumoto19].
APPENDIX
The WinBUGS code used to fit the Bayesian model is given below: 
In this WinBUGS code, a 1, a 2, a 3 are the hyperparameters of the uniform prior distribution of α1, α2 and α3, respectively; b 1, b 2, b 3 are the hyperparameters of the uniform prior distribution of β1, β2 and β3, respectively [see equation (4)]; and c 1, d 1, c 2 and d 2 are the hyperparameters in equation (5). These hyperparameters are declared in the last line of the WinBugs code.
DECLARATION OF INTEREST
None.
