R Notebook

# Authors: Rahul Bamotra, Maximilian Kantsler, Ryan Ansett, Anson Abraham, Maxim Skliarov

# For proper execution of this code the following files need to be located in the Working Directory
#
# BMI Index.csv
# hapind2020.csv
# HIV-AIDS.csv
# Malaria.csv
# Tuberculosis.csv

# Packages used: countrycode; readxl; wpp2019; wbstats; dplyr; rvest; ggplot2; ggfortify; car; pedometrics; plotly

library(countrycode)
library(readxl)
library(wpp2019)
library(wbstats)
library(dplyr)
library(rvest)
library(ggplot2)
library(ggfortify)
library(car)
library(pedometrics)
library(plotly)


rm(list = ls())

# Scraping from Global Residence Website
# Scraping code taken from: https://beanumber.github.io/sds192/lab-import.html

src <- "https://globalresidenceindex.com/wp-content/uploads/2017/10/STC-Climate-18.xlsx"
lcl <- basename(src)
download.file(url = src, destfile = lcl,mode='wb')
weather <- as.data.frame(read_excel(lcl))
c.code<-countrycode(weather[,2],"country.name","iso3c")
weather<-cbind(weather,c.code)

src <- "https://globalresidenceindex.com/wp-content/uploads/2017/10/STC-Health-Index-2018.xlsx"
lcl <- basename(src)
download.file(url = src, destfile = lcl,mode='wb')
health<-as.data.frame(read_excel(lcl))
c.code<-countrycode(health[,2],"country.name","iso3c")
health<-cbind(health,c.code)

# Getting Life Expectancy from World Bank
inf <- wb_search("Life Expectancy")
life_exp<-wb_data(country="countries_only",indicator = c("SP.DYN.LE00.FE.IN","SP.DYN.LE00.IN","SP.DYN.LE00.MA.IN"),start_date = 2016,end_date = 2016)
life_exp<-na.exclude(life_exp)
colnames(life_exp)[2] <- "c.code"

happiness<-read.csv(file='hapind2020.csv')
hapdata<-happiness[c(1:3,7)]
c.code<-countrycode(hapdata[,1],"country.name","iso3c")
hapdata<-cbind(hapdata,c.code)

rm(c.code,inf,lcl,src,happiness)

data(pop) #Call for data of population
#remove data that is not required
rm(popF)
rm(popM)
rm(popFT)
rm(popMT)

#Get a column of country code data in the population data frame

pop$country.code <- countrycode(
  pop[,2],
  "country.name",
  "iso3c",
  
)

#Data cleanup - only get the columns required for the code
pop <- pop[, c(2,16,18)]
pop <- pop[complete.cases(pop),]
pop <- pop[-1,]
colnames(pop) <- c("country.name", "population","c.code")

variables <- weather[,c(16,4,6,10,11,13,15)]
cor(variables[-1])

# Check for correlations and patterns within the variables (same for every regression set)
plot(variables, pch = 15, col = "blue")

# Merging the variables with the response variables
data.lif <- right_join(life_exp[,c(2,6)], variables, by ="c.code")
data.lif <- data.lif[complete.cases(data.lif),]
attr(data.lif$SP.DYN.LE00.IN, "label") <- NULL # Removing attribute from column name
colnames(data.lif) <- c("c.code", "Response", "Index", "Days32", 
                        "Precipitation.Annual", "Precipitation.Days",
                        "Heating.Days", "Sunshine.Hours")
data.hap <- right_join(hapdata[,c(5,3)], variables, by ="c.code")
data.hap <- data.hap[complete.cases(data.hap),]
colnames(data.hap) <- c("c.code", "Response", "Index", "Days32", 
                        "Precipitation.Annual", "Precipitation.Days",
                        "Heating.Days", "Sunshine.Hours")
lif <- lm(Response ~ . -c.code, data = data.lif) # Problem with country codes
summary(lif)

# Check for relationships of variables to response variable (normalized value)
plot(data.lif, pch = 15, col = "blue")

# After checking for model accuracy we decided to include higher order terms of variables to improve the model
lif <- lm(Response ~ Precipitation.Days + poly(Heating.Days, 4) + poly(Index, 4) + poly(Precipitation.Annual, 4), data = data.lif)
summary(lif)

predict_data.lif <- data.frame(data.lif[,-c(2)])
predict_data.lif$predict_respone <- predict.lm(lif,predict_data.lif)



t1 <- list(family = "Arial, sans-serif", size = 13, color = "black")
t2 <- list(family = "Arial, sans-serif", size = 12, color = "black")

line_45 <- data.frame("Prediction"= c(70,82),"Actual"=c(70,82))


