Econometrics in R

The notebook is used exclusivelly by the students within the seminars of Econometrics at the University of János Selye. We have kindly used the WHO Life Expactancy Data database.

Below, I introduced all the necessary packages needed for the work.

library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
rm(list=ls())

I set set the Session in the Main menu to Sorce File Location. I can to it also in the interactive way (or include to the script) as

setwd(“/cloud/project/ExternalStudy/final”)

But, organization of the folders and subfolders can change from one project to another, that is, why I do not put it into the Chunk.

The Life expectancy data are organized in the csv file, columns separated by “,” and using the decimal comma. I created in my working folde a subfolder called adot, to divide the data and rest of the project. My frame is also called adot.

Not all the data will be used, that is, why I outfiltered just som of the columns for the later use.

As some data are missing, I imputed them by the median values of the variable in the consideration.

adot <- read.csv("adot/Life_Expectancy_Data.csv",dec=".",sep=",",header = TRUE)
# select just the record from 2015
adot.2015 <- adot[adot$Year==2015,c("Life.expectancy","Adult.Mortality","infant.deaths","Hepatitis.B","percentage.expenditure","Measles","thinness..1.19.years","BMI","Income.composition.of.resources","Polio","Diphtheria","HIV.AIDS","GDP","Schooling")]

# data imputation

# Compute column medians
#column_medians <- sapply(adot.2015, median, na.rm = TRUE)

# Impute missing values with column medians
# Compute column medians
column_medians <- sapply(adot.2015, median, na.rm = TRUE)

# Impute missing values with column medians
adot_imputed <- adot.2015
for (col in names(adot.2015)) {
  adot_imputed[[col]][is.na(adot_imputed[[col]])] <- column_medians[col]
}

adot.2015 <- adot_imputed

Now, we want to see the shape of the data (if there are not some irregularities - 0 values for example)

# Create histograms for each variable in the data frame
for (col in names(adot.2015)) {
  boxplot(adot.2015[[col]], main = col, xlab = "Value")
}

The variables infant.deaths, percentage.expenditure, and Measles are not well distributed. It can be the cosequence of the data error and that is, why we exclude these data.

# Exclude specified columns from the data frame
adot_excluded <- adot.2015[, !names(adot.2015) %in% c("infant.deaths", "percentage.expenditure", "Measles")]
adot.2015 <- adot_excluded

At first, I estimated the overparametrized model having the life expectancy as the dependent variable, and all the others as the independet variables.

model <- lm(Life.expectancy ~ .,data=adot.2015)
summary(model)

Call:
lm(formula = Life.expectancy ~ ., data = adot.2015)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.124  -1.467  -0.146   1.709   8.327 

Coefficients:
                                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      5.379e+01  2.116e+00  25.418  < 2e-16 ***
Adult.Mortality                 -2.921e-02  3.484e-03  -8.383 1.78e-14 ***
Hepatitis.B                      3.647e-02  2.525e-02   1.444   0.1506    
thinness..1.19.years            -1.683e-01  6.846e-02  -2.458   0.0149 *  
BMI                              3.161e-03  1.428e-02   0.221   0.8250    
Income.composition.of.resources  2.538e+01  4.705e+00   5.395 2.24e-07 ***
Polio                            1.308e-02  1.332e-02   0.982   0.3274    
Diphtheria                       2.198e-03  2.891e-02   0.076   0.9395    
HIV.AIDS                        -4.656e-01  2.347e-01  -1.984   0.0489 *  
GDP                              2.881e-06  2.502e-05   0.115   0.9085    
Schooling                        1.008e-01  2.116e-01   0.477   0.6343    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.116 on 172 degrees of freedom
Multiple R-squared:  0.8609,    Adjusted R-squared:  0.8528 
F-statistic: 106.5 on 10 and 172 DF,  p-value: < 2.2e-16

The summary of the estimated model provides us with a set of the estimated regression coefficients, the signs of which will be discussed later. If speaking about the properties of the model as a whole, let us check firstly the following Figures. Based on the Q-Q graph, we get the impression about possible problems of the violation of the residuals normality.

plot(model)

At the same time, we run the consequent tests, where we reject the hypothesis of the normal disributed model errors, as well as the presence of the outliers. The violation of the homoskedasticity and improper model specification has been also proved.

# normality tests
residuals <- residuals(model)
jb_test <- jarque.bera.test(residuals)
jb_test

    Jarque Bera Test

