Assignment

The dataset for this assignment contains real-world data from 2008. Using the dataset answer the following questions:

  1. Provide a scatterplot of LifeExp~TotExp, and run simple linear regression. Do not transform the variables. Provide and interpret the F statistics, R^2, standard error,and p-values only. Discuss whether the assumptions of simple linear regression met.
  2. Raise life expectancy to the 4.6 power (i.e., LifeExp^4.6). Raise total expenditures to the 0.06 power (nearly a log transform, TotExp^.06). Plot LifeExp^4.6 as a function of TotExp^.06, and r re-run the simple regression model using the transformed variables. Provide and interpret the F statistics, R^2, standard error, and p-values. Which model is “better?”
  3. Using the results from 3, forecast life expectancy when TotExp^.06 =1.5. Then forecast life expectancy when TotExp^.06=2.5.
  4. Build the following multiple regression model and interpret the F Statistics, R^2, standard error, and p-values. How good is the model? LifeExp = b0+b1 x PropMd + b2 x TotExp +b3 x PropMD x TotExp
  5. Forecast LifeExp when PropMD=.03 and TotExp = 14. Does this forecast seem realistic? Why or why not?

Dataset

Variable Description
Country name of the country
LifeExp average life expectancy for the country in years
InfantSurvival proportion of those surviving to one year or more
Under5Survival proportion of those surviving to five years or more
TBFree proportion of the population without TB
PropMD proportion of the population who are MDs
PropRN proportion of the population who are RNs
PersExp mean personal expenditures on healthcare in US dollars at average exchange rate
GovtExp mean government expenditures per capita on healthcare, US dollars at average exchange rate
TotExp sum of personal and government expenditures 
library(DT)
data <- read.csv("https://raw.githubusercontent.com/suswong/DATA-605/main/DATA%20605%20Week%2012%20dataset.csv")
datatable(data, options = list(scrollX = TRUE))

Packages

library(openintro)
## Loading required package: airports
## Loading required package: cherryblossom
## Loading required package: usdata
library(psych)
library(PerformanceAnalytics)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
library(Hmisc)
## 
## Attaching package: 'Hmisc'
## The following object is masked from 'package:psych':
## 
##     describe
## The following objects are masked from 'package:base':
## 
##     format.pval, units
library(ggplot2)
## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha

Exploratory Analysis

Correlation

The average life expectancy has a strong relationship with the following variables:

suppressWarnings({
chart.Correlation(data[2:10], histogram = TRUE, method = "pearson") })

Average Life Expectancy

It is not a normal distribution. The distribtion is skewed to the left. This make sense as most people like to around 70-80 years old, which is near the median and mean of average life expectancy.

#library(Hmisc)
#hist.data.frame(data[2:10]) # This code will plot the histogram for variables in columns 2 to 10. However, it is a bit difficult to see the distribution because the scales are the same

#for (col in 2:ncol(data)){
#  hist(data[,col], main = colnames(data)[col], xlab = colnames(data)[col])
#}

ggplot(data = data, aes(x = LifeExp)) + geom_histogram(bins = 10)  +  
  labs(x = "Average Life Expectancy" , y = "count") + 
  ggtitle("Average Life Expectancy Histogram") + theme(plot.title=element_text(hjust=0.5))+
  geom_vline(xintercept = mean(data$LifeExp), 
             color = "indianred") +
  geom_vline(xintercept = median(data$LifeExp), 
             color = "cornflowerblue")

Infant Survival

hist(data$InfantSurvival)

Linear Regression

Correlation

Although there is a correlation between the the total expenditures and average life expectancy, it does not show a liner relationship.

ggplot(data = data, aes(x = TotExp, y = LifeExp)) + geom_point()  +  geom_smooth() +
  labs(x = "Sum of personal and government expenditures" , y = "Average Life Expectancy") + 
  ggtitle("Total Expenditures v. Average Life Expectancy") + theme(plot.title=element_text(hjust=0.5))
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Model 1

\(LifeExp = b_{0}+b_{1} \cdot TotExp\)

\(LifeExp = 64.75 + 6.297\cdot10^{-5} \cdot TotExp\)

