Exercice

In a public health study investigating the impact of environmental factors, socio-economic conditions and healthcare interventions on malaria incidence, data were collected from several health districts during a peak transmission period. Each observation describes local environmental and living conditions, a specific medical intervention and the resulting malaria incidence rate, measured in cases per 1,000 inhabitants. This practical exercise aims to model and explain malaria incidence using a multiple linear regression approach based on available environmental, socioeconomic and prevention variables.

The dataset includes the following variables:

  1. Rainfall_mm: Total rainfall in the district (mm)

  2. Average_temp_C: Average monthly temperature (°C)

  3. Humidity_%: Average relative humidity (%)

  4. Stagnant_water_index: Index of stagnant water surface area (scale from 0 to 10)

  5. Bednet_coverage_%: Percentage of the population regularly using insecticide-treated bed nets (%)

  6. Healthcare_distance_km: Average distance to the nearest health center (km)

  7. Poverty_rate_%: Percentage of the population living below the poverty line (%)

  8. Malaria_incidence_rate: Number of confirmed malaria cases per 1,000 inhabitants

Production

Packages and data importation

Variables recoding

Here is the part where we will rename variables if need and add the labels for these variables

## 'data.frame':    200 obs. of  8 variables:
##  $ Rainfall_mm           : num  179 130 153 177 171 ...
##  $ Average_temp_C        : num  32.2 26.2 30.1 30.5 28.1 ...
##  $ Humidity_pct          : num  73.8 67.1 64 64.7 71.6 ...
##  $ Stagnant_water_index  : num  7.17 5.81 6.54 5.68 4.39 6.28 5.54 6.89 8.17 3.38 ...
##  $ Bednet_coverage_pct   : num  21.7 37.5 39.5 68.3 25.3 ...
##  $ Healthcare_distance_km: num  20.33 6.92 19.61 7.25 18.16 ...
##  $ Poverty_rate_pct      : num  25.5 45.9 67.8 41 40 ...
##  $ Malaria_incidence_rate: num  107.7 72.7 110 83.8 87.7 ...

Exploratory data analysis

Visualisation of all the variables

The following figure shows the histograms of all the independant variables.

From this visualisation, we can notice that, the rainfall, average temperature, humidity and stagnant water index seams to follows normal distribution. But for the others, by using shapiro-test we notice that the p-value was less that 0.05, the we rejected the null hypothesis and confirm that, they don’t follow normal distribution.

hist(Malaria[-8], main = "Histograms for all the variables in the dataset")

Here you can find the results for the shapiro test

(Rain <- shapiro.test(Malaria$Rainfall))
## 
##  Shapiro-Wilk normality test
## 
## data:  Malaria$Rainfall
## W = 0.99453, p-value = 0.6794
(Temp <- shapiro.test(Malaria$temperature))
## 
##  Shapiro-Wilk normality test
## 
## data:  Malaria$temperature
## W = 0.99622, p-value = 0.9064
(Hum <- shapiro.test(Malaria$Humidity))
## 
##  Shapiro-Wilk normality test
## 
## data:  Malaria$Humidity
## W = 0.99337, p-value = 0.5102
(SWI <- shapiro.test(Malaria$Stagnant_water_index))
## 
##  Shapiro-Wilk normality test
## 
## data:  Malaria$Stagnant_water_index
## W = 0.98922, p-value = 0.1369
(Bnet <- shapiro.test(Malaria$Bednet))
## 
##  Shapiro-Wilk normality test
## 
## data:  Malaria$Bednet
## W = 0.92625, p-value = 1.717e-08
(Dist <- shapiro.test(Malaria$distance))
## 
##  Shapiro-Wilk normality test
## 
## data:  Malaria$distance
## W = 0.9322, p-value = 5.073e-08
(Povty <- shapiro.test(Malaria$Poverty))
## 
##  Shapiro-Wilk normality test
## 
## data:  Malaria$Poverty
## W = 0.95242, p-value = 3.253e-06
Particular attention for the dependant variable
hist(Malaria$incidence_rate, col = "skyblue", probability = TRUE, ylim = c(0,0.025), main = "Histogram of malaria incidence rate with the bellshaped curve", xlab="Malaria incidence ratio")
curve(dnorm(x, mean = mean(Malaria$incidence_rate), sd(Malaria$incidence_rate)), col = "red", lwd = 2, add = TRUE)
legend("topright", legend = c("Bellshaped curve", "p-value = 0.3066"), col=c("red","skyblue"), lwd =2)

