1 Intro

World Happiness Report adalah survei penting mengenai keadaan kebahagiaan global. Laporan pertama diterbitkan pada tahun 2012, yang kedua pada tahun 2013, yang ketiga pada tahun 2015, dan yang keempat pada Pembaruan 2016. Laporan Kebahagiaan Dunia 2017, yang menilai 155 negara berdasarkan tingkat kebahagiaan mereka, dirilis di Perserikatan Bangsa-Bangsa dalam acara peringatan Hari Internasional Kebahagiaan pada tanggal 20 Maret. Laporan ini terus mendapatkan pengakuan global karena pemerintah, organisasi, dan masyarakat sipil semakin menggunakan indikator kebahagiaan untuk menginformasikan keputusan pembuatan kebijakan mereka. Para ahli terkemuka di berbagai bidang - ekonomi, psikologi, analisis survei, statistik nasional, kesehatan, kebijakan publik, dan lain-lain - menjelaskan bagaimana pengukuran kesejahteraan dapat digunakan secara efektif untuk menilai kemajuan suatu negara. Laporan-laporan ini meninjau keadaan kebahagiaan di dunia saat ini dan menunjukkan bagaimana ilmu kebahagiaan yang baru menjelaskan variasi kebahagiaan secara pribadi dan nasional.

1.1 Project’s Goal

Dalam projek ini, data yang digunakan adalah laporan index kebahagiaan pada tahun 2019 pada 154 negara di dunia. Projek ini bertujuan untuk melihat seberapa besar GDP, Social Support, Angka Harapan hidup, freedon of choices, generocity, dan persepsi terhadap korupsi berpengaruh terhadap index kebahagiaan di suatu negara. Oleh karena itu, metode yang digunakan dalam projek ini adalah metode linear regression untuk menilai korelasi hubungan antar beberapa variabel prediktor terhadap variabel target.

2 Data Preparation

#Import library
library(dplyr)
library(ggplot2)
library(lubridate)
library(MLmetrics)

#read data
happiness <- read.csv("2019.csv")
head(happiness)

Deskripsi Kolom:

GDP.per.capita (GDP per capita): Kolom ini mencerminkan nilai-nilai PDB (Produk Domestik Bruto) per kapita di berbagai negara.
Social.support (Dukungan sosial): Kolom ini menggambarkan tingkat dukungan sosial yang diterima oleh individu dalam suatu negara.
Healthy.life.expectancy (Harapan hidup sehat): Kolom ini mencerminkan harapan hidup sehat di suatu negara.
Freedom.to.make.life.choices (Kebebasan dalam membuat pilihan hidup): Kolom ini mencerminkan tingkat kebebasan yang dimiliki oleh individu dalam membuat keputusan hidup mereka.
Generosity (Kegenerosan): Kolom ini mencerminkan tingkat kebaikan hati dan kecenderungan untuk memberi individu dalam suatu negara.
Perceptions.of.corruption (Persepsi tentang korupsi): Kolom ini mencerminkan persepsi individu tentang tingkat korupsi di suatu negara.

#Check datatypes
glimpse(happiness)

## Rows: 156
## Columns: 9
## $ Overall.rank                 <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13…
## $ Country.or.region            <chr> "Finland", "Denmark", "Norway", "Iceland"…
## $ Score                        <dbl> 7.769, 7.600, 7.554, 7.494, 7.488, 7.480,…
## $ GDP.per.capita               <dbl> 1.340, 1.383, 1.488, 1.380, 1.396, 1.452,…
## $ Social.support               <dbl> 1.587, 1.573, 1.582, 1.624, 1.522, 1.526,…
## $ Healthy.life.expectancy      <dbl> 0.986, 0.996, 1.028, 1.026, 0.999, 1.052,…
## $ Freedom.to.make.life.choices <dbl> 0.596, 0.592, 0.603, 0.591, 0.557, 0.572,…
## $ Generosity                   <dbl> 0.153, 0.252, 0.271, 0.354, 0.322, 0.263,…
## $ Perceptions.of.corruption    <dbl> 0.393, 0.410, 0.341, 0.118, 0.298, 0.343,…