data:  residuals
X-squared = 22.873, df = 2, p-value = 1.079e-05

# outlier test (see p-value for Bonferroni correction)
outlier_test <- outlierTest(model)
outlier_test[3]

$bonf.p
       433 
0.02850533 

# heteroskedasticity tests
bp_test <- bptest(model)
bp_test

    studentized Breusch-Pagan test

data:  model
BP = 27.605, df = 10, p-value = 0.002088

# Ramsey Reset test
reset_test <- resettest(model,power=2, type="fitted") # further modification - see help(resettest)
reset_test

    RESET test

data:  model
RESET = 7.448, df1 = 1, df2 = 171, p-value = 0.007015

# multicollinearity
vif(model)

                Adult.Mortality                     Hepatitis.B            thinness..1.19.years 
                       2.164814                        6.610713                        1.461006 
                            BMI Income.composition.of.resources                           Polio 
                       1.617980                        9.246130                        1.953644 
                     Diphtheria                        HIV.AIDS                             GDP 
                       7.692364                        1.866332                        1.377535 
                      Schooling 
                       6.726593 

Based on the test results, we assume there is problem with the residuals normality as well as with the heteroskedasticity. Ramsey Reset test indicates the problem with the bad specification of the model. VIF also does not indicate some signifficant problem with multicollinearity, Variable Income.composition.of.resources seems to be highly correlated with the others (however not exceeding the tresh value 10). As we have a lot of the other variables, we exclude it.

After exclusion of the Income.composition.of.resources variable, we got the following results:

model2 <- lm(Life.expectancy ~ . - Income.composition.of.resources ,data=adot.2015)
summary(model2)

Call:
lm(formula = Life.expectancy ~ . - Income.composition.of.resources, 
    data = adot.2015)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.8044 -1.8035 -0.0944  2.3544  8.2056 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           5.865e+01  2.065e+00  28.394  < 2e-16 ***
Adult.Mortality      -3.435e-02  3.613e-03  -9.508  < 2e-16 ***
Hepatitis.B           1.855e-02  2.699e-02   0.687  0.49278    
thinness..1.19.years -2.027e-01  7.350e-02  -2.758  0.00643 ** 
BMI                   1.587e-02  1.518e-02   1.045  0.29727    
Polio                 2.357e-02  1.420e-02   1.659  0.09884 .  
Diphtheria            1.895e-02  3.099e-02   0.612  0.54166    
HIV.AIDS             -5.976e-01  2.517e-01  -2.375  0.01867 *  
GDP                   3.613e-05  2.614e-05   1.382  0.16868    
Schooling             1.045e+00  1.283e-01   8.144 7.31e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.36 on 173 degrees of freedom
Multiple R-squared:  0.8374,    Adjusted R-squared:  0.8289 
F-statistic: 98.99 on 9 and 173 DF,  p-value: < 2.2e-16

The signs of the statistically significant coeficients are expected but besides the signe of Polio, which has a positive sign. If checking the following figures, we got impression on the signifficant improving of the model statistical wquality

plot(model2)

Our impressions were corroborated as introduced in the following set of tests. The only remaining problem is the problem of heteroskedasticity.

# normality tests
residuals <- residuals(model2)
jb_test <- jarque.bera.test(residuals)
jb_test

    Jarque Bera Test

data:  residuals
X-squared = 1.609, df = 2, p-value = 0.4473

# outlier test (see p-value for Bonferroni correction)
outlier_test <- outlierTest(model2)
outlier_test[3]

$bonf.p
      401 
0.4061314 

# heteroskedasticity tests
bp_test <- bptest(model2)
bp_test

    studentized Breusch-Pagan test

data:  model2
BP = 29.948, df = 9, p-value = 0.0004478

# Ramsey Reset test
reset_test <- resettest(model2,power=2, type="fitted") # further modification - see help(resettest)
reset_test

    RESET test

data:  model2
RESET = 3.7825, df1 = 1, df2 = 172, p-value = 0.05342

# multicollinearity
vif(model2)

     Adult.Mortality          Hepatitis.B thinness..1.19.years                  BMI                Polio 
            2.002605             6.496412             1.448322             1.573912             1.911980 
          Diphtheria             HIV.AIDS                  GDP            Schooling 
            7.603608             1.846048             1.293914             2.127115 

Diphteria is not statistically signifficant and is highly correlated with other variables (see VIF). Let us exclude it.

