Introduction

HIV / AIDS infection is an epidemic that has affected the whole world in a relevant way since the early eighties [1]. This disease has led to high morbidity and mortality in infected patients. [2]. The epidemiology of infection has changed over the years in relation to different factors and especially with the introduction of the combination antiretroviral therapy (cART) from 1997. The most relevant change found with the introduction of cART has been the decrease in opportunistic infections (OI)

Toxoplasmosis is the commonest central nervous system (CNS) infection in patients with HIV/AIDS who are not on appropriate prophylaxis. The epidemiology of Toxoplasmosis encephalitis in patients with cART hs changes in Europe.

Trends for toxoplasmosis in HIV patients admission have not been examined in Spain with nationwide population-based from of the Spanish Ministry of Health (SMH). The Spanish National Hospital Discharge Database (SNHDD) belongs to the SMH and includes data from all patients discharged from the public and privates hospital [11],[12].

The objective of this manuscript is to describe the change in the of HIV hospital admissions due to toxoplasmosis and association with other OI and mortality.

Goals

The objective of this paper is to exploit the database of HIV hospital admissions in Spain during the period 1997-2015. The main goals are:

Database

We conducted an observational retrospective study using the SNHDD [11],[12]. The SNHDD includes information on the sex, age, dates of admission and discharge, up to 14 discharge diagnoses, and up to 20 procedures performed during the hospitalization. Data were collected from January 1, 1997, to December 31, 2015 (19 years). Country of origin is not available in the SNHDD. The SMH conducts regular audits to assess the accuracy of the SNHDD [11],[13]. The database can be accessed upon request and data are only anonymously given.

The International Classification of Diseases-Ninth Revision, Clinical Modification (ICD-9-CM) is used in the SNHDD. which is the one used by the Spanish SNHDD. We selected hospital admissions for patients with a diagnosis within the ICD-9-CM code 042 [Human immunodeficiency virus [HIV] disease] and V08 [Asymptomatic human immunodeficiency virus [HIV] infection status] from 1997 to 2015 at any position in the diagnostic list for each episode of hospital admission.

Material and methods

The available data consists on dessagregated data from VIH individuals being admitted to hospitals during the period 1997-2015. The database does not contain any identifier for each person, so it is impossible to assume independency between data: the same person probably will have several entries to the hospital during a year and also along the years during the period of interest.

Given the high variability in available data, coming from individuals and also political regions in Spain, in order to study trends along the years, data will be aggregated by years, in terms of counts (for categorical variables) or means (for continuous variables). For visualization, plots of means for continuous variables and ratios for count data are used.

We’ll focus on studying differences due to gender, and sometimes region (acting as main effects), as well as its evolution along the years under the period of interest. Evolution of the observed trend for different groups (Males/Females) along the years will be studied throug interaction effects. When the response variable is continuous, linear regression models will be fitted. When the response variable consists on counts or proportions, generalized linear models will be fitted, assuming a binomial likelihood (and link logit: \(log(p/(1-p))\)) is the response is dicotomous, or poisson (and link \(log\)) when it is not associated to a dicotomous variable. The independency assumptions of linear and generalized linear models will lack irreversibly, as commented before, but will be undertaken.

Variable Year will be treated as covariable in order to provide the effect of change in time along the period of interest. It will be differenced to start in t=1 for 1997 and finish at t=15 por year 2015, in order to provide a clear effect of a unit change in time on the response (it will be referred to as NYear.from.1996=number of years from 1996). When the observed trend time-response is non-linear, the year variable will entry the model through a cubic spline transformation. A cubic spline model is a nonparametric regression technique defined piecewise by polynomials, and used in the literature to address longitudinal data fitting problems, like those we are concerned in this paper. It has been profusely used in applied statistic studies, even in Medicine (see Durrleman and Simon, 1989, or Desquilbet and Mariotti, 2010, among many others). Splines models are adaptive in the sense that fit the data in a continuous way, by dividing data in segments and assessing the fit by taking into acount the data from the closest entries (in our case from the closest years). Knots defining segments are defined in terms of the quantiles of the data. The number of knots, or equivalently the degrees of freedom, must be chosen dependent on the data. A fit where the predictor x enters through a spline model with df=1 is equivalent to the usual linear model; in this case, for simplicity of interpretations, a usual linear model will be fitted. A model with df=2 implies dividing the x values in 2 segments (through the median as a knot) and in some way similar to fit a polynomial of degree 2 at each segment; a model with df=3 implies dividing the x values in 3 segments and fitting something similar to a polynomial of degree 3 at each segment. See details in De Boor (1978) and Welham et al (2007). The natural splines models are implemented with the library splines in R (ns).