Nampaknya semua data telah tersimpan ke dalam tipe yang sudah tepat. Namun, karena kita ingin melakukan prediksi terhadap Score sehingga kolom Country.or.region dan Overall.rank akan dihapus karena tidak akan dimasukkan sebagai variabel prediktor dalam analisis ini.

#Select data
happiness_data <- happiness %>% 
  select(-Overall.rank, -Country.or.region)
head(happiness_data)

Langkah selanjutnya adalah mengecek apakah terdapat kolom yang memiliki nilai NULL sehingga kita bisa menentukan tindakan yang akan dilakukan jika hal sedemikin terdapat pada data yang kita miliki.

#Check if there is NULL columns
colSums(is.na(happiness_data))

##                        Score               GDP.per.capita 
##                            0                            0 
##               Social.support      Healthy.life.expectancy 
##                            0                            0 
## Freedom.to.make.life.choices                   Generosity 
##                            0                            0 
##    Perceptions.of.corruption 
##                            0

Tidak ada data yang NULL sehingga kita bisa melanjutkan ke tahap selanjutnya yaitu Exploratory Data Analysis (EDA).

3 Exploratory Data Analysis (EDA)

Cek Persebaran data dengan menggunakan fungsi boxplot():

#Check data distrbution
boxplot(happiness_data)

Meskipun terdapat beberapa outlier pada data diatas, namun jumlahnya sangat kecil sehingga kita bisa mengabaikannya dan lanjut melakukan pengecekan pada korelasi antara variabel prediktor dengan variabel target menggunakan fungsi ggcorr() dari library GGALLY:

library(GGally)
ggcorr(happiness_data, label = T, label_size = 2.9, hjust = 1, layout.exp = 2)

Dari hasil visualisasi diatas, kita bisa mengasumsikan bahwa variable prediktor yang memiliki korelasi yang kuat dengan variabel target adalah variabel dengan nilai korelasi >0.05, sehingga dapat diambil beberapa variabel, yaitu GDP per Kapita, Social Support, Angka harapan hidup, dan Freedom to make choices.

4 Train-test split

Sebelum kita membuat model, kita perlu membagi data menjadi Data Train dan Data Test. Kita akan menggunakan Data Train untuk melatih model regresi linear. Data Test akan digunakan sebagai pembanding dan melihat apakah model terlalu overfit dan tidak dapat memprediksi data baru yang belum pernah dilihat selama fase pelatihan. Kita akan menggunakan 80% dari happiness_data sebagai Data Train dan sisanya sebagai Data Test.

set.seed(123)
#splitting data
index <- sample(x = nrow(happiness_data),
                size = nrow(happiness_data)*0.8)
happiness_train <- happiness_data[index, ]
happiness_test <- happiness_data[-index, ]

5 Model Fitting and Future Selection

Sekarang kita akan membuat sebuah model regresi dengan variabel target adalah Score:

5.1 Model dengan semua prediktor

#Building model with all predictor variable
happiness_all <- lm(formula = Score ~.,
                    data = happiness_train)
summary(happiness_all)

## 
## Call:
## lm(formula = Score ~ ., data = happiness_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.76924 -0.36065  0.07753  0.38675  1.12317 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   1.89423    0.23099   8.200 3.49e-13 ***
## GDP.per.capita                0.92691    0.24917   3.720 0.000307 ***
## Social.support                1.01779    0.26753   3.804 0.000228 ***
## Healthy.life.expectancy       0.90900    0.36597   2.484 0.014415 *  
## Freedom.to.make.life.choices  1.57301    0.44380   3.544 0.000567 ***
## Generosity                    0.80114    0.57719   1.388 0.167775    
## Perceptions.of.corruption     0.07922    0.65034   0.122 0.903254    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5512 on 117 degrees of freedom
## Multiple R-squared:  0.7633, Adjusted R-squared:  0.7511 
## F-statistic: 62.87 on 6 and 117 DF,  p-value: < 2.2e-16

