Berikut adalah analisis mengenai data World Happiness Report 2019 yang disediakan oleh Prompt Cloud. Data hasil survey ini yang menunjukan pengaruh beragam aspek kehidupan terhadap tingkat kebahagiaan masyarakat.

Menggunakan linear regression analysis, akan dibuat suatu model untuk memprediksi faktor-faktor yang menjadikan nilai harapan hidup menjadi suatu aspek yang penting dalam kebahagiaan masyarakat.

1. Load Dataset

happy <- read.csv("data_input/happiness.csv")

2. Explanatory Data Analysis

str(happy)

## 'data.frame':    156 obs. of  9 variables:
##  $ Country..region.: Factor w/ 156 levels "Afghanistan",..: 44 37 106 58 100 134 133 101 24 7 ...
##  $ Ladder          : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Positive.affect : int  41 24 16 3 12 44 34 22 18 64 ...
##  $ Social.support  : int  2 4 3 1 15 13 25 5 20 31 ...
##  $ Freedom         : int  5 6 3 7 19 11 10 8 9 26 ...
##  $ Corruption      : int  4 3 8 45 12 7 6 5 11 19 ...
##  $ Generosity      : int  47 22 11 3 7 16 17 8 14 25 ...
##  $ GDP             : int  22 14 7 15 12 8 13 26 19 16 ...
##  $ Life.exp        : int  27 23 12 13 18 4 17 14 8 15 ...

Berikut adalah keterangan dari setiap kolom:

-Country: negara

-Ladder: peringkat tingkat kebahagiaan

-Positive.affect: tingkat emosi positif

-Social.support: pengaruh dukungan sosial terhadap kebahagiaan

-Freedom: pengaruh kebebasan individu terhadap kebahagiaan

-Corruption: pengaruh kejadian korupsi terhadap kebahagiaan

-Generosity: pengaruh kondisi saling membantu antar individu terhadap kebahagiaan

-GDP: pengaruh Gross Domestic Product per Kapita terhadap kebahagiaan

-Life.exp: pengaruh nilai harapan hidup terhadap kebahagiaan

# subsetting data
h <- happy[,-c(1:2)]
# menguji keberadaan missing value
anyNA(h)

## [1] TRUE

table(is.na(h))

## 
## FALSE  TRUE 
##  1070    22

Karena proporsi missing value (NA) dari data hanya sebanyak 2% maka baris yang mengandung missing value dapat dihapus.

# menghapus baris yang terdapat 'NA'
h <- na.omit(h)
# melihat korelasi antar variabel
library(GGally)
ggcorr(h, label = TRUE, label_size = 2.9, hjust = 1, layout.exp = 2)

Pada grafik korelasi, terlihat bahwa semua variabel memiliki pengaruh positif terhadap Life.exp dimana faktor GDP memiliki korelasi positif yang paling tinggi dibandingkan faktor-faktor lain.

Berikut adalah distribusi nilai dari masing-masing variabel.

boxplot(h)

Berdasarkan visualisasi boxplot, tidak ditemukan outlier pada masing-masing kolom variabel dan data berdistribusi normal, sehingga data tersebut sudah dapat dianalisis lebih lanjut.

3. Pembuatan Model Regresi Linear

Selanjutnya dapat dibuat model regresi linear dengan variabel prediktor GDP karena variabel tersebut memiliki korelasi positif tertinggi terhadap variabel target Life.exp.

m <- lm(Life.exp ~ GDP, h)
m

## 
## Call:
## lm(formula = Life.exp ~ GDP, data = h)
## 
## Coefficients:
## (Intercept)          GDP  
##      5.4565       0.8862

summary(m)

## 
## Call:
## lm(formula = Life.exp ~ GDP, data = h)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -49.507 -16.586  -0.902  13.416  53.868 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.45652    3.78381   1.442    0.152    
## GDP          0.88619    0.04202  21.091   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.48 on 138 degrees of freedom
## Multiple R-squared:  0.7632, Adjusted R-squared:  0.7615 
## F-statistic: 444.8 on 1 and 138 DF,  p-value: < 2.2e-16

