RPubs link information

• Rpubs link: https://rpubs.com/Pasmans/681096

Introduction

• Life expectancy is a metric used world-wide to establish the expected age of death of a person born in a specific year according to the mortality rates of that year.

• While people born in some countries can expect to live well into their 80’s those born in areas amidst political conflicts or with less access to healthcare can be unlikely to live past 60.

• Life Expectancy in a country can be influenced by many factors as it is essentially a measure of all causes of mortality combined. By investigating and understanding which variables influence life expectancy, plans to increase the overall life expectancy can be established.

• This study focuses on the relationship between life expectancy and GNI per Capita.

• GNI (gross national income) per capita is the sum of value added by all resident producers of a country divided by the midyear population of that country. It is generally considered a good measure of a countries wealth.

• Since wealthier countries generally have better access to food, water and healthcare it is expected that there will be a positive linear relationship between life expectancy and GNI per Capita.

Problem Statement

• Is there a linear relationship between life expectancy and GNI per Capita?

• By utilizing linear regression and correlation testing this investigation seeks to establish whether there is a statistically significant linear relationship between GNI per Capita and life expectancy.

• Summary statistics and visualizations of the relevant variables will establish if any transformations are needed to increase normality.

• A linear regression model will be calculated and it’s results tested at the 95% significance level to determine if they are statistically significant.

• The slope, intercept and multiple r squared together with the pearson correlation will determine strength and characteristics of the relationship between GNI per Capita and life expectancy.

Data

• The Dataset was taken from the https://databank.worldbank.org/source/world-development-indicators#

• The World Bank is a financial institution which offers loans and grants to lower income countries. They have a substantial open database that sources data from a number of institutions.

• They list the source for their life expectancy data as United Nations Population Division, Census reports and other statistical publications from national statistical offices. GNI data was taken from their World Bank national accounts data, and OECD National Accounts data files.

• The database allows the user to select the variables and layout of a dataset. This specific dataset was taken from their World Development Indicators database.

• All individual countries were selected excluding country groups such as the European Union. The variables chosen were “Life expectancy at birth, total (years)” and “GNI per capita, Atlas method (current US$)”. Finally, the year 2018 was selected as that was the most recent year that had information for life expectancy in the database.

Variables.

• Country: (chr) The name of the country or region the observations information is taken from.
• Code: (chr) A 3 character unique identifier for that country
• GNI per Capita USD: (dbl) GNI is calculated as the total amount of money earned by a nation’s people and businesses. Including income the nations gross domestic product and all income received from foreign sources. This figure is then converted to current USD using the average exchange rate for that nations local currency in 2018 and the two preceding years. This is then divided by that nations mid year population in 2018.
• Life Expectancy: (dbl) An estimate of the expected age (in years) members of a population will be when they die if current mortality rates for each age remained the same.
• Income Group: (ordered Factor) Income groups uses 2018 gross national income (GNI) per capita of each country to separate them into the following levels:
  • low income: $1,035 or less
  • lower middle income: $1,036 - 4,045
  • upper middle income: $4,046 - 12,535
  • high income: $12,536 or more
• Region: (Factor) The region that the given nation resides in.
  • East Asia & Pacific
  • Europe & Central Asia
  • Latin America & Caribbean
  • Middle East & North Africa
  • North America
  • South Asia
  • Sub-Saharan Africa

Processing data

• Column names were changed to remove redundant information
• As the dataset contained “..” for missing data, these values were imputed with NA’s so that the relevant columns could be converted to numerics
• Rows with NA’s were excluded from the dataset. as they made up a small percentage of the data they could likely be excluded without introducing a bias.
• Income group and region were converted to an ordered factor and an un-ordered factor respectively.

# Reading the Data
gniLife <- read_excel("GNILifeExpectancy.xlsx", sheet = "Data")

# Excluding unwanted rows and columns
gniLife <- gniLife[1:217, 3:6]