model3 <- lm(Life.expectancy ~ . - Income.composition.of.resources - Diphtheria,data=adot.2015)
summary(model3)

Call:
lm(formula = Life.expectancy ~ . - Income.composition.of.resources - 
    Diphtheria, data = adot.2015)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.7979 -1.8515 -0.0687  2.3003  8.2266 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           5.882e+01  2.042e+00  28.805  < 2e-16 ***
Adult.Mortality      -3.444e-02  3.603e-03  -9.559  < 2e-16 ***
Hepatitis.B           3.273e-02  1.378e-02   2.375  0.01863 *  
thinness..1.19.years -2.022e-01  7.336e-02  -2.756  0.00648 ** 
BMI                   1.504e-02  1.509e-02   0.996  0.32052    
Polio                 2.593e-02  1.364e-02   1.901  0.05900 .  
HIV.AIDS             -6.122e-01  2.501e-01  -2.448  0.01536 *  
GDP                   3.814e-05  2.588e-05   1.474  0.14238    
Schooling             1.052e+00  1.274e-01   8.259 3.56e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.354 on 174 degrees of freedom
Multiple R-squared:  0.837, Adjusted R-squared:  0.8295 
F-statistic: 111.7 on 8 and 174 DF,  p-value: < 2.2e-16

The only statistically insignifficant variables remained GDP and BMI. All the statistically signifficant coefficients have expected signs, besides Polio.

To get some imagine about the statistical properties of the model as a whole, let us look at the following figures. All the graphical analysis are rather optimistic if speaking about the statistical quality of the model.

plot(model3)

# normality tests
residuals <- residuals(model3)
jb_test <- jarque.bera.test(residuals)
jb_test

    Jarque Bera Test

data:  residuals
X-squared = 1.7349, df = 2, p-value = 0.42

# outlier test (see p-value for Bonferroni correction)
outlier_test <- outlierTest(model3)
outlier_test[3]

$bonf.p
    401 
0.40155 

# heteroskedasticity tests
bp_test <- bptest(model3)
bp_test

    studentized Breusch-Pagan test

data:  model3
BP = 30.266, df = 8, p-value = 0.0001897

# Ramsey Reset test
reset_test <- resettest(model3,power=2, type="fitted") # further modification - see help(resettest)
reset_test

    RESET test

data:  model3
RESET = 3.7192, df1 = 1, df2 = 173, p-value = 0.05543

# multicollinearity
vif(model3)

     Adult.Mortality          Hepatitis.B thinness..1.19.years                  BMI                Polio 
            1.999196             1.700070             1.448099             1.561153             1.770545 
            HIV.AIDS                  GDP            Schooling 
            1.829413             1.273366             2.106336 

If putting out GDP variable, the results worsened. that is, why the model3 is considered as the best one and we are going to solve the heteroskedasticity problem by the introduction of the t-tests using the White heteroskedasticity consistent matrix as follows

# Load libraries
library("lmtest")
library("sandwich")

# Robust t test
coeftest(model3, vcov = vcovHC(model3, type = "HC0"))

t test of coefficients:

                        Estimate  Std. Error t value  Pr(>|t|)    
(Intercept)           5.8820e+01  2.6650e+00 22.0714 < 2.2e-16 ***
Adult.Mortality      -3.4441e-02  5.8121e-03 -5.9258 1.627e-08 ***
Hepatitis.B           3.2734e-02  1.4013e-02  2.3360   0.02063 *  
thinness..1.19.years -2.0217e-01  7.8692e-02 -2.5691   0.01104 *  
BMI                   1.5036e-02  1.7635e-02  0.8526   0.39504    
Polio                 2.5933e-02  1.5622e-02  1.6600   0.09872 .  
HIV.AIDS             -6.1216e-01  3.4428e-01 -1.7781   0.07714 .  
GDP                   3.8144e-05  1.8118e-05  2.1054   0.03669 *  
Schooling             1.0525e+00  1.5447e-01  6.8135 1.496e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

##Conclusion##

Variables Adult.mortality, thinness and HIV.AIDS shorten the expected length of life. On the other side, Polio and Hepatitis.B seem to extend the length of life, what is the fact, we can not interpret. Probably the prevent has some long term impacts, or some vaccination treatment coming after the period of higher amount of these diseases was helping? Schooling, percentage.expenditure and GDP are the economic factors significantly improving the expected length of life.