Membuat Regresi Linier Berganda Data Covid DKI Jakarta dan Google Mobility Index Pada Bulan April 2020 Khusus Pasien Yang Meninggal
Hafidz Nurrois
2022-03-28
1. Regresi Linier Berganda
Regresi Linear Berganda adalah model regresi linear dengan melibatkan lebih dari satu variable bebas atau predictor. Dalam bahasa inggris, istilah ini disebut dengan multiple linear regression.
1.1. Pengambilan Dan Hasil Data
Disini data diambil dari Excel. Excel adalah salah satu jenis file eksternal yang sering digunakan untuk menyimpan data. Kita dapat menggunakan package {readxl} dengan fungsi read_excel() untuk import data dari file Excel. Argumen path = adalah lokasi dan nama file Excel yang akan kita gunakan.
library(readxl)
dataKematianApril <- read_excel(path = "KematianCovidApril.xlsx")
dataKematianApril## # A tibble: 30 x 8
## Tanggal Meninggal retail_and_recreation_percent~ grocery_and_pha~
## <dttm> <dbl> <dbl> <dbl>
## 1 2020-04-01 00:00:00 1219 -24 -2
## 2 2020-04-02 00:00:00 1237 -25 -6
## 3 2020-04-03 00:00:00 1253 -25 -6
## 4 2020-04-04 00:00:00 1260 -22 -5
## 5 2020-04-05 00:00:00 1277 -25 -4
## 6 2020-04-06 00:00:00 1289 -30 -10
## 7 2020-04-07 00:00:00 1318 -25 -8
## 8 2020-04-08 00:00:00 1330 -26 -9
## 9 2020-04-04 00:00:00 1347 -26 -10
## 10 2020-04-10 00:00:00 1365 -26 -10
## # ... with 20 more rows, and 4 more variables:
## # parks_percent_change_from_baseline <dbl>,
## # transit_stations_percent_change_from_baseline <dbl>,
## # workplaces_percent_change_from_baseline <dbl>,
## # residential_percent_change_from_baseline <dbl>1.2. Summary Data
summary() atau ringkasan data digunakan untuk mencari nilai statistik lima serangkai (minimum, Q1, Q2 atau median, Q3, maksimum).
summary(datacovidmobility
summary(dataKematianApril)## Tanggal Meninggal
## Min. :2020-04-01 00:00:00 Min. :1219
## 1st Qu.:2020-04-07 06:00:00 1st Qu.:1334
## Median :2020-04-15 12:00:00 Median :1483
## Mean :2020-04-15 08:00:00 Mean :1479
## 3rd Qu.:2020-04-22 18:00:00 3rd Qu.:1644
## Max. :2020-04-30 00:00:00 Max. :1731
## retail_and_recreation_percent_change_from_baseline
## Min. :-53.00
## 1st Qu.:-41.75
## Median :-38.50
## Mean :-35.00
## 3rd Qu.:-26.00
## Max. :-22.00
## grocery_and_pharmacy_percent_change_from_baseline
## Min. :-30.00
## 1st Qu.:-18.75
## Median :-12.00
## Mean :-13.47
## 3rd Qu.: -8.00
## Max. : -2.00
## parks_percent_change_from_baseline
## Min. :-76.00
## 1st Qu.:-65.00
## Median :-61.50
## Mean :-61.83
## 3rd Qu.:-59.00
## Max. :-51.00
## transit_stations_percent_change_from_baseline
## Min. :-55.0
## 1st Qu.:-52.0
## Median :-47.0
## Mean :-45.6
## 3rd Qu.:-41.0
## Max. :-34.0
## workplaces_percent_change_from_baseline
## Min. :-39.00
## 1st Qu.:-37.00
## Median :-32.00
## Mean :-30.27
## 3rd Qu.:-23.25
## Max. :-13.00
## residential_percent_change_from_baseline
## Min. : 8.0
## 1st Qu.:13.0
## Median :14.0
## Mean :14.9
## 3rd Qu.:18.0
## Max. :19.01.3. Matriks Dengan Fungsi Pairs
Untuk membuat matriks scatterplot kita hanya perlu memasukkan objek datacovidmobility kedalam fungsi pairs(). Berikut adalah sintaks yang digunakan dan output yang dihasilkan.
pairs(dataKematianApril)
Kita juga dapat melakukan drop terhadap panel bawah grafik tersebut, yaitu dengan memasukkan argumen lower.panel=NULL.
pairs(dataKematianApril, lower.panel=NULL)
1.4. Visualisasi Data Dengan Fungsi Plot
plot(dataKematianApril$Meninggal ~ dataKematianApril$retail_and_recreation_percent_change_from_baseline, data = dataKematianApril)
Memvisualisasikan Data kematian sebagai variabel Y dan Google Mobility Index sebagai variabel X
plot(dataKematianApril$Meninggal ~ dataKematianApril$
retail_and_recreation_percent_change_from_baseline+ dataKematianApril$
grocery_and_pharmacy_percent_change_from_baseline+dataKematianApril$
parks_percent_change_from_baseline+dataKematianApril$
transit_stations_percent_change_from_baseline+dataKematianApril$
workplaces_percent_change_from_baseline+dataKematianApril$
residential_percent_change_from_baseline, data = dataKematianApril)
1.5. Korelasi Antar Variabel
Korelasi merupakan keterhubungan antar variabel. Untuk mengukur seberapa jauh hubungan antara satu variabel dengan variabel yang lain kita dapat menggunakan fungsi cor(). ##### a) Korelasi Antara Variabel Y Dengan X1
cor(dataKematianApril$Meninggal,dataKematianApril$
retail_and_recreation_percent_change_from_baseline)## [1] -0.8508943b) Korelasi Antara Variabel Y Dengan X2
cor(dataKematianApril$Meninggal,dataKematianApril$
grocery_and_pharmacy_percent_change_from_baseline)## [1] -0.5516343c) Korelasi Antara Variabel Y Dengan X3
cor(dataKematianApril$Meninggal,dataKematianApril$
parks_percent_change_from_baseline)## [1] -0.5532132d) Korelasi Antara Variabel Y Dengan X4
cor(dataKematianApril$Meninggal,dataKematianApril$
transit_stations_percent_change_from_baseline)## [1] -0.6831169e) Korelasi Antara Variabel Y Dengan X5
cor(dataKematianApril$Meninggal,dataKematianApril$
workplaces_percent_change_from_baseline)## [1] -0.2632221f) Korelasi Antara Variabel Y Dengan X6
cor(dataKematianApril$Meninggal,dataKematianApril$
residential_percent_change_from_baseline)## [1] 0.5879545Dari hasil seluruh output di atas dapat disimpulkan bahwa untuk tingkat keterhubungan antara variabel y dengan x1, x2, x3, x4, dan x5 memiliki hubungan karena nilai yang dihasilkan berjumlah lebih dari 0.3. Sedangkan untuk tingkat keterhubungan antara variabel y dan x6 tidak sangat terhubung karena nilai yang dihasilkan berjumlah kurang dari 0.
1.6. Melakukan Permodelan
Berikut cara melakukan permodelan regresi linier berganda.
modelData <- lm(dataKematianApril$Meninggal ~ dataKematianApril$Tanggal, data = dataKematianApril)Selanjutnya dengan model yang telah dibuat di atas, kita akan menggunakan fungsi summary() untuk menjelaskan atau mereview hasil dari model tersebut. Dan dengan ringkasan summary(model) kita dapat melihat informasi terperinci tentang kinerja dan koefisian model.
summary(modelData)##
## Call:
## lm(formula = dataKematianApril$Meninggal ~ dataKematianApril$Tanggal,
## data = dataKematianApril)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.273 -13.785 -3.993 5.850 77.201
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.383e+05 7.717e+03 -43.83 <2e-16 ***
## dataKematianApril$Tanggal 2.141e-04 4.863e-06 44.02 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.31 on 28 degrees of freedom
## Multiple R-squared: 0.9858, Adjusted R-squared: 0.9853
## F-statistic: 1938 on 1 and 28 DF, p-value: < 2.2e-161.7. Penjelasan Model Dengan Fungsi Summary()
Di atas merupakan rincian dari model yang telah dibuat.
Di posisi paling atas terdapat lm formula adalah dataKematianApril dataKematianApril Tanggal, data = dataKematianApril.
Lalu di bawahnya terdapat 5 nilai residual, sebelumnya kita perlu tahu bahwa Residual adalah perbedaan antara nilai nyata dan nilai prediksi. Yang mana semakin kecil nilai residual maka semakin baik atau benar model yang kita buat. Berikut nilai-nilai residual yang dihasilkan:
Nilai minimum = -26.273
Nilai maximum = 77.201
Nilai median = -3.993
Nilai quartil 1 = -13.785
Nilai quartil 3 = 5.850
Dari nilai-nilai tersebut dapat kita lihat bahwa dalam konteks ini berupa nilai minimum, maximum, median, quartil 1 dan quartil 3. Dapat kita simpulkan bahwa model yang telah kita buat belum bisa dikatakan baik atau benar karena nilai-nilai yang dihasilkan tidak mendekati nol.
Di bawah nilai residual terdapat koefisien, yang mana dalam koefisien tersebut terdapat nilai intersep, dan tanggal. Selain itu juga terdapat nilai-p dari koefisien
Selanjutnya terdapat dua R2 yaitu:
Multiple R-squared (pelipatan daripada variabel):0.9858. Hal ini menunjukkan bahwa 0.009858% variasi variabel respon, y, dapat dijelaskan oleh variabel prediktor, x.
Multiple R-squared tidak dapat berkurang saat kita menambahkan lebih banyak variabel independen ke model yang kita buat.
Adjusted R-squared (tidak pelipatan pada variabel): 0.9853.
Adjusted R-squared lebih baik ada penambahan variabel. Jadi jika kita menambahkan lebih dari satu variabel ke model, itu hanya meningkat jika itu mengurangi kesalahan prediksi secara keseluruhan
1.8. Aplikasi Fungsi anova()
ANOVA (analysis of variance) adalah pengujian yang dilakukan dengan membandingkan varians. Dengan membandingkan varians tersebut, dapat diketahui ada tidaknya perbedaan rata-rata dari tiga atau lebih kelompok. Asumsi normalitas pada ANOVA adalah pada residual yaitu selisih antara Y Prediksi dengan Y Aktual. Tepatnya residual dapat dihitung sebagai berikut: Y Aktual – Y Prediksi. Dimana Y Aktual adalah Y sesungguhnya atau kenyataan. Sedangkan Y prediksi adalah Y hasil persamaan ANOVA.
anova(modelData)## Analysis of Variance Table
##
## Response: dataKematianApril$Meninggal
## Df Sum Sq Mean Sq F value Pr(>F)
## dataKematianApril$Tanggal 1 799478 799478 1938.2 < 2.2e-16 ***
## Residuals 28 11550 412
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11.9. Membuat plot() Model Data Real Dan Data Prediksi
plot(dataKematianApril$Meninggal ~ dataKematianApril$
Tanggal,
data = dataKematianApril, col = "red", pch = 20, cex = 1.5,
main = "Data Orang Meninggal Covid-19 di DKI Jakarta dan Google Mobility Index")
abline(modelData) 
Titik-titik merah pada grafik merupakan data real Garis hitam yang
membelah diagonal adalah data prediksi.
plot(cooks.distance(modelData), pch = 16, col = "blue") 
plot(modelData)
1.10. Penggunaan AIC dan BIC
AIC dan BIC banyak digunakan dalam kriteria pemilihan model. AIC berarti Kriteria Informasi Akaike dan BIC berarti Kriteria Informasi Bayesian. Meskipun kedua istilah ini membahas pemilihan model, keduanya tidak sama. Seseorang dapat menemukan perbedaan antara dua pendekatan pemilihan model.
AIC(modelData)## [1] 269.7332BIC(modelData)## [1] 273.93681.11. Menampilkan Dan Memvisualisasikan Nilai Predicted
head(predict(modelData), n = 20)## 1 2 3 4 5 6 7 8
## 1214.308 1232.805 1251.302 1269.799 1288.296 1306.793 1325.290 1343.788
## 9 10 11 12 13 14 15 16
## 1269.799 1380.782 1399.279 1417.776 1436.273 1454.770 1473.268 1491.765
## 17 18 19 20
## 1510.262 1528.759 1547.256 1565.753plot(head(predict(modelData), n = 20))
1.12. Menampilkan Nilai Residuals
head(resid(modelData), n = 15)## 1 2 3 4 5 6 7
## 4.692439 4.195292 1.698146 -9.799001 -11.296148 -17.793295 -7.290442
## 8 9 10 11 12 13 14
## -13.787589 77.200999 -15.781883 -17.279030 -13.776177 -26.273324 -14.770471
## 15
## -5.267618coef(modelData)## (Intercept) dataKematianApril$Tanggal
## -3.382638e+05 2.140873e-041.13. Membuat Tabel Untuk Menambah Data Residuals dan Data Predicted
Tabel di bawah merupakan semua proses yang telah dilakukan, sehingga terbuat tabel yang ada nilai residuals dan nilai protected.
dataKematianApril$residuals <- modelData$residualsdataKematianApril$predicted <- modelData$fitted.values
dataKematianApril## # A tibble: 30 x 10
## Tanggal Meninggal retail_and_recreation_percent~ grocery_and_pha~
## <dttm> <dbl> <dbl> <dbl>
## 1 2020-04-01 00:00:00 1219 -24 -2
## 2 2020-04-02 00:00:00 1237 -25 -6
## 3 2020-04-03 00:00:00 1253 -25 -6
## 4 2020-04-04 00:00:00 1260 -22 -5
## 5 2020-04-05 00:00:00 1277 -25 -4
## 6 2020-04-06 00:00:00 1289 -30 -10
## 7 2020-04-07 00:00:00 1318 -25 -8
## 8 2020-04-08 00:00:00 1330 -26 -9
## 9 2020-04-04 00:00:00 1347 -26 -10
## 10 2020-04-10 00:00:00 1365 -26 -10
## # ... with 20 more rows, and 6 more variables:
## # parks_percent_change_from_baseline <dbl>,
## # transit_stations_percent_change_from_baseline <dbl>,
## # workplaces_percent_change_from_baseline <dbl>,
## # residential_percent_change_from_baseline <dbl>, residuals <dbl>,
## # predicted <dbl>1.14. Visualisasi Data Menggunakan scatter.smooth, boxplot dan plot
scatter.smooth(x=dataKematianApril$Tanggal, y=dataKematianApril$Meninggal,
main="Tanggal ~ MENINGGAL")
boxplot(dataKematianApril$Meninggal, main="Kasus Meninggal Covid",
boxplot.stats(dataKematianApril$Kematian)$out)## Warning: Unknown or uninitialised column: `Kematian`.
plot(density(dataKematianApril$Meninggal), main="Google Mobility Index: Kasus Covid Meninggal",
ylab="Frequency")
coefs <- coef(modelData)
plot(Meninggal ~ Tanggal, data = dataKematianApril)
abline(coefs)
text(x = 12, y = 10, paste('expression = ', round(coefs[1], 2), '+',
round(coefs[2], 2), '*Meninggal'))
1.15. Uji Korelasi Antar Variabel
Keterkaitan atau keterhubungan antar variabel bisa dilakukan dengan cara melakukan visualisasi menggunakan scatterplot dengan hitungan matematis yang menggunakan metode Pearson. Sedangkan metode parametrik, metode rangking Spearman dan Kendall digunakan untuk metode non-parametrik. Pada R uji keterhubungan dilakukan dengan menggunakan fungsi cor.test().
a) Uji Korelasi Antara Variabel Y Dengan X1
cor.test(dataKematianApril$
retail_and_recreation_percent_change_from_baseline,
dataKematianApril$Meninggal)##
## Pearson's product-moment correlation
##
## data: dataKematianApril$retail_and_recreation_percent_change_from_baseline and dataKematianApril$Meninggal
## t = -8.5707, df = 28, p-value = 2.579e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9269929 -0.7075143
## sample estimates:
## cor
## -0.8508943b) Uji Korelasi Antara Variabel Y Dengan X2
cor.test(dataKematianApril$
grocery_and_pharmacy_percent_change_from_baseline,
dataKematianApril$Meninggal)##
## Pearson's product-moment correlation
##
## data: dataKematianApril$grocery_and_pharmacy_percent_change_from_baseline and dataKematianApril$Meninggal
## t = -3.4996, df = 28, p-value = 0.001578
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.7607204 -0.2388293
## sample estimates:
## cor
## -0.5516343c) Uji Korelasi Antara Variabel Y Dengan X3
cor.test(dataKematianApril$
parks_percent_change_from_baseline,
dataKematianApril$Meninggal)##
## Pearson's product-moment correlation
##
## data: dataKematianApril$parks_percent_change_from_baseline and dataKematianApril$Meninggal
## t = -3.514, df = 28, p-value = 0.00152
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.7616761 -0.2409709
## sample estimates:
## cor
## -0.5532132d) Uji Korelasi Antara Variabel Y Dengan X4
cor.test(dataKematianApril$
transit_stations_percent_change_from_baseline,
dataKematianApril$Meninggal)##
## Pearson's product-moment correlation
##
## data: dataKematianApril$transit_stations_percent_change_from_baseline and dataKematianApril$Meninggal
## t = -4.9496, df = 28, p-value = 3.183e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.8373173 -0.4282402
## sample estimates:
## cor
## -0.6831169e) Uji Korelasi Antara Variabel Y Dengan X5
cor.test(dataKematianApril$
workplaces_percent_change_from_baseline,
dataKematianApril$Meninggal)##
## Pearson's product-moment correlation
##
## data: dataKematianApril$workplaces_percent_change_from_baseline and dataKematianApril$Meninggal
## t = -1.4438, df = 28, p-value = 0.1599
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5694865 0.1072143
## sample estimates:
## cor
## -0.2632221f) Uji Korelasi Antara Variabel Y Dengan X6
cor.test(dataKematianApril$
residential_percent_change_from_baseline,
dataKematianApril$Meninggal)##
## Pearson's product-moment correlation
##
## data: dataKematianApril$residential_percent_change_from_baseline and dataKematianApril$Meninggal
## t = 3.8462, df = 28, p-value = 0.0006337
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2888754 0.7824778
## sample estimates:
## cor
## 0.5879545Berdasarkan seluruh output yang dihasilkan, metode Pearson menghasilkan output berupa nilai t uji, derajat kebebasan, nilai p-value, rentang estimasi nilai korelasi berdasarkan tingkat kepercayaan, dan estimasi nilai korelasi.