5.2 Model dengan beberapa prediktor (Selection)

#Building model with Predictors based on correlation

happiness_selection <- lm(formula = Score ~ GDP.per.capita + Social.support + Healthy.life.expectancy + Freedom.to.make.life.choices,
                          data = happiness_train)
summary(happiness_selection)

## 
## Call:
## lm(formula = Score ~ GDP.per.capita + Social.support + Healthy.life.expectancy + 
##     Freedom.to.make.life.choices, data = happiness_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.91855 -0.30804  0.04522  0.42198  1.09834 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    2.0082     0.2179   9.215 1.35e-15 ***
## GDP.per.capita                 0.9257     0.2445   3.787 0.000241 ***
## Social.support                 0.9858     0.2595   3.799 0.000230 ***
## Healthy.life.expectancy        0.8972     0.3658   2.453 0.015620 *  
## Freedom.to.make.life.choices   1.8094     0.3916   4.621 9.77e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5517 on 119 degrees of freedom
## Multiple R-squared:  0.7588, Adjusted R-squared:  0.7507 
## F-statistic: 93.62 on 4 and 119 DF,  p-value: < 2.2e-16

5.3 step-wise (Backward)

happiness_backward <- step(object = happiness_all,
                          direction = "backward")

## Start:  AIC=-140.92
## Score ~ GDP.per.capita + Social.support + Healthy.life.expectancy + 
##     Freedom.to.make.life.choices + Generosity + Perceptions.of.corruption
## 
##                                Df Sum of Sq    RSS     AIC
## - Perceptions.of.corruption     1    0.0045 35.554 -142.90
## <none>                                      35.550 -140.92
## - Generosity                    1    0.5854 36.135 -140.89
## - Healthy.life.expectancy       1    1.8745 37.424 -136.55
## - Freedom.to.make.life.choices  1    3.8172 39.367 -130.27
## - GDP.per.capita                1    4.2049 39.755 -129.06
## - Social.support                1    4.3977 39.948 -128.46
## 
## Step:  AIC=-142.9
## Score ~ GDP.per.capita + Social.support + Healthy.life.expectancy + 
##     Freedom.to.make.life.choices + Generosity
## 
##                                Df Sum of Sq    RSS     AIC
## <none>                                      35.554 -142.90
## - Generosity                    1    0.6599 36.214 -142.62
## - Healthy.life.expectancy       1    1.8864 37.441 -138.49
## - Freedom.to.make.life.choices  1    4.3941 39.948 -130.45
## - GDP.per.capita                1    4.4297 39.984 -130.34
## - Social.support                1    4.5926 40.147 -129.84

summary(happiness_backward)

## 
## Call:
## lm(formula = Score ~ GDP.per.capita + Social.support + Healthy.life.expectancy + 
##     Freedom.to.make.life.choices + Generosity, data = happiness_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.76896 -0.35863  0.07999  0.39185  1.11975 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    1.8948     0.2300   8.239 2.71e-13 ***
## GDP.per.capita                 0.9329     0.2433   3.834 0.000204 ***
## Social.support                 1.0100     0.2587   3.904 0.000158 ***
## Healthy.life.expectancy        0.9110     0.3641   2.502 0.013713 *  
## Freedom.to.make.life.choices   1.5910     0.4166   3.819 0.000215 ***
## Generosity                     0.8199     0.5540   1.480 0.141556    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5489 on 118 degrees of freedom
## Multiple R-squared:  0.7632, Adjusted R-squared:  0.7532 
## F-statistic: 76.08 on 5 and 118 DF,  p-value: < 2.2e-16

5.4 Membandingkan Model

