Masoud Amiri
email: masoud.amiri@yahoo.com

Fanny Janssen

Anton E Kunst

The decline in stomach cancer mortality: exploration of future trends in seven European countries

Keywords: Stomach cancer, mortality, trends, Europe, projections

Mortality from stomach cancer has fallen steadily during the past decades. The aim of this paper is to assess the implication of a possible continuation of the decline in stomach cancer mortality until the year 2030. Annual rates of decline in stomach cancer mortality from 1980 to 2005 were determined for the Netherlands, United Kingdom, France, and four Nordic countries on the basis of regression analysis. Mortality rates were extrapolated until 2030, assuming the same rate of decline as in the past, using three possible scenarios. The absolute numbers of deaths were projected taking into account data on the ageing of national populations. Stomach cancer mortality rates declined between 1980 and 2005 at about the same rate (3.6-4.9 percent per year) for both men and women in all countries. The rate of decline did not level off in recent years, and it was not smaller in countries with lower overall mortality rates in 1980. If this decline were to continue into the future, stomach cancer mortality rates would decline with about 66 percent between 2005 and 2030 in most populations, while the absolute number of stomach cancer deaths would diminish by about 50 percent. Thus, in view of the strong, stable and consistent mortality declines in recent decades, and despite population ageing, stomach cancer is likely to become far less important as a cause of death in Europe in the future.

Introduction

Stomach cancer is the fourth most frequent cancer [START_REF] Crew | Epidemiology of gastric cancer[END_REF][START_REF] Prinz | H pylori and gastric cancer: shifting the global burden[END_REF] and the world's second leading cause of cancer mortality [START_REF] Crew | Epidemiology of gastric cancer[END_REF] . However, it has been fallen throughout Europe during the past decades in terms of both incidence and mortality rates 3 4 , mainly as a result of remarkable improvement of living conditions in European societies [START_REF] Bosetti | Cancer mortality in the European Union, 1970-2003, with a joinpoint analysis[END_REF][START_REF] Lepage | Trends in incidence of digestive cancers in France[END_REF][START_REF] Levi | Continuing declines in cancer mortality in the European Union[END_REF][START_REF] Jemal | Cancer statistics[END_REF] . Efforts to reduce global cancer disparities begin with an understanding of geographic patterns in cancer incidence, mortality and prevalence, by studies such as GLOBOCAN [START_REF] Parkin | Estimating the world cancer burden: Globocan 2000[END_REF][START_REF] Parkin | Global cancer statistics[END_REF] , EUROCARE [START_REF] Francisci | The cure of cancer: a European perspective[END_REF] , and Five Continents databases [START_REF] Kamangar | Patterns of cancer incidence, mortality, and prevalence across five continents: defining priorities to reduce cancer disparities in different geographic regions of the world[END_REF] . Survival increased and mortality decreased through the combination of earlier detection, better access to care and improved treatment [START_REF] He | Recent trends of cancer in Europe: a combined approach of incidence, survival and mortality for 17 cancer sites since the 1990s[END_REF] . There has also been a concomitant change in lifestyle and environmental exposures over successive generations [START_REF] Fitzsimmons | Trends in stomach and pancreatic cancer incidence and mortality in England and Wales, 1951-2000[END_REF] , including changes in exposures to risk factors in early life [START_REF] Aragones | Gastric cancer in the European Union (1968-1992): mortality trends and cohort effect[END_REF] .

Several studies conducted projections of future trends in stomach cancer mortality in European countries [START_REF] Levi | Continuing declines in cancer mortality in the European Union[END_REF][START_REF] Kamangar | Patterns of cancer incidence, mortality, and prevalence across five continents: defining priorities to reduce cancer disparities in different geographic regions of the world[END_REF][START_REF] Fitzsimmons | Trends in stomach and pancreatic cancer incidence and mortality in England and Wales, 1951-2000[END_REF][START_REF] Aragones | Gastric cancer in the European Union (1968-1992): mortality trends and cohort effect[END_REF][START_REF] Quinn | Cancer mortality trends in the EU and acceding countries up to 2015[END_REF][START_REF] Coebergh | Kanker in Nederland[END_REF] . Mortality rates were generally expected to decline further. Based on the cancer mortality trends in the European Union until 2000, one study projected a further fall by 11% in age-standardized cancer mortality from 2000 to 2015 [START_REF] Quinn | Cancer mortality trends in the EU and acceding countries up to 2015[END_REF] . A Dutch study also projected a substantial decline in stomach cancer mortality until 2015, based on trends until 2000 [START_REF] Coebergh | Kanker in Nederland[END_REF] . Similarly, an Irish study predicted mortality from stomach cancer to fall further in the near future, although at a slower rate than in the recent past [START_REF] O'lorcain | Mortality predictions for esophageal, stomach, and pancreatic cancer, Ireland, up to 2015[END_REF] .

European studies included projections of the absolute number of cancer deaths, thereby taking into account the rapid ageing of European populations in the forthcoming decades [START_REF] Fitzsimmons | Trends in stomach and pancreatic cancer incidence and mortality in England and Wales, 1951-2000[END_REF][START_REF] Aragones | Gastric cancer in the European Union (1968-1992): mortality trends and cohort effect[END_REF][START_REF] Quinn | Cancer mortality trends in the EU and acceding countries up to 2015[END_REF] . Population ageing alone may induce a strong increase in the absolute number of incident cases and deaths for many types of cancer. However, in the case of stomach cancer, this demographic effect may be offset by an opposite epidemiological effect if the recent decline in age-specific mortality rates were to continue in the future.

The aim of this study was to assess, for seven European countries, the implications of a possible continuation of the recent declines in stomach cancer rates. Future trends were evaluated in terms of not only mortality rates, but also absolute numbers of death, thereby taking into account the ageing of national populations.

Projections of future mortality trends should be based on a careful analysis of trends in the past [START_REF] Janssen | The choice among past trends as a basis for the prediction of future trends in old-age mortality[END_REF][START_REF] Wilmoth | Demography of longevity: past, present, and future trends[END_REF] . We therefore started the analysis with a description of trends in stomach cancer mortality over a long period, from 1980 up to 2005, in order to check whether mortality declines continued undiminished until recent years. In addition, we assessed whether the rate of decline was similar for all seven countries and both sexes, despite differences in overall levels of stomach mortality rates. We thus aimed to identify a stable and consistent long-term mortality trend that could serve as a basis for projections into the future. In addition, we formulated three scenarios to express uncertainties regarding future mortality trends.

Material and Methods

historical national data for England and Wales and France were transformed into data for United Kingdom and France metropolitan, by applying age and sex-specific correction factors, based on the comparison of stomach cancer mortality rates in 1995.

To describe trends in stomach cancer mortality between 1980 and 2005, sexspecific age-standardized mortality rates were calculated using direct standardization, taking the average national male and female populations by five year age groups over 2000-2005 (2000-2001 for Denmark) as the standard populations. In addition, we estimated sex-and country-specific annual mortality changes (%) over the period 1980-2005 by means of age-period log-linear regression analysis. The dependent variable was the number of deaths, with the person-years at risk (estimated by the midyear population) as offset variable. We used age by five year age groups (from 0 to 80+) (categorical) and single calendar year (continuous) as independent variables. Annual mortality changes (%) were calculated by this formula: 100*(exp(b)-1), in which b is the parameter estimate of the calendar year variable.

To project age, sex, and country-specific mortality rates until 2030, we assumed that the mortality trends in the past 25 years (from 1980 to 2005) would continue for the coming 25 years. Therefore, we applied the estimated annual changes in mortality in the period 1980-2005 for all age groups combined to average national age-specific mortality rates in 2000-2005.

To express uncertainties regarding future mortality trends, we formulated three scenarios. In the reference scenario, the future rate of change was assumed to be equal for all countries and therefore was calculated as the unweighted average of the rates of change of the seven countries. In the pessimistic scenario, we assumed the future rate of change to be equal to the slowest observed rate of change observed among all countries. In the optimistic scenario, we took the highest rate of change observed in the past. All three scenarios distinguished men and women.

The formula used to make projections, per country and sex, was:

t x n t x n t t x n c M M ) 1 ( * 0 0
, where M denotes the mortality rate in the age interval x to x+n. De year t is expressed as the difference with the baseline year t 0 (=2005); thus t=25 for the year 2030. The change factor c * is the observed annual decline (exp(b)-1) for all age groups combined, multiplied by a ratio expressing the annual decline in the age groups 40-59, 60-79, and 80+ relative to the annual decline for all age groups (unweighted average over all countries). In the reference scenario, the applied ratios for men were 0.92, 1.03, and 0.99 for the age groups 40-59, 60-79, and 80 and more, respectively; and the ratios for women were 0.85, 0.97, and 1.09, respectively.

Future stomach cancer mortality numbers were estimated by applying the future sex and age-specific rates to the projected age-specific population sizes for the different countries up to 2030. These future population numbers have been obtained from Eurostat (baseline variant), except for Norway for which the data have been obtained from Statistics Norway (medium national Growth variant) (http://statbank.ssb.no/statistikkbanken). Looking into the projected mortality rates and absolute number of deaths in 2030, it can be measured if stomach cancer is a less important cause of death in the future.

Results

Figure 1 shows trends in age-standardized mortality rates between 1980 and 2005. In all countries, stomach cancer mortality rates declined steadily. The rate of decline is within a narrow range, from 3.61 percent per year for Norwegian men to 4.90 percent for Danish and Finnish women. The overall average rate of decline was 4.2 percent per year. The rate of decline remained similar over time for both sexes in most of these countries. There was no indication that the rate of decline (expressed in percent decline per year) had been levelling off in recent years. A main exception is the strong decline in mortality among Finnish women in the 1980's, when Finnish women were catching up with other countries.

In Figure 2, the rate of declines in stomach cancer mortality over the period 1980 to 2005 is plotted against the age-standardized mortality rates in 1980. The aim of this scatter plot is to assess whether the decline in stomach cancer mortality was higher in populations with higher initial mortality rates. However, both among men and women, no such positive association was found.

Table 1 presents projections for 2030 of stomach cancer mortality rates. Values projected for 2030 are compared to values observed for 2005 by means of ratios with 2005 as reference. In all seven European countries, age-standardized mortality rates would have declined in 2030 compared to 2005 if past trends were to continue. This decline would be substantial in all populations. The percentages of declines in 2030 compared to 2005 (calculated by one minus rate ratios or 1-RR) for men are -on average-0.66, 0.59, and 0.75 for the reference, pessimistic, and optimistic scenarios, respectively. These figures for women are 0.67, 0.59, and 0.77 respectively.

Table 2 presents projections for 2030 of the absolute number of stomach cancer deaths. Values projected for 2030 are compared to values observed for 2005 by means of ratios with 2005 as reference. The decline in the absolute number of deaths would be substantial as well, although with more variability between countries. The percentages of declines (in 2030 compared to 2005) for men are -on average-0.43, 0.32, and 0.59 for the reference, pessimistic, and optimistic scenarios, respectively. These figures for women are 0.54, 0.42, and 0.67 respectively.

Stomach Cancer Men

Discussion

The strong decline in stomach cancer mortality was found to continue at an undiminished rate until 2005 in each of the seven European countries. If this decline were to continue into the future, stomach cancer mortality rates would decline with about 66 percent between 2005 and 2030. The absolute number of stomach cancer deaths would diminish by about 50 percent despite the ageing of national populations. Even though weaker declines were predicted in the pessimistic scenario, this scenario as well implied that stomach cancer would become far less important as a cause of death in Europe in the next decades.

A key factor for the past declines in stomach cancer mortality is a decrease in the exposure to Helicobacter pylori (H. pylori) infection [START_REF] Ferlay | Estimates of the cancer incidence and mortality in Europe in 2006[END_REF] . H. pylori prevalence had been found to decline in parallel with the decrease of incidence of stomach cancer [START_REF] Axon | Helicobacter pylori: what do we still need to know?[END_REF] . However, H. Pylori infection may not be the only factor. It has been estimated to account for only half of the global total burden of stomach cancer because of synergistic or antagonistic roles of dietary and other exogenous factors [START_REF] Key | Diet, nutrition and the prevention of cancer[END_REF] . Complementary explanations for the observed declines in stomach cancer mortality may include better food preservation and refrigeration, improvement in environmental conditions and lifestyle changes [START_REF] Verdecchia | Explaining gastric cancer survival differences among European countries[END_REF][START_REF] Corley | Influence of site classification on cancer incidence rates: an analysis of gastric cardia carcinomas[END_REF] . Similarly, declines in the near future may still be largely driven by recent improvements in the living conditions that current generations had experienced during their childhood and adulthood.

Improved medical treatment may also have contributed to the decline in mortality from stomach cancer [START_REF] Corley | Influence of site classification on cancer incidence rates: an analysis of gastric cardia carcinomas[END_REF][START_REF] Verdecchia | Explaining gastric cancer survival differences among European countries[END_REF] like it contributed to the observed differences in stomach cancer mortality between European countries [START_REF] Verdecchia | Explaining gastric cancer survival differences among European countries[END_REF] . Therefore, further advances in the treatment for stomach cancer and H. pylori infections may be important to maintain a strong decline in stomach cancer mortality rates in the future. This role is supported by the observation that, while the incidence of stomach cancer had declined as well, mortality rates seemed to have declined at faster rates, indicating a decline in case fatality rates [START_REF] Fitzsimmons | Trends in stomach and pancreatic cancer incidence and mortality in England and Wales, 1951-2000[END_REF] . However, we should note that different trends in incidence and mortality of cancer are found for most solid tumours, and may be a consequence of increased detection of early tumours with often good prognosis [START_REF] Ferlay | Estimates of cancer incidence and mortality in Europe in 2008[END_REF] .

Empirical support for the expectation that the decline in mortality from stomach cancer will continue in the future comes from the trends that were observed for the past 25 years. Firstly, there is a strong consistency in the recent trends in stomach cancer mortality among both sexes and among each of the seven European countries. Second, these declines had persisted up to recent years in each of these populations, including those with the lowest initial mortality rates. Furthermore, steady declines in stomach cancer mortality were observed in the middle-aged and the young population as well, suggesting that they are likely to continue in the near future 4, 5, 16 . The latter observation is consistent with observations of cohort-wise patterns of decline in stomach cancer mortality in European countries 29 14, 30 , which may reflect the effects of improvement in living conditions in childhood [START_REF] Amiri | Trends in stomach cancer mortality in relation to living conditions in childhood. A study among cohorts born between 1860 and 1939 in seven European countries[END_REF] .

It should be emphasized that, even though future declines may seem likely, our study was primarily aimed at exploring possible future trends by extrapolating past trends. This extrapolation provides a baseline scenario against which new studies may formulate more specific scenarios of future trends. For example, policy-based scenarios may focus on the potential effects of specific preventive policies or advancement in the treatment of stomach cancer. As stomach cancer may become ever less important in terms of mortality, scenario studies will need to also include measures of incidence, prognosis and prevalence of stomach cancer.

Figure 1 -Figure 2 -

 12 Figure 1-Age-standardized stomach cancer mortality rates (per 100,000) in seven European countries by sex: 1980-2005 *

Table 1 -

 1 Stomach cancer mortality rates (age standardized, per 100,000), by scenario # , country and sex, for 2005 (observed) and 2030 (projected) The figure for 2030 divided by the related figure for 2005. For Denmark, the values for 2005 are not observed but projected. # Reference scenario = Equal future rate of change for all countries; Pessimistic scenario = Slowest observed rate of change for all countries; Optimistic scenario = Highest rate of change for all countries.

	Year

*

Table 2 -

 2 Stomach cancer absolute number of deaths, by scenario # , country and sex, for 2005 (observed) and 2030 (projected)

	Year

Data for Denmark, Finland, France metropolitan, the Netherlands, Norway, Sweden and England and Wales on the number of deaths from stomach cancer (ICD code 151 for ICD-8 and ICD-9, C16 for ICD-10 21 ) and corresponding numbers of population at risk, by sex and five-year age groups for the years 1980 to 2005 have been obtained from Eurostat (Denmark from 1999-2001, France from 1998-2005, and other countries from 2000 to 2005) and national data sources (1980-1999) [START_REF] Janssen | Trends in old-age mortality in seven European countries, 1950-1999[END_REF] . The


S Shoma

Nelianne J Verkaik
email: n.j.verkaik@erasmusmc.nl

C P Vogel

Peter W M Hermans

S Selm

Tim J Mitchell

M Roosmalen

Shahadat Hossain

Mahbubur Rahman

Hubert Ph Endtz

Shereen Shoma #

Corné P De Vogel

Saskia Van Selm

Maarten Van Roosmalen

Willem J B Van Wamel

Alex Van Belkum

Development of a multiplexed bead-based immunoassay for the

Keywords: Streptococcus pneumoniae, vaccine, proteins, antibody detection, Luminex

come

Introduction

Streptococcus pneumoniae (S. pneumoniae, pneumococcus) is an important human pathogen that causes life-threatening diseases such as pneumonia and meningitis, as well as less serious but highly prevalent diseases such as otitis media and sinusitis. Between 2002 and 2003, pneumonia accounted for 19% of the 10.6 million deaths per annum among children younger than 5 years [START_REF] Bryce | WHO estimates of the causes of death in children[END_REF]. Ninety percent of childhood deaths occurred in developing countries [START_REF] Williams | Estimates of world-wide distribution of child deaths from acute respiratory infections[END_REF], including Bangladesh [START_REF] Arifeen | Causes of death in children under five years of age Bangladesh Demographic and Health Survey[END_REF][START_REF] Baqui | Causes of childhood deaths in Bangladesh: an update[END_REF]. The global incidence of pneumococcal meningitis in children is 17 cases per 100,000. The case-fatality rate (CFR) for pneumococcal meningitis is high. In 2000, the global pneumococcal meningitis CFR was 59%, ranging from 29% in the Western Pacific to 73% in Africa [START_REF] O'brien | Burden of disease caused by Streptococcus pneumoniae in children younger than 5 years: global estimates[END_REF].

Fortunately, the availability of the 23-valent capsular polysaccharide and 7-valent pneumococcal conjugate vaccines (PCV23 and PCV7) has resulted in a dramatic reduction in morbidity and mortality of pneumococcal diseases. However, the high costs and still limited vaccine-mediated protection, which is restricted to the included serotypes, have prevented its implementation in large-scale immunization programs in developing countries. For these reasons, there is considerable interest in designing alternative and more cost-effective strategies.

Presently, several pneumococcal proteins are being evaluated as potential vaccine candidates.

Optimally, these would provide protection against pneumococcal infection regardless of serotype [START_REF] Black | Efficacy, safety and immunogenicity of heptavalent pneumococcal conjugate vaccine in children. Northern California Kaiser Permanente Vaccine Study Center Group[END_REF]. To date, the most promising protein vaccine candidates include pneumolysin (PLY), pneumococcal surface adhesin A (PsaA), pneumococcal surface protein A (PspA) and Pneumococcal surface protein C (PspC) [START_REF] Simell | Pneumococcal carriage and acute otitis media induce serum antibodies to pneumococcal surface proteins PspC and PhtD in children[END_REF]. These proteins are produced by virtually all clinical isolates of the pneumococcus. PLY is a 53-kDa protein that causes cytolysis, induces complement activation and the production of cytokines and nitric oxide [START_REF] Berry | Comparative virulence of Streptococcus pneumoniae strains with insertion-duplication, point, and deletion mutations in the pneumolysin gene[END_REF][START_REF] Braun | Pneumolysin, a protein toxin of Streptococcus pneumoniae, induces nitric oxide production from macrophages[END_REF][START_REF] Jedrzejas | Pneumococcal virulence factors: structure and function[END_REF][START_REF] Mitchell | Complement activation and antibody binding by pneumolysin via a region of the toxin homologous to a human acute-phase protein[END_REF][START_REF] Paton | The contribution of pneumolysin to the pathogenicity of Streptococcus pneumoniae[END_REF]. In addition, PLY has been assigned several functions with respect to modification of the immune response. PLY has recently been shown to interact with Toll-like receptor 4 (TLR-4) [START_REF] Mitchell | Streptococcus pneumoniae: virulence factors and variation[END_REF]. PsaA is a surfaceexposed 37-kDa lipoprotein that plays a major role in pneumococcal attachment to the host cell and virulence [START_REF] Rajam | Pneumococcal surface adhesin A (PsaA): a review[END_REF]. PspA is a choline-binding surface protein, which inhibits complementmediated phagocytosis, binds to lactoferrin and as such prevents lactoferrin-mediated killing [START_REF] Hammerschmidt | Identification of pneumococcal surface protein A as a lactoferrin-binding protein of Streptococcus pneumoniae[END_REF].

Antibodies to pneumococcal proteins PspA, PsaA and PLY have been shown to develop early in life [START_REF] Adrian | Development of antibodies against pneumococcal proteins alpha-enolase, immunoglobulin A1 protease, streptococcal lipoprotein rotamase A, and putative proteinase maturation protein A in relation to pneumococcal carriage and Otitis Media[END_REF][START_REF] Rapola | Natural development of antibodies to pneumococcal surface protein A, pneumococcal surface adhesin A, and pneumolysin in relation to pneumococcal carriage and acute otitis media[END_REF]. PspC (also known as CbpA or SpsA) acts as an adhesin and binds the complement regulatory protein factor H, to provide resistance to complement [START_REF] Mitchell | Streptococcus pneumoniae: virulence factors and variation[END_REF][START_REF] Dave | PspC, a pneumococcal surface protein, binds human factor H[END_REF].

We developed a multiplex bead-based immunoassay using Luminex xMAP ® Technology to gather novel insights into the immunogenicity of PLY, PsaA, PspA and PspC and 13 other pneumococcal proteins. With this assay, antibodies to these 17 pneumococcal proteins can be quantified simultaneously. Thus far, for the pneumococcus, this technology was only used for measurement of antibodies directed to different pneumococcal capsular polysaccharides [START_REF] Biagini | Method for simultaneous measurement of antibodies to 23 pneumococcal capsular polysaccharides[END_REF][START_REF] Lal | Development and validation of a nonaplex assay for the simultaneous quantitation of antibodies to nine Streptococcus pneumoniae serotypes[END_REF][START_REF] Pickering | A multiplexed fluorescent microsphere immunoassay for antibodies to pneumococcal capsular polysaccharides[END_REF].

Materials and methods

Antigens. The pneumococcal proteins PLY, PsaA, PspA, PspC, neuraminidase A (NanA), Hyaluronidase (Hyl), putative proteinase maturation protein A (PmpA), streptococcal lipoprotein rotamase A (SlrA), -enolase (Eno), immunoglobulin A1 protease (IgA1-protease), PdBD and BVH-3, SP1003, SP1633, SP1651, SP0189, SP0376 were used. NanA plays an important role in biofilm formation and promotes pneumococcal brain endothelial cell invasion [START_REF] Trappetti | Sialic acid: a preventable signal for pneumococcal biofilm formation, colonization, and invasion of the host[END_REF][START_REF] Uchiyama | The surface-anchored NanA protein promotes pneumococcal brain endothelial cell invasion[END_REF]. Hyl is present on the majority of strains. The enzyme degrades essential components of the host's extracellular matrix and as such it may facilitate bacterial spread in host tissues [START_REF] Mitchell | Streptococcus pneumoniae: virulence factors and variation[END_REF][START_REF] Humphrey | Hyaluronidase production by pneumococci[END_REF]. PpmA and SlrA, two surface-exposed lipoproteins have been shown to play a role in virulence and colonization [START_REF] Hermans | The streptococcal lipoprotein rotamase A (SlrA) is a functional peptidyl-prolyl isomerase involved in pneumococcal colonization[END_REF][START_REF] Cron | Surface-associated lipoprotein PpmA of Streptococcus pneumoniae is involved in colonization in a strain-specific manner[END_REF]. Eno has been identified both in the cytoplasm and on the surface of pneumococci and is involved in virulence by the activation of plasmin into plasminogen, a host collagen-degrading enzyme [START_REF] Bergmann | Alpha-Enolase of Streptococcus pneumoniae is a plasmin(ogen)-binding protein displayed on the bacterial cell surface[END_REF]. IgA1-protease plays an important role in pneumococcal colonization. IgA1-protease modifies IgA1 antibodies so that it promotes rather than inhibits pneumococcal adherence to epithelial cells in a model of colonization [START_REF] Weiser | Antibodyenhanced pneumococcal adherence requires IgA1-protease[END_REF]. PdBD is a mutant version of PLY with less cytolytic activity and reduced ability to activate complement [START_REF] Rijneveld | Roles of interleukin-6 and macrophage inflammatory protein-2 in pneumolysin-induced lung inflammation in mice[END_REF]. BVH-3 is capable of conferring protection against lethal experimental S. pneumoniae infection. This protein family was also designated Pht (for pneumococcal histidine triad) [START_REF] Adamou | Identification and characterization of a novel family of pneumococcal proteins that are protective against sepsis[END_REF][START_REF] Hamel | Prevention of pneumococcal disease in mice immunized with conserved surface-accessible proteins[END_REF]. SP1003, SP1633, SP1651, SP0189, SP0376 are pneumococcal proteins with currently unknown functions.

Recently, a comprehensive review summarizing the surface-exposed virulence factors and their functions was published [START_REF] Mitchell | Streptococcus pneumoniae: virulence factors and variation[END_REF].

The genes encoding the proteins NanA, PsaA, PspA, PspC, SP1633, SP1651, SP0189, SP0376, Hyl, PLY and PpmA were isolated from S. pneumoniae strain TIGR4 chromosomal DNA and cloned in the vector pOPINF using In-Fusion Technology. Cell extracts were made from the recombinant Escherichia coli (Rosetta) strains and the recombinant proteins purified by immobilized metal affinity chromatography using the poly-His tag added to the N-terminal end of the protein during the cloning procedure. The genes used for the production of the recombinant antigens Eno, SlrA and PpmA were amplified by PCR from S. pneumoniae D39. IgA1-protease was amplified from S. pneumoniae TIGR4. The amplified DNA was cloned into a pET11a expression vector (Stratagene) and electrotransformed into E. coli BL21(DE3). The expression of recombinant protein was induced by the addition of IPTG and the recombinant proteins were purified by Ni + affinity chromatography, as described before [START_REF] Adrian | Development of antibodies against pneumococcal proteins alpha-enolase, immunoglobulin A1 protease, streptococcal lipoprotein rotamase A, and putative proteinase maturation protein A in relation to pneumococcal carriage and Otitis Media[END_REF][START_REF] Hermans | The streptococcal lipoprotein rotamase A (SlrA) is a functional peptidyl-prolyl isomerase involved in pneumococcal colonization[END_REF][START_REF] Mccool | Serum immunoglobulin G response to candidate vaccine antigens during experimental human pneumococcal colonization[END_REF]. The genes encoding BVH-3 and PdBD were cloned into plasmid pPA195 and pPA180, respectively, and transformed into Lactococcus lactis PA1001. Overexpression of BVH-3 and PdBD was induced by nisin, essentially as described previously [START_REF] Van Roosmalen | Mucosal vaccine delivery of antigens tightly bound to an adjuvant particle made from food-grade bacteria[END_REF]. The purity of the recombinant proteins was determined by SDS-PAGE gel electrophoresis. The amino acid sequence of the proteins was confirmed with mass spectrometry (Ultraflex MALDI-ToF, Bruker Daltonics).

Coupling methods. To quantify antibodies directed against the 17 S. pneumoniae proteins simultaneously, the bead-based flow cytometry technique (xMAP ® ; Luminex Corporation) was applied. The purified proteins were coupled to fluorescent SeroMAP beads. The coupling procedure was performed as described elsewhere [START_REF] Verkaik | Comparison of carboxylated and Penta-His microspheres for semi-quantitative measurement of antibody responses to His-tagged proteins[END_REF][START_REF] Verkaik | Antistaphylococcal humoral immune response in persistent nasal carriers and non-carriers of Staphylococcus aureus[END_REF]. In each experiment, control beads were included to determine nonspecific binding. For control beads the coupling procedure was followed, only no S. pneumoniae protein was added. In case of nonspecific binding, the median fluorescence intensity (MFI) values were subtracted from the antigen-specific results. As a negative control, PBS-BN was included. Human pooled serum was used as a standard.

Multiplex S. pneumoniae antibody assay. The multiplex assay (serum incubated with the different fluorescence-colored antigen coupled beads mixed in 1 well) was validated by

comparing the MFI values for HPS obtained with this multiplex assay with the results for HPS obtained with singleplex assays (serum incubated with individual single-color antigen-coupled beads in separate wells). After validation, the different antigen-coupled microspheres were mixed to a working concentration of 4000 beads per color per well. The procedure used was the same as described elsewhere [START_REF] Verkaik | Comparison of carboxylated and Penta-His microspheres for semi-quantitative measurement of antibody responses to His-tagged proteins[END_REF][START_REF] Verkaik | Antistaphylococcal humoral immune response in persistent nasal carriers and non-carriers of Staphylococcus aureus[END_REF][START_REF] Verkaik | Induction of antibodies by Staphylococcus aureus nasal colonization in young children[END_REF]. To optimize dilutions, serum samples of children were diluted 1:25, Statistical analysis. The Mann-Whitney U test was used to compare differences in antipneumococcal antibody levels. Differences were considered to be statistically significant when 2sided P-values were <0.05.

Results and Discussion

Validation of the multiplex assay and reproducibility. The MFI values obtained for HPS with the multiplex assay were between 82% and 111% (median 94%) of those obtained with the singleplex assays. Therefore, it was considered legitimate to use the multiplex assay. Serum incubated with control beads resulted in median MFI values for IgG, IgA and IgM of 77 (range, 21-7210), 92 (range, 15-11549) and 1282 (range, 171-7966), respectively. This indicates that there was a quite variable level of non-specific binding that was consistently subtracted from the antigen-specific results. The negative control (PBS-BN) incubated with protein-coupled beads resulted in low MFI values (<10). Inter-assay variation was calculated from MFI values obtained from serum samples (n=54) run in 2 separate assays and averaged per protein and antibody isotype. For IgG, the median CV was 5%, range 2% (PsaA) to 9% (SlrA). For IgA, the median CV was 7%, range 3% (PspC and PsaA) to 14% (SlrA). For IgM, the median CV was 13%, range 11% (SP1633) to 15% (PspC, PpmA and PspA). Earlier studies found equal CVs for inter-assay variation using other types of antigens [START_REF] Lal | Development and validation of a nonaplex assay for the simultaneous quantitation of antibodies to nine Streptococcus pneumoniae serotypes[END_REF][START_REF] Verkaik | Comparison of carboxylated and Penta-His microspheres for semi-quantitative measurement of antibody responses to His-tagged proteins[END_REF][START_REF] Ray | Development, validation, and implementation of a multiplex immunoassay for the simultaneous determination of five cytokines in human serum[END_REF]. Therefore, this new multiplexed antipneumocococcal immunoassay is considered reproducible from experiment to experiment.

Anti-pneumococcal antibodies. Anti-pneumococcal antibody levels showed extensive interindividual variability (Figure 1), probably due to the variable number of previous encounters with different S. pneumoniae strains (colonization, subclinical infection) as well as inter-individual differences in the ability to mount a humoral immune response. The levels of IgM were not significantly different between children suffering from a pneumonia/meningitis caused by S. pneumoniae or by a different bacterial species. The levels of IgG directed against NanA and the levels of IgA directed against Hyl and PpmA were higher in children suffering from a pneumonia caused by S. pneumoniae than in children suffering from a pneumonia caused by other bacterial pathogens (P<0.05, Figure 1a). In contrast, significantly higher levels of antibodies directed against IgA1-protease were detected in serum samples of patients suffering from a non-S. pneumoniae pneumonia. This might be explained by the fact that IgA1-protease is also present in Haemophilus influenzae [START_REF] Kilian | Antigenic heterogeneity of immunoglobulin A1 proteases from encapsulated and non-encapsulated Haemophilus influenzae[END_REF], another frequent cause of pneumonia in non-vaccinated children. In meningitis patients, the levels of IgG directed against PspC, PdBD and Ply were significantly higher in children suffering from meningitis caused by a bacterial pathogen other than S. pneumoniae (P<0.05, Figure 1b). Furthermore, the levels of IgA directed to PdBD, PLY, SP0189 and SP1003 were higher in this group (P<0.05). We speculate that the lower S. pneumoniae antibody levels in the children suffering from S. pneumoniae meningitis might be due to consumption of antibodies. However, such data can be obscured by the fact that we do not know the exact patient delay (time between the onset of disease and the arrival in the hospital and the collection of serum). Furthermore, sera from children with different ages were included.

In conclusion, using our novel non-polysaccharide-based high-throughput immunological test we were able to determine the level of antibodies to 17 pneumococcal proteins simultaneously in small volume serum samples. The requirement for small sample volumes only is a great advantage as compared to ELISA, especially when limited amounts of samples are available, which is usually the case in young children. Novel antigens can be easily included and combinations of protein and polysaccharide antigens may help profiling the complete immune response of both patients and healthy individuals to a common opportunistic pathogen, especially relevant in children in developing countries.

1 :

 1 50 and 1:100 in PBS-BN. The secondary antibodies were diluted 1:50, 1:100 and 1:200 in PBS-BN. Checkerboard titrations were performed. Considering results in MFI values and amounts of serum and secondary antibody needed, optimal serum dilutions were 1:100 for measurement of IgG and 1:50 for measurement of antigen-specific IgA and IgM. The optimal secondary antibody dilutions were 1:200 for IgG and 1:100 for IgA and IgM. Measurements were performed on the Luminex 100 instrument (BMD) using Luminex IS software (version 2.2). Tests were performed in duplicate, and the MFI values, reflecting semi-quantitative antibody levels, were averaged. Anti-pneumococcal antibodies. The multiplexed immunoassay was used to compare differences in anti-pneumococcal antibodies in serum samples obtained from 54 children under the age of 5 years with clinical suspicion of pneumonia (n=26) or meningitis (n=28). Half of these children suffered from pneumonia/meningitis caused by S. pneumoniae, the other half suffered from pneumonia/meningitis caused by a bacterial species other than S. pneumoniae. The serum samples were stored at -80°C until use. The research protocol was approved by the Research Review and Research Ethical Committee of ICDDR,B and written informed consent from the parents was obtained.

Figure caption


Tobias ; Rydén

Søren Asmussen

A note on skewness in regenerative simulation

Keywords: batch means, busy cycle, coverage probability, regenerative process, regenerative simulation, skewness AMS 1991 classification. Primary 00A72, 65C05. Secondary 68M20, 90B05, 90B22 ABC interval, batch means, busy cycle, busy period, coverage probability, Fieller test, jackknife, M/G/1 queue, regenerative process, regenerative simulation, skewness

The purpose of this paper is to show, empirically and theoretically, that performance evaluation by means of regenerative simulation often involves random variables with distributions that are heavytailed and heavily skewed. This, in turn, leads to the variance of estimators being poorly estimated, and confidence intervals having actual coverage quite different from (typically lower than) the nominal. We illustrate these general ideas by estimating the mean occupancy and tail probabilities in M/G/1 queues, comparing confidence intervals computed from batch means to various intervals computed from regenerative cycles. In addition we provide theoretical results on skewness to support the empirical findings.

Introduction

Let X = (X(t)) t≥0 be a stochastic process which is regenerative with independent cycles. That is, there exists a renewal process (S k ) ∞ k=0 such that the regenerative cycles (X(t)) S n-1 ≤t<Sn , n = 1, 2, 3, . . . , are i.i.d., and not dependent on the process up to S 0 (cf. Asmussen, 2003, p. 169). If S 0 = 0 the renewal process is pure, while it is delayed if S 0 > 0. The renewal epochs S n are referred to as the regeneration times of X. Typical examples of regenerative processes are found in queueing and storage theory, with the canonical example being the GI/GI/1 queue which regenerates every time a customer arrives to an empty system. Now write E 0 for the conditional expectation given S 0 = 0 (i.e., given that the process is initialised by a regeneration), τ = S 1 for the first subsequent regeneration time and µ = E 0 τ for the mean inter-renewal time. Then, provided that µ is finite, the distribution of τ is non-lattice, and that X has right-continuous paths, it holds that X(t) has an asymptotic distribution, P e say, defined by

E e [f (X(•))] = 1 µ E 0 τ 0 f (X(s)) ds (1.1)
for any real measurable bounded function f of X (Asmussen, 2003, Theorem VI.1.2). Regenerative simulation amounts to simulating i.i.d. copies of (τ, C), with C being the integral in (1.1), and using these for estimating the ratio of means on the right-hand side of (1.1). A difficulty with regenerative simulation however, that we want to highlight in this paper, is that such methods may suffer from distributions being heavy-tailed and heavily skewed. Our canonical example will be that of estimating the mean number of customers in an M/GI/1 queue; then X(t) is the number of customers in the system at t, f is the identity function, and C is the area under X during a busy cycle. It has been established that for the M/M/1 queue with unit service rate, this variable has tail probabilities of the form P(C > x) ∼ (1 -ρ)/(ρ √ 2πγ) × x -1/4 e -γ √ x , where ρ is the system load and γ = 2 -2(1 -ρ) -(1 + ρ) log ρ; see Guillemin and Pinchon (1998, Eq. (5.20)) and Kearney (2004, Eq. (24)). This tail is subexponential, implying that although a central limit theorem (CLT) for the sample mean of i.i.d. copies of C holds, convergence will be slow. For extensions of these results to GI/GI/1 queues and regularly varying service times, see [START_REF] Kulik | Tail behaviour of the area under the queue length process of the single server queue with regularly varying service times[END_REF]. We will also show that as ρ ↑ 1, the distributions of τ and C become totally unbalanced in the sense that asymptotically the probability that they exceed their means tend to zero. From a simulation perspective this means that simulating the expectation of C for instance is difficult for large ρ, since one must take a sufficient number of replications to pick up a few very large values of C before the sample mean becomes a good estimate of its expectation.

The outline of the paper is as follows. First in Section 2 we review the estimators that will be used in the paper (for recent surveys about of the general state-of-the-art in the areas, see e.g. [START_REF] Alexopoulos | Statistical analysis of simulation output: State of the art[END_REF][START_REF] Law | Statistical analysis of simulation output data: The practical state of the art[END_REF]. Then in Section 3 we give some numerical examples, in which the method of batch means performs better than regenerative simulation. We also illustrate, by means of simulation, the skewness discussed above, and Section 4 provides proofs of this asymptotic unbalancedness. for some large T (TA stands for time average). If the simulated X is not stationary the distribution of X(t) will however be different from P e for any t > 0, although approaching P e as t → ∞. Such a transient will induce a bias in θTA . The standard way of addressing this problem is to first simulate X(t) during a burn-in period of length T 0 , and then use the sample path over [T 0 , T 0 + T ] to compute the empirical average. However, deciding on what is a sufficiently large value for T 0 is not always easy.

Estimators of steady state expectations

A straightforward approach to estimating θ = E e [f (X(•))], is through the empirical average θTA = T -1 T 0 f (X(t)) dt (2.
Another difficulty with time averages concerns the computation of confidence intervals. Often a CLT with rate T 1/2 holds for θTA , i.e. T 1/2 ( θTA -θ) → d N(0, σ 2 ). The asymptotic variance σ 2 , often referred to as the time average variance constant (TAVC), is then given by

σ 2 = ∞ -∞ Cov e [X(0), X(s)] ds, (2.2)
with Cov e denoting covariance under the stationary regime. Estimating the TAVC is a non-trivial task, so that computing a confidence interval for θ is not straightforward. See Asmussen and Glynn (2007, Section IV.3) for further reading. An alternative method is that of batch means. This amounts to splitting [0, T ] into B equally long subintervals, and letting θ(k) = (T /B) -1 kT /B (k-1)T /B f (X(t)) dt be the average over the k-th subinterval. If the process X is sufficiently weakly dependent, then the θ(k) will be roughly Normally distributed and independent. Thus a (two-sided) confidence interval for θ with approximate coverage 1 -α, is given by θTA ± t 1-α/2 (B -1)s B / √ B with s 2 B the sample variance of the θ(k) and t 1-α/2 (B -1) being the (1 -α/2)-quantile of the t-distribution with B -1 degrees of freedom.

We now turn to regenerative simulation, and let ((

τ k , C k )) n k=1 be i.i.d. copies of (τ, C). The standard estimator of θ = EC/Eτ is θREG = n k=1 C k n k=1 τ k = C n τ n , (2.3)
where

τ n = n -1 n 1 τ k and C n = n -1 n 1 C k .
That this estimator is consistent is evident from (1.1). It is not unbiased however, because θREG is a ratio estimator. Provided that τ and C both have finite second moment, the CLT n 1/2 (τ n -Eτ, C n -EC) ⊤ → d N(0, Σ) holds for some covariance matrix Σ. Using the delta method one finds that θREG is asymptotically Normal at rate n 1/2 , with asymptotic variance η 2 that can be estimated as

η2 = 1 n -1 n k=1 (C k -θREG τ k ) 2 τ 2 n ;
(2.4) see e.g. Asmussen and Glynn (2007, Proposition IV.4.1). Thus a (two-sided) confidence interval can be constructed as θREG ± z 1-α/2 η/ √ n, where z 1-α/2 is the (1 -α/2)-quantile of the standard Normal distribution.

Another standard approach to constructing confidence intervals is through the jackknife (e.g. Ripley, 1987, pp. 158-160). This estimator is

θJACK = n θREG -(n -1) θ(-) with θ(-) = n -1 n k=1 θ(k) and θ(k) = i =k C i i =k τ i ,
and the corresponding variance estimator is s 2 JACK = (n -1)/n × n k=1 ( θ(-) -θ(k) ) 2 . Confidence intervals are usually computed as θJACK ± z 1-α/2 s JACK .

The appealing features of regenerative simulation are thus that the method avoids dealing with transient periods (burn-in) and dependence. The fact that θREG is a ratio estimator however hints that even if (τ n , C n ) is Normal distributed with reasonable accuracy, it is not obvious that the same conclusion may be drawn about θREG itself. One purpose of this paper is indeed to compare confidence intervals computed using batch means, to intervals computed using regenerative simulation. In addition to the standard delta method intervals and jackknife intervals, we will also study three other methods described below.

One way of circumventing the ratio estimator is to construct a linear test procedure based on a t-test. Forming Z k (θ) = C k -θτ k it is clear from (1.1) that EZ k (θ) = 0; testing θ = θ 0 can thus be done by testing if EZ k (θ 0 ) = 0. By inverting a t-test, a (two-sided) confidence interval for θ with approximate coverage 1 -α is obtained as

{θ : |Z n (θ)/s Z(θ) | ≤ t 1-α/2 (n -1)}, (2.5)
where Z n (θ) and s 2 Z(θ) are the sample mean and variance, respectively, of {Z k (θ)} n k=1 . This interval was originally derived by Fieller; see e.g. [START_REF] Fieller | Some problems in interval estimation[END_REF] or Roy and Potthoff (1958, Section 4), who both considered it as an exact test for a bivariate Normal distribution. In the appendix we show that this procedure is also equivalent to several tests in multivariate analysis of variance (MANOVA). A point estimate can be obtained as the mid-point of the interval, or as the θ that maximises the p-value of the test. The latter is the θ such that the test statistic has minimum absolute value, which occurs when Z n (θ) = 0 and is attained for θREG . The interval is in general not symmetric around θREG however.

Further refinements of the above methods will be done by means of ABC (approximate bootstrap confidence) intervals, introduced by DiCiccio and [START_REF] Diciccio | More accurate confidence intervals in exponential families[END_REF]. Their basic purpose is to provide confidence intervals whose end-points are correct up to order O(n -1 ), rather than O(n -1/2 ) as for the standard intervals, and this is done by including third-order properties (i.e., skewness) of the observed variables. In the present setting the ratio estimator θREG = C n /τ n and the Fieller statistic Z(θ) = C n -θτ n are non-parametric estimators, and they are both functions of the sample means τ n and C n . The influence function for the i-th cycle

(τ i , C i ) is given by ℓ i = C i /τ n -C n /τ 2
n × τ i for the ratio estimator and by ℓ i = Z i (θ) -Z n (θ) for the Fieller statistic, where Z i (θ) = C i -θτ i as above. The remaining steps of the ABC procedure are given e.g. in Davison and Hinkley (1997, Eq. (5.46)). For the Fieller statistic the expressions become particularly simple; the upper end-point of the one-sided ABC interval for the mean of Z(θ) with approximate coverage 1 - and a two-sided confidence interval is obtained by computing end-points for approximate coverage 1 -α/2 and α/2, respectively.

α is Z n + (a + z 1-α )/(1 -a(a + z 1-α )) 2 × n -1/2 s Z(θ) with a = skew(Z(θ))/(6 √

Numerical examples

This section contains numerical comparisons of confidence intervals computed using the methods discussed above. For the numerical results we used the M/GI/1 queue with unit mean service time and X being the number of customers in the system. The function f was chosen either as the identity, so that θ becomes the mean number of customers in stationarity, or as the indicator I(X(t) ≥ b), so that θ becomes the probability P e (X(•) ≥ b). In the following we chose b = 3. We studied three types of queues: the M/M/1 queue, the M/D/1 queue, and the queue with service time density

b(x) = √ 3 π 1 1 + x 6 /27 for x > 0. (3.1)
The corresponding distribution has unit mean and variance 1/2, so that its variance is intermediate between that of the exponential distribution, with unit variance, and that of the constant unity with no variance. We will denote this queue as M/HT/1 (HT for heavy tail).

The three cases are meant to represent light tail, no tail, and heavy tail, respectively, of the service time distribution. For each case we studied the two loads ρ = 0.5 and ρ = 0.8. For the case of the number of customers in the system (i.e., f is the identity function), the Pollaczek-Khintchine formula (e.g. King, 1990, Eq (5.6)) yields the mean number of customers ρ/(1 -ρ), ρ(1 -ρ/2)/(1 -ρ), and ρ(1 -ρ/4)/(1 -ρ) respectively, for the three queues. We also remark that the TAVC σ 2 in (2.2) is finite even for the M/HT/1 queue. This follows as σ 2 is finite provided that C has finite second moment (e.g. Asmussen and Glynn, 2007, Proposition 4.2), and Daley and Jacobs (1969, Eq. (8.7)) show that so is the case if the service time distribution has finite fourth moment. Obviously this holds true for the density (3.1).

Regarding the tail probability P e (X(•) ≥ b), the stationary distribution of the number of customers in the M/M/1 queue is geometric with parameter ρ, so that P e (X(•) ≥ b) = ρ b for this queue. For the M/D/1 and M/HT/1 queues, let π k be the stationary probability of k customers in the system. We have π 0 = 1 -ρ, and we computed π 1 and π 2 numerically by solving the first two global balance equations of the Markov chain embedded after departures (e.g. Asmussen, 2003, Eq. (5.5)); for the M/HT/1 queue the integrals corresponding to the probabilities of 0 or 1 arrivals during a service time were evaluated using Maple. In this way we obtained P e (X(•) ≥ 3) = 0.0530 and P e (X(•) ≥ 3) = 0.3655 respectively for ρ = 0.5 and ρ = 0.8 for the M/D/1 queue, and P e (X(•) ≥ 3) = 0.0914 and P e (X(•) ≥ 3) = 0.4564 respectively for ρ = 0.5 and ρ = 0.8 for the M/HT/1 queue.

Tables 1234show empirical coverage probabilities obtained from regenerative simulation and time-average simulation. For the regenerative simulation, regeneration epochs were defined as arrivals to an empty queue, and confidence intervals were computed using (i) the estimator θREG with variance estimate (2.4) and Normal quantiles, (ii) the jackknife estimator with jackknife variance estimator and Normal quantiles, (iii) the Fieller intervals (2.5), (iv) ABC intervals derived from the ratio estimator θREG (denoted ratio-ABC in the tables), and (v) ABC intervals derived from the Fieller statistic Z n (θ) (denoted F-ABC in the tables). Neither the Fieller intervals nor the Fieller-ABC intervals were inverted into confidence intervals for θ however, since we were only interested in whether these tests rejected the true θ or not.

For the time-average simulation, confidence intervals were computed using the estimator θTA , variance estimates obtained by splitting the total simulation time into B = 5, 10 or 25 batches, and t-quantiles (Asmussen and Glynn, 2007, Section IV.5, recommend 5-30 batches). [START_REF] Alexopoulos | Statistical analysis of simulation output: State of the art[END_REF] describes elaborate algorithms for adaptively choosing the number of batches; our aim here is however not to use the most sophisticated algorithms, but rather to show that even with a simple choice for the number of batches, this approach performs better than intervals based on regenerative simulation.

Sample sizes were n = 500, 1,000, and 2,000 regenerative cycles for the estimates based on regenerative simulation. Since a regenerative cycle is a busy period plus a waiting time until another arrival, its mean value is Eτ = 1/(1 -ρ) + 1/ρ = 1/ρ(1 -ρ) in a queue with unit mean service time (see e.g. King, 1990, Eq. (8.1)). For the time-average estimator the total simulation time was chosen as T = nEτ , preceded by a burn-in period of duration T 0 = T /5 for the maximal value of T (corresponding to n = 2,000). All simulations were done in Matlab.

From Tables 1-2, reporting empirical coverages for the mean number of customers, we can draw the general conclusions that performance is best, i.e. estimation of θ is easiest, for the M/D/1 queue and, maybe surprisingly, worst for the M/M/1 queue with the M/HT/1 queue in between. A possible explanation for this is that the M/M/1 queue has the largest variance of the service time distribution. Of the five methods based on regenerative simulation, we see that the jackknife and ABC intervals perform best; this finding is in agreement with those of Iglehart (1975, Section 5), who found the jackknife intervals to be superior to several other ones, including the Fieller and the standard delta method intervals. A closer analysis reveals that the differences towards the other intervals are statistically significant at the 5%level (one-sided test) in most cases for ρ = 0.8, but in only a few cases for ρ = 0.5. The Fieller-ABC intervals perform the best in all cases for ρ = 0.8. For the M/HT/1 queue this is peculiar, as the cumulative variable C does not have finite third moment for the service time density (3.1). The ABC intervals involve the empirical skewness, whence these procedures are highly questionable from a theoretical perspective for the M/HT/1 queue, despite their better performance in the simulations.

The overall best intervals are however those computed from batch means with 5 batches. As indicated by the tables, this method has consistently better coverage than the other methods, and these differences are significant (5%-level, one-sided) in most cases. (We did however not attempt to correct for multiple testing.) However, also for this method the empirical coverage is significantly below the nominal 95% for all queues when ρ = 0. Table 1: Empirical coverage probabilities for different simulation strategies and confidence intervals applied to the mean number of customers in an M/GI/1 queue with unit mean service time and load ρ = 0.5. Top part: M/M/1 queue; middle part: M/D/1 queue; bottom part: M/GI/1 queue with service time density (3.1). For the time-average estimator, the simulation periods were T = n/ρ(1-ρ), preceded by a burn-in period of length T 0 = T /5 for the maximal T . All confidence intervals and tests were two-sided with nominal coverage probability 95% (size 5%). The number of replications was 8,000, giving two-sided confidence intervals for the actual coverage probabilities with half-widths of about 0.007. Figures in italics mark coverages that are significantly smaller, at the 5% level (one-sided test for difference of probabilities in two binomial distributions), than the coverage (on the same row) for the time-average estimator and 5 batches. for the M/M/1 and M/HT/1 queues also when ρ = 0.5, even for the largest sample size, although the the differences are small for the M/D/1 queue and for the M/M/1 queue when ρ = 0.5. The results for estimating the tail probability P e (X(•) ≥ 3), reported in Tables 34, are in many respects similar to those regarding the mean number of customers, although the differences between the two loads are less pronounced. The two types of ABC intervals are indistinguishable for ρ = 0.5, and for ρ = 0.8 there are no notable differences at all between these and the jackknife intervals. Again batch means with 5 batches performs best, and often significantly better than the other methods. The coverage probabilities of all methods are better than for the mean number of customers, however.

Obviously the Fieller intervals (2.5) do not perform better than the other tests based on regenerative simulation, which are ratio estimators. This suggests that these methods perform worse than batch means not because they are ratio estimators, but because the average (τ n , C n ) is not approximately Normal with sufficient accuracy. A Normal probability plot (not shown) with simulated replications of the C n for the number of customers in an M/M/1 queue with ρ = 0.8 and n = 2,000, does indeed reveal some deviations from Normality. From a sample of 250,000 simulated replications of C, we obtained the sample mean 25.00 (the true value is 25), and the sample standard deviation 155.2. Moreover, the sample skewness was 20.7, and the fraction of replications larger than 25 was 0.112; these are both numbers that indicate a large degree of skewness of the distribution of C. An estimate of the density of C is plotted in Figure 1.

We also simulated, for the same setting, a sample of 250,000 replications of (τ, C) and the corresponding Z(θ) = C -θτ with the true θ = ρ/(1 -ρ) = 4. The mean of this sample was 0.074 (the true value is 0) and its sample standard deviation was 104.2. Moreover the sample skewness was 29.3, and the fraction of replications larger than 0 was 0.046. An estimate of its density is found in Figure 1, and one should note the very sharp decrease of the density just left of zero. Our conclusion here is that splitting a sample trajectory as is done in regenerative simulation, implies a way of estimating the TAVC that involves distributions which are both heavy-tailed and heavily skewed, and that this explains the inferior performance of estimators built on regenerative simulation. Even the ABC intervals, which do adjust for skewness, are unable to fully compensate for it and achieve coverage better than for the batch means. We did also compute the BC a (bias-corrected accelerated) intervals for the Fieller statistic Z n (θ), which is another method of constructing second-order correct confidence intervals for the mean of this estimate (see Efron, 1987, Example 2). These intervals did however in fact perform worse than the Fieller intervals, and we do not report these results.

Regarding the variable τ , its distribution is not heavy-tailed in the M/M/1-setting, because the busy cycle distribution has an exponential tail; this follows as its Laplace transform (e.g. Asmussen, 2003, Proposition III.8.10) is finite to the left of the origin. Its distribution is however asymptotically unbalanced in a sense detailed in the next section; for the sample used here, a fraction 0.200 of the replications were larger than the true mean 1/ρ(1 -ρ) = 6.25. These observations may explain why the confidence intervals built on regenerative simulation 9 perform worse than those built on batch means also for the tail probability P e (X(•) ≥ 3), for which the function f (x) = I(x ≥ 3) is bounded, but with smaller differences than for estimating the mean number of customers.

For the time-average estimator, the distribution of concern is the stationary distribution P e (f (X(•) ∈ •)). For the M/M/1 and M/D/1 queues, this is not heavy-tailed for the examples studied here. The concern is then rather about the dependence over time. This dependence is however light-tailed (exponentially decaying) both for the M/M/1 queue and for the M/D/1 queue. For these queues there are thus no heavy-tailed phenomena disturbing the timeaverage estimator, while there are for the estimators built on regenerative simulation. For the M/HT/1 queue we however expect heavy-tailed dependence over time and also a heavy-tailed distribution for the number of customers (although a Bernoulli distribution when using the indicator function for the tail probability). This however does not have a notable degrading effect on coverage in the simulations reported here. As noted above, for 5 batches the batch means intervals perform much better than those built on regenerative simulation. We did not explore further why coverage is better with 5 batches than with 10 or 25, but two possible explanations are that (i) with the smaller batch sizes the batch averages are correlated to an extent that the variance estimate (based on the assumption of no correlation) becomes too small, or that (ii) with the smaller batch sizes the distribution of the batch means becomes too far from Normal, so that the distribution of the t-statistic is not well approximated by the t-distribution. Of course, also both of (i) and (ii) may occur. In this section we address from a theoretical perspective the empirical observations about skewness made in the previous section. First we make the following definition.

Definition 1 Let (ν s ) be a family of distributions on R, for s in some index set of integers or reals, and let m s = R x ν s (dx) be the mean of ν s . We say that this family is asymptotically totally unbalanced as s ↑ a, where a is a limit point, finite or infinite, if the mean m s is finite for all s, and ν s ({x : x > m s }) tends to either 0 or 1 as s ↑ a. For s ↓ a the definition is analogous.

Consider first the variable C when X is the number of customers in an M/M/1 queue and f is the identity function. In the Introduction we mentioned that C then has a subexponential distribution. Write P ρ for the distribution of (X(t)) t≥0 , when stationary with traffic load ρ, and E ρ for the corresponding expectation.

Proposition 1 For X = (X(t)) t≥0 being the number of customers in an M/M/1 queue with load ρ, P ρ (C ∈ •) and P ρ (τ ∈ •) are asymptotically totally unbalanced as ρ ↑ 1, namely

P ρ (C > E ρ C) → 0 as ρ ↑ 1 and P ρ (τ > E ρ τ ) → 0 as ρ ↑ 1.
Proof. We do the proof for C only, as the proof for τ is completely analogous. We first notice that C is finite P 1 -a.s., because even for ρ = 1 a busy cycle is of finite duration a.s.; this is since the embedded X-process, i.e. X observed at arrivals and departures only, is then a nullrecurrent random walk. Second, for ρ < ρ ′ ≤ 1 the P ρ ′ -distribution of C is stochastically larger than the P ρ -distribution of C. This can be established, for instance, by a coupling argument: run two queues with identical sequences of service times, with unit mean, and inter-arrival times constructed as U k /ρ and U k /ρ ′ respectively, where (U k ) is an i.i.d. sequence of standard negative exponential random variables. Letting C and C ′ correspond to the first busy cycle of the respective queues, it then holds that C ′ ≥ C a.s.

From this stochastic domination result we conclude that

P ρ (C > E ρ C) ≤ P 1 (C > E ρ C) for all 0 ≤ ρ ≤ 1.
Since C is finite P 1 -a.s., the proof is now completed by letting ρ ↑ 1 and noting that then

E ρ C → ∞. 2
In the previous section we saw that even for a moderate value like ρ = 0.8, skewness can be severe.

The above argument also applies to M/GI/1 queues, and also to more general systems like the GI/GI/1 queue. For the M/GI/1 queue it holds that (again assuming unit mean service time) E ρ C = 1/(1-ρ)+(ρb 2 /2)/(1-ρ) 2 with b 2 being the second moment of the waiting time distribution (Daley and Jacobs, 1969, Eq. (8.6)). Moreover, P 1 (C > x) ∼ Kx for some constant K (Kearney, 2004, Eq. (17)). These results are however unimportant for the above argument.

We now turn to the Fieller statistic Z(θ) = C -θτ in (2.5), where we again let X be the number of customers in an M/M/1 queue. Then θ = E e X(•) = ρ/(1 -ρ), and EZ(θ) = 0.

Proposition 2 For X = (X(t)) t≥0 being the number of customers in an M/M/1 queue with load ρ, P ρ (Z(θ(ρ)) ∈ •) is asymptotically totally unbalanced as ρ ↑ 1, namely

P ρ (Z(θ(ρ)) > 0) → 0 as ρ ↑ 1.
Proof. Consider a busy cycle during which in total d ≥ 1 customers are served. Clearly C is then maximal when all customers arrive at the very start of the busy cycle, and the maximal value of C is S 1 + 2S 2 + . . . + dS d , where S k is the service time of the k-th customer served during the busy cycle. The length of the regenerative cycle is S 1 + . . . + S d + A, where A is the time from the end of the busy cycle to the next arrival (then to an empty queue). Thus, for such a busy cycle we have Z

(θ) ≤ (1 -θ)S 1 + (2 -θ)S 2 + . . . + (d -θ)S d . If ρ is large enough that θ = θ(ρ) ≥ d, say ρ ≥ ρ(d)
, the right-hand side of this inequality will be non-positive. Since θ(ρ) → ∞ as ρ ↑ 1, such a ρ(d) can always be found. Now write M for the number of customers served during the busy period. Since the above argument does not involve the actual traffic load, we conclude that for any ρ such that θ(ρ) ≥ d we have

P ρ (Z(θ(ρ)) > 0) = P ρ (Z(θ(ρ)) > 0 | M ≤ d)P ρ (M ≤ d) + P ρ (Z(θ(ρ)) > 0 | M > d)P ρ (M > d) ≤ 0 + P ρ (M > d) ≤ P 1 (M > d).
When ρ ↑ 1 we can let d → ∞, and hence, because M is finite P 1 -a.s. (cf. the proof of Proposition 1), the result follows.

2

This result explains why the Fieller intervals work poorly for ρ close to 1. We however remark that, again, the numerical results indicate that also for moderate values of ρ these intervals can suffer substantially from skewness, and that the Fieller-ABC intervals are not able to fully compensate for this skewness. We also remark that again the above argument is valid for more general systems, like M/GI/1 and GI/GI/1 queues.

A The Fieller test as a MANOVA test

Consider an 2 × n-matrix Y of random variables, where the columns of Y are i.i.d. random vectors with some unknown mean vector ξ and unknown covariance matrix. We shall think of each column as being an independent copy of (τ, C) ⊤ from a regenerative cycle. To test the hypothesis H 0 : θξ 1 -ξ 2 = 0, we compute the quantities where C = [θ, -1] and A is an 1 × n-vector of ones; see [START_REF] Khatri | A note of a MANOVA model applied to problems in growth curve[END_REF] or Srivastava and Carter (1983, Section 6.3.2). These expressions follow from Khatri (1966, Eqs. ( 17)-( 18)) after some simple algebra (in Khatri's notation, B = I 2×2 , V = 1, Y 1 = Y, A * = A). Under the assumption that each column of Y is bivariate Normal, the likelihood ratio tests rejects H 0 if Q/(Q + P ) is small. Another test, known as Pillai's trace, rejects H 0 if P/(P + Q) is large. Since P and Q are scalars in this problem and P/(P + Q) = 1 -Q/(P + Q), these tests are equivalent. Pillai's trace test is known to possess robustness against non-normality (Kotz and Johnson, 1985, p. 25), which is appealing as in regenerative simulation (τ, C) is in general not bivariate Normal. It is also immediate that the LR test and the Pillai trace test are both equivalent to rejecting H 0 for large values of P/Q. However, using the notation of (2.5), we have P = nZ 2 n (θ) and Q = (n -1)s 2 Z(θ) , and we find that P/Q is, up to a multiplicative constant, equal to the Fieller statistic (2.5) squared. Thus the both of the MANOVA test statistics discussed here, and the Fieller t-test, are all equivalent.

P = CYn -1 A ⊤ AY ⊤ C ⊤ , Q = CY(I n -n -1 A ⊤ A)Y ⊤ C ⊤ ,

  n) and skew(Z(θ)) the sample skewness of the Z i (θ

5

  Page 6 of 16 URL: http://mc.manuscriptcentral.com/lssp E-mail: comstat@univmail.cis.mcmaster.

Figure 1 :

 1 Figure 1: Estimates of the densities of the random variables C (left) and Z(θ) (right), for the number of customers in the M/M/1 queue with ρ = 0.8, and θ set to its true value ρ/(1-ρ) = 4. The estimates are renormalised histograms computed from 250,000 replications of C and Z(θ) respectively.

  Estimates of the densities of the random variable C, for the number of customers in the M/M/1 queue with ρ = 0.8, and θ set to its true value ρ/(1-ρ) = 4. The estimate is a renormalised histograms computed from 250,000 replications of C. 187x146mm (600 x 600 DPI)Estimates of the densities of the random variable Z(θ), for the number of customers in the M/M/1 queue with ρ = 0.8, and θ set to its true value ρ/(1-ρ) = 4. The estimate is a renormalised histograms computed from 250,000 replications of Z(θ). 189x146mm (600 x 600 DPI)

Table 2 :

 2 Empirical coverage probabilities for mean number of customers as in Table1, but with load

								Page 8 of 16
			regenerative simulation			batch means
		delta jackknife Fieller ratio-ABC F-ABC	5	10	25
	n = 500	0.912	0.914	0.908	0.918	0.918	0.934 0.925 0.911
	n = 1,000 0.920	0.923	0.918	0.926	0.927	0.945 0.941 0.934
	n = 2,000 0.934	0.935	0.933	0.932	0.932	0.944 0.940 0.940
	n = 500	0.928	0.928	0.928	0.928	0.926	0.946 0.942 0.938
	n = 1,000 0.941	0.941	0.941	0.940	0.939	0.944 0.944 0.940
	n = 2,000 0.945	0.946	0.944	0.942	0.942	0.947 0.946 0.944
	n = 500	0.920	0.922	0.918	0.925	0.923	0.942 0.932 0.923
	n = 1,000 0.932	0.935	0.932	0.932	0.931	0.941 0.937 0.933
	n = 2,000 0.945	0.947	0.944	0.943	0.943	0.945 0.940 0.941
	ρ = 0.8.
			regenerative simulation			batch means
		delta jackknife Fieller ratio-ABC F-ABC	5	10	25
	n = 500	0.832	0.847	0.808	0.857	0.864	0.899 0.871 0.841
	n = 1,000 0.876	0.886	0.861	0.893	0.897	0.916 0.899 0.887
	n = 2,000 0.903	0.908	0.892	0.912	0.916	0.929 0.918 0.908
	n = 500	0.879	0.888	0.860	0.897	0.903	0.914 0.902 0.880
	n = 1,000 0.903	0.910	0.893	0.917	0.921	0.933 0.917 0.904
	n = 2,000 0.914	0.917	0.907	0.923	0.925	0.940 0.931 0.921
	n = 500	0.855	0.868	0.832	0.878	0.885	0.904 0.884 0.862
	n = 1,000 0.883	0.891	0.870	0.898	0.904	0.926 0.911 0.894
	n = 2,000 0.907	0.913	0.898	0.916	0.920	0.935 0.927 0.917
					7

URL: http://mc.manuscriptcentral.com/lssp E-mail: comstat@univmail.cis.mcmaster.

ca Communications in Statistics -Simulation and Computation

	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

Table 3 :

 3 Empirical coverage probabilities for tail probability P e (X(•) ≥ 3) and load ρ = 0.5; simulation set-up is otherwise as described in Table1.

			regenerative simulation			batch means
		delta jackknife Fieller ratio-ABC F-ABC	5	10	25
	n = 500	0.922	0.926	0.921	0.934	0.934	0.938 0.932 0.925
	n = 1,000 0.932	0.934	0.930	0.938	0.938	0.948 0.948 0.942
	n = 2,000 0.938	0.938	0.938	0.942	0.942	0.950 0.951 0.948
	n = 500	0.916	0.917	0.916	0.932	0.932	0.935 0.932 0.926
	n = 1,000 0.936	0.937	0.936	0.942	0.942	0.942 0.942 0.939
	n = 2,000 0.943	0.943	0.944	0.945	0.945	0.948 0.946 0.939
	n = 500	0.922	0.923	0.922	0.939	0.938	0.938 0.932 0.927
	n = 1,000 0.939	0.940	0.937	0.940	0.940	0.943 0.937 0.936
	n = 2,000 0.946	0.946	0.946	0.949	0.949	0.946 0.944 0.941

Table 4 :

 4 Empirical coverage probabilities for tail probability P

			regenerative simulation			batch means
		delta jackknife Fieller ratio-ABC F-ABC	5	10	25
	n = 500	0.917	0.928	0.910	0.926	0.926	0.944 0.935 0.914
	n = 1,000 0.935	0.939	0.930	0.936	0.936	0.945 0.945 0.938
	n = 2,000 0.943	0.946	0.939	0.944	0.945	0.942 0.944 0.940
	n = 500	0.929	0.934	0.919	0.933	0.934	0.939 0.939 0.931
	n = 1,000 0.940	0.943	0.935	0.944	0.944	0.947 0.946 0.940
	n = 2,000 0.935	0.938	0.935	0.941	0.941	0.948 0.948 0.944
	n = 500	0.929	0.934	0.921	0.933	0.933	0.943 0.935 0.921
	n = 1,000 0.939	0.943	0.930	0.941	0.941	0.949 0.945 0.936
	n = 2,000 0.945	0.948	0.942	0.947	0.947	0.950 0.951 0.948
					8

e (X(•) ≥ 3) as in Table

3

, but with load ρ = 0.8.