Descriptive statistics
Malaria  %>%  tbl_summary(
  include = everything(-c(5,6,7)),
  statistic = all_continuous()~ "{mean} ± {sd}"
)
Characteristic N = 2001
Total rainfall in the district (mm) 183 ± 49
Average monthly temperature (°C) 27.79 ± 2.55
Average relative humidity (%) 73 ± 10
1 Mean ± SD
Table 1 : Descriptive statistics for rainfall, temparature, humidity and stagnant water.
Variable Minimum Mean Standard deviation Maximun Normality test
Total rainfall in the district (mm) 40.91 182.64 48.57 302.62 0.6794
Average monthly temperature (°C) 20.17 27.79 2.55 35.1 0.9064
Average relative humidity (%) 45.58 72.88 10.04 95 0.5102
Index of stagnant water surface area (scale from 0 to 10) 1.51 6.47 1.79 10 0.1369
Table 2: Descriptive statistics for the use of bednets, distance to health and poverty.
Variable Minimum 1st quatile median 2nd quartile Normality test
Percentage of the population regularly using insecticide-treated bed nets (%) 20.74 32.23 51.72 73.36 1.72^{-8}
Average distance to the nearest health center (km) 2.12 7.1 10.45 14.8 5.07^{-8}
Percentage of the population living below the poverty line (%) 15.38 28.26 45.59 57.82 3.2526159^{-6}
Box-plots

The following diagrams shows the boxplots for all the variables inside the dataset. One can notice that there are some outliers. We observe them for the variables rainfall, incidence rate, temparature and distance. After investigation, we can see that they are not errorneous. So we will keep it and analyse them at the end.

par(mfrow =c(2,2))
boxplot(Malaria[,1])
boxplot(Malaria[,c(3,5,7,8)])
boxplot(Malaria[,c(2,4,6)])

Part for linear regression

In this part, we are going to model and malaria incidence using a multiple linear regression approach based on available environmental, socioeconomic and prevention variables

The model will be

incidence_rate = \(\beta_0\) + \(\beta_1\) × Rainfall + \(\beta_2\) × temperature + \(\beta_3\) × Humidity+ \(\beta_4\) × Stagnant_water_index + \(\beta_5\) × Bednet + \(\beta_6\) × distance + \(\beta_7\) × Poverty + \(\epsilon\)

Verification of the correlation between variables

The figure above present the correlation between all the variables in the dataset with the scatterplots. Pearson’s correlation analysis revealed strong positive correlations among the environmental variables. Rainfall was highly correlated with humidity (r = 0.88) and the stagnant water index (r = 0.78), whereas humidity was also strongly correlated with the stagnant water index (r = 0.73), suggesting potential multicollinearity among these predictors. Temperature showed negligible correlations with the remaining explanatory variables. Malaria incidence was moderately and positively correlated with rainfall (r = 0.37), temperature (r = 0.37), stagnant water index (r = 0.37), humidity (r = 0.32), poverty (r = 0.29), and distance (r = 0.23), but bednet coverage exhibited a weak negative correlation with malaria incidence (r = −0.27). Overall, these findings indicate that climatic factors are moderately associated with malaria incidence, while strong intercorrelations among rainfall, humidity, and stagnant water warrant assessment of multicollinearity before multivariable modeling.