# Reading Country metadata and joining Country metadata
metadata <- read_xlsx("Data_Extract_From_World_Development_Indicators_Metadata.xlsx", sheet = "Country - Metadata")
gniLife <- left_join(gniLife, metadata[c(1, 3, 4)], by = c("Country Code" = "Code")) 

# Simplifying column names
colnames(gniLife) <- c("Country", "Code", "GNI.per.capita.USD", "Life.Expectancy", "Income.Group", "Region")

# Changing missing values to NA
gniLife$`Life.Expectancy`[gniLife$`Life.Expectancy` == ".."] <- NA
gniLife$`GNI.per.capita.USD`[gniLife$`GNI.per.capita.USD` == ".."] <- NA

# converting to numerics
gniLife$`Life.Expectancy` <- as.numeric(gniLife$`Life.Expectancy`)
gniLife$`GNI.per.capita.USD` <- as.numeric(gniLife$`GNI.per.capita.USD`)

# Excluding NA's
gniLife <- na.omit(gniLife) 

# Converting to Factors
gniLife$Region <- factor(gniLife$Region)
gniLife$`Income.Group` <- factor(gniLife$`Income.Group`,
                                           levels = c("Low income", "Lower middle income", "Upper middle income", "High income"),
                                           labels = c("Low", "Lower Middle", "Upper Middle", "High"),
                                           ordered = TRUE)

# Head of complete dataset
head(gniLife)

Descriptive Statistics

• Observing GNI summary statistics the Median (5340) is considerably lower than the mean (13916.94), this suggests a right skewed distribution. There is also a large amount of variability in the data with a standard deviation of 18771.78 and large range between the minimum (280) and maximum (84430) value.

• The median of Life Expectancy (73.82) was slightly higher than the mean (72.48) suggesting a slightly left-skewed distribution but overall the life expectancy had a much smaller variation. standard deviation = 7.51.

# GNI Descriptive Statistics
summaryGNI <- gniLife %>% summarise(Variable = "GNI per Capita",
                                        Min = min(`GNI.per.capita.USD`),
                                        Q1 = quantile(`GNI.per.capita.USD`,probs = .25),
                                        Median = median(`GNI.per.capita.USD`),
                                        Q3 = quantile(`GNI.per.capita.USD`,probs = .75),
                                        Max = max(`GNI.per.capita.USD`),
                                        Mean = mean(`GNI.per.capita.USD`),
                                        IQR = IQR(`GNI.per.capita.USD`),
                                        SD = sd(`GNI.per.capita.USD`),
                                        n = n())
                                      

# Life Expectancy Descriptive Statistics
summaryLife <- gniLife %>% summarise(Variable = "Life.Expectancy",
                                        Min = min(`Life.Expectancy`),
                                        Q1 = quantile(`Life.Expectancy`,probs = .25),
                                        Median = median(`Life.Expectancy`),
                                        Q3 = quantile(`Life.Expectancy`,probs = .75),
                                        Max = max(`Life.Expectancy`),
                                        Mean = mean(`Life.Expectancy`),
                                        IQR = IQR(`Life.Expectancy`),
                                        SD = sd(`Life.Expectancy`),
                                        n = n())

summary <- rbind(summaryGNI, summaryLife)
knitr::kable(summary, align = "c", digits = 2, caption = "Summary Statistics")

life expectancy Descriptive Statistics grouped by Income Group

• Grouping life expectancy summary statistics by Income Group, a pattern can be seen.
• The low Income Group had the lowest mean (62.04) and median (61.88) while the high Income Group had a much higher mean (79.79) and median (80.89).
• Using summary statististics alone it can be seen that as income group increases life expectancy also tends to increase