plot(h$GDP, h$Life.exp)
abline(m$coefficients[1],m$coefficients[2])

Dapat dilihat bahwa adjusted R-squared memiliki nilai 0.7615

Selanjutnya akan dicoba pemilihan variabel prediktor secara automatis menggunakan step-wise regression dengan metode backward elimination.

m1 <- lm(Life.exp ~ ., h)
step(m1, direction = "backward")

## Start:  AIC=861.07
## Life.exp ~ Positive.affect + Social.support + Freedom + Corruption + 
##     Generosity + GDP
## 
##                   Df Sum of Sq    RSS    AIC
## - Positive.affect  1         6  59412 859.08
## - Generosity       1       176  59581 859.48
## - Corruption       1       589  59995 860.45
## <none>                          59406 861.07
## - Social.support   1       918  60324 861.22
## - Freedom          1       958  60364 861.31
## - GDP              1     50845 110251 945.64
## 
## Step:  AIC=859.08
## Life.exp ~ Social.support + Freedom + Corruption + Generosity + 
##     GDP
## 
##                  Df Sum of Sq    RSS    AIC
## - Generosity      1       171  59582 857.49
## - Corruption      1       583  59995 858.45
## <none>                         59412 859.08
## - Social.support  1       917  60329 859.23
## - Freedom         1      1289  60701 860.09
## - GDP             1     51457 110869 944.42
## 
## Step:  AIC=857.49
## Life.exp ~ Social.support + Freedom + Corruption + GDP
## 
##                  Df Sum of Sq    RSS    AIC
## - Corruption      1       528  60111 856.72
## <none>                         59582 857.49
## - Social.support  1       956  60538 857.71
## - Freedom         1      1769  61351 859.58
## - GDP             1     54513 114096 946.44
## 
## Step:  AIC=856.72
## Life.exp ~ Social.support + Freedom + GDP
## 
##                  Df Sum of Sq    RSS    AIC
## <none>                         60111 856.72
## - Social.support  1      1168  61279 857.42
## - Freedom         1      1292  61402 857.70
## - GDP             1     54013 114124 944.48

## 
## Call:
## lm(formula = Life.exp ~ Social.support + Freedom + GDP, data = h)
## 
## Coefficients:
##    (Intercept)  Social.support         Freedom             GDP  
##        0.65753         0.11127         0.07697         0.76100

summary(lm(formula = Life.exp ~ Social.support + Freedom + GDP, data = h))

## 
## Call:
## lm(formula = Life.exp ~ Social.support + Freedom + GDP, data = h)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -43.794 -15.351   0.688  11.391  56.272 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     0.65753    4.21863   0.156   0.8764    
## Social.support  0.11127    0.06845   1.626   0.1063    
## Freedom         0.07697    0.04502   1.710   0.0896 .  
## GDP             0.76100    0.06884  11.055   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.02 on 136 degrees of freedom
## Multiple R-squared:  0.7764, Adjusted R-squared:  0.7715 
## F-statistic: 157.4 on 3 and 136 DF,  p-value: < 2.2e-16

# membuat formula model regresi yang baru menjadi objek
m1  <- lm(formula = Life.exp ~ Social.support + Freedom + GDP, data = h)
m1

## 
## Call:
## lm(formula = Life.exp ~ Social.support + Freedom + GDP, data = h)
## 
## Coefficients:
##    (Intercept)  Social.support         Freedom             GDP  
##        0.65753         0.11127         0.07697         0.76100

Metode step-wise regression ini akan menghasilkan formula optimum berdasarkan nilai AIC yang terendah, dimana semakin rendah nilai AIC tersebut, maka nilai observasi yang tidak tertangkap semakin kecil.

Bila dibandingkan dengan model awal yang hanya menggunakan variabel GDP, model regresi yang menggunakan variabel prediktor Social.support dan Freedom dan GDP memiliki adjusted R-squared 0.7715. lebih tinggi dibandingkan model sebelumnya yaitu 0.7615.