cormat <- cor(Malaria)
pval <- cor.mtest(Malaria)
par(mfrow=c(1,2))
corrplot(cormat,method="number", type = "full")
corrplot(cormat, method="color",type = "upper", p.mat = pval$p, insig = "p-value", number.digits = 3, sig.level = 0.05)

Estimation of the parameters for model with all the variables inside the dataset
Model1 <- lm(incidence_rate~.,data = Malaria)

4. Verification of the assumption to validate the model

4.1. Global significance test of the model

anova(Model1)
## Analysis of Variance Table
## 
## Response: incidence_rate
##                       Df Sum Sq Mean Sq F value    Pr(>F)    
## Rainfall               1  15297 15297.1 50.0879 2.686e-11 ***
## temperature            1  13799 13798.7 45.1817 1.995e-10 ***
## Humidity               1     89    89.3  0.2924 0.5892852    
## Stagnant_water_index   1   1649  1648.5  5.3979 0.0212087 *  
## Bednet                 1   7022  7022.0 22.9925 3.254e-06 ***
## distance               1   4812  4811.8 15.7554 0.0001018 ***
## Poverty                1   8343  8343.5 27.3193 4.478e-07 ***
## Residuals            192  58638   305.4                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.2. Significance test of regression coefficients

summary(Model1)
## 
## Call:
## lm(formula = incidence_rate ~ ., data = Malaria)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -72.554  -8.504   1.042   9.341  53.755 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -48.12435   19.66504  -2.447  0.01530 *  
## Rainfall               0.11129    0.05913   1.882  0.06135 .  
## temperature            3.11819    0.49166   6.342 1.59e-09 ***
## Humidity              -0.21322    0.26746  -0.797  0.42631    
## Stagnant_water_index   3.03431    1.11948   2.710  0.00733 ** 
## Bednet                -0.24267    0.05683  -4.270 3.07e-05 ***
## distance               0.99401    0.21930   4.533 1.02e-05 ***
## Poverty                0.37611    0.07196   5.227 4.48e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.48 on 192 degrees of freedom
## Multiple R-squared:  0.4652, Adjusted R-squared:  0.4457 
## F-statistic: 23.86 on 7 and 192 DF,  p-value: < 2.2e-16

4.3. Multicollinearity assessment among explanatory variables

Variance Inflation Factors (VIFs) were computed to assess multicollinearity among the explanatory variables. VIF values ranged from 1.02 to 5.37, indicating generally low multicollinearity. Rainfall (VIF = 5.37) and humidity (VIF = 4.69) exhibited moderate multicollinearity, reflecting their strong pairwise correlation (r = 0.88). However, all VIF values were well below the commonly accepted threshold of 10, suggesting that multicollinearity was not severe and that all predictors were retained in the final model.

vif(Model1)
##             Rainfall          temperature             Humidity 
##             5.374762             1.021030             4.694208 
## Stagnant_water_index               Bednet             distance 
##             2.614735             1.039260             1.022684 
##              Poverty 
##             1.018634

When we used `step()` to variable selection, the process excluded the humidity but this wasn’t have any difference with the model with all the variable (p.value=0.4263). So, we decided to keep all the variables for prédiction.

