Modeling Notebook

Introduction

The dataset constructed in previous scripts is a panel-data set and is the type of data which would typically analyzed in a social science. As such we will use the linear models typical of those fields, however for this project we will also be using modeling techniques typically found in predictive analytics, i.e. crossvalidation, with our primary focus being prediction. This notebook will be focused on these models, see previous notebooks and scripts for data-sourcing, wrangling, and EDA.

Loading Data and Packages

library(zoo)
library(MASS)
library(tidyverse)
library(lme4)
library(caTools)
library(caret)
final_data <- read.csv("final_data.csv")

Setting Up the Data for Modeling

First we need to handle some lingering issues with the dataset, i.e. missing values so we can approprately model our data. First we interpolate missing values for GDP_Per_Capita by country. This should only leave us with countries that have no infromation on GDP_Per_capita, which we will omit.

model_data <- final_data %>% 
  group_by(Entity) %>% 
  mutate(GDP_Per_Capita = na.approx(GDP_Per_Capita, maxgap = 25, rule = 2))

model_data <- model_data %>% 
  na.omit()

First Model

We’ll begin with a basic two-way fixed effects model with year and country as the fixed effects. We can evaluate how well this model does by examining our RMSE and residual plots. Here we run into one of our first obstacles in the modeling process, since our data is panel-data when we segment our data with validation techniques, K-fold CV, data segmentation, etc. we run into issues with stability in our regression. There is however a technique called LOOCV which will be helpful for the modeling here. LOOCV fits a model with the ith observation removed and then calculates the residual at that data point. This is alternatively known as the PRESS residual. I write a function which returns a dataframe with the actual value, predicted value, residual value, standardized residual value, PRESS residual value, and studentized residual value. We will save this data frame as well as one containing the PRESS statistic (PRESS RMSE), and train test error for comparison later on.

The following plots depict these residuals in multiple ways. The first plots our actual values against our predicted values. The interpretation of this plot is that the closer the points match the diagonal line the more appropriately the model fits the data. The further the points stray from the line we may find issues such as heteroskedasticity, or lack of fit. The next four depict the residual vs predicted plots. These plots ideally should resemble a random cloud. Additionally, for studentized and standardized residuals we may depict outliers, those points which lie beyond 4.

formula <- (cardiovasc_deaths ~ mean_bmi +percent_obesity + percent_severe_obesity  + GDP_Per_Capita +factor(Entity) + factor(Year))

model1  <- lm(formula, data=model_data)

PRESS <- function(linear.model) {
  #' calculate the predictive residuals
  pr <- residuals(linear.model)/(1-lm.influence(linear.model)$hat)
  #' calculate the PRESS
  PRESS <- sum(pr^2)
  
  return(pr)
}

model_eval<- function(model, data){

  stanresid<-model %>% 
    rstandard()
  
  yhat <-model %>% 
    predict()
  
  resid <- model$residuals
  
  pressres <- model %>% 
    PRESS()
  
  stures <-  model %>% 
  studres()
  
  stanres <- model %>% 
    rstandard()
  fulldf<-cbind(data$cardiovasc_deaths, yhat, resid, stanres, pressres, stures)
  fulldf = as.data.frame(fulldf)
  names(fulldf) <- c("Y","Y_Hat", "residuals", "standardized_residuals", "press_residuals", "studentized_residuals")
  return(fulldf)
}

residuals_model1<- model_eval(model1, model_data)

errors<- data.frame("Full_Test" = numeric(),"LOOCV" = numeric(),"Segment_Test" = numeric(),"Segment_Train" = numeric())
errors[1,1] <- sqrt(sum((residuals_model1$Y - residuals_model1$Y_Hat)^2)/length(residuals_model1[,1]))
up <-sum(residuals_model1$press_residuals^2)/length(residuals_model1[,1])   
errors[1,2]<-sqrt(up)

