Linear Regression - Mathematics Transformation

Overview

Here, we will demonstrate the power of mathematics transformation when performing linear regression.

We are going to use 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.

Data Exploration

Let first load the data, and take a look on how it looks like:

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

LifeExp~TotExp: Simple linear regression

It is good to visualize the regression first

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

Here, we are going to build the model and have a look on its summary

# Linear model

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

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

From the summary above, we realize that:

R square = 0.2577; The model explains 25.77% of variability of the response data (LifeExp) around its mean. rse = 9.371; spread is lower F-stat = 65.26 on 1 and 188 DF p-value = 7.714e-14; Reject the null hypothesis

Let now check on linear model assumptions and verify if a linear model is suitable for this data set.

Assumptions

Let take a look on each variable and its residuals…

1. Normality

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

Residuals are not nearly normal.

2. Linearity

The scatter plot will have to check the linearity…

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

Looking at the above plot, we can say at first glace that the data don’t look linear.

3. Constant variability

ggplot(data = df_lm, aes(x = .fitted, y = .resid)) +
  geom_point()+
  
  geom_hline(yintercept = 0, linetype = "dashed") +
  xlab("Fitted values") +
  ylab("Residuals")

The variability on the residuals is not constant.

At the end, we can say that the assumptions of simple linear regression are not met. That said, it is a good time to do some transformation. We are going to add power on the variables to transfrom the data.

Model transformed

# Transform the data

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

# Re-run the model

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

## 
## 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_lm2$residuals)
qqline(df_lm2$residuals)

Take Away

It is important to think of mathematics transformations when a linear model can’t fit the data. Here, I choose some value arbitrary to power the variables but you can use PowerTransform function of “car” library for not guessing on which power you need for your variables.

As we can see, the transformed variables has significantly improve the model. The residuals plot shows there is more constant variability and also comparing different parameters(r square, F-Stat), that’s it, the transformed model fits the data better.