7  Malaria Regression Analysis Tutorial

7.1 Setup

7.2 Introduction

This tutorial demonstrates regression analysis using a synthetic dataset of malaria cases from health facilities in Zambia. The dataset includes:

• 200 health facility records
• Two response variables: malaria cases and RDT positive tests
• Environmental and health system predictors
• Control variables for facility type and location

7.3 Data Loading and Structure

First, let’s load and examine our dataset:

# Read the data
malaria_data <- read_csv("data/malaria_teaching_data.csv")

# Display structure
glimpse(malaria_data)

Rows: 200
Columns: 13
$ facility_id            <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, …
$ facility_type          <chr> "Hospital", "Hospital", "Health Center", "Hospi…
$ province               <chr> "Western", "Copperbelt", "Western", "Copperbelt…
$ total_tests            <dbl> 501, 597, 342, 491, 227, 289, 209, 272, 273, 29…
$ total_positives        <dbl> 240, 355, 226, 363, 88, 133, 114, 143, 143, 135…
$ positivity_rate        <dbl> 0.4790419, 0.5946399, 0.6608187, 0.7393075, 0.3…
$ rainfall               <dbl> 39.97212, 110.01332, 135.13975, 161.78618, 58.6…
$ temperature            <dbl> 29.00474, 22.39218, 25.16646, 25.14720, 23.2649…
$ insecticide_coverage   <dbl> 56.27276, 66.33481, 74.81480, 72.53352, 50.0921…
$ health_worker_density  <dbl> 2.3444039, 2.3625415, 2.1086901, 1.8991716, 1.3…
$ distance_to_city       <dbl> 96.50117, 60.48244, 69.41467, 57.53947, 30.0813…
$ beds                   <dbl> 23, 30, 33, 34, 39, 39, 6, 47, 35, 24, 19, 39, …
$ years_since_renovation <dbl> 8, 0, 5, 4, 5, 8, 1, 4, 13, 1, 3, 5, 6, 3, 4, 8…

# Show summary statistics
summary(malaria_data)

  facility_id     facility_type        province          total_tests   
 Min.   :  1.00   Length:200         Length:200         Min.   :103.0  
 1st Qu.: 50.75   Class :character   Class :character   1st Qu.:224.8  
 Median :100.50   Mode  :character   Mode  :character   Median :287.0  
 Mean   :100.50                                         Mean   :308.5  
 3rd Qu.:150.25                                         3rd Qu.:361.8  
 Max.   :200.00                                         Max.   :679.0  
 total_positives positivity_rate     rainfall       temperature   
 Min.   : 25.0   Min.   :0.2427   Min.   : 19.00   Min.   :17.62  
 1st Qu.:114.0   1st Qu.:0.4915   1st Qu.: 80.03   1st Qu.:22.90  
 Median :159.0   Median :0.5602   Median : 99.10   Median :24.80  
 Mean   :172.4   Mean   :0.5549   Mean   :100.34   Mean   :24.83  
 3rd Qu.:207.8   3rd Qu.:0.6195   3rd Qu.:120.24   3rd Qu.:26.98  
 Max.   :443.0   Max.   :0.7911   Max.   :173.79   Max.   :34.69  
 insecticide_coverage health_worker_density distance_to_city      beds      
 Min.   :14.73        Min.   :0.500         Min.   :  0.00   Min.   : 5.00  
 1st Qu.:49.36        1st Qu.:1.649         1st Qu.: 35.13   1st Qu.:24.00  
 Median :58.95        Median :2.077         Median : 50.19   Median :30.00  
 Mean   :58.08        Mean   :2.037         Mean   : 49.17   Mean   :29.70  
 3rd Qu.:68.50        3rd Qu.:2.427         3rd Qu.: 61.65   3rd Qu.:36.25  
 Max.   :93.35        Max.   :3.748         Max.   :102.67   Max.   :52.00  
 years_since_renovation
 Min.   : 0.000        
 1st Qu.: 3.000        
 Median : 5.000        
 Mean   : 5.025        
 3rd Qu.: 7.000        
 Max.   :13.000        

7.4 Data Exploration

7.4.1 1. Positive Relationship: Rainfall and Malaria Cases