Model2 <- step(Model1, direction = "both")
## Start:  AIC=1152.16
## incidence_rate ~ Rainfall + temperature + Humidity + Stagnant_water_index + 
##     Bednet + distance + Poverty
## 
##                        Df Sum of Sq   RSS    AIC
## - Humidity              1     194.1 58832 1150.8
## <none>                              58638 1152.2
## - Rainfall              1    1081.7 59720 1153.8
## - Stagnant_water_index  1    2243.7 60882 1157.7
## - Bednet                1    5568.0 64206 1168.3
## - distance              1    6274.7 64913 1170.5
## - Poverty               1    8343.5 66981 1176.8
## - temperature           1   12284.3 70922 1188.2
## 
## Step:  AIC=1150.82
## incidence_rate ~ Rainfall + temperature + Stagnant_water_index + 
##     Bednet + distance + Poverty
## 
##                        Df Sum of Sq   RSS    AIC
## <none>                              58832 1150.8
## + Humidity              1     194.1 58638 1152.2
## - Rainfall              1    1095.9 59928 1152.5
## - Stagnant_water_index  1    2082.8 60915 1155.8
## - Bednet                1    5376.6 64208 1166.3
## - distance              1    6217.6 65050 1168.9
## - Poverty               1    8327.8 67160 1175.3
## - temperature           1   12869.4 71701 1188.4
summary(Model2) 
## 
## Call:
## lm(formula = incidence_rate ~ Rainfall + temperature + Stagnant_water_index + 
##     Bednet + distance + Poverty, data = Malaria)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -72.777  -9.213   0.338   9.302  52.399 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -58.09070   15.16544  -3.830 0.000173 ***
## Rainfall               0.07709    0.04065   1.896 0.059439 .  
## temperature            3.16688    0.48739   6.498 6.78e-10 ***
## Stagnant_water_index   2.87877    1.10131   2.614 0.009656 ** 
## Bednet                -0.23519    0.05600  -4.200 4.08e-05 ***
## distance               0.98908    0.21900   4.516 1.09e-05 ***
## Poverty                0.37575    0.07189   5.227 4.46e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.46 on 193 degrees of freedom
## Multiple R-squared:  0.4635, Adjusted R-squared:  0.4468 
## F-statistic: 27.78 on 6 and 193 DF,  p-value: < 2.2e-16
vif(Model2) 
##             Rainfall          temperature Stagnant_water_index 
##             2.545299             1.005276             2.535319 
##               Bednet             distance              Poverty 
##             1.010973             1.021871             1.018593
anova(Model1, Model2)
## Analysis of Variance Table
## 
## Model 1: incidence_rate ~ Rainfall + temperature + Humidity + Stagnant_water_index + 
##     Bednet + distance + Poverty
## Model 2: incidence_rate ~ Rainfall + temperature + Stagnant_water_index + 
##     Bednet + distance + Poverty
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1    192 58638                           
## 2    193 58832 -1    -194.1 0.6355 0.4263

4.4. Residual analysis

Residuals <- residuals(Model1) 
Fitted <- fitted(Model1)

4.4.1. Residual normality test checking

H0: Residuals are normally distributed

H1: Residuals are not normally distributed

shap <- shapiro.test(Residuals)
par(mfrow = c(1,2))
hist(Residuals, main = "Histogram", probability = TRUE, ylim=c(0,0.03), lwd=1.5)
curve(dnorm(x, mean(Residuals),sd(Residuals)), col="red", add=TRUE,lwd=2)
legend("topright", legend = c("Bellshaped curve", "p-value = 0.0004942"), col=c("red","gray"), lwd =2)
qqnorm(Residuals) 
qqline(Residuals, lty = 2, col = "blue",lwd=2)

4.4.3. Residual homoscedasticity checking

H0: Homoscedasticity (constant variance of residuals)

H1: Heteroscedasticity

By inspecting the residuals vs the fitted values, we see that the values present a structure. Meaning the non constant variance. After checking for homoscedasticity test of Breusch-Pagan, we notice that the p-value is less that 0.05, this lead to an absence of homoscédasticity. So this condition didn’t holds.

plot(Fitted, Residuals, xlab = "Fitted values", ylab = "Residuals")
abline(h = 0, lty = 2)
text(Fitted, Residuals, labels = rownames(Malaria), pos = 4)

bptest(Model1)
## 
##  studentized Breusch-Pagan test
## 
## data:  Model1
## BP = 15.988, df = 7, p-value = 0.02522

4.4.4. Residual autocorrelation

H0: Autocorrelation = 0 (independent residuals)

H1: Positive autocorrelation (non-independent residuals)