# Life Expectancy grouped by Income Group descriptive statistics
summaryLifeByIncome <- gniLife %>% group_by(gniLife$`Income.Group`) %>% 
                                        summarise(Min = min(`Life.Expectancy`),
                                        Q1 = quantile(`Life.Expectancy`,probs = .25,na.rm = TRUE),
                                        Median = median(`Life.Expectancy`, na.rm = TRUE),
                                        Q3 = quantile(`Life.Expectancy`,probs = .75,na.rm = TRUE),
                                        Max = max(`Life.Expectancy`,na.rm = TRUE),
                                        Mean = mean(`Life.Expectancy`, na.rm = TRUE),
                                        IQR = IQR(`Life.Expectancy`),
                                        SD = sd(`Life.Expectancy`, na.rm = TRUE),
                                        n = n())

knitr::kable(summaryLifeByIncome, align = "c", digits = 2, caption = "Life Expectancy Summary Statistics by Income Group")

Boxplots

• Plotting a boxplot for both variables shows a strong right skew for GNI and a moderate left skew for life expectancy.
• For the purpose of this study to reduce the variation caused by outliers and increase normality in order to apply a linear regression test, a transformation will be applied to both variables

# GNI Histograms
par(mfrow=c(1,2))  
boxplot(gniLife$`GNI.per.capita.USD`,
      main = "GNI per Capita", 
      ylab = "GNI per Capita USD",
      col = I("thistle3"))

boxplot(gniLife$`Life.Expectancy`,
      main = "Life Expectancy", 
      ylab = "Life Expectancy at Birth (Years)",
      col = I("paleturquoise3"))

Histograms

• Histograms are used to visualize the distribution of the data before and after transformation.
• The transformed variables are overlayed with a normal distribution curve to see if they approximate normality

# GNI Histograms
par(mfrow=c(2,2))
plot1 <- ggplot() + geom_histogram(aes(gniLife$`GNI.per.capita.USD`, y = ..density..), bins = 20, fill = "thistle3", color = "Black" ) +
    labs(title = "GNI per Capita Histogram", x = "GNI.per.capita.USD")

log10GNI <- log10(gniLife$`GNI.per.capita.USD`)
plot2 <- ggplot() + geom_histogram(aes(log10GNI, y = ..density..), bins = 20, fill = "thistle3", color = "Black" ) +
    labs(title = "log10(GNI per Capita) Histogram", x = "log10(GNI.per.capita.USD)") +
    stat_function(fun = dnorm, args = list(mean = mean(log10GNI), sd = sd(log10GNI)))

# Life Expectancy Histograms
plot3 <- ggplot() + geom_histogram(aes(gniLife$`Life.Expectancy`, y = ..density..), bins = 20, fill = "paleturquoise3", color = "Black" ) +
    labs(title = "Life Expectancy", x = "Life.Expectancy")

lifeCubed <- gniLife$`Life.Expectancy`^ 3
plot4 <- ggplot() + geom_histogram(aes(lifeCubed, y = ..density..), bins = 20, fill = "paleturquoise3", color = "Black" ) +
    labs(title = "Life Expectancy^3 Histogram", x = "Life.Expectancy^3") +
    stat_function(fun = dnorm, args = list(mean = mean(lifeCubed), sd = sd(lifeCubed)))

grid.arrange(plot1, plot2, plot3, plot4, ncol=2)

Q-Q plots of transformed variables

• Testing the transformed variables for normality using qqplots it is clear that the transformed GNI mostly falls within the 95% CI for the normal quantiles. There are however a few points edging out at the tails of the distribution suggesting the tails are slightly more dense than they would be with a normal distribution
• Life Expectancy^3 also appears to approximate normality with a few values falling outside the lines of the 95% CI.
• While we should be cautious in assuming normality completely. Both transformed variables appear to approximate a normal distribution

# qqplots to test normality
par(mfrow=c(1,2))
log10GNI %>% qqPlot(dist="norm",
                    main = "log10(GNI per capita)",
                    col.lines = ("blue"))

## [1]  27 100

lifeCubed %>% qqPlot(dist="norm",
                      main = "Life Expectancy^3", 
                      col.lines = ("red"))

## [1] 32 93

Hypothesis Testing