#Split data into train and test sets to approximate error for new obs.
set.seed(789)
training_samples <- model_data$cardiovasc_deaths %>% 
    createDataPartition(p= .9, list= FALSE)

training_samples <- training_samples %>% 
  as.vector()

train_data <- model_data[training_samples, ]
test_data <- model_data[-training_samples, ]

train_model  <- lm(formula, data=train_data)

predictions <- train_model %>% 
  predict(test_data)

errors[1,3] <- RMSE(predictions, test_data$cardiovasc_deaths)

predictions <- train_model %>% 
  predict(train_data)

errors[1,4] <- RMSE(predictions, train_data$cardiovasc_deaths)

ggplot(residuals_model1, aes( Y, Y_Hat))+
       labs(x='Actual', y ='Predicted')+
       geom_point()+
      geom_abline(intercept =0, slope =  1, color='red', size =2)

ggplot(residuals_model1, aes( Y_Hat, residuals))+
  labs(x='Predicted', y ='Residual')+
  geom_point()+
  geom_abline(intercept =0, slope =  0, color='blue', size =1)

ggplot(residuals_model1, aes( Y_Hat, standardized_residuals))+
  labs(x='Predicted', y ='Standardized Residual')+
  geom_point()+
  geom_abline(intercept =0, slope =  0, color='blue', size =1)

ggplot(residuals_model1, aes( Y_Hat, press_residuals))+
  labs(x='Predicted', y ='PRESS Residual')+
  geom_point()+
  geom_abline(intercept =0, slope =  0, color='blue', size =1)

ggplot(residuals_model1, aes( Y_Hat, studentized_residuals))+
  labs(x='Predicted', y ='Studentized Residual')+
  geom_point()+
  geom_abline(intercept =0, slope =  0, color='blue', size =1)

Second Model

As we see our model appears to fit well. There does appear to be a degree of heteroskedasticity in the higher ranges of our data, this can be seen in all four plots. We will try to overcome this by introducing a transformation for our response variable.

formula <- (log(cardiovasc_deaths) ~ mean_bmi +percent_obesity + percent_severe_obesity  + GDP_Per_Capita +factor(Entity) + factor(Year))

model2  <- lm(formula, data=model_data)


residuals_model2 <- model_eval(model2, model_data)

residuals_model2$Y <- log(residuals_model2$Y)

We can do the same analysis on our log-transformed response model. The plots for these are provided below as well as the previously defined function which saves our residuals, both LOOCV, and train-test split based errors are saved.

formula <- (log(cardiovasc_deaths) ~ mean_bmi +percent_obesity + percent_severe_obesity  + GDP_Per_Capita +factor(Entity) + factor(Year))

model2  <- lm(formula, data=model_data)


residuals_model2 <- model_eval(model2, model_data)

residuals_model2$Y <-  log(residuals_model2$Y)

errors[2,1] <- RMSE(exp(residuals_model2$Y), exp(residuals_model2$Y_Hat))

Yresid <- exp(residuals_model2$Y) - exp(residuals_model2$Y_Hat)
pr <- Yresid/(1-lm.influence(model2)$hat)
errors[2,2] <- sqrt(sum(pr^2)/length(pr))

train_model  <- lm(formula, data=train_data)

predictions <- train_model %>% 
  predict(test_data) %>% 
  exp()

errors[2,3] <- RMSE(predictions, test_data$cardiovasc_deaths)

predictions <- train_model %>% 
  predict(train_data) %>% 
  exp()
errors[2,4] <- RMSE(predictions, train_data$cardiovasc_deaths)

#Same plots as before
ggplot(residuals_model2, aes( Y, Y_Hat))+
  labs(x='Actual', y ='Predicted')+
  geom_point()+
  geom_abline(intercept =0, slope =  1, color='red', size =2)

ggplot(residuals_model2, aes( Y_Hat, residuals))+
  labs(x='Predicted', y ='Residual')+
  geom_point()+
  geom_abline(intercept =0, slope =  0, color='blue', size =1)