dwtest(Model1)
## 
##  Durbin-Watson test
## 
## data:  Model1
## DW = 1.9127, p-value = 0.2688
## alternative hypothesis: true autocorrelation is greater than 0

4.4.5. Test for zero mean residuals

H0: Mean of residuals = 0

H1: Mean of residuals ≠ 0

t.test(Residuals, mu = 0)
## 
##  One Sample t-test
## 
## data:  Residuals
## t = -5.748e-16, df = 199, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -2.393562  2.393562
## sample estimates:
##     mean of x 
## -6.976936e-16

5. Linearity assessment (Ramsey RESET test)

H0: Linearity (the model is correctly specified)

H1: Non-linearity (model misspecification)

resettest(Model1, power = 2:3, type = "regressor")
## 
##  RESET test
## 
## data:  Model1
## RESET = 0.97411, df1 = 14, df2 = 178, p-value = 0.4817

Analysis of extreme values

In this section, we evaluated the extreme values and their influences in the regression line. The following plot represents the influence plot. Which combines Leverage, studentized residuals, and Cook’s distance (the size of the circles). As we can see in this figure, the observations 171, 92, 42 and 50 have the biggest values for Leverage and Cook’s distances. But only the 171 and 50 values are out of the acceptance region for studentized residuals. The critical value of the Leverage was calculated using the following formula :

Leverage critique = 2(p+1)​ /n

Where :

  • p = 7 the number of explanatory variables

  • n = 200 the sample size

  • \(\approx\) 200

influencePlot(Model1,
              id=TRUE,
              main="Influence Plot", fill.col = carPalette()[2], fill.alpha = 0.5)
##       StudRes        Hat      CookD
## 42  -1.085378 0.09227491 0.01495542
## 50  -4.445311 0.04252868 0.09994977
## 92   1.450012 0.10060282 0.02922974
## 171  3.292099 0.08225477 0.11550298
abline(v=0.08, col="red", lwd=2, lty=2)
abline(h=-2, col="blue", lwd=2, lty=2)
abline(h=2, col="blue", lwd=2, lty=2)

Variables transformation

bcmodel <- boxCox(Model1,lambda = seq(-2,2,1/10), plotit=TRUE)

The Box–Cox profile log-likelihood indicated an optimal transformation parameter close to \(\lambda\) = 1, with the 95% confidence interval including 1. Therefore, no transformation of the response variable was considered necessary, and the original scale of the malaria incidence was retained for subsequent analyses.

For independant variables

The overall Box-Tidwell test indicated no evidence against the linearity assumption (χ² = 5.07, df = 7, p = 0.652). At the individual level, only the variable distance showed evidence of departure from linearity according to the Score test (p = 0.040), whereas the remaining predictors satisfied the linearity assumption.

btmodel <- bocotir::boxTidwell(Model1)
btmodel$boxtid$main$Overall
## Chi2(df=7) Pr(> Chi2) 
##  5.0656723  0.6519491
btmodel$boxtid$main$Score
##                        Estimate    z-value  Pr(> |z|)
## Rainfall              0.2126226 -0.1012924 0.91931831
## temperature           1.5423375  0.1158761 0.90775073
## Humidity             -7.9867857 -0.3869323 0.69880629
## Stagnant_water_index -0.7696041 -0.2715343 0.78598016
## Bednet                0.2039661 -0.5575119 0.57717773
## distance              2.2849314  2.0533993 0.04003387
## Poverty               0.7482078 -0.1994503 0.84191052

Correction of the model taking into account the power tranrform

