Multinomial Logistic Model for Coinfection Diagnosis Between Arbovirus and Malaria in Kedougou

3.1 Multinomial logit model

We recall that Y is the response variable indicating the class of the disease: “Other febrile illnesses” (Y = 0), “arboviral monoinfection” (Y = 1), “malaria monoinfection” (Y = 2) and “coinfection” (Y = 3). Let X=(1,X1,…,Xp) be the vector of the p covariates. For an individual with covariates X = x, we want to predict the probability of belonging to the class k given x,

The multinomial logit model assumes the existence of β₁, β₂, β3∈Rp+1 such that, for each k = 1, 2, 3 and each vector of covariates x,

where

and x₀ = 1 to include the intercept parameters β_k0, k = 1, 2, 3. The reference modality is class 0.

Consequently, for each k = 1, 2, 3 and each vector of covariates x,

and

From the computation of the maximum likelihood estimates βˆk, we derive for k = 1, 2, 3,

3.2 Fitting strategy for handling imbalanced IgM data

The IgM data set contains 18 arboviral monoinfection cases, 21 coinfection cases, 5 180 other febrile illness cases and 7 069 malaria monoinfection cases. Trained on the original IgM data set, the fitted logit model only predicted classes 0 and 2, which means it ignores the two minority classes 1 and 3 in favour of the majority classes. Applying resampling strategies to obtain a more balanced data sample is an effective solution to the imbalance problem (see [16] for a survey of existing methods). Two of the most simple resampling approaches are undersampling and oversampling. Since the IgM is highly imbalanced with a large number of observations in the two majority classes, we used a random undersampling strategy that removes observations and reduces the sample size. We sampled without replacement 50 cases from each of the two majority classes to create a balanced sub-sample of size 18+21+50+50=139. Trained on a sub-sample, the model predicted four classes.

Undersampling results in loss of information and the risk of removing relevant observations is present. To overcome this problem, we repeated the sampling step a thousand times and worked with 1 000 balanced sub-samples of the IgM data set. The multinomial model was fitted to each sub-sample and a stepwise covariate selection was performed (see Figure 5 in the supplementary material). The observed variability of the 1 000 covariate selections raised robustness questions. To answer this point, we conducted a nonparametric analysis based on the RF algorithm. In recent years, several methods involving the combination of resampling and ensemble learning have appeared in the imbalanced distributions literature ([16]). We found that the importance score based on random forests yielded a convenient way to summarize the information obtained from the 1 000 sub-samples.

3.3 Variable selection using random forests

A random forest is an ensemble of unpruned trees, induced from bootstrap samples of the training data, that uses random covariate selection in the tree construction process. Prediction is made by aggregating the predictions of the ensemble, using the majority vote rule.

One of the most widely used RF score of importance of a given variable is the Mean Decrease of Accuracy (MDA) in predictions. It is based on the out-of-bag (OOB) error. For each tree t of the forest, consider the associated OOB_t sample (data not included in the bootstrap sample used to construct t). Denote by errOOB_t the misclassification rate of tree t computed on this OOB_t sample. Then, randomly permute the observed values of covariate X_j in OOB_t to get a perturbed sample and compute errOOBtj, the error of t on the perturbed sample. Variable importance of X_j is then given by

where ntree denotes the number of trees of the RF. The higher the MDA, the more important the variable is. Several variable selection procedures using RF are based on this quantification of variable importance.

Using R packages, we made the following implementation choices: randomForest for RF fitting and MDA calculation, VSURF for selecting the important variables. The main parameters of randomForest were calibrated and set to their default values, ntree=500 and mtry=p=3 (number of variables tried at each split of a tree of the RF). The variable selection strategy of VSURF is based on a two-stage procedure [17]: 1. the covariates are ranked by sorting their variable importance measures in descending order and the covariates whose importance is less than a threshold (the minimum value of the standard deviations of the importance measures) are eliminated; 2. a sequence of nested models starting from the one with only the most important variable and ending with the one involving all important variables kept previously is considered; the variables of the model leading to the smallest OOB error are selected. An advantage of using VSURF is that this procedure does not require the choice of tuning parameters.

Figure 2 ranks the variable importances (MDA) of the 15 covariates across the 1 000 sub-samples. First, rainfall is the most important covariate; a second group of less important covariates is formed by cough, age and joint pain; then comes a group of five covariates: number of sick days, temperature, nausea or vomiting, eye pain and nasal congestion; finally, six unimportant covariates are displayed: muscle pain, chills, cephalalgia, jaudice, diarrhea and sex . The boundary between the two last groups is not clear and we used the VSURF procedure to separate the important covariates from the other ones. We can notice on the plot that both MDA level and variability are larger for relevant variables; as explained by [14], this is expected and the VSURF threshold value is based on MDA standard deviation estimation. Figure 3 summarizes the results of the VSURF selection procedure based on the 1 000 sub-samples. The covariate rainfall (95.2%) is almost always selected. Next, the more often selected variables are cough (29.1%), age (28.3%), joint pain (19.8%), nausea or vomiting (16.4%), number of sick days (16.1%), temperature (16.1%) and nasal congestion (11%), in decreasing order. The other covariates are selected in less than 10% of the samples.