• Plotting GNI and life expectancy on a scatterplot exhibits a clear positive linear relationship.

• As GNI increases, life expectancy also tends to increase.

• As the data appears to show signs of a positive linear relationship we can move on to fitting a linear regression model

• The linear regression equation \(y = α + βx + ϵ\) was used to calculate the line of best fit.

  • Predictor (x) = log10 GNI per Capita
  • Dependent (y) = Life Expectancy^3

• This line was then applied to the scatterplot.

# Preparing line of best fit
sum_x <- sum(log10GNI)
sum_y <- sum(lifeCubed)
sum_x_sq <- sum(log10GNI^2)
sum_y_sq <- sum(lifeCubed^2)
sum_xy <- sum(log10GNI*lifeCubed)
n <- length(log10GNI)
Lxx <-  sum_x_sq-((sum_x^2)/n)
Lyy <- sum_y_sq-((sum_y^2)/n)
Lxy = sum_xy - (((sum_x)*(sum_y))/n)
b = Lxy/Lxx
a = mean(lifeCubed - b*mean(log10GNI))

# Plotting the line
plot(lifeCubed ~ log10GNI, main = "Line of Best Fit", 
     xlab = "log10 GNI per Capita", 
     ylab = " Life Expectancy^3")
abline(a = a, b = b, col = "red")

Linear regression - F test and Multiple R-Squared

• To test whether the data fits a linear regression model the F test is used with a \(α = 0.05\) level of significance.
  • \(H_0\): The data does not fit the linear regression model
  • \(H_A\): The data fits the linear regression model.
• The F-test reported a p-value < 0.001 using an F-statistic of 553 on 1 and 178 DF.
• This indicates that the overall regression model is statistically significant.
• The Multiple R-Squared value of 0.756 indicates 75.6% of the variability in life expectancy^3 can be explained by log10 GNI.

gniLifeModel <- lm(lifeCubed ~ log10GNI)
gniLifeModel %>% summary()

## 
## Call:
## lm(formula = lifeCubed ~ log10GNI)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -203643  -29571    4764   36552  102128 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -214909      26178   -8.21 4.36e-14 ***
## log10GNI      161595       6872   23.52  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 56240 on 178 degrees of freedom
## Multiple R-squared:  0.7565, Adjusted R-squared:  0.7551 
## F-statistic:   553 on 1 and 178 DF,  p-value: < 2.2e-16

Linear regression - Intercept and slope

• The regression model reported a statistically significant intercept of \(α = -214908.7\), \(p < 0.001\).
• When log10 GNI is 0, the expected Life Expectancy^3 is -214908.7
• The slope was reported as \(β = 161594.9\) and was also statistically significant \(p < 0.001\).
• Every one unit increase in log10 GNI corresponds to an increase of 161594.9 in life expectancy^3

# Slope And intercept
gniLifeModel %>% summary() %>% coef()

##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) -214908.7  26177.520 -8.209664 4.358542e-14
## log10GNI     161594.9   6871.854 23.515475 1.721335e-56

• Testing the confidence intervals for intercept and slope at a 0.05 significance level it is found that \(H_0:α=0\) and \(H_0:β=0\) are not captured in the confidence interval and therefore the null hypothesis is rejected.
• The slope and intercept are statistically significant.

# Testing confidence intervals
gniLifeModel %>% confint()

##                 2.5 %    97.5 %
## (Intercept) -266566.9 -163250.4
## log10GNI     148034.1  175155.7

Linear regression - Diagnostics

• The following assumptions have been confirmed.
  • Independence: All observations of countries are separate and no country is measured multiple times.
  • Linearity: A bivariate linear relationship between log10 GNI and Life Expectancy^3 was confirmed using a scatterplot