Since the p-value (7.714e-14) associated with the F-statistic (65.26 with degrees of freedom = 1) is less than 0.05, it means that at least 1 independent value is related to the average life expectancy.

Since our model contains one independent variable, the p-value associated with the F-statistic is the same as the p-value of the independent variable, TotExp.

The adjusted R-squared is 0.2537, which means this model can explain 25.37% of the variability in the average life expectancy.

m1 <- lm(LifeExp~TotExp , data = data) 
summary(m1)
## 
## Call:
## lm(formula = LifeExp ~ TotExp, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.764  -4.778   3.154   7.116  13.292 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.475e+01  7.535e-01  85.933  < 2e-16 ***
## TotExp      6.297e-05  7.795e-06   8.079 7.71e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.371 on 188 degrees of freedom
## Multiple R-squared:  0.2577, Adjusted R-squared:  0.2537 
## F-statistic: 65.26 on 1 and 188 DF,  p-value: 7.714e-14

Model Diagnostic

This model does not meet all the conditions.

Linearity and Constant Variability: Residuals indicates the difference between the actual and predicted value. We want the the mean of residuals to be close to 0. Based on the plot below, the points are distributed around 0 in no apparent pattern.

Both linearity and constant variability are not met. The points are not distributed randomly around 0.

ggplot(data = m1, aes(x = .fitted, y = .resid)) + geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") + xlab("Fitted values") +
  ylab("Residuals")  + 
  ggtitle("Residuals v. Fitted") + theme(plot.title=element_text(hjust=0.5))

Normality of Residuals: All of the points should lie roughly on the line. The assumption of normality is not met. The histogram of residual does not show a normal distribution.

ggplot(data = m1, aes(x = .resid)) + geom_histogram() + xlab("Residuals")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(data = m1, aes(sample = .resid)) + stat_qq() +stat_qq_line(col = "red")  + 
  xlab("Theoretical Quantities") + ylab("Standardized Resididuals")+ggtitle("Normal Q~Q ") + theme(plot.title=element_text(hjust=0.5))

Model 2

We want to transform the varaible to achieve approximately symmetry and homoscedasticity of the residuals. Thus, raise life expectancy to the 4.6 power (i.e., LifeExp^4.6) and raise total expenditures to the 0.06 power (nearly a log transform, TotExp^.06)

data$LifeExp_4.6 <- (data$LifeExp)^4.6
data$TotExp_0.06 <- (data$TotExp)^0.06

Correlation

After the transformations, the the total expenditures and average life expectancy shows a more linear relationship than before transformations.

ggplot(data = data, aes(x = TotExp_0.06, y = LifeExp_4.6)) + geom_point()  +  geom_smooth() +
  labs(x = "Sum of personal and government expenditures" , y = "Average Life Expectancy") + 
  ggtitle("Total Expenditures v. Average Life Expectancy") + theme(plot.title=element_text(hjust=0.5))
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Model 1 v. Model 2

\(LifeExp^{4.6} = 620060216 \cdot TotExp^{0.06} - 736527910\)

Since the p-value (< 2.2e-16) associated with the F-statistic (507.7 with degrees of freedom = 1) is less than 0.05, it means that at least 1 independent value is related to the average life expectancy.

Since our model contains one independent variable, the p-value associated with the F-statistic is the same as the p-value of the independent variable, TotExp.

The adjusted R-squared is 0.7283, which means this model can explain 72.83% of the variability in the average life expectancy. The adjusted R-squared is higher than the last model. This can indicate that this is a better model than the first model.

m2 <- lm(LifeExp_4.6 ~ TotExp_0.06 , data = data) 
summary(m2)
## 
## Call:
## lm(formula = LifeExp_4.6 ~ TotExp_0.06, data = data)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -308616089  -53978977   13697187   59139231  211951764 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -736527910   46817945  -15.73   <2e-16 ***
## TotExp_0.06  620060216   27518940   22.53   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 90490000 on 188 degrees of freedom
## Multiple R-squared:  0.7298, Adjusted R-squared:  0.7283 
## F-statistic: 507.7 on 1 and 188 DF,  p-value: < 2.2e-16