compare_df <- data.frame(predict_data.lif$predict_respone,data.lif$Response)
colnames(compare_df) <- c("Prediction","Actual")
fig <- plot_ly(compare_df)
fig <- fig %>% add_trace(x=~Actual,y=~Prediction,type="scatter",name="Life Expectancy Prediction",marker = list(color = 'cornflowerblue'))
fig <- fig %>% add_trace(x=line_45$Prediction,y=line_45$Actual,name="1:1 45°Line",type="scatter",mode="line",line=list(width=3,dash='dash',color='red'))
fig <- fig %>% layout(title="Life Expectancy Comparison",
                      yaxis = list(title = 'Predicted Life Expectancy',titlefont=t1,tickfont=t2,showgrid=F,showline = T,range=c(65,85)), 
                      xaxis=list(autotick=FALSE, dtick=5,title="Actual Life Expectancy",titlefont=t1,tickfont=t2,showgrid=F,showline = T,range=c(65,85)),
                      legend = list(orientation = 'v',x=0.1,y=0.8),
                      showlegend=T)
fig

# Check normality of errors - Histogram for Life Expectancy

fig1 <- plot_ly(x=~lif$residuals, type = "histogram",marker=list(color="cornflowerblue"))
fig1 <- fig1 %>% layout(title="Histogram of Residuals - Life Expectancy",
                        yaxis=list(title="Frequency",titlefont=t1,tickfont=t2,showgrid=F,showline=T),
                        xaxis=list(title="Residuals - Life Expectancy",titlefont=t1,tickfont=t2,showgrid=F,showline=T))
fig1
# Check normality of errors - Q-Q Plot for Life Expectancy
autoplot(lif, which = 2, ncol = 1,colour="cornflowerblue") + theme_bw() + labs(title="Q-Q Plot for Life Expectancy")

# Check for Heteroskedasticity
for (i in 2:length(data.lif)){
  col.names <- colnames(data.lif)
  p <- qplot(data.lif[[i]], lif$residuals,color =I("cornflowerblue")) +
    xlab(col.names[i]) + ylab('Residuals') + 
    labs(title="Life Expectancy Residual vs. Predictor") + theme_bw() 
  print(p)
  Sys.sleep(2)
}

# Checking for Multicollinearity
vif(lif)
lif <- stepVIF(lif, threshold = 5, verbose  = TRUE)
summary(lif)

# Checking for outliers
p <- lif$rank - 1 # Calculate number of predictors in the model
n <- nrow(data.lif) # Number of observations in the dataset
autoplot(lif, which = 4, ncol = 1,colour="cornflowerblue") + theme_bw() + labs(title="Cook's Distance - Life Expectancy") 
lif.cooks <-  as.vector(cooks.distance(lif))
lif.outl1 <- which(lif.cooks > 1)
lif.outl1
lif.outl2 <- which(lif.cooks > qf(0.50, p + 1, n - p - 1))
lif.outl2

hap <- lm(Response ~ . -c.code, data = data.hap)
summary(hap)

# Check for relationships of variables to response variable (normalized value)
plot(data.hap, pch = 15, col = "blue")

# After checking for model accuracy we have decided to include higher order terms
# in order to generate a better fit of the model
hap <- lm(Response ~ Precipitation.Days + poly(Heating.Days, 4) + poly(Index, 4) + poly(Precipitation.Annual, 4), data = data.hap)
summary(hap)


# Check normality of errors
fig7 <- plot_ly(x=~hap$residuals, type = "histogram",marker=list(color="aquamarine"))
fig7 <- fig7 %>% layout(title="Histogram of Residuals - Happiness Scores",
                        yaxis=list(title="Frequency",titlefont=t1,tickfont=t2,showgrid=F,showline=T),
                        xaxis=list(title="Residuals - Happiness Scores",titlefont=t1,tickfont=t2,showgrid=F,showline=T))
fig7

autoplot(hap, which = 2, ncol = 1,colour="aquamarine") + theme_bw() + labs(title="Q-Q Plot for Happiness Scores")

# Check for Heteroskedasticity
for (i in 2:length(data.hap)){
  col.names <- colnames(data.hap)
  p <- qplot(data.hap[[i]], hap$residuals,color=I("aquamarine")) +
    xlab(col.names[i]) + ylab('Residuals') +
    labs(title="Happiness Scores Residual vs. Predictor") + theme_bw()
  print(p)
  Sys.sleep(2)
}

# Checking for Multicollinearity
vif(hap)
hap <- stepVIF(hap, threshold = 5, verbose  = TRUE)
summary(hap)

# Checking for outliers
p <- hap$rank - 1 # Calculate number of predictors in the model
n <- nrow(data.hap) # Number of observations in the dataset
autoplot(hap, which = 4, ncol = 1,colour="aquamarine") + theme_bw() + labs(title="Cook's Distance - Happiness")
hap.cooks <-  as.vector(cooks.distance(hap))
hap.outl1 <- which(hap.cooks > 1)
hap.outl1 # No such outliers exist
hap.outl2 <- which(hap.cooks > qf(0.50, p + 1, n - p - 1))
hap.outl2 # No such outliers exist