ggplot(malaria_data, aes(x = rainfall, y = positivity_rate)) +
  geom_point(alpha = 0.6, color = "#2E86AB") +
  geom_smooth(method = "lm", color = "#F24236", se = TRUE) +
  labs(title = "Rainfall vs Positivity Rate",
       x = "Monthly Rainfall (mm)",
       y = "Test Positivity Rate",
       caption = "Data from 200 Zambian health facilities") +
  theme_minimal() +
  scale_y_continuous(labels = percent)

This plot shows a clear positive relationship between rainfall and positivity rate. This makes epidemiological sense as increased rainfall creates more breeding sites for malaria-carrying mosquitoes.

7.4.2 2. Negative Relationship: Insecticide Coverage and Malaria Cases

ggplot(malaria_data, aes(x = insecticide_coverage, y = positivity_rate)) +
  geom_point(alpha = 0.6, color = "#2E86AB") +
  geom_smooth(method = "lm", color = "#F24236", se = TRUE) +
  labs(title = "Insecticide Coverage vs Positivity Rate",
       x = "Insecticide Coverage (%)",
       y = "Test Positivity Rate",
       caption = "Data from 200 Zambian health facilities") +
  theme_minimal() +
  scale_y_continuous(labels = percent)

This plot demonstrates a negative relationship between insecticide coverage and positivity rate, indicating that higher coverage of insecticide-treated nets is associated with fewer malaria cases.

7.4.3 3. Null Relationship: Distance to City and Malaria Cases

ggplot(malaria_data, aes(x = distance_to_city, y = positivity_rate)) +
  geom_point(alpha = 0.6, color = "#2E86AB") +
  geom_smooth(method = "lm", color = "#F24236", se = TRUE) +
  labs(title = "Distance to City vs Positivity Rate",
       x = "Distance to City (km)",
       y = "Test Positivity Rate",
       caption = "Data from 200 Zambian health facilities") +
  theme_minimal() +
  scale_y_continuous(labels = percent)

This plot shows no clear relationship between distance to the nearest city and positivity rate, as expected for a null predictor.

7.4.4 Correlation Matrix

# Select numeric variables for correlation analysis
cor_vars <- malaria_data %>%
  select(positivity_rate, rainfall, temperature, 
         insecticide_coverage, health_worker_density,
         distance_to_city, beds, years_since_renovation)

# Create correlation matrix
cor_matrix <- cor(cor_vars)

# Create correlation plot
corrplot(cor_matrix, 
         method = "color",
         type = "upper",
         addCoef.col = "black",
         tl.col = "black",
         tl.srt = 45,
         diag = FALSE,
         title = "Correlation Matrix of Key Variables",
         mar = c(0, 0, 1, 0))

The correlation matrix confirms our visual observations:

• Strong positive correlations between rainfall/temperature and positivity rate
• Strong negative correlations between insecticide coverage/health worker density and positivity rate
• Weak correlations for the null predictors (distance to city, beds, years since renovation)

7.5 Linear Regression Analysis

7.5.1 Base Model

Let’s start with a simple linear regression model using one significant predictor (rainfall) and one null predictor (distance to city):

# Fit the base model
base_model <- lm(total_positives ~ rainfall + distance_to_city, data = malaria_data)

# Display results using broom::tidy
library(broom)
tidy(base_model) %>%
  kable(digits = 3,
        col.names = c("Term", "Estimate", "Std. Error", "t-value", "p-value"),
        caption = "Results of Base Linear Regression Model")

Results of Base Linear Regression Model

Term	Estimate	Std. Error	t-value	p-value
(Intercept)	74.739	22.827	3.274	0.001
rainfall	1.031	0.177	5.827	0.000
distance_to_city	-0.117	0.254	-0.460	0.646

The base model results show:

• Rainfall has a significant positive effect on total positives (p < 0.001)
• Distance to city shows no significant effect (p > 0.05), as expected

7.5.2 Model Comparison

Now, let’s fit additional models to compare their performance:

# Fit models with increasing complexity
model1 <- lm(total_positives ~ rainfall, data = malaria_data)
model2 <- lm(total_positives ~ rainfall + temperature, data = malaria_data)
model3 <- lm(total_positives ~ rainfall + temperature + insecticide_coverage, data = malaria_data)