Model selection -selection of significative effects and also the number of knots for the splines- will be resolved in terms of step searchs based on the AIC criteria. The AIC, Akaike Information Criterion, (Akaike, 1974), is used as a model selection index that takes into account the likelihood of the proposed model for the fitted model, and also the number of parameters used for the fit, so promoting simpler models that provide enough likelihood. Goodness of fit is reported in terms or R-squared indexes for linear models and significance of the full model based on the F-statistic for the residuals. For generalized linear models, the percentage of deviance explained by the model, 1-Residual.Deviance/Null.Deviance, could reproduce the interpretation of R-squared in the linear model, and it will be reported. Significance of effects will also be reported in this paper in order to justify differences by the variables of interest.

The R statistical package version 3.6.3 (The R Foundation) was used to perform the statistical analysis.

Results

In this study 66.451.094 hospital admissions were recorded in Spain during 19 years (1997-2015). A total of 472.269 HIV hospital admissions were obtained, and 9006 had toxoplasmosis. The prevalence of toxoplasmosis has decreased from 4.2% of HIV admission in 1997 to 0.8% in 2015 (p<0.001). The most common opportunistic co-infections was: tuberculosis (27.5%), esophageal candidiasis (9-.7%), cytomegalovirus diseases (8.5%), kaposii diseases (4.6%). The 70.1% with toxoplasmosis were in drugs usurers, 44.4 had hepatitis c virus infection and 10.4% had hepatitis B infection.

Table 1. Number of cases per year (and percentages in blue).
Year Toxoplasmosis of the
brain 		Total
1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015
No 25162
5.4 % 		25026
5.4 % 		25100
5.4 % 		25728
5.6 % 		26151
5.6 % 		25907
5.6 % 		26390
5.7 % 		26675
5.8 % 		26436
5.7 % 		25071
5.4 % 		25908
5.6 % 		25855
5.6 % 		25110
5.4 % 		22886
4.9 % 		22261
4.8 % 		21094
4.6 % 		21600
4.7 % 		21032
4.5 % 		19871
4.3 % 		463263
100 %
Yes 1115
12.4 % 		836
9.3 % 		713
7.9 % 		645
7.2 % 		602
6.7 % 		517
5.7 % 		506
5.6 % 		500
5.6 % 		480
5.3 % 		453
5 % 		421
4.7 % 		395
4.4 % 		372
4.1 % 		287
3.2 % 		266
3 % 		270
3 % 		237
2.6 % 		234
2.6 % 		157
1.7 % 		9006
100 %
Total 26277
5.6 % 		25862
5.5 % 		25813
5.5 % 		26373
5.6 % 		26753
5.7 % 		26424
5.6 % 		26896
5.7 % 		27175
5.8 % 		26916
5.7 % 		25524
5.4 % 		26329
5.6 % 		26250
5.6 % 		25482
5.4 % 		23173
4.9 % 		22527
4.8 % 		21364
4.5 % 		21837
4.6 % 		21266
4.5 % 		20028
4.2 % 		472269
100 %
χ2=1691.247 · df=18 · Cramer’s V=0.060 · p=0.000

Fig.1 displays the total number of VIH admissions with Toxoplasmosis of the brain, differentiated by gender. In 1997 the number of males was four times bigger than the number of males, and in 2015, it was double the females one. A continuous decrease in absolute numbers is perceived along the period. This reduction is not only due to a decrease in the total number of VIH admissions along the period. Table 1 displays the counts and relative frequencies of Toxoplasmosis diagnoses with respect to the total VIH admissions per year. In this table the descendant trend of Toxoplasmosis prevalence. Also these relative frequencies are displayed in Fig.2, and the descendant trend is observed both for male and female patients.

Fig.3 displays the counts of Toxoplasmosis diagnoses in VIH patients, in every political region in Spain and attending trienial periods from 1997. The descendant trend is common to all regions.