ggplot(residuals_model2, aes( Y_Hat, standardized_residuals))+
  labs(x='Predicted', y ='Standardized Residual')+
  geom_point()+
  geom_abline(intercept =0, slope =  0, color='blue', size =1)

ggplot(residuals_model2, aes( Y_Hat, press_residuals))+
  labs(x='Predicted', y ='PRESS Residual')+
  geom_point()+
  geom_abline(intercept =0, slope =  0, color='blue', size =1)

ggplot(residuals_model2, aes( Y_Hat, studentized_residuals))+
  labs(x='Predicted', y ='Studentized Residual')+
  geom_point()+
  geom_abline(intercept =0, slope =  0, color='blue', size =1)

Final Model

In the above plots we see our second model does appear to fit the data better. However, there are a couple of serious problems with this form of our model. First, looking at our cross-validation methods we see that the error for this term is higher in all metrics than the previous model. Additionally, there is difficulty in the interpretation of our coefficients when we introduce transformations into the model.

With this in mind let’s try a final different model. Our two previous models have both been fixed effects, we may be able to improve our model by introducing what are known as random effects by country for mean_bmi. These random effects allow the slope, for mean_bmi, and the intercept to vary depending on the entity. LOOCV are not as easy to derive computationally for this type of model, so we will have to omit them in our model evaluation. The evaluation plots for this model are provided below.

library(lme4)
formula <- (cardiovasc_deaths ~ percent_obesity + percent_severe_obesity + GDP_Per_Capita + 1+(1 + mean_bmi|Entity) + factor(Year))
model_3 <- lmer(formula, data=model_data, REML=FALSE)

Some predictor variables are on very different scales: consider rescalingModel failed to converge with max|grad| = 0.00319404 (tol = 0.002, component 1)

predictions <- model_3 %>% 
  predict(model_data)

errors[3,1]<- RMSE( model_data$cardiovasc_deaths, predictions)

train_model <- lmer(formula, data=train_data, REML=FALSE)

Some predictor variables are on very different scales: consider rescaling

predictions <- train_model %>% 
  predict(test_data)
errors[3,3]<- RMSE( test_data$cardiovasc_deaths, predictions)

predictions <- train_model %>% 
  predict(train_data)
errors[3,4]<- RMSE( train_data$cardiovasc_deaths, predictions)

predictions <- model_3 %>% 
  predict(model_data)
residuals_model3 <- as.data.frame(cbind(predictions, model_data$cardiovasc_deaths))
residuals_model3$residuals <- residuals_model3$V2 - residuals_model3$predictions 
names(residuals_model3)[2] <- "actual"
residuals_model3$Entity <- model_data$Entity

ggplot(residuals_model3, aes( actual, predictions))+
  labs(x='Actual', y ='Predicted')+
  geom_point()+
  geom_abline(intercept =0, slope =  1, color='red', size =2)

ggplot(residuals_model3, aes( predictions, residuals))+
  labs(x='Predicted', y ='Residual')+
  geom_point()+
  geom_abline(intercept =0, slope =  0, color='blue', size =1)

We see from the above plots that our new model fits the data well. From our residual plot however it does appear as though we have reintroduced heteroskedasticity. Generally heteroskedasticity is considered to be less of a severe violation of our linear regression assumptions.

errors[1, 5] <- "model_1"
errors[2, 5] <- "model_2"
errors[3, 5] <- "model_3"
names(errors)[5] <- "Model"
  
ggplot(errors, aes(x=Model, y=Segment_Test ))+
  labs( x= "Model", y= "Test RMSE") +
  geom_bar(stat="identity", width=.009)

Model Comparison and Conclusion

Plotted above is the test RMSE for all three models. We see from this that our third model clearly fits the data the best, additionally none of the models seem to seriously violate linear assumptions. Using this information, we can now predict the decrease in deaths which would result from a country decreasing their mean_bmi. This is an extremely useful insight for public health officials, and though the models presented here are relatively simple, should not be underestimated.