Regresi Linear pada Kasus Covid-19

Model Prediksi Jumlah Kasus COVID-19 di Provinsi DKI Jakarta

knitr::include_graphics("Assets/covid.jpg")

Introduction

Background and Goals

Berdasarkan World Health Organization (WHO), Corona Virus Disease 19 (Covid-19) sampai saat ini telah tersebar lebih dari 122 negara, termasuk Indonesia. Data menunjukan bahwa DKI Jakarta sebagai Ibu Kota Negara Indonesia menjadi salah satu provinsi dengan jumlah kasus terbanyak dengan persentase sebesar 21,4% per 31 Juli 2021 menurut website https://covid19.go.id . Tingginya kasus yang terjadi di DKI Jakarta, menjadikan DKI Jakarta sebagai provinsi penyumbang terbanyak atas kasus konfirmasi positif Covid-19 di Indonesia.

Akan digunakan analisis regresi dalam mengetahui faktor-faktor yang mempengaruhi jumlah kasus positif Covid-19 per kecamatan di Provinsi DKI Jakarta. Hasil analisis yang kami lakukan dapat menjadi acuan dan perhatian untuk masyarakat luas, dan diharapkan angka pertumbuhan kasus positif Covid-19 di wilayah DKI Jakarta dapat diminimalisir.

Dataset