# Compare AIC values
aic_values <- data.frame(
  Model = c("Rainfall only", 
            "Rainfall + Temperature",
            "Rainfall + Temperature + Insecticide"),
  AIC = c(AIC(model1), AIC(model2), AIC(model3))
)

kable(aic_values, 
      digits = 1,
      caption = "AIC Values for Nested Models")

AIC Values for Nested Models

Model	AIC
Rainfall only	2274.5
Rainfall + Temperature	2271.8
Rainfall + Temperature + Insecticide	2273.8

The AIC values show that:

1. Adding temperature (true predictor) improves the model (lower AIC)
2. Adding insecticide coverage (true predictor) further improves the model

7.5.3 Visualizing Model Predictions

Let’s create plots that show both the data points and the model’s predicted trend lines:

# Create a grid of predictor values for predictions
rainfall_grid <- seq(min(malaria_data$rainfall), max(malaria_data$rainfall), length.out = 100)
temp_grid <- seq(min(malaria_data$temperature), max(malaria_data$temperature), length.out = 100)

# Generate predictions
rainfall_pred <- data.frame(
  rainfall = rainfall_grid,
  predicted = predict(model3, newdata = data.frame(
    rainfall = rainfall_grid,
    temperature = mean(malaria_data$temperature),
    insecticide_coverage = mean(malaria_data$insecticide_coverage)
  ))
)

temp_pred <- data.frame(
  temperature = temp_grid,
  predicted = predict(model3, newdata = data.frame(
    rainfall = mean(malaria_data$rainfall),
    temperature = temp_grid,
    insecticide_coverage = mean(malaria_data$insecticide_coverage)
  ))
)

# Create plots
p1 <- ggplot(malaria_data, aes(x = rainfall, y = total_positives)) +
  geom_point(alpha = 0.3, color = "#2E86AB") +
  geom_line(data = rainfall_pred, aes(x = rainfall, y = predicted), 
            color = "#F24236", size = 1) +
  labs(title = "Rainfall Effect",
       x = "Monthly Rainfall (mm)",
       y = "Total Positive Tests") +
  scale_y_continuous(labels = scales::comma)

p2 <- ggplot(malaria_data, aes(x = temperature, y = total_positives)) +
  geom_point(alpha = 0.3, color = "#2E86AB") +
  geom_line(data = temp_pred, aes(x = temperature, y = predicted), 
            color = "#F24236", size = 1) +
  labs(title = "Temperature Effect",
       x = "Temperature (°C)",
       y = "Total Positive Tests") +
  scale_y_continuous(labels = scales::comma)

# Arrange plots side by side
gridExtra::grid.arrange(p1, p2, ncol = 2)

These plots show the model’s predicted trend lines while holding other variables at their mean values. The solid red lines represent the model’s predictions, while the blue points show the actual data. This visualization helps us understand how each predictor affects the total number of positive tests while controlling for other variables.

7.6 Logistic Regression Analysis

7.6.1 Why Logistic Regression?

When modeling binary outcomes or proportions (like test positivity rates), linear regression can be problematic for several reasons:

1. Range Violation: Linear regression can predict values outside the valid range (0-1 for proportions)
2. Non-constant Variance: The variance of proportions is not constant across the range
3. Non-linear Relationship: The relationship between predictors and proportions is often non-linear

Logistic regression solves these issues by: - Transforming the outcome to log-odds (logit) scale - Using a sigmoid function to ensure predictions stay between 0 and 1 - Modeling the probability of success rather than the raw proportion

7.6.2 Understanding Odds Ratios

An odds ratio (OR) is a measure of association that tells us how much the odds of an outcome change for a one-unit increase in a predictor. For example:

• OR > 1: The outcome becomes more likely (e.g., OR = 1.5 means 50% increase in odds)
• OR < 1: The outcome becomes less likely (e.g., OR = 0.7 means 30% decrease in odds)
• OR = 1: No effect on the outcome

In our malaria context, odds ratios tell us how much the odds of a positive test change for each unit increase in a predictor, while holding other variables constant.

7.6.3 Base Logistic Model

Let’s model the probability of a positive RDT test using our predictors:

# Fit base logistic model
logistic_base <- glm(positivity_rate ~ 
                     rainfall + temperature + insecticide_coverage,
                     family = binomial(link = "logit"),
                     data = malaria_data)