We set different random seeds and we found that, for our purpose of selecting significant covariates, aggregation of 1 000 RF classifiers learned from 1 000 randomly balanced sub-samples yielded stable selected variable sets.

Figure 2:

A variable importance plot for the IgM data set. Each boxplot summarizes the distribution of the variable importance among 1000 IgM sub-samples.

Figure 3:

Ranking by VSURF: for each variable, the length of the bar corresponds to the empirical probability to be selected by VSURF among 1000 IgM sub-samples.

3.4 Influence of selected covariates on disease status

In the previous sections, the RF variable importance results on the IgM sub-samples produced a robust ranking of the covariates. From these results, we decided to fit multinomial model with eight covariates (age, temperature, number of sick days, rainfall, nausea or vomiting, cough, nasal congestion and joint pain) to the data set of our analysis and to further quantify the effects of the covariates in this model.

Within the multinomial logit model, we can quantify the effect of a variable in terms of an odds ratio or its logarithm. The odds that Y = k occurs for an individual with covariates X = x is the ratio of P(Y=k|X=x) divided by P(Y=0|X=x), k = 1, 2, 3. Then, the log odds of category k is given by eq. (1):

Thus the multinomial logit model is a linear regression model in the log odds. The parameter component β_kj can be interpreted as the change in the log odds per unit change in the continuous covariate X_j, if all other covariates are held constant. The odds ratio (OR) of category k for a d units increase of X_j, all other covariates remaining constant, is defined as

Once β is estimated, one can estimate any odds or odds ratios. An OR equal to one means that a change in covariate X_j has no effect on the odds of category k; if OR_k(d) > 1 (OR_k(d) < 1), the effect of an increase of X_j is to increase (decrease) the odds of category k. The risk ratio P(Y=k|Xj+d)/P(Y=k|Xj), which could be more interpretable in terms of predicted probabilities instead of odds, depends on the values of all other covariates. ORs are similar to risk ratios if the risk is small, otherwise ORs overestimate risk ratios.

For each covariate, we computed the odds ratios OR_k, k = 1, 2, 3 and their confidence intervals for each disease. Figure 4 display the OR by which the odds increases for a certain change in a covariate, holding all other covariates constant. The ORs associated with binary variables (nausea/vomiting, cough, nasal congestion and joint pain) were computed by comparing the two modalities: 0 for absence and 1 for presence of the symptom. We computed the ORs resulting from increasing temperature from 38 to 40 degrees Celsius (d = 2) and from increasing Number of sick days from 2 to 6 days (d = 4). The outer quartiles of Age are 8 and 28 years (d = 20), so we computed the half-sample OR for age. Similarly, we computed the half-sample OR for a rainfall of 14 mm compared to a rainfall of 370 mm (d = 356).

Figure 4:

IgM data: boxplots of 1000 odds ratios with respect to the reference category (grey full line boxplot); boxplots of the associated confidence intervals are shown in dotted line (green for upper and blue for lower bound); (a) Arbovirus (b) Coinfection (c) Malaria.

The ORs defined previously are relative to the reference category Y = 0. We also computed the ORs between two diseases Y = k and Y = l in order to differentiate the effect of each covariable between the three clinical groups, arbovirus vs malaria, coinfection vs arbovirus and coinfection vs malaria (Figure 5):

Figure 5:

IgM data: boxplots of 1000 odds ratios between two categories (grey full line boxplot); boxplots of the associated confidence intervals are shown in dotted line (green for upper and blue for lower bound); (a) Arbovirus vs Malaria (b) Coinfection vs Arbovirus (c) Coinfection vs Malaria.

The confidence intervals are derived from the fitted multinomial logit model and their accuracy is based on the parametric assumption that the true data generating distribution does fall in the model.

Figure 4 and Figure 5 display the sampling distribution of ORs based on the fitting of the 1000 sub-samples of the IgM data set. According to Figure 4, we can say that rainfall and vomiting symptoms are hightly correlated with malaria monoinfections whereas joint pain is correlated with arboviral monoinfections. The odds of coinfection increases with high fever. It corroborates the conclusion of the paper Sow et al. [13]. Figure 5(a), (b) and (c) can be interpreted in the same way. They show that a high temperature and the presence of nausea or vomiting symptoms are mostly indicative of malaria parasite infections whereas an increase of age and of number of sick days are indicative of arboviral infections. The effects of nasal congestion and joint pain symptoms on the disease status are not clear enough to be interpreted. The main question of the study was to identify risk factors that can help doctors to diagnose a concurrent malaria and arbovirus infection. From these results, temperature is the only risk factor that differentiates between coinfection and single infections.

4.1 Testing independence between arbovirus and malaria

In the multinomial model given by (1) in Section 3.1, we can test the independence between arboviral and malaria infections.

