Linear Regression Application - Forecasting

Introduction

Prediction is one most the main goal of building a regression model. In one my blogs ( read here), I spoke about mathematics transformation for linear regression. Today, we are going to use the same data set and go over a transformed model to predict some values.

Just to recall, it ’s a real-world data from 2008 collected by WHO, and can be accessed on my github repository here

The variables included follow. 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.

We are going to focus on predicting “LifeExp” based on “TotExp”.

Data Exploration

First, let load the data and have a look on that.

df <- read.csv('https://raw.githubusercontent.com/jnataky/Computational_Mathematics/main/Assignments/who.csv')

glimpse(df)

## Rows: 190
## Columns: 10
## $ Country        <chr> "Afghanistan", "Albania", "Algeria", "Andorra", "Angola~
## $ LifeExp        <int> 42, 71, 71, 82, 41, 73, 75, 69, 82, 80, 64, 74, 75, 63,~
## $ InfantSurvival <dbl> 0.835, 0.985, 0.967, 0.997, 0.846, 0.990, 0.986, 0.979,~
## $ Under5Survival <dbl> 0.743, 0.983, 0.962, 0.996, 0.740, 0.989, 0.983, 0.976,~
## $ TBFree         <dbl> 0.99769, 0.99974, 0.99944, 0.99983, 0.99656, 0.99991, 0~
## $ PropMD         <dbl> 0.000228841, 0.001143127, 0.001060478, 0.003297297, 0.0~
## $ PropRN         <dbl> 0.000572294, 0.004614439, 0.002091362, 0.003500000, 0.0~
## $ PersExp        <int> 20, 169, 108, 2589, 36, 503, 484, 88, 3181, 3788, 62, 1~
## $ GovtExp        <int> 92, 3128, 5184, 169725, 1620, 12543, 19170, 1856, 18761~
## $ TotExp         <int> 112, 3297, 5292, 172314, 1656, 13046, 19654, 1944, 1907~

We are going to build a simple linear regression of LifeExp as a function of TotExp

The Model

We are going here to apply directly a transformation to our data knowing that a simple model not transformed is not fitting the data well (See the blog link mentioned above).

Let transform the variable by applying some power and then build the linear model.

# Transform the data

df$LifeExp <- df$LifeExp**4.6
df$TotExp <- df$TotExp**0.06

df_lm <- lm(data = df, LifeExp~TotExp)
summary(df_lm)

## 
## Call:
## lm(formula = LifeExp ~ TotExp, data = df)
## 
## 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       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

Here’s what we’ve got by transforming the data:

R square = 0.7298; The model explains 72.98% of variability of the response data (LifeExp) around its mean. rse = 90490000; Larger deviation from the true regression line F-stat = 507.7 on 1 and 188 DF p-value = 2.2e-16; Reject the null hypothesis

Visualize the new model…

# Scatter plot
df %>%
  ggplot(aes(x = TotExp, y = LifeExp)) +
  geom_point()

Residuals plot…

qqnorm(df_lm$residuals)
qqline(df_lm$residuals)

The residuals plot shows there is more constant variability. Looking at different parameters(r square, F-Stat), the transformed model fits the data better.

Forecasting

Now we are going to forest life expectancy when sum of personal and government expenditures TotExp = 14

We know from the model that:

LifeExp = -736527910 + (620060216)*(TotExp)

Let do the predictions:

TotExp2 <- 1.5
LifeExp11 <- -736527910 + (620060216)*(TotExp2)

paste0("Life expectancy when  TotExp^.06 = 1.5 is ", LifeExp11)

## [1] "Life expectancy when  TotExp^.06 = 1.5 is 193562414"

TotExp4 <- 2.5
LifeExp22 = -736527910 + (620060216)*(TotExp4)

paste0("Life expectancy when  TotExp^.06 = 2.5 is ", LifeExp22)

## [1] "Life expectancy when  TotExp^.06 = 2.5 is 813622630"

Transforming back LifeExp = LifeExp^(1/4.6)…

LifeExp21 <- round(LifeExp11**(1/4.6), 0)

paste0("Life expectancy when  TotExp^.06 = 1.5 and  transorming back LifeExp is ", LifeExp21)

## [1] "Life expectancy when  TotExp^.06 = 1.5 and  transorming back LifeExp is 63"

LifeExp41 <- round(LifeExp22**(1/4.6), 0)

paste0("Life expectancy when  TotExp^.06 = 1.5 and  transorming back LifeExp is ", LifeExp41)

## [1] "Life expectancy when  TotExp^.06 = 1.5 and  transorming back LifeExp is 87"

And these predictions seems realistic based on WHO data.