Kita telah membuat beberapa model diatas, oleh karena itu pada tahap kali ini kita akan coba membandingkan ketiga model tersebut dengan tujuan untuk mengambil model yang terbaik. Ada dua pembandingan yang dilakukan pada projek ini yaitu:

Berdasarkan Goodness of Fit
Berdasarkan nilai Root Mean Squared Error (RMSE)

5.4.1 Goodness of Fit

Tujuan: untuk menentukan seberapa baik model dalam menjelaskan variansi dari target variabel.

Bandingkan nilai R-Squared dari model-model yang telah dibuat

summary(happiness_all)$adj.r.squared

## [1] 0.7511315

summary(happiness_selection)$adj.r.squared

## [1] 0.7507409

summary(happiness_backward)$adj.r.squared

## [1] 0.7532092

💡 Kesimpulan: Berdasarkan nilai R-Squared, dapat diambil kesimpulan bahwa model terbaik yang merepresentasikan data yang kita gunakan adalah model happiness_backward yaitu sebesar 0.7532092 atau 75.32%.

5.5 RMSE

Tujuan : Mengukur sejauh mana perbedaan antara nilai-nilai yang diobservasi dan nilai-nilai yang diprediksi oleh model regresi.

Lakukan prediksi dari model yang telah dibuat terhadap Data Train yang kita miliki:

happiness_train$predict_all <- predict(object = happiness_all, newdata = happiness_train)
happiness_train$predict_selection <- predict(object = happiness_selection, newdata = happiness_train)
happiness_train$predict_backward <- predict(object = happiness_backward, newdata = happiness_train)
head(happiness_train)

Lakukan perhitungan RMSE menggunakan fungsi RMSE() dari package MLmetrics:

#RMSE for happiness_all
RMSE(y_pred = happiness_train$predict_all,
     y_true = happiness_train$Score)

## [1] 0.5354361

#RMSE FOR happiness_selection
RMSE(y_pred = happiness_train$predict_selection,
     y_true = happiness_train$Score)

## [1] 0.5404166

#RMSE for Happiness_backward
RMSE(y_pred = happiness_train$predict_backward,
     y_true = happiness_train$Score)

## [1] 0.53547

💡 Kesimpulan: Dari hasil perhitungan RMSE diatas, dapat dilihat bahwa model dengan seluruh prediktor (happiness_all) dengan nilai RMSE sebesar 0.5354361 merupakan paling kecil dibandingkan dengan kedua model lainnya. Oleh sebab itu, meskipun nilai Goodness of Fit model ini lebih kecil dibandingkan dengan model happiness_backward, kita akan menganggap model ini sebagai model terbaik karena dalam kasus ini nilai error lebih dikonsiderasikan dibandingkan dengan nilai Goodness of Fit dari model. Oleh sebab itu, model ini akan digunakan sebagai final model untuk melakukan prediksi terhadap data test untuk melihat kondisi apakah model final kita overfit atau underfit.

happiness_test$predict_score <- predict(object = happiness_all, newdata = happiness_test)
#RMSE
RMSE(y_pred = happiness_test$predict_score,
    y_true = happiness_test$Score)

## [1] 0.4961612

Kesimpulan: Berdasarkan hasil RMSE antara data train yaitu 0.5354361 dan pada data test yaitu 0.4961612, perbedaannya cukup kecil dan keduanya memiliki nilai yang rendah. Ini menunjukkan bahwa model yang telah dibuat relatif baik dalam memprediksi data train maupun data test. Sehingga dapat diambil kesimpulan bahwa model happiness_all optimum.

6 Model Interpretation

summary(happiness_all)