Kandidat Model:

Life.exp = 0.11127(Social.support) + 0.07697(Freedom) + 0.76100(GDP) + 0.65753

Life.exp = 0.8862(GDP) + 5.4565

4. Prediksi Model & Error

Akan dicoba prediksi nilai GDP berdasarkan nilai variabel prediktor, dan hasilnya akan dibandingkan dengan data aktual yang kita miliki.

# prediksi nilai GDP berdasarkan model `m`
predict(m, data.frame(GDP = 22), interval = "confidence", level = 0.95)

##        fit      lwr      upr
## 1 24.95281 19.00976 30.89586

# prediksi nilai GDP berdasarkan model `m1`
predict(m1, data.frame(Social.support = 2, Freedom = 5, GDP = 22), interval = "confidence", level = 0.95)

##       fit      lwr      upr
## 1 18.0069 10.40975 25.60405

# menghitung error m
sqrt((27-24.95281)^2)

## [1] 2.04719

# menghitung error m1
sqrt((27-18.0069)^2)

## [1] 8.9931

Analisis penghitungan RSE menghasilkan model m memiliki nilai RSE yang lebih kecil sehingga lebih baik dibandingkan model m1.

5. Evaluasi Model

Normalitas

hist(m$residuals, breaks = 20)

shapiro.test(m$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  m$residuals
## W = 0.98511, p-value = 0.1338

hist(m1$residuals, breaks = 20)

shapiro.test(m1$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  m1$residuals
## W = 0.9884, p-value = 0.293

Untuk kedua model, P-value > 0.05 sehingga H0 diterima. Hal ini juga berarti residual menyebar normal agar model kita memiliki error disekitar mean-nya.

Heteroscedasticity

library(lmtest)
plot(h$GDP, m$residuals)
abline(h = 0, col = "red")

bptest(m)

## 
##  studentized Breusch-Pagan test
## 
## data:  m
## BP = 0.060444, df = 1, p-value = 0.8058

plot(h$GDP, m1$residuals)
abline(h = 0, col = "red")

bptest(m1)

## 
##  studentized Breusch-Pagan test
## 
## data:  m1
## BP = 3.2444, df = 3, p-value = 0.3555

Untuk kedua model, P-value > 0.05 sehingga H0 diterima. Hal ini juga berarti residual tidak memiliki pola (Heteroscedasticity) dimana semua pola yang ada sudah berhasil ditangkap oleh model yang dibuat.

Variance Inflation Factor (Multicollinearity)

library(car)

## Loading required package: carData

vif(m1)

## Social.support        Freedom            GDP 
##       3.092686       1.297157       2.801478

# untuk model m yang hanya menggunakan satu variabel prediktor, tidak dapat digunakan `vif()` untuk analisis Multicollinearity

Tidak ada nilai sama dengan atau lebih dari 10 sehingga tidak ditemukan Multicollinearity antar variabel (antar variabel prediktor saling independen).

Berdasarkan hasil analisis, kedua model memiliki kriteria yang baik sebagai model linear regression. Kemudian, bila dibandingkan RSE antara kedua model, model m memberikan nilai RSE yang lebih rendah. Oleh karena itu model m dipilih sebagai model yang lebih baik.

Simpulan dan Saran

Model m yang didapatkan memiliki R-square 0.7615 dan memiliki RSE sebesar 2.04719. Selain itu setelah dilakukan uji analisis, model memiliki kriteria yang sudah baik.

Berdasarkan model ini, nilai harapan hidup berkorelasi positif dengan nilai GDP. Dalam kata lain, masyarakat yang menganggap besarnya nilai GDP akan meningkatkan kebahagiaan mereka, juga cenderung menganggap nilai harapan hidup yang tinggi akan meningkatkan kebahagiaan mereka.

World Happiness Report 2019

Nabiilah Ardini Fauziyyah

7/7/2019