Attaching package: 'dplyr'The following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, unionlinear_regression
Implementing Linear Regression Model in R
Today we will be implementing a simple linear regression model on the World Bank and Our World in Data data set of GDP per capita and Happiness report across countries around the world.
Install and load packages
Install the packages if not yet installed.
install.packages("ggplot2")
install.packages("dplyr")
install.packages("ggpubbr")Load the packages so that we can use the needed functions
Load the data set
We will be using the data from the World Happiness Report (2012-2024) – with major processing by Our World in Data. “Self-reported life satisfaction – WHR” [dataset]. Wellbeing Research Centre, “World Happiness Report 2024”; Various sources, “Population” [original data].
You can download the csv online into your project’s directory, and import the data set it in your current project for it to be directly accessible with the following call.
data <- read_csv("gdp_happiness.csv", show_col_types = FALSE)New names:
• `time` -> `time...7`
• `time` -> `time...8`
• `time` -> `time...9`#data <- gdp_happiness
happiness <- data$`Cantril ladder score`
gdp <- data$`GDP per capita, PPP (constant 2017 international $)`
# Display the first few rows
head(data)# A tibble: 6 × 9
Entity Code Year `Cantril ladder score` GDP per capita, PPP (constant…¹
<chr> <chr> <dbl> <dbl> <dbl>
1 Afghanistan AFG 2022 2.40 1516.
2 Albania ALB 2022 5.28 15492.
3 Algeria DZA 2022 5.33 11198.
4 Angola AGO 2022 3.79 7216.
5 Argentina ARG 2022 6.02 22461.
6 Armenia ARM 2022 5.34 16057.
# ℹ abbreviated name: ¹`GDP per capita, PPP (constant 2017 international $)`
# ℹ 4 more variables: `World regions according to OWID` <chr>, time...7 <dbl>,
# time...8 <dbl>, time...9 <dbl># Get summary statistics
summary(data) Entity Code Year Cantril ladder score
Length:150 Length:150 Min. :2022 Min. :2.404
Class :character Class :character 1st Qu.:2022 1st Qu.:4.568
Mode :character Mode :character Median :2022 Median :5.524
Mean :2022 Mean :5.417
3rd Qu.:2022 3rd Qu.:6.187
Max. :2022 Max. :7.804
NA's :1
GDP per capita, PPP (constant 2017 international $)
Min. : 711.4
1st Qu.: 5135.2
Median :13393.6
Mean :20400.1
3rd Qu.:33169.5
Max. :74917.7
World regions according to OWID time...7 time...8 time...9
Length:150 Min. :2023 Min. :2017 Min. :2017
Class :character 1st Qu.:2023 1st Qu.:2022 1st Qu.:2022
Mode :character Median :2023 Median :2022 Median :2022
Mean :2023 Mean :2022 Mean :2022
3rd Qu.:2023 3rd Qu.:2022 3rd Qu.:2022
Max. :2023 Max. :2022 Max. :2022
Check assumptions
Normality of the dependent variable
The Happiness report, measured by the Cantrill ladder score, is the dependent variable. Plotting it as a histogram we can see that it is roughly normally distributed
hist(happiness)
Linearity
We visually check if our data, the GDP per country and the happiness report, is linear through a scatter plot.
The scatter plot suggests a positive linear relationship between GDP per capita and happiness.
scatter_plot <- ggplot(data, aes(x = `GDP per capita, PPP (constant 2017 international $)`, y = `Cantril ladder score`)) +
geom_point() +
labs(title = "Scatter Plot of Happiness Score vs. GDP per Capita",
x = "GDP per Capita (PPP)",
y = "Happiness Score")
scatter_plot
Check for outliers
We can check for outliers for our GDP per capita data by plotting a box plot
boxplot(gdp)
boxplot(happiness)
Linear regression
Linear regression is one of the simplest algorithms for doing supervised learning. We will proceed with a simple linear regression since we only have a single independent variable (GDP per capita), predicting our dependent variable (Happiness, measured through the Cantril score).
Simple linear regression assumes that the statistical relationship between two continuous variables is approximately linear, for our variables \(GDP\) and \(H\), each data point must satisfy
\[ GDP_i = \beta_0 + \beta_1H_1 + \epsilon_i \]
for \(i=1,...,n\).
What this means is that if we want to find a linear model function for the \(GDP\) and \(H\), then we must minimize the errors from the data points and this would give us the linear model
\[ GDP = \alpha + \beta H \]
where \(\alpha\) and \(\beta\) are the constants that would would minimize the error for between our variables.
Measuring the strength of association
Pearson’s correlation coefficient is often used to quantify the strength of the linear association between two continuous variable.
Minimizing the square of errors to find the best fit line
Performing the linear regression
model <- lm(`Cantril ladder score` ~ `GDP per capita, PPP (constant 2017 international $)`, data = data)
summary(model)
Call:
lm(formula = `Cantril ladder score` ~ `GDP per capita, PPP (constant 2017 international $)`,
data = data)
Residuals:
Min 1Q Median 3Q Max
-2.12312 -0.40606 0.07336 0.45329 1.52890
Coefficients:
Estimate Std. Error
(Intercept) 4.455e+00 8.499e-02
`GDP per capita, PPP (constant 2017 international $)` 4.714e-05 3.072e-06
t value Pr(>|t|)
(Intercept) 52.42 <2e-16 ***
`GDP per capita, PPP (constant 2017 international $)` 15.34 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7031 on 148 degrees of freedom
Multiple R-squared: 0.614, Adjusted R-squared: 0.6114
F-statistic: 235.4 on 1 and 148 DF, p-value: < 2.2e-16The given estimates for the intercept and \(GDP\) are the constants for our linear model \(\alpha\) and \(\beta\) are. The provided p-value is also near-zero indicating the model fits the data well.
Visualize the linear model
scatter_plot <- scatter_plot + geom_smooth(method = "lm") + stat_regline_equation()
scatter_plot`geom_smooth()` using formula = 'y ~ x'
Checking for homoscedasticity
One assumption of linear regression is homoscedasticity. It is when the error is constant across the values of the dependent variable.
Homoscedasticity occurs when the variance in a data set is constant, making it easier to estimate the standard deviation and variance of a data set.
hist(model$residuals, main = "Histogram of Residuals", xlab = "Residuals")
qqnorm(model$residuals)
qqline(model$residuals)
plot(model, which = 1)
Exploration of the data set
References
boehmke, b., & greenwell, b.m. (2019). hands-on machine learning with r. chapman and hall/crc.
Understanding Diagnostic Plots for Linear Regression Analysis | UVA Library