Fig.4 displays the points with the rates of toxoplasmosis diagnosis, assessed with respect to the total VIH admissions each year. A cuadratic glm binomial model is fitted to the data points in order to predict the rates, and it is represented by the line. A continuous descendant trend is observed along the studied period, showing that Toxoplasmosis in the brain is a disease that has considerably evolved in VIH patients, in the sense that its prevalence in 2015 is much smaller than in 1997.

In Table 2 the model fitted to capture the descendant trend in Toxoplasmosis prevalence is displayed, with coefficients (odds ratios), confidence intervales (CI) and p-values (\(<0.001\)), as well as goodness of fit statistics as Deviance and AIC.

Table 2. Predicted model for explaining Toxoplasmosis decrease in the period 1997-2015.
  Toxoplasmosis proportion
Predictors Odds Ratios CI p
(Intercept) 0.02 0.02 – 0.02 <0.001
Year.linear 0.18 0.16 – 0.20 <0.001
Year.cuadratic 1.39 1.26 – 1.53 <0.001
Observations 19
Deviance 74.841
AIC 230.336

Next we study possible differences due to gender in this observed trend along the period of interest.

Figure 5 displays the points with the rates of toxoplasmosis diagnosis, assessed with respect to the total VIH admissions each year, and differentiated by gender A smoothed-glm binary model (of order 3) is fitted to the data points in order to predict the rates, and it is represented by the lines. A continuous descendant trend is observed along the studied period for males. For females, we notice some stationarity between 2003 and 2010; after 2010 a descendant trend is also evident for female VIH patients.

Next we study the evolution in mean age of VIH patients admitted to hospitals, and diagnosed with Toxoplasmosis at admission. We fit a linear regression to variable Age, including the sex effect and also differetiating by Toxoplasmosis diagnosis. In mean terms we can appreciate (see Fig.6), that both for males and females, the mean age of HIV patients with a Toxoplasmosis diagnosis at admission has decreased with respect to the mean age of patients with other diagnoses. The fitted linear regression provides significative differences at 95% (see Table 3) for the slope of the prediction lines (which show the trend with years) among the four groups emerged from crosses: male/Female and with/without Toxoplasmosis diagnosis. The \(R^2\) statistic from regression is so small due to the high variability of data. However this variability, the observed trend for the means is clear and clearly fitted by the model.

Table 3. Fitted model to predict mean age for HIV patients admitted in hospitals along 1997-2015.
  Mean Age
Predictors Estimates CI p
(Intercept) 34.01 33.94 – 34.08 <0.001
Year 0.82 0.81 – 0.82 <0.001
SexFem -3.99 -4.12 – -3.86 <0.001
Toxoplas 1.40 0.99 – 1.81 <0.001
Year:SexFem 0.04 0.03 – 0.05 <0.001
Year:Toxoplas -0.30 -0.35 – -0.26 <0.001
SexFem:Toxoplas 2.28 1.44 – 3.12 <0.001
Year:SexFem:Toxoplas -0.10 -0.19 – -0.01 0.023
Observations 472269
R2 / R2 adjusted 0.190 / 0.190

Evolution of Re-entry

Next we study the proportion of VIH patients which enter the hospital due to a re-entry within a period of 30 days from the last entry, and try to stablish some trend/behavior for the patients diagnosed with Toxoplasmosis versus those diagnosed with other diseases. The gender effect is also considered. Fig.7 displays the percentage of re-entries with respect to total VIH admissions, for male and female patients, as well as Toxoplasmosis diagnosed and others. While a stable behavior in re-entries is perceived for VIH patients with other diagnoses than Toxoplasmosis (red lines), the behavior in patients with Toxoplasmosis is really variable along the period. A glm binomial model has been fitted to capture the trend. However the high variability, estimates of re-entries for Toxo-patients are higher than for non-Toxo-patients.

Table 4 shows the fit. Just the main effects and the interaction Year:Sex is significant. This justify that lines fitting re-entry incidence are parallels for males and also for females, and do not differentiate the effect of Year on Toxo-patients and non-Toxo-patients. The difference in height of these regression lines explains a higher re-entry incidence in Toxo-patients.

Table 4. Fitted model to predict Re-Entry for HIV patients admitted in hospitals along 1997-2015.
  Re-Entry incidence