Data yang dipakai ialah data sekunder. * Untuk variabel respon didapatkan melalui website Open Data Covid-19 Provinsi DKI Jakarta (https://riwayat-file-covid-19-dki-jakarta-jakartagis.hub.arcgis.com/). * Untuk variabel jumlah rumah sakit rujukan penanggulangan Covid-19 per kecamatan, jumlah penduduk per kecamatan, jumlah penduduk usia produktif, dan jumlah kedatangan penduduk dari luar wilayah Provinsi DKI Jakarta didapatkan melalui website Open Data Pemerintah DKI Jakarta (https://data.jakarta.go.id/). * Data jumlah rumah sakit rujukan penanggulangan Covid-19 per kecamatan merupakan data per tahun 2020, data jumlah penduduk per kecamatan dan jumlah penduduk usia produktif merupakan data per tahun 2019, data jumlah kedatangan penduduk dari luar wilayah Provinsi DKI Jakarta merupakan data bulan September 2020. Untuk variabel jumlah pelanggaran kasus pada PSBB didapatkan melalui tautan (https://datastudio.google.com/)

Load Packages

Berikut merupakan basic packages yang digunakan selama pengerjaan analisis:

library(readr) # untuk membaca data
library(dplyr) # untuk data manipulation
library(GGally) # untuk membuat matriks korelasi
library(lmtest) # untuk pengujian Breusch-Pagan/heteroscedasticity
library(car) # untuk pengujian multikolinearitas
library(caret)
library(performance) # for model comparison and assumption check)
library(tidyr)
library(rmdformats)

Load Dataset

covid <- readxl::read_xlsx("dataset/data_covid.xlsx")
rmarkdown::paged_table(covid)

Data Wrangling

glimpse(covid)

#> Rows: 45
#> Columns: 10
#> $ kota           <chr> "JAKARTA BARAT", "JAKARTA TIMUR", "JAKARTA BARAT", "JAK…
#> $ kecamatan      <chr> "CENGKARENG", "CAKUNG", "KALI DERES", "CILINCING", "DUR…
#> $ jumlah         <dbl> 1104, 993, 525, 1309, 1278, 1613, 877, 690, 1356, 713, …
#> $ rumah_sakit    <dbl> 2, 2, 1, 2, 9, 8, 2, 7, 3, 3, 5, 4, 3, 1, 4, 1, 10, 2, …
#> $ penduduk       <dbl> 406806, 399252, 315170, 302248, 322504, 305222, 262859,…
#> $ pelanggaran    <dbl> 89, 114, 43, 68, 100, 140, 77, 139, 60, 95, 33, 86, 119…
#> $ usia_produktif <dbl> 246606, 243331, 193692, 188296, 181266, 176350, 151794,…
#> $ kedatangan     <dbl> 616, 709, 482, 526, 535, 401, 342, 306, 294, 297, 364, …
#> $ Latitude       <dbl> -6.146312, -6.181968, -6.133928, -6.121266, -6.232147, …
#> $ Longitude      <dbl> 106.7349, 106.9474, 106.7058, 106.9476, 106.9152, 106.8…

Deskripsi untuk setiap kolomnya antara lain:

kota : Kota administratif pada Provinsi DKI Jakarta
kecamatan : Nama kecamatan pada masing-masing kota administratif pada Provinsi DKI Jakarta
jumlah : Jumlah kasus positif Covid-19 per kecamatan
rumah_sakit : Jumlah rumah sakit rujukan penanggulangan Covid-19 per kecamatan
penduduk : Jumlah penduduk per kecamatan
pelanggaran : Jumlah pelanggaran kasus pada PSBB
usia_produktif : Jumlah penduduk usia produktif
kedatangan : Jumlah kedatangan penduduk dari luar wilayah Provinsi DKI Jakarta
Latitude : Titik Latitude untuk suatu Kecamatan
Longitude : Titik Longitude untuk suatu kecamatan

Selanjutnya, kita akan melakukan data transformasi dengan mengubah tipe data kota dan kecamatan menjadi factor serta pemilihan variabel-varibel yang digunakan untuk kebutuhan analisis EDA, di mana akan dihilangkan variabel Latitude dan Longitude.

covid <- covid %>%
  mutate(
    kota = as.factor(kota),
    kecamatan = as.factor(kecamatan)) %>% 
  select(-c('Latitude', 'Longitude'))
glimpse(covid)

#> Rows: 45
#> Columns: 8
#> $ kota           <fct> JAKARTA BARAT, JAKARTA TIMUR, JAKARTA BARAT, JAKARTA UT…
#> $ kecamatan      <fct> CENGKARENG, CAKUNG, KALI DERES, CILINCING, DUREN SAWIT,…
#> $ jumlah         <dbl> 1104, 993, 525, 1309, 1278, 1613, 877, 690, 1356, 713, …
#> $ rumah_sakit    <dbl> 2, 2, 1, 2, 9, 8, 2, 7, 3, 3, 5, 4, 3, 1, 4, 1, 10, 2, …
#> $ penduduk       <dbl> 406806, 399252, 315170, 302248, 322504, 305222, 262859,…
#> $ pelanggaran    <dbl> 89, 114, 43, 68, 100, 140, 77, 139, 60, 95, 33, 86, 119…
#> $ usia_produktif <dbl> 246606, 243331, 193692, 188296, 181266, 176350, 151794,…
#> $ kedatangan     <dbl> 616, 709, 482, 526, 535, 401, 342, 306, 294, 297, 364, …

Exploratory Data Analysis

Exploratory Data Analysis mencakup tentang proses kritis uji investigasi awal pada sebuah data untuk melakukan pengecekan missing value, membuat ringkasan data, dan representasi grafis (visual).

1. Melakukan pengecekan missing value

colSums(is.na(covid))

#>           kota      kecamatan         jumlah    rumah_sakit       penduduk 
#>              1              1              0              1              0 
#>    pelanggaran usia_produktif     kedatangan 
#>              1              1              1

Karena terdapat missing value pada data, akan dilakukan pembuangan

covid <- covid %>%
  drop_na()

Selanjutnya akan dilakukan pengecekan kembali untuk kemunculan missing value pada data

colSums(is.na(covid))

#>           kota      kecamatan         jumlah    rumah_sakit       penduduk 
#>              0              0              0              0              0 
#>    pelanggaran usia_produktif     kedatangan 
#>              0              0              0

2. Correlation plot untuk masing-masing variabel

Pada bagian ini akan dicari nilai korelasi antar masing-masing variabel.

ggcorr(covid, hjust = 1, layout.exp = 1, label = TRUE)

Insight:

variabel penduduk dan usia_produktif memiliki korelasi yang cukup kuat dengan variabel jumlah
variabel penduduk dan kedatangan serta variabel usia_produktif dan kedatangan memiliki hubungan korelasi positif yang kuat sebesar 0.9.

covidx <- covid %>%  
  select(-c(kota, kecamatan))

rmarkdown::paged_table(covidx)

3. Melihat visualisasi korelasi antar 2 variabel

Jumlah kasus positif vs Jumlah ketersedian RS pada suatu kecamatan

covid %>% 
  ggplot(data = covid, mapping = aes(x = rumah_sakit, y = jumlah)) +
  geom_point(aes(colour = rumah_sakit)) +
  theme_minimal() +
  labs(title = "Jumlah kasus positif vs Jumlah ketersedian RS pada suatu kecamatan",
       x = "Jumlah ketersedian rumah sakit",
       y = "Jumlah kasus positif")

Insight:

Jumlah ketersediaan rumah sakit per kecamatan lebih banyak berada pada kisaran rentang nilai < 6 RS. Tidak terlihat adanya korelasi antara kedua nilai tersebut.

Jumlah kasus postif vs Jumlah penduduk

covid %>% 
  ggplot(data = covid, mapping = aes(x = penduduk, y = jumlah)) +
  geom_point(aes(colour = penduduk)) +
  theme_minimal() +
  labs(title = "Jumlah kasus positif vs Jumlah Penduduk",
       x = "Jumlah penduduk",
       y = "Jumlah kasus positif")

Insight:

Semakin banyak jumlah penduduk, semakin banyak pula jumlah kasus positif covid per masing-masing kecamatannya

Jumlah kasus positif vs Kasus pelanggaran PSBB

covid %>% 
  ggplot(data = covid, mapping = aes(x = pelanggaran, y = jumlah)) +
  geom_point(aes(colour = pelanggaran)) +
  theme_minimal() +
  labs(title = "Jumlah kasus positif vs Kasus pelanggaran PSBB",
       x = "Kasus pelanggaran PSBB",
       y = "Jumlah kasus positif")

Insight:

Jumlah kasus pelanggaran PSBB sedikit banyak menunjukkan hubungan signifikan antara jumlah kasus positif. Terdapat 1 observasi dengan kasus pelanggaran PSBB terbanyak menunjukkan kasus positif covid-19 yang tinggi.

Jumlah kasus positif vs Jumlah kedatangan penduduk dari luar Jakarta

covid %>% 
  ggplot(data = covid, mapping = aes(x = kedatangan, y = jumlah)) +
  geom_point(aes(colour = kedatangan)) +
  theme_minimal() +
  labs(title = "Jumlah kasus positif vs Jumlah kedatangan penduduk dari luar Jakarta",
       x = "Jumlah kedatangan penduduk dari luar Jakarta",
       y = "Jumlah kasus positif")

Insight :

Jumlah kedatangan penduduk dari luar jakarta ke jakarta sedikit banyak menunjukkan hubungan signifikan antara jumlah kasus positif covid.

Model Building

Simple Linear Regression

Akan dibuat model regresi linear dengan variabel prediktor penduduk karena variabel tersebut memiliki korelasi positif tertinggi di antara variabel lainnya terhadap variabel target jumlah.

model_slr <- lm(jumlah ~ penduduk, data = covid)
summary(model_slr)

#> 
#> Call:
#> lm(formula = jumlah ~ penduduk, data = covid)
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#> -502.0 -168.6  -41.2  123.4  606.5 
#> 
#> Coefficients:
#>                Estimate  Std. Error t value  Pr(>|t|)    
#> (Intercept) 377.1098770  86.3623934   4.367 0.0000808 ***
#> penduduk      0.0020620   0.0004216   4.890 0.0000152 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 251.9 on 42 degrees of freedom
#> Multiple R-squared:  0.3628, Adjusted R-squared:  0.3477 
#> F-statistic: 23.92 on 1 and 42 DF,  p-value: 0.00001519

ggplot(data = covid, mapping=aes(x = penduduk, y = jumlah))+
  geom_point()+
  scale_color_brewer(palette="Accent")+
  geom_smooth(method="lm", se=FALSE)+
  labs(title="Plot Regresi Penduduk vs Jumlah Positif Covid")+
  theme_minimal()

Insight:

Dengan Menguji Signifikansi Parameter, jumlah penduduk memiliki nilai p-value 0.0000152 < 0.05, sehingga disimpulkan bahwa jumlah penduduk memiliki pengaruh yang signifikan terhadap peningkatan kasus positif covid-19 per kecamatan di DKI Jakarta
Adj R-Squared = 0.3477, dengan kata lain hanya sebesar 34,77% dari variabel jumlah penduduk dapat menjelaskan model (pengaruh jumlah kasus positif).
Walaupun secara signifikansi parameter, variabel penduduk berpengaruh signifikan, akan tetapi dari nilai $adj R^2$, nilai yang dihasilkan masih sangat kecil, sehingga perlu penambahan faktor variabel lainnya.

Multiple Linear Regression

Akan dibuat model dengan keseluruhan variabel, dengan variabel target adalah jumlah dan variabel prediktor antara lain rumah_sakit, penduduk, pelanggaran, usia_produktif, dan kedatangan.

mlr_model <- lm(jumlah ~ rumah_sakit + penduduk + pelanggaran + usia_produktif + kedatangan, data = covid)
summary(mlr_model)

#> 
#> Call:
#> lm(formula = jumlah ~ rumah_sakit + penduduk + pelanggaran + 
#>     usia_produktif + kedatangan, data = covid)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -401.18 -163.58  -58.09  142.61  503.63 
#> 
#> Coefficients:
#>                  Estimate Std. Error t value Pr(>|t|)  
#> (Intercept)    252.761094 113.873377   2.220   0.0325 *
#> rumah_sakit     15.426457  16.685515   0.925   0.3610  
#> penduduk         0.002215   0.006385   0.347   0.7305  
#> pelanggaran      1.182860   1.573550   0.752   0.4569  
#> usia_produktif   0.001298   0.011106   0.117   0.9076  
#> kedatangan      -0.710305   0.787732  -0.902   0.3729  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 252.4 on 38 degrees of freedom
#> Multiple R-squared:  0.4212, Adjusted R-squared:  0.345 
#> F-statistic:  5.53 on 5 and 38 DF,  p-value: 0.000638

Insight:

Dapat dilihat dari summary() di atas bahwa untuk mlr_model hanya nilai intercept yang bernilai signifikan terhadap model.
Nilai $Adj R^2$ sebesar 0.345 yang mana hanya sebesar 34,5% dari variabel-variabel dapat dijelaskan oleh model.

Evaluasi Model Keseluruhan

Performa Model

Kinerja model (seberapa baik model kami memprediksi variabel target) dapat dihitung dengan menggunakan nilai RMSE (Root Mean Square Error) dengan formula berikut:

\[ RMSE = \sqrt{\frac{\sum (prediction_i - actual_i)^2}{n}}\] RMSE merupakan suatu nilai pengujian performa yang lebih baik dibandingkan MAE (Mean Absolute Error). Hal ini disebabkan nilai RMSE mengkuadratkan selisih antara nilai aktual dan nilai prediksi, artinya prediksi dengan error yang lebih tinggi akan dikenakan sanksi yang besar. Kita dapat menggunakan fungsi RMSE () dari package caret

# make prediction
covidx$positif_pred <- predict(mlr_model, covidx)

# calculating RMSE
sqrt(mean((covidx$positif_pred - covidx$jumlah)^2))

#> [1] 234.5659

range(covid$jumlah)

#> [1]   24 1613

Berdasarkan hasil ini, model mlr_model masih memberikan kesalahan yang cukup tinggi untuk prediksi kami dan oleh karena itu masih perlu perbaikan.

Pengujian Asumsi

Terdapat 4 pengujian asumsi residual yang harus dilakukan antara lain:

Linearitas
Normalitas
Autokorelasi
Heteroskedastisitas
Multikolinearitas

Linearitas

Pengecekan asumsi linearitas sudah diperiksa sebelum membangun model menggunakan matriks korelasi pada langkah EDA dan juga scatter plot yang dibuat menggunakan variabel prediktor dan variabel target. Akan tetapi akan dilakukan pemeriksaan cor.test() pada salah satu variabel prediktor dengan target.

Hipotesis

H0 : tidak ada korelasi yang signifikan (p=0)
H1 : terdapat korelasi signifikan (p=!0)

$\alpha$ = 0.05

Statistik Uji: Pearson dengan nilai p-value

cor.test(covidx$pelanggaran, covidx$jumlah)

#> 
#>  Pearson's product-moment correlation
#> 
#> data:  covidx$pelanggaran and covidx$jumlah
#> t = 3.0788, df = 42, p-value = 0.003655
#> alternative hypothesis: true correlation is not equal to 0
#> 95 percent confidence interval:
#>  0.1515307 0.6439521
#> sample estimates:
#>      cor 
#> 0.429107

Kesimpulan: nilai p-value = 0.003655 < alpha = 0.95, sehingga H0 ditolak. Sehingga terdapat korelasi yang signifikan antara variabel pelanggaran dan jumlah.

Selain itu, akan dilakukan pengecekan asumsi linearitas dengan membuat plot residual vs fitted untuk mlr_model.

resfit <- data.frame(residual = mlr_model$residuals, fitted = mlr_model$fitted.values)

resfit %>% ggplot(aes(fitted, residual)) + 
  geom_point() + 
  geom_hline(aes(yintercept = 0, colour= "red")) +
  theme_minimal() + 
  labs(title = "Residual vs Fitted Plot") + 
  theme(legend.position = "none")

Insight:

Dari plot di atas dapat dilihat bahwa residual/error tidak memantul secara acak di sekitar 0 (tidak seragam) sehingga akan menghasilkan rata-rata tidak sama dengan 0. Maka syarat $E(\epsilon) = 0$ tidak terpenuhi. Disimpulkan bahwa data kita mungkin tidak linear. Asumsi tidak terpenuhi.

Normalitas Residual

Akan dibuat plot histogram dan juga normal probability plot (normal qq plot) untuk melihat distribusi dari eror.

par(mfrow=c(1,2))
hist(mlr_model$residuals, breaks = 10)
qqnorm(residuals(mlr_model), main="Normal Q-Q Plot")
qqline(residuals(mlr_model))

Insight:

Dapat dilihat dari plot, bahwa terdapat beberapa residual yang tidak berada pada garis. Hal itu mengindikasikan bahwa mungkin saja distribusi dari eror tidak normal. Untuk lebih meyakinkan, akan dilakukan pengujian normalitas dengan menggunakan shapiro wilk test.

Hipotesis:

H0: error berdistribusi normal
H1: error TIDAK berdistribusi normal

$\alpha$ = 0.05

Statistik Uji: Shapiro Wilk dengan nilai p-value

#Shapiro wilk test
shapiro.test(mlr_model$residuals)

#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  mlr_model$residuals
#> W = 0.94619, p-value = 0.03977

Aturan Keputusan p-value < 0.05; H0 ditolak.

Kesimpulan Maka dari itu, error/residual tidak berdistribusi normal. Asumi normalitas residual tidak terpenuhi.

Autokorelasi

Hipotesis

H0 : tidak ada autokorelasi
H1 : ada auatokorelasi

$\alpha$ = 0.05

Statistik Uji: Durbin Watson dengan nilai p-value

durbinWatsonTest(mlr_model)

#>  lag Autocorrelation D-W Statistic p-value
#>    1     -0.08709257      2.142334   0.826
#>  Alternative hypothesis: rho != 0

Aturan Keputusan Nilai p-value = 0.826 > 0.05 sehingga H0 gagal ditolak

Kesimpulan Karena H0 gagal ditolak, sehingga tidak terdapat autokorelasi pada model. Asumsi terpenuhi

Homoskedastisitas

Hipotesis

H0: variansi residual seragam (homoskedastisitas)
H1: variansi residual tidak seragam (heteroskedastisitas)

$\alpha$ = 0.05

Statistik Uji: Breusch Pagan dengan nilai p-value

bptest(mlr_model)

#> 
#>  studentized Breusch-Pagan test
#> 
#> data:  mlr_model
#> BP = 7.8203, df = 5, p-value = 0.1664

Aturan Keputusan Nilai p-value = 0.1664 > 0.05 sehingga H0 gagal ditolak

Kesimpulan Karena H0 gagal ditolak, sehingga residual model yang kita miliki memiliki variansi yang seragam (homoskedastisitas). Asumsi terpenuhi.

Multikolinearitas

Asumsi multikolinearitas perlu diuji untuk mengetahui apakah terdapat variabel bebas yang saling berkorelasi secara kuat. Variabel bebas yang berkorelasi secara kuat dapat memengaruhi kekuatan prediksi model yang diajukan. Suatu variabel tidak menunjukan gejala multikolinearitas apabila memiliki nilai VIF < 10. Akan dianalisis menggunakan nilai VIF masing-masing variabel

vif(mlr_model)

#>    rumah_sakit       penduduk    pelanggaran usia_produktif     kedatangan 
#>       1.227968     228.357484       1.633895     259.699084      10.844641

Dari output diketahui bahwa variabel kedatangan, penduduk, dan usia_produktif memiliki nilai VIF > 10, sehingga model terdapat gejala multikolinaritas.

Model Improvement

Dari analisis residual yang dilakukan, disimpulkan bahwa residual dari mlr_model (full model) menunjukan ketidaklinearan, tidak normal, dan multikolinearitas. Dengan demikian akan dilakukan transformasi data untuk variabel prediktor dan respon. Salah satu jenis transformasi yang dapat dilakukan adalah transformasi logaritmik.

Berdasarkan hasil summary mlr_model diketahui bahwa nilai p-value terbesar terdapat pada variabel usia_produktif , dengan demikian akan dipertimbangkan untuk melakukan transformasi logaritmik pada variabel jumlah dan variabel usia_produktif.

Tuning Model

Transformasi dengan Komponen Logaritmik

trans_model <- lm(log(jumlah) ~ rumah_sakit + penduduk + pelanggaran + log(usia_produktif) + kedatangan, data = covid)
summary(trans_model)

#> 
#> Call:
#> lm(formula = log(jumlah) ~ rumah_sakit + penduduk + pelanggaran + 
#>     log(usia_produktif) + kedatangan, data = covid)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -0.52036 -0.25318 -0.04665  0.26545  0.77290 
#> 
#> Coefficients:
#>                         Estimate   Std. Error t value       Pr(>|t|)    
#> (Intercept)         -8.240166046  1.562640467  -5.273 0.000005638696 ***
#> rumah_sakit          0.011424252  0.022585758   0.506         0.6159    
#> penduduk            -0.000005279  0.000002375  -2.222         0.0323 *  
#> pelanggaran          0.002865193  0.002095211   1.367         0.1795    
#> log(usia_produktif)  1.368993562  0.158156707   8.656 0.000000000161 ***
#> kedatangan          -0.000438848  0.000998781  -0.439         0.6629    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.3532 on 38 degrees of freedom
#> Multiple R-squared:  0.8049, Adjusted R-squared:  0.7793 
#> F-statistic: 31.36 on 5 and 38 DF,  p-value: 0.000000000001642

Dari trans_model didapatkan persamaan regresi sebagai berikut: $ln(jumlah)=−8.240 + 0.01142rumahsakit − 0.00000527penduduk + 0.00286pelanggaran + 1.369ln(usiaproduktif) − 0.0004388kedatangan$

Berdasarkan summary model trans_model, nilai p-value variabel jumlah penduduk per kecamatan dan variabel jumlah penduduk usia produktif lebih kecil dari 0.1 sehingga dapat disimpulkan bahwa kedua variabel bebas tersebut terdapat pengaruh terhadap variabel target.

Model Backward Regression pada Model Transformasi

Selanjutnya, akan dilakukan backward regression dengan tujuan menemukan model terbaik pada trans_model. Kriteria AIC dilakukan untuk mencari model terbaik yang digunakan pada kasus ini. Nilai AIC terkecil menjadi kriteria pemilihan model yang paling baik. Dilakukan metode backward regression untuk pemilihan model menggunakan kriteria AIC. Backward regression ialah metode yang akan mengeluarkan variabel bebas yang tidak berpengaruh signifikan terhadap variabel respon secara satu per satu dilihat dari nilai p-value masing-masing variabel bebas, variabel bebas yang pertama dikeluarkan ialah yang memiliki p-value paling besar dan melebihi taraf signifikansi

set.seed(123)
step(trans_model, direction = "backward")

#> Start:  AIC=-86.02
#> log(jumlah) ~ rumah_sakit + penduduk + pelanggaran + log(usia_produktif) + 
#>     kedatangan
#> 
#>                       Df Sum of Sq     RSS     AIC
#> - kedatangan           1    0.0241  4.7658 -87.800
#> - rumah_sakit          1    0.0319  4.7736 -87.728
#> <none>                              4.7417 -86.023
#> - pelanggaran          1    0.2333  4.9750 -85.910
#> - penduduk             1    0.6163  5.3580 -82.646
#> - log(usia_produktif)  1    9.3492 14.0909 -40.101
#> 
#> Step:  AIC=-87.8
#> log(jumlah) ~ rumah_sakit + penduduk + pelanggaran + log(usia_produktif)
#> 
#>                       Df Sum of Sq     RSS     AIC
#> - rumah_sakit          1    0.0373  4.8031 -89.457
#> <none>                              4.7658 -87.800
#> - pelanggaran          1    0.2633  5.0290 -87.435
#> - penduduk             1    2.8269  7.5927 -69.308
#> - log(usia_produktif)  1   10.1885 14.9543 -39.484
#> 
#> Step:  AIC=-89.46
#> log(jumlah) ~ penduduk + pelanggaran + log(usia_produktif)
#> 
#>                       Df Sum of Sq     RSS     AIC
#> <none>                              4.8031 -89.457
#> - pelanggaran          1    0.3603  5.1633 -88.275
#> - penduduk             1    2.9930  7.7961 -70.145
#> - log(usia_produktif)  1   10.5117 15.3147 -40.436

#> 
#> Call:
#> lm(formula = log(jumlah) ~ penduduk + pelanggaran + log(usia_produktif), 
#>     data = covid)
#> 
#> Coefficients:
#>         (Intercept)             penduduk          pelanggaran  
#>         -8.47624797          -0.00000626           0.00335114  
#> log(usia_produktif)  
#>          1.39652025

# memasukan hasil backward stepwise ke dalam suatu object
step_model <- lm(log(jumlah) ~ penduduk + pelanggaran + log(usia_produktif), data = covid)
summary(step_model)

#> 
#> Call:
#> lm(formula = log(jumlah) ~ penduduk + pelanggaran + log(usia_produktif), 
#>     data = covid)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -0.59641 -0.25521 -0.03329  0.27984  0.76748 
#> 
#> Coefficients:
#>                         Estimate   Std. Error t value        Pr(>|t|)    
#> (Intercept)         -8.476247972  1.486487199  -5.702 0.0000012389731 ***
#> penduduk            -0.000006260  0.000001254  -4.993 0.0000121229114 ***
#> pelanggaran          0.003351142  0.001934680   1.732           0.091 .  
#> log(usia_produktif)  1.396520248  0.149258874   9.356 0.0000000000127 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.3465 on 40 degrees of freedom
#> Multiple R-squared:  0.8024, Adjusted R-squared:  0.7876 
#> F-statistic: 54.15 on 3 and 40 DF,  p-value: 0.00000000000003807

Dari step_model didapatkan persamaan regresi sebagai berikut: $ln(jumlah)= −8.476 + − 0.0000062 penduduk + 0.00335pelanggaran + 1.39652ln(usiaproduktif)$

Berdasarkan summary model step_model, nilai p-value variabel jumlah penduduk per kecamatan, pelanggaran pada peraturan PSBB, dan variabel jumlah penduduk usia produktif lebih kecil dari 0.1 sehingga dapat disimpulkan bahwa ketiga variabel bebas tersebut terdapat pengaruh terhadap variabel target.

Performa Model

Untuk kedua model kandidat yaitu trans_model dan step_model, dilakukan perbandingan performa untuk kedua model tersebut. Tujuannya adalah mencari model regresi terbaik untuk memodelkan data

model_comparison <- compare_performance(trans_model, step_model)
                                         
as.data.frame(model_comparison) %>% 
  select("Name","AIC","BIC","R2","R2_adjusted","RMSE")

#>          Name      AIC      BIC        R2 R2_adjusted      RMSE
#> 1 trans_model 610.4879 622.9772 0.8049452   0.7792801 0.3282759
#> 2  step_model 607.0540 615.9749 0.8024195   0.7876010 0.3303944

Dapat dilihat bahwa model yang terbaik ialah step_model karena memiliki tiga performa kriteria yang lebih unggul dibandingkan performa kriteria pada trans_model, yaitu nilai Adj R-Square terbesar, nilai AIC terkecil, dan nilai BIC terkecil. Sehingga disimpulkan bahwa model regresi terbaik untuk memodelkan data pada kasus ini ialah step_model.

Pengujian Asumsi

1. Linearitas

resfit1 <- data.frame(residual = step_model$residuals, fitted = step_model$fitted.values)

resfit1 %>% ggplot(aes(fitted, residual)) + 
  geom_point() + 
  geom_hline(aes(yintercept = 0, colour= "red")) +
  theme_minimal() + 
  labs(title = "Residual vs Fitted Plot") + 
  theme(legend.position = "none")

Insight: Dari plot di atas dapat dilihat bahwa residual/error cukup memantul secara acak di sekitar 0 (seragam) sehingga akan menghasilkan rata-rata dengan 0. Maka syarat $E(\epsilon) = 0$ terpenuhi. Disimpulkan bahwa model kita linear. Asumsi terpenuhi.

2. Normalitas Residual

Hipotesis:

H0: error berdistribusi normal
H1: error TIDAK berdistribusi normal

$\alpha$ = 0.05

Statistik Uji: Shapiro Wilk dengan nilai p-value

#Shapiro wilk test
shapiro.test(step_model$residuals)

#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  step_model$residuals
#> W = 0.96357, p-value = 0.1767

Aturan Keputusan p-value > 0.05; H0 tidak ditolak.

Kesimpulan Maka dari itu, error/residual berdistribusi normal. Asumsi normalitas residual terpenuhi.

3. Autokorelasi

Hipotesis

H0 : tidak ada autokorelasi
H1 : ada auatokorelasi

$\alpha$ = 0.05

Statistik Uji: Durbin Watson dengan nilai p-value

durbinWatsonTest(step_model)

#>  lag Autocorrelation D-W Statistic p-value
#>    1      0.04982623      1.855166   0.432
#>  Alternative hypothesis: rho != 0

Aturan Keputusan Nilai p-value = 0.432 > 0.05 sehingga H0 gagal ditolak

Kesimpulan Karena H0 gagal ditolak, sehingga tidak terdapat autokorelasi pada model. Asumsi terpenuhi

4. Homoskedastisitas

Hipotesis

H0: variansi residual seragam (homoskedastisitas)
H1: variansi residual tidak seragam (heteroskedastisitas)

$\alpha$ = 0.05

Statistik Uji: Breusch Pagan dengan nilai p-value

bptest(step_model)

#> 
#>  studentized Breusch-Pagan test
#> 
#> data:  step_model
#> BP = 0.19602, df = 3, p-value = 0.9782

Aturan Keputusan Nilai p-value = 0.9782 > 0.05 sehingga H0 gagal ditolak

Kesimpulan Karena H0 gagal ditolak, sehingga residual model yang kita miliki memiliki variansi yang seragam (homoskedastisitas). Asumsi terpenuhi.

5. Multikolinearity

vif(step_model)

#>            penduduk         pelanggaran log(usia_produktif) 
#>            4.673168            1.310454            4.698178

Dari output diketahui bahwa semua nilai VIF < 10. Dengan demikian tidak terdapat gejala multikolinearitas pada masing-masing variabel bebas yang terdapat pada step_model.

Interpretasi Model

Dari pengujian asumsi step_model diketahui bahwa keseluruhan asumsi terpenuhi. Selanjutnya akan dilakukan interpretasi model berikut:

$ln(jumlah)= −8.476 + − 0.0000062 penduduk + 0.00335pelanggaran + 1.39652ln(usiaproduktif)$

Apabila variabel jumlah penduduk per kecamatan bertambah 1 satuan, maka rata-rata banyaknya jumlah kasus positif Covid-19 per kecamatan di Provinsi DKI Jakarta naik sebesar e^(−0.0000062)$ = 0.99999374 kali dengan variabel lain tetap.
Apabila variabel jumlah pelanggaran kasus pada PSBB bertambah 1 satuan, maka rata-rata banyaknya jumlah kasus positif Covid-19 per kecamatan di Provinsi DKI Jakarta naik sebesar e^(0.00335) = 1.0033556175 kali dengan variabel lain tetap.
Apabila variabel logaritma natural jumlah penduduk usia produktif bertambah 1 satuan, maka rata-rata banyaknya jumlah kasus positif Covid-19 per kecamatan di Provinsi DKI Jakarta naik 1.39652 kali dengan variabel lain tetap.

Kesimpulan

1. Berdasarkan proses analisis faktor-faktor pada data jumlah kasus positif Covid-19 per kecamatan di DKI Jakarta, didapatkan hasil bahwa faktor yang berpengaruh secara signifikan terhadap jumlah kasus positif Covid-19 per kecamatan di Provinsi DKI Jakarta yaitu jumlah penduduk per kecamatan, jumlah pelanggaran kasus PSBB per kecamatan, dan logaritma natural jumlah penduduk usia produktif per kecamatan.
1. Jumlah penduduk per kecamatan yang tinggi berdampak pada tingginya jumlah kasus positif Covid-19 per kecamatan di Provinsi DKI Jakarta. Hal ini dapat menjadi bahan pertimbangan dan juga saran untuk masyarakat setempat maupun pemerintah terkait untuk melakukan pengawasan yang jauh lebih komprehensif di daerah padat penduduk. Misalnya, dengan aktif mensosialisasikan protokol kesehatan, 3M, maupun menghadirkan lebih banyak posko-posko Satgas Covid-19 di kecamatan tersebut.
1. Jumlah penduduk usia produktif per kecamatan yang tinggi berdampak pada tingginya kasus positif Covid-19 per kecamatan di Provinsi DKI Jakarta. Hal ini mengindikasikan bahwa penduduk dengan usia produktif memiliki mobilitas yang jauh lebih tinggi dibandingkan dengan penduduk usia muda atau lansia. Kami melihat perlu diberlakukannya 3M yang lebih ketat untuk penduduk usia produktif baik di luar rumah maupun di dalam rumah. Harapannya, hal itu dapat menjaga diri sendiri dan lingkungan dari terpaparnya Covid-19.
1. Jumlah pelanggaran kasus PSBB per kecamatan yang tinggi juga berdampak pada tingginya jumlah kasus positif Covid-19 per kecamatan di Provinsi DKI Jakarta. Kebijakan dan upaya yang dilakukan pemerintah DKI Jakarta terkait dengan PSBB sepertinya belum mampu menekan persebaran Covid-19 per kecamatan di Provinsi DKI Jakarta. Hambatan ini berasal dari ketidakpatuhan masyarakat terhadap protokol kesehatan dan kebijakan PSBB yang diterapkan pemerintah. Kami melihat perlu diberlakukan hukuman yang lebih tegas dan nyata bagi pelanggar berupa denda ataupun sanksi sosial untuk menimbulkan efek jera bagi masyarakat yang melakukan pelanggaran PSBB, sehingga jumlah kasus postif Covid-19 dapat diminimalisir.