Model Diagnostic

This model meets the conditions.

Linearity and Constant Variability: Residuals indicates the difference between the actual and predicted value. We want the the mean of residuals to be close to 0. Based on the plot below, the points are distributed around 0 in no apparent pattern.

Both linearity and constant variability are met.

ggplot(data = m2, aes(x = .fitted, y = .resid)) + geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") + xlab("Fitted values") +
  ylab("Residuals")  + 
  ggtitle("Residuals v. Fitted") + theme(plot.title=element_text(hjust=0.5))

Normality of Residuals: All of the points should lie roughly on the line. The assumption of normality is met.

ggplot(data = m2, aes(x = .resid)) + geom_histogram() + xlab("Residuals")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(data = m2, aes(sample = .resid)) + stat_qq() +stat_qq_line(col = "red")  + 
  xlab("Theoretical Quantities") + ylab("Standardized Resididuals")+ggtitle("Normal Q~Q ") + theme(plot.title=element_text(hjust=0.5))

Forecast Life Expectancy Using Model 2

Using the results from model 2, forecast life expectancy when TotExp^.06 =1.5. Then forecast life expectancy when TotExp^.06=2.5.

\(LifeExp^{4.6} = 620060216 \cdot TotExp^{0.06} - 736527910\)

When TotExp^.06 = 1.5, the average life expectancy is 19,3562,414.

x<- 1.5
620060216*x - 736527910
## [1] 193562414

When TotExp^.06 = 1.5, the average life expectancy is 813,622,630.

x<- 2.5
620060216*x - 736527910
## [1] 813622630

Multiple Regression Model

Model 3

Build the following multiple regression model and interpret the F Statistics, R^2, standard error, and p-values. How good is the model?

\(LifeExp = b_{0}+b_{1} \cdot PropMd + b_{2} \cdot TotExp + b_{3} \cdot PropMD \cdot TotExp\)

\(LifeExp = 6.277e+01 +1.497e+03 \cdot PropMd + 7.233e-05 \cdot TotExp + -6.026e-03 \cdot PropMD \cdot TotExp\)

Since the p-value (< 2.2e-16) associated with the F-statistic (34.49 with degrees of freedom = 3) is less than 0.05, it means that at least 1 independent value is related to the average life expectancy.

All of the p-values of the independent variables is less than 0.05. This indicates they are statistically signnifiacnt predictor of our dependent variable.

The adjusted R-squared is 0.3471, which means this model can explain 34.71% of the variability in the average life expectancy. The adjusted R-squared is lower than the adjusted R-squared in model 2.

m3 <- lm(LifeExp ~ PropMD + TotExp + (PropMD * TotExp), data = data)
summary(m3)
## 
## Call:
## lm(formula = LifeExp ~ PropMD + TotExp + (PropMD * TotExp), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -27.320  -4.132   2.098   6.540  13.074 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    6.277e+01  7.956e-01  78.899  < 2e-16 ***
## PropMD         1.497e+03  2.788e+02   5.371 2.32e-07 ***
## TotExp         7.233e-05  8.982e-06   8.053 9.39e-14 ***
## PropMD:TotExp -6.026e-03  1.472e-03  -4.093 6.35e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.765 on 186 degrees of freedom
## Multiple R-squared:  0.3574, Adjusted R-squared:  0.3471 
## F-statistic: 34.49 on 3 and 186 DF,  p-value: < 2.2e-16

Forecast Life Expectancy Using Model 3

Forecast LifeExp when PropMD=.03 and TotExp = 14. Does this forecast seem realistic? Why or why not?

\(LifeExp = 62.77 + 1,497 \cdot PropMd + 0.00007233 \cdot TotExp - 0.006026 \cdot PropMD \cdot TotExp\)

When PropMD is \(0.03\) and TotExp is \(14\), the forecasted average life expectancy is 107.6785. This is not realistic as the highest average life expectancy in a country (Monaco) is 87 in the world.

PropMd <-.03

TotExp <- 14

LifeExp <-  62.77 + 1497 * PropMd + 0.00007233 * TotExp - 0.006026 * PropMd * TotExp

LifeExp
## [1] 107.6785