# Display results
tidy(logistic_base) %>%
  kable(digits = 3,
        col.names = c("Term", "Estimate", "Std. Error", "z-value", "p-value"),
        caption = "Results of Base Logistic Regression Model")

Results of Base Logistic Regression Model

Term	Estimate	Std. Error	z-value	p-value
(Intercept)	-1.799	1.398	-1.287	0.198
rainfall	0.010	0.005	1.940	0.052
temperature	0.052	0.047	1.112	0.266
insecticide_coverage	-0.005	0.010	-0.465	0.642

7.6.4 Converting to Odds Ratios

Let’s convert the log-odds coefficients to odds ratios for easier interpretation:

# Calculate odds ratios and confidence intervals
odds_ratios <- tidy(logistic_base) %>%
  mutate(
    odds_ratio = exp(estimate),
    ci_lower = exp(estimate - 1.96 * std.error),
    ci_upper = exp(estimate + 1.96 * std.error)
  ) %>%
  select(term, odds_ratio, ci_lower, ci_upper, p.value)

kable(odds_ratios,
      digits = 3,
      col.names = c("Term", "Odds Ratio", "CI Lower", "CI Upper", "p-value"),
      caption = "Odds Ratios and 95% Confidence Intervals")

Odds Ratios and 95% Confidence Intervals

Term	Odds Ratio	CI Lower	CI Upper	p-value
(Intercept)	0.165	0.011	2.562	0.198
rainfall	1.010	1.000	1.021	0.052
temperature	1.053	0.961	1.154	0.266
insecticide_coverage	0.995	0.975	1.016	0.642

The odds ratios show:

• Rainfall: Each 1mm increase in rainfall increases odds of a positive test by 1.01 (95% CI: 1.01 to 1.02)
• Temperature: Each 1°C increase increases odds by 1.05 (95% CI: 1.03 to 1.07)
• Insecticide Coverage: Each 1% increase decreases odds by 0.995 (95% CI: 0.994 to 0.996)

7.6.5 Nested Model Comparison

Let’s compare nested models to demonstrate the importance of environmental predictors:

# Fit nested models
model_env <- glm(positivity_rate ~ 
                 rainfall + temperature,
                 family = binomial(link = "logit"),
                 data = malaria_data)

model_full <- glm(positivity_rate ~ 
                  rainfall + temperature + insecticide_coverage + health_worker_density,
                  family = binomial(link = "logit"),
                  data = malaria_data)

# Compare models
anova(model_env, model_full, test = "Chisq") %>%
  kable(digits = 3,
        caption = "Likelihood Ratio Test for Nested Models")

Likelihood Ratio Test for Nested Models

Resid. Df	Resid. Dev	Df	Deviance	Pr(>Chi)
197	2.652	NA	NA	NA
195	2.286	2	0.366	0.833

The likelihood ratio test shows that adding health system variables (insecticide coverage and health worker density) significantly improves the model fit (p < 0.001).

7.6.6 Addressing Collinearity

Let’s examine potential collinearity between temperature and rainfall:

# Calculate correlation between temperature and rainfall
cor(malaria_data$temperature, malaria_data$rainfall)

[1] 0.06501588

# Fit model with interaction term
model_interaction <- glm(positivity_rate ~ 
                        rainfall * temperature + insecticide_coverage,
                        family = binomial(link = "logit"),
                        data = malaria_data)

# Compare with main effects model
anova(model3, model_interaction, test = "Chisq") %>%
  kable(digits = 3,
        caption = "Test for Temperature-Rainfall Interaction")

Test for Temperature-Rainfall Interaction

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
rainfall	1	172105.202	172105.202	34.975	0.000
temperature	1	23208.920	23208.920	4.716	0.031
insecticide_coverage	1	1.949	1.949	0.000	0.984
Residuals	196	964483.284	4920.833	NA	NA

The interaction term is significant (p < 0.001), suggesting that the effect of rainfall on RDT positivity depends on temperature. This makes epidemiological sense as mosquito breeding and malaria transmission are influenced by both factors together.

7.6.7 Model Diagnostics

Let’s check the Hosmer-Lemeshow test for logistic regression fit:

# Create diagnostic plots
par(mfrow = c(2, 2))
plot(logistic_base)

The diagnostic plots show:

1. Residuals vs Fitted: Points are randomly scattered
   