The joint statistical distribution of arboviral infection (Y∈{1,3}) and malaria infection (Y∈{2,3}), is given in Table 2. Independence between arboviral and malaria infections means that for all (a,m)∈{0,1},

which is equivalent to

for all (a,m)∈{0,1}, where PY∈{1,3}|X=x corresponds to the probability to belonging of categories 1 or 3, PY∈{2,3}|X=x corresponds to the probability of categories 2 or 3. The independence hypothesis can be written in terms of parameters as:

The Wald statistic to test H₀ against its two-sided alternative is computed as

with h(βˆ)=βˆ3−βˆ1−βˆ2 and Σ = DVD^T where D=(−Idp+1,−Idp+1,Idp+1); Id_p is the p × p identity matrix and V is an estimator of the variance of βˆ=(βˆ1,βˆ2,βˆ3)T. Under H₀, W is asymptotically distributed as a chi-square variable with (p + 1) degrees of freedom. Under H₁, W converges to infinity as the sample size goes to infinity.

4.2 Diagnosis of arboviral disease

In absence of rapid arbovirus detection tests, the aim is to provide a decision support tool to determine if an arbovirus could be responsible for the clinical symptoms of the patient coinfection. We propose to base the diagnosis on the conditional probability to be coinfected q(x)=PY=3|Y∈{2,3},X=x given that malaria infection is observed. This probability is the quantity of interest because arboviral infections are considered by healthcare workers only if malaria tests are negative.

In the previous section, it is shown that we can test the independence hypothesis between malaria and arboviral infections. If this test is rejected, then we can derive the probability q to be coinfected given that malaria infection is observed. This probability can be computed in function of the π_k probabilities estimated from the multinomial logit model. For X = x,

This probability can be used to differentiate whether the illness to be treated should be arbovirus or malaria. We propose a binary classification rule and we predict an arbovirus illness if the estimated coinfection probability is greater than a threshold value γ:

The evaluation of the classification is based on the confusion matrix and the overall classification accuracy. The confusion matrix is used to compute true arbovirus positives (TP), false arbovirus positives (FP), true negatives (TN) and false negatives (FN). A global performance measure is the miss-classification rate (MCR) defined as:

with N=TP+FP+TN+FN.

Our analysis is based on a real-life medical data set. In the original IgM data set, arbovirus positive individuals are identified as individuals likely to be in the early stages of arbovirus illness. It is the relevant data set for the classification problem. However, the positive cases constitute only a very small minority class of the data (39 positive cases over 12288 individuals). Based on these data, the computation of the independence test is very sensitive to the fluctuations of the sub-sampling procedure and the classification procedure could not be implemented. Instead, we propose a simulation study based on a balanced data set to illustrate our classification procedure.

4.3 Simulation study

Simulated data. Taking advantage of the previous influencing factors analysis (Section 3.4), we simulated data using a multinomial model similar to the one previously estimated. The eigth covariates were generated from distributions similar to those observed in the real data set. The beta parameters values are given by Table 7 in supplementary material. They were chosen according to the conclusions of the statistical analysis of the influential factors. The larger values emphasize the influence of the associated covariates that are positively correlated to each disease category. For example, the parameter value associated with the number of sick days covariate is larger for the arbovirus category than for the malaria category. Based on this generative model, we computed the probabilities of belonging to each category and generated the Y response to be the modality with the greatest probability. We used this procedure to simulate a data set of size n = 5000 which is summarized in Table 6 and Table 8 in the supplementary material.

Independence between abovirus and malaria. We fitted the multinomial model to the simulated data and tested the independence between malaria and arbovirus. We obtained that the independence hypothesis was rejected with a p-value equal to 1.13 × 10^– 4. Then we derived the probability to be coinfected given that malaria infection is observed and performed the classification procedure.

Diagnosis of arboviral disease. We randomly divided the simulated data set into two part, a training data set of size 3333 and a test data set of size 1667. The classification was applied only to individuals infected with malaria parasites, namely 1925 individuals in the training data set and 626 individuals in the test set. We computed the five-fold cross-validation estimator of the MCR and we chose the classification threshold value γ as the minimizer of the MCR. We can see on Figure 6 that the optimal value of this threshold is γ = 0.45. Five-fold cross-validation was run several times and the optimal value of γ was found to be quite stable. Based on this γ value, we performed the classification to predict the type of illness that has affected the patient. Predicted and actual arbovirus cases were compared using the test set, as presented in Table 3. The rows of the matrix are actual classes and the columns are the predicted classes. We observe that the corresponding test MCR is 7.83%. The ROC curve of the classification is presented in Figure 7 of the supplementary material. Based on the simulated data set, the accuracy of the classification is quite good (92.17%). This suggests that this predictive analysis can be medically valuable to identify arboviral cases among coinfection cases.

Figure 6:

Cross-validation miss-classification rate. The MCR is shown in black as full line. Increasing γ increases the number of FN (green line) and decreases the FP (red line).