## 
## Call:
## lm(formula = Score ~ ., data = happiness_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.76924 -0.36065  0.07753  0.38675  1.12317 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   1.89423    0.23099   8.200 3.49e-13 ***
## GDP.per.capita                0.92691    0.24917   3.720 0.000307 ***
## Social.support                1.01779    0.26753   3.804 0.000228 ***
## Healthy.life.expectancy       0.90900    0.36597   2.484 0.014415 *  
## Freedom.to.make.life.choices  1.57301    0.44380   3.544 0.000567 ***
## Generosity                    0.80114    0.57719   1.388 0.167775    
## Perceptions.of.corruption     0.07922    0.65034   0.122 0.903254    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5512 on 117 degrees of freedom
## Multiple R-squared:  0.7633, Adjusted R-squared:  0.7511 
## F-statistic: 62.87 on 6 and 117 DF,  p-value: < 2.2e-16

Interpretation:

\[Score = 1.89423+ 0.92691*GDP.per.capita+ 1.01779*Social.support +\\ 0.90900*Healthy.life.expectancy+ 1.57301*Freedom.to.make.life.choices + \\ \\0.80114* Generosity+ 0.07922 * Perceptions.of.corruption\]

Intercept (Intersep): Estimasi intercept adalah 1.89423. Ini mengindikasikan bahwa jika semua variabel independen diatur ke nilai nol, prediksi Happiness Index akan memiliki nilai sekitar 1.89423. Dalam konteks ini, intercept mewakili pengaruh faktor-faktor lain yang tidak termasuk dalam model terhadap Happiness Index.
GDP per capita (PDB per kapita): Estimasi koefisien untuk variabel GDP per capita adalah 0.92691. Ini berarti bahwa setiap peningkatan satu unit dalam GDP per capita akan menyebabkan peningkatan sekitar 0.92691 dalam prediksi Happiness Index, dengan asumsi variabel lain tetap konstan.
Social support (Dukungan sosial): Estimasi koefisien untuk variabel social support adalah 1.01779. Artinya, setiap peningkatan satu unit dalam tingkat dukungan sosial akan berkontribusi sekitar 1.01779 dalam meningkatkan prediksi Happiness Index, dengan asumsi variabel lain tetap konstan.
Healthy life expectancy (Harapan hidup sehat): Estimasi koefisien untuk variabel healthy life expectancy adalah 0.90900. Ini menunjukkan bahwa setiap peningkatan satu unit dalam harapan hidup sehat akan berdampak sekitar 0.90900 dalam meningkatkan prediksi Happiness Index, dengan asumsi variabel lain tetap konstan.
Freedom to make life choices (Kebebasan dalam membuat pilihan hidup): Estimasi koefisien untuk variabel freedom to make life choices adalah 1.57301. Ini berarti bahwa setiap peningkatan satu unit dalam tingkat kebebasan dalam membuat pilihan hidup akan berdampak sekitar 1.57301 dalam meningkatkan prediksi Happiness Index, dengan asumsi variabel lain tetap konstan.
Generosity (Kegenerosan): Estimasi koefisien untuk variabel generosity adalah 0.80114. Namun, karena koefisien ini tidak signifikan secara statistik (p-value > 0.05), kita tidak dapat membuat kesimpulan yang kuat tentang pengaruh variabel ini terhadap prediksi Happiness Index.
Perceptions of corruption (Persepsi tentang korupsi): Estimasi koefisien untuk variabel perceptions of corruption adalah 0.07922. Karena koefisien ini juga tidak signifikan secara statistik (p-value > 0.05), kita tidak dapat menyimpulkan adanya pengaruh yang kuat dari variabel ini terhadap prediksi Happiness Index.

7 Asumsi Linear Regression

Sebagai salah satu model statistik, linear regression adalah model yang ketat asumsi. Berikut beberapa asumsi yang harus dicek untuk memastikan apakah model yang kita buat dianggap sebagai Best Linear Unbiased Estimator (BLUE) model, yaitu model yang dapat memprediksi data baru secara konsisten.

Asumsi model linear regression:

Linearity
Normality of Residuals
Homoscedasticity of Residuals
No Multicollinearity

7.1 1. Linearity

plot(happiness_all, which = 1)

💡 Kesimpulan: Nilai residual tersebar secara acak diantara -0.5 dan 0.5, artinya model yang kita miliki memenuhi asumsi linieriats.