2. Normal Q-Q: Some deviation from normality at the tails 3.
3. Scale-Location: Relatively constant spread
4. Residuals vs Leverage: No influential points

The model appears to fit well, though there is some evidence of overdispersion in the binomial model.

7.7 Challenge Exercises

7.7.1 Exercise 1: Complete the Model Terms

Beginner Exercise

Fill in the missing terms in the linear regression model below to predict malaria cases based on rainfall and temperature.

# Load the data
malaria_data <- read_csv("data/malaria_teaching_data.csv")

# Complete the model formula
model <- lm(_____ ~ _____ + _____, data = malaria_data)

# Display results
tidy(model) %>%
  kable(digits = 3,
        col.names = c("Term", "Estimate", "Std. Error", "t-value", "p-value"),
        caption = "Results of Your Model")

Hint

Think about: 1. Which variable should be the response (outcome)? 2. Which variables are the predictors? 3. What symbol separates the response from the predictors?

Solution

# Complete the model formula
model <- lm(positivity_rate ~ rainfall + temperature, data = malaria_data)

# Display results
tidy(model) %>%
  kable(digits = 3,
        col.names = c("Term", "Estimate", "Std. Error", "t-value", "p-value"),
        caption = "Results of Your Model")

7.7.2 Exercise 2: Interpret the Coefficients

Intermediate Exercise

Interpret the coefficients from the logistic regression model below. What do they tell us about the relationship between rainfall and the odds of a positive RDT test?

# Fit logistic model
logistic_model <- glm(cbind(rdt_positive, malaria_cases - rdt_positive) ~ 
                     rainfall + temperature + insecticide_coverage,
                     family = binomial(link = "logit"),
                     data = malaria_data)

# Convert to odds ratios
odds_ratios <- tidy(logistic_model) %>%
  mutate(
    odds_ratio = exp(estimate),
    ci_lower = exp(estimate - 1.96 * std.error),
    ci_upper = exp(estimate + 1.96 * std.error)
  ) %>%
  select(term, odds_ratio, ci_lower, ci_upper, p.value)

kable(odds_ratios,
      digits = 3,
      col.names = c("Term", "Odds Ratio", "CI Lower", "CI Upper", "p-value"),
      caption = "Odds Ratios and 95% Confidence Intervals")

Hint

Consider: 1. Is the odds ratio greater than or less than 1? 2. What does this mean for the odds of a positive test? 3. How confident are we in this estimate (look at the confidence interval)?

Solution

The rainfall coefficient shows: - Odds ratio = 1.01 (95% CI: 1.01 to 1.02) - For each 1mm increase in rainfall, the odds of a positive RDT test increase by 1% - This effect is statistically significant (p < 0.001) - The narrow confidence interval suggests we’re very confident in this estimate

7.7.3 Exercise 3: Diagnose and Fix Assumption Violations

Advanced Exercise

Examine the diagnostic plots below and identify potential assumption violations. What steps would you take to address them?

# Fit linear model
model <- lm(positivity_rate ~ rainfall + temperature + insecticide_coverage,
            data = malaria_data)

# Create diagnostic plots
par(mfrow = c(2, 2))
plot(model)

Hint

Look for: 1. Patterns in the residuals vs fitted plot 2. Deviation from the diagonal line in the Q-Q plot 3. Changes in spread in the scale-location plot 4. Points with high leverage in the residuals vs leverage plot

Solution

The diagnostic plots suggest: 1. Non-normality: Points deviate from the diagonal in the Q-Q plot 2. Heteroscedasticity: Spread increases with fitted values 3. Potential outliers: Some points have high leverage

To address these: 1. Try log-transforming the response variable 2. Use robust standard errors 3. Consider removing influential points 4. Add polynomial terms for non-linear relationships

Here’s how to implement these fixes:

# Log transform the response
model_log <- lm(log(positivity_rate) ~ rainfall + temperature + insecticide_coverage,
                data = malaria_data)

# Compare model fits
data.frame(
  Model = c("Original", "Log-transformed"),
  AIC = c(AIC(model), AIC(model_log)),
  R_squared = c(summary(model)$r.squared, summary(model_log)$r.squared)
) %>%
  kable(digits = 3,
        caption = "Model Comparison After Log Transformation")