• Normality of residuals and Homoscedasticity are tested using the following diagnostics.
  • Residuals vs Fitted: The variance of residuals remains approximately constant. The generally flat trend line has a slight bump in the middle but is mostly a good indication of homoscedasticity.
  • Normal Q-Q: While most residuals fall close to the line of normality there are some values falling off towards the tails. values 51, 121 and 53 could potentially be removed from the model.
  • Scale-Location: The trend line is close to flat and the standardized residuals are generally consistent further indicating homoscedasticity.
  • Residuals vs Leverage: All points fall well within red bands. There are no values unduly influencing the fit of the regression model.

# Diagnostics of fitted regression model
par(mfrow = c(2, 2))
gniLifeModel %>% plot(which=1)
gniLifeModel %>% plot(which=2)
gniLifeModel %>% plot(which=3)
gniLifeModel %>% plot(which=5)

Pearson Correlation

• The Pearson correlation shows a statistically significant strong positive linear relationship \(r = .87\), \(p < .001\).
• The 95% confidence interval for the correlation is tested using the hypothesis test… \[H_0:r=0 \] \[H_A:r≠0 \]
• The Confidence interval [0.8289567, 0.9013615] did not capture 0 so the null hypothesis was rejected.
• The positive correlation was statistically significant.

# Pearson Correlation
bivariate <- as.matrix(cbind(lifeCubed,log10GNI))
rcorr(bivariate, type = "pearson")

##           lifeCubed log10GNI
## lifeCubed      1.00     0.87
## log10GNI       0.87     1.00
## 
## n= 180 
## 
## 
## P
##           lifeCubed log10GNI
## lifeCubed            0      
## log10GNI   0

# Testing Confidence interval
r <- cor(lifeCubed,log10GNI)
CIr(r = r, n = 180, level = .95)

## [1] 0.8289567 0.9013615

Discussion

• The dependent variable (life expectancy) and predictor variable (GNI per Capita) were transformed in order to increase Normality and reduce spread.
• A scatterplot confirmed a positive linear relationship between the two variables and a linear regression model was established.
• The overall regression model was confirmed to be statistically significant \(F(1,178) = 553\), \(p < 0.001\) and reported that 75.6% of the variability in life expectancy^3 can be explained by a positive linear relationship with log10(GNI). \(R^2 = 0.756\).
• The equation for the regression was estimated to be Life Expectancy^3 = -214909 + 161594.9 * log10(GNI). Both the slope and intercept were found to be statistically significant.
  • \(α = -214908.7\), \(p < 0.001\), 95% CI [-266566.9 -163250.4].
  • \(β = 161594.9\), \(p<.001\), 95% CI [148034.1 175155.7].
• The strength of the positive linear relationship was measured by calculating the Pearson correlation. A statistically significant strong positive correlation of \(r = 0.87\) was found \(p < .001\), 95% CI[0.8289567 0.9013615]
• The linear regression assumptions of normality of residuals and homoscedasticity were validated by inspecting the distribution of residuals. While there were some minor deviations from normality in the residuals the overall distribution was approximately normal. In future tests the offending points could be considered for removal from the regression model.

Discussion Continued

• The linear regression test shows conclusively that in 2018, wealthier countries (according to GNI per Capita) tend to have higher life expectancy.

• The largest limitation in this investigation is that GNI per Capita and life expectancy are both quite broad variables that are influenced by and have influence on many other variables.

• Assumptions could be made that wealthier countries are more likely to have lower crime rates and better access to healthcare and that those factors are the main influence on life expectancy. However, those assumptions would be outside the scope of this investigation and while the data shows there is a strong positive correlation between GNI per Capita and life expectancy, the cause of this correlation would require further study.

• Future investigations could focus on more specific aspects of countries, for example, access to healthcare or food and water. Grouping countries by income group to investigate which factors influence Life expectancy within that specific group could also provide further insight.

• While further studies would be needed to establish the underlying cause of the correlation, the results of this study showed a strong, statistically significant positive linear relationship between GNI and Life Expectancy during the year 2018.

References

Data was taken from the https://databank.worldbank.org/source/world-development-indicators#