Predictors Odds Ratios CI p
(Intercept) 0.20 0.20 – 0.20 <0.001
Year 1.00 1.00 – 1.00 0.016
SexFem 0.94 0.91 – 0.98 0.001
Toxoplas 1.30 1.23 – 1.37 <0.001
Year:SexFem 1.00 0.99 – 1.00 0.010
Observations 76
R2 Tjur 0.039
Deviance 159.668
AIC 714.536

Evolution of Exitus rate

Next we study the exitus rate for VIH patients diagnosed with and without Toxoplasmosis. In Fig.8 we can observe that the exitus rate is considerably higher for VIH diagnosed with Toxoplasmosis of the brain, that globally for all the other diagnoses, both for males and females. However it is remarkable the high variability in the exitus rate for patients with Toxoplasmosis diagnosis, during all the period under study. A glm binomial model has been fitted to capture the trend. However the high variability of the exitus rates for Toxo-patients, these are higher than for non-Toxo-patients, as shown by the fit, represented by lines.

Table 5 shows the fit. Just the triple interaction and double one Year:Toxoplas dissapear from the fitted model. The year effect vary (in slope) with Sex, and also Sex with presence or absence of Toxoplasmosis. The difference in height of these regression lines explains a higher exitus incidence in Toxo-patients.

Table 5. Fitted model to predict Exitus incidence for HIV patients admitted in hospitals along 1997-2015.
  Exitus rate
Predictors Odds Ratios CI p
(Intercept) 0.08 0.08 – 0.08 <0.001
Year 0.98 0.98 – 0.99 <0.001
SexFem 0.59 0.55 – 0.62 <0.001
Toxoplas 2.03 1.88 – 2.19 <0.001
Year:SexFem 1.01 1.00 – 1.02 0.001
SexFem:Toxoplas 1.43 1.22 – 1.66 <0.001
Observations 76
R2 Tjur 0.064
Deviance 122.064
AIC 620.312

Association Toxoplasmosis

In order to investigate which diseases are more correlated with Toxoplasmosis of the brain as diagnostiscs in hospital admissions for the VIH patients, we assess the number of cases that two diseases are diagnosed jointly for any patient admitted during the period of interest.

Fig.9 shows the relative frequencies for associations with all diagnoses considered. The most common for patients diagnosed with Toxoplasmosis is “Use of drugs by parentals” (pardrug), followed by hepathitis and pneumonia.

Table.2. Cases per 10.000 where Toxoplasmosis of the brain was jointly diagnosed with other disease at admission. Computed on the VIH admissions 1997-2015.
Disease Cases Cases/10.000 Disease Cases Cases/10.000
pardrug 3310 70.09 mycobac 133 2.82
hepac 2099 44.45 enceph 100 2.12
tuberc 1300 27.53 crytococ 81 1.72
pneumo 640 13.55 leishman 40 0.85
hepab 491 10.40 crytospor 26 0.55
candieso 459 9.72 lymburkit 18 0.38
cytomega 401 8.49 isospor 17 0.36
pneumocys 317 6.71 candiresp 12 0.25
progleukoen 283 5.99 salmon 9 0.19
wasting 280 5.93 lymbrain 6 0.13
kaposi 216 4.57 histoplas 1 0.02
herpes 163 3.45 cervican 0 0.00

References

De Boor, C. (1978). A practical guide to splines. New York : Springer. https://www.springer.com/gp/book/9780387953663

Welham, S.J., Cullis, B.R., Kenward, M.G., Thompson, R. (2007). A comparison of mixed model splines for curve fitting. Australian and New Zealand Journal of Statistics, 49,1, 1-23. DOI https://doi.org/10.1111/j.1467-842X.2006.00454.x

Desquilbet, L. and Mariotti, F. Dose‐response analyses using restricted cubic spline functions in public health research, Statistics in Medicine, 29, 9, (1037-1057), (2010). DOI: https://doi.org/10.1002/sim.3841 Durrleman, S, and Simon, R. Flexible regression models with cubic splines. Statistics in Medicine, 8, 5, (551-561), (1989). DOi: https://doi.org/10.1002/sim.4780080504

Akaike, Hirotugu (1974), «A new look at the statistical model identification», IEEE Transactions on Automatic Control 19 (6): 716-723, MR 0423716, doi:10.1109/TAC.1974.1100705.

The R foundation. Available on https://www.r-project.org/