7.2 2. Normality of Residuals

Visualisasi histogram residual menggunakan fungsi hist()

hist(happiness_all$residuals)

Uji statistik dengan shapiro.test()

Shapiro-Wilk hypothesis test:

H0: error berdistribusi normal
H1: error TIDAK berdistribusi normal

H0 ditolak jika p-values < 0.05 (alpha)

shapiro.test(happiness_all$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  happiness_all$residuals
## W = 0.97999, p-value = 0.06273

💡 Kesimpulan: p-value = 0.06273 > 0.05, artinya residual data berdistribusi normal.

Membuat plot plot

plot(happiness_all, which = 2)

Berdasarkan pada plot diatas, nilai residual dapat dikatakn berdistribusi normal sebab titik data terdapat pada strip garis.

7.3 3. Homoscedasticity of Residuals

Diharapkan error yang dihasilkan oleh model menyebar secara acak atau dengan variasi konstan. Apabila divisualisasikan maka error tidak berpola. Kondisi ini disebut juga sebagai homoscedasticity.

Visualisasi scatter plot: fitted.values vs residuals

plot(x = happiness_all$fitted.values,
     y = happiness_all$residuals)
abline(h = 0, col = "red")

Berdasarkan hasil plot diatas, dapat diambil kesimpulan bahwa terjadi homoscedasticity pada data sebab titik nilai residual dan fitted.value dari model tidak menunjukkan pola tertentu. Untuk lebih membuktikan asumsi ini maka akan dilakukan bptest() dari library lmtest sebagai berikut:

Breusch-Pagan hypothesis test:

H0: error menyebar konstan atau homoscedasticity
H1: error menyebar TIDAK konstan atau heteroscedasticity

tolak H0 jika nilai p-value < 0.05 (alpha)

library(lmtest)
bptest(happiness_all)

## 
##  studentized Breusch-Pagan test
## 
## data:  happiness_all
## BP = 12.006, df = 6, p-value = 0.06183

💡 Kesimpulan: p-value = 0.06183 > 0.05 artinya error menyebar konstan atau homoscedasticity.

7.4 4. No multicolinearity

Uji VIF (Variance Inflation Factor) dengan fungsi vif() dari package car:

nilai VIF > 10: terjadi multicollinearity pada model
nilai VIF < 10: tidak terjadi multicollinearity pada model

library(car)
vif(happiness_all)

##               GDP.per.capita               Social.support 
##                     4.034052                     2.751096 
##      Healthy.life.expectancy Freedom.to.make.life.choices 
##                     3.436382                     1.601491 
##                   Generosity    Perceptions.of.corruption 
##                     1.232681                     1.427068

💡 Kesimpulan: Berdasarkan hasil uji VIF diatas, tidak terdapat satupun variabel prediktor yang memiliki nilai VIF > 10, artinya antara satu variabel independen dengan yang lainnya tidak terdapat korelasi yang kuat.

8 Kesimpulan

Variabel prediktor yang paling baik dalam menjelaskan distribusi dari score happiness index adalah GDP per Kapita, Social Support, life expectancy, freedom to make choices, Generosity, dan percepstion of corruption. Meskipun secara statistik Generosity, dan percepstion of corruption tidak menunjukkan pengaruh yang signifikan terhadap score, namun variabel ini tetap menjadi salah satu konsiderasi dalam menjelaskan score dari happiness index di setiap negara.
Nilai R-Squared dari model yang dimiliki cukup tinggi dengan 75.11% dari variabel bisa menjelaskan variasi dari score happiness index.
Nilai RMSE pada data train adalah sebesar 0.5354361 dan pada data test yaitu 0.4961612, hal ini menunjukkan angka serta range yang cukup kecil menandakan bahwa model kita optimum.

Happiness Index Score Prediction Using Linear Regression

Miftahul Labiib Syam

2023-06-10