Model3 <- lm(incidence_rate~Rainfall+temperature+Humidity+Stagnant_water_index+Bednet+I(distance^2)+Poverty, data=Malaria)
summary(Model3)
## 
## Call:
## lm(formula = incidence_rate ~ Rainfall + temperature + Humidity + 
##     Stagnant_water_index + Bednet + I(distance^2) + Poverty, 
##     data = Malaria)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -71.001  -9.018   1.085   9.374  49.058 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -42.430911  19.351432  -2.193   0.0295 *  
## Rainfall               0.111971   0.058626   1.910   0.0576 .  
## temperature            3.100998   0.487437   6.362 1.43e-09 ***
## Humidity              -0.208484   0.265145  -0.786   0.4327    
## Stagnant_water_index   3.054151   1.109983   2.752   0.0065 ** 
## Bednet                -0.242894   0.056328  -4.312 2.58e-05 ***
## I(distance^2)          0.035623   0.007239   4.921 1.85e-06 ***
## Poverty                0.370814   0.071179   5.210 4.86e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.33 on 192 degrees of freedom
## Multiple R-squared:  0.4743, Adjusted R-squared:  0.4551 
## F-statistic: 24.75 on 7 and 192 DF,  p-value: < 2.2e-16

Prédictions tests

New data

new_data <- read_excel("C:/Users/hp/Downloads/new_data.xlsx")
# 2. Vérifier les noms des variables
names(new_data)
##  [1] "Year"                   "Health_zone_ID"         "Department"            
##  [4] "Health_zone"            "Rainfall_mm"            "Average_temp_C"        
##  [7] "Humidity_%"             "Stagnant_water_index"   "Bednet_coverage_%"     
## [10] "Healthcare_distance_km" "Poverty_rate_%"         "Malaria_incidence_rate"
## [13] "Data_status"            "Scientific_use_note"
new_data <- new_data %>% 
  rename(
    "Rainfall" = "Rainfall_mm",
    "temperature" = "Average_temp_C",
    "Humidity"="Humidity_%",
    "Bednet"="Bednet_coverage_%",
    "distance" = "Healthcare_distance_km",
    "Poverty" = "Poverty_rate_%",
    "incidence_rate" = "Malaria_incidence_rate"
  ) %>% 
  set_variable_labels(
    Rainfall = "Total rainfall in the district (mm)",
    temperature = "Average monthly temperature (°C)",
    Humidity = "Average relative humidity (%)",
    Stagnant_water_index = "Index of stagnant water surface area (scale from 0 to 10)",
    Bednet = "Percentage of the population regularly using insecticide-treated bed nets (%)",
    distance = "Average distance to the nearest health center (km)",
    Poverty = "Percentage of the population living below the poverty line (%)",
    incidence_rate = "Number of confirmed malaria cases per 1,000 inhabitants"
  )

prediction_data <- new_data[5:12,]

Prédictions

Here it’s clear that we don’t have enougth eveidence for predicting Malaria incidence. This model is not very correct event if it have R²a=0.4551. There remain others factor non mesured.

pred <- predict(
  Model3,
  newdata = prediction_data,
  interval = "prediction",
  level = 0.95
)


prediction_results <- prediction_data %>%
  mutate(
    predicted_incidence = pred[, "fit"],
    lower_95 = pred[, "lwr"],
    upper_95 = pred[, "upr"],
    residual_prediction = incidence_rate - predicted_incidence
  )

head(prediction_results[,c(1,4,15:17)],8)
## # A tibble: 8 × 5
##    Year Health_zone                  predicted_incidence lower_95 upper_95
##   <dbl> <chr>                                      <dbl>    <dbl>    <dbl>
## 1  2023 Kouande-Pehunco-Kerou                       168.     67.4     269.
## 2  2023 Tanguiéta-Materi-Cobly                      136.     47.5     225.
## 3  2023 Abomey-Calavi-So-Ava                        189.     52.2     325.
## 4  2023 Allada-Toffo-Ze                             197.     59.3     334.
## 5  2023 Ouidah-Kpomasse-Tori-Bossito                193.     47.5     339.
## 6  2023 Bembereke-Sinende                           156.     62.8     249.
## 7  2023 Nikki-Kalale-Perere                         164.     65.0     264.
## 8  2023 Parakou-N'Dali                              140.     55.8     225.