Regresi Linear 3 model (Smoothing Spline Model)

~ untuk Kasus Dirawat Covid-19 DKI Jakarta dibandingkan dengan Mobilitas di Terminal/Stasiun ~


Oleh Muhammad Aji Permana
Dosen Pengampu Prof. Dr M. Suhartono, M.Kom
Program Magister Informatika
Tanggal 10 Mei 2022

1. Latar Belakang

Smoothing merupakan salah satu metode yang digunakan dalam analisis data non parametrik. Tujuan dari smoothing adalah untuk meminimalkan keragaman karakteristik data dari data yang tidak memiliki pengaruh sehingga ciri-ciri dari data akan tampak lebih jelas.

Smooth Spline Model adalah Salah satu model regresi dengan pendekatan non parametrik yang dapat digunakan untuk menduga kurva regresi adalah regresi spline.

2. Import dataset

library(readxl)
## Warning: package 'readxl' was built under R version 4.1.3
mobilityjakarta <- read_excel(path = "Datagab.xlsx")
mobilityjakarta

3. Aktivasi Library

library(npreg)
## Warning: package 'npreg' was built under R version 4.1.3
library(jmv)
## Warning: package 'jmv' was built under R version 4.1.3
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.1.3
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.6     v dplyr   1.0.8
## v tidyr   1.2.0     v stringr 1.4.0
## v readr   2.1.2     v forcats 0.5.1
## Warning: package 'ggplot2' was built under R version 4.1.3
## Warning: package 'tibble' was built under R version 4.1.3
## Warning: package 'tidyr' was built under R version 4.1.3
## Warning: package 'readr' was built under R version 4.1.3
## Warning: package 'purrr' was built under R version 4.1.3
## Warning: package 'dplyr' was built under R version 4.1.3
## Warning: package 'stringr' was built under R version 4.1.3
## Warning: package 'forcats' was built under R version 4.1.3
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(skimr)
## Warning: package 'skimr' was built under R version 4.1.3
library(stats)
library(npreg)
descriptives(mobilityjakarta, vars = vars(Dirawat, Station), freq = TRUE)
## 
##  DESCRIPTIVES
## 
##  Descriptives                                    
##  ----------------------------------------------- 
##                          Dirawat     Station     
##  ----------------------------------------------- 
##    N                           30           30   
##    Missing                      0            0   
##    Mean                  1155.800    -37.53333   
##    Median                1044.000    -38.00000   
##    Standard deviation    499.9613     3.234868   
##    Minimum               510.0000    -44.00000   
##    Maximum               2172.000    -29.00000   
##  -----------------------------------------------

4. Histrogram plot

hist(mobilityjakarta$Dirawat, main="Distribusi Dirawat di RS untuk Covid 19", xlab="Dirawat di RS Untuk Covid 19", xlim=c(0,2500),ylim=c(0,15)) # histogram variabel wt
mean=mean(mobilityjakarta$Dirawat) # Menghitung nilai mean variabel wt
abline(v=mean, col="red",lwd=2) # Menambahkan garis vertical pada plot

5. Density plot

d<-density(mobilityjakarta$Dirawat) # menghitung density variabel wt
plot(d,  xlab="Dirawat di RS untuk Covid 19", main="Distribusi Dirawat di RS untuk Covid 19")

6. Pemodelan Mobilitas

Pada pemodelan ini dilakukan dengan linier dan smooth spline

6.1 Linier Model Mobilitas di Terminal/Stasiun

# linear / tren
model.linear <- lm(mobilityjakarta$Station ~ mobilityjakarta$Tanggal, data = mobilityjakarta)
model.linear
## 
## Call:
## lm(formula = mobilityjakarta$Station ~ mobilityjakarta$Tanggal, 
##     data = mobilityjakarta)
## 
## Coefficients:
##             (Intercept)  mobilityjakarta$Tanggal  
##              -3.701e+03                2.245e-06

VieW Data Prediksi linier model Mobilitas di Terminal/stasiun

mobilityjakarta$prediksi_model.linear <- model.linear$fitted.values
mobilityjakarta$prediksi_model.linear
##         1         2         3         4         5         6         7         8 
## -40.34624 -40.15224 -39.95825 -39.76426 -39.57026 -39.37627 -39.18228 -38.98828 
##         9        10        11        12        13        14        15        16 
## -38.79429 -38.60030 -38.40630 -38.21231 -38.01832 -37.82432 -37.63033 -37.43634 
##        17        18        19        20        21        22        23        24 
## -37.24234 -37.04835 -36.85436 -36.66036 -36.46637 -36.27238 -36.07838 -35.88439 
##        25        26        27        28        29        30 
## -35.69040 -35.49640 -35.30241 -35.10842 -34.91442 -34.72043

6.2 Smooth Spline Model Mobilitas di Terminal/Stasiun

# fit using smooth.spline
model.sspline <- smooth.spline(mobilityjakarta$Tanggal, mobilityjakarta$Station, nknots = 10)
model.sspline
## Call:
## smooth.spline(x = mobilityjakarta$Tanggal, y = mobilityjakarta$Station, 
##     nknots = 10)
## 
## Smoothing Parameter  spar= 0.03974373  lambda= 4.743395e-06 (16 iterations)
## Equivalent Degrees of Freedom (Df): 11.24592
## Penalized Criterion (RSS): 108.0067
## GCV: 9.21256

VieW Data Prediksi Smooth Spline Model Kasus Dirawat Covid-19

mobilityjakarta$prediksi_model.sspline <- model.sspline$y
mobilityjakarta$prediksi_model.sspline
##  [1] -41.46553 -41.76956 -39.20704 -37.10966 -37.86240 -40.06354 -41.36470
##  [8] -40.14184 -37.66834 -35.94193 -36.38456 -38.11505 -39.67645 -39.94961
## [15] -39.16659 -37.89725 -36.71146 -36.07978 -36.07560 -36.67298 -37.69930
## [22] -38.39501 -37.85385 -35.63092 -33.12673 -32.20317 -33.97623 -36.57839
## [29] -37.39621 -33.81630

7. Plot Regresi Linear Mobility Index DKI Jakarta di Terminal/Stasiun

Mobility Index di Terminal/stasiun (station) akan di bandingkan dengan Data Real, Linear Model dan Smooth Spline Model seperti pada gambar dibawah ini

Dari visualisasi diatas dengan menggunakan Smooth Spline Model prediksi mobilitas di tempat kerja antara -45 dan -35 dengan kurva naik yang lebih stabil. Prediksi ini dapat menggambarkan bahwa mobilitas di tempat kerja mempunyai tren kenaikan selama bulan September 2020

8. Hasil Regresi antara Kasus Dirawat Covid-19 dengan Mobility Index Terminal/Stasiun

# fit using ss
mod.ss <- ss( mobilityjakarta$Tanggal,mobilityjakarta$Dirawat, nknots = 10)
mod.ss
## 
## Call:
## ss(x = mobilityjakarta$Tanggal, y = mobilityjakarta$Dirawat, 
##     nknots = 10)
## 
## Smoothing Parameter  spar = 0.1836522   lambda = 1.265074e-06
## Equivalent Degrees of Freedom (Df) 9.294061
## Penalized Criterion (RSS) 18187.95
## Generalized Cross-Validation (GCV) 1272.668
# fit using smooth.spline
mod.smsp <- smooth.spline( mobilityjakarta$Tanggal,mobilityjakarta$Dirawat, nknots = 10)
mod.smsp
## Call:
## smooth.spline(x = mobilityjakarta$Tanggal, y = mobilityjakarta$Dirawat, 
##     nknots = 10)
## 
## Smoothing Parameter  spar= 0.1909029  lambda= 5.863705e-05 (14 iterations)
## Equivalent Degrees of Freedom (Df): 8.959171
## Penalized Criterion (RSS): 19714.47
## GCV: 1335.921
# rmse between solutions
sqrt(mean(( mod.ss$y - mod.smsp$y )^2))
## [1] 3.347557
# rmse between solutions and f(x)
sqrt(mean(( mobilityjakarta$Dirawat - mod.ss$y )^2))
## [1] 24.62245
sqrt(mean(( mobilityjakarta$Dirawat - mod.smsp$y )^2))
## [1] 25.63492
mobilityjakarta$prediksi_ss <- mod.ss$y
mobilityjakarta$prediksi_smsp <- mod.smsp$y
mobilityjakarta
# plot results
plot(mobilityjakarta$Tanggal, mobilityjakarta$Dirawat)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_ss , lty = 2, col = 'black', lwd = 4)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_smsp, lty = 4, col = 'red', lwd = 2)

# plot method
plot(mobilityjakarta$Tanggal, mobilityjakarta$Dirawat, lty = 10, col = 'black', lwd =5)
#plot(mod.ss)
# add lm fit
abline(coef(lm( mobilityjakarta$Dirawat ~ mobilityjakarta$Tanggal  , data = mobilityjakarta)), lty = 10, col = 'blue', lwd =5)
rug(mobilityjakarta$Tanggal)  # add rug to plot
#points(mobilityjakarta$Tanggal, mobilityjakarta$DIRAWAT, lty = 2, col = 2, lwd = 2)
points(mobilityjakarta$Tanggal,mod.ss$y , lty = 2, col = 'red', lwd = 4)
#points(mobilityjakarta$Tanggal,mod.smsp$y , lty = 2, col = 4, lwd = 2)
legend("topright", 
       legend = c("Real", "Model SS", "Trends"), 
       lty = 1:3, col = c("Black","Red","Blue"), lwd = 3, bty = "p")

Dari Visualisasi diatas Smooth Spline Model hampir mendekati Data Real, sedangkan Linear Model menunjukkan prediksi negatif atau nilai kasus semakin turun.

# Function to add histograms
panel.hist <- function(x, ...) {
    usr <- par("usr")
    on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5))
    his <- hist(x, plot = FALSE)
    breaks <- his$breaks
    nB <- length(breaks)
    y <- his$counts
    y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col = rgb(0, 1, 1, alpha = 0.5), ...)
    # lines(density(x), col = 2, lwd = 2) # Uncomment to add density lines
}
# Equivalent with a formula
pairs(mobilityjakarta$Dirawat~mobilityjakarta$Station, data = mobilityjakarta,
      upper.panel = NULL,         # Disabling the upper panel
      diag.panel = panel.hist)    # Adding the histograms

model <- lm(mobilityjakarta$Dirawat ~ mobilityjakarta$Station, data = mobilityjakarta)
summary(model)
## 
## Call:
## lm(formula = mobilityjakarta$Dirawat ~ mobilityjakarta$Station, 
##     data = mobilityjakarta)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -574.56 -390.81  -58.22  249.36  853.11 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)   
## (Intercept)             -1746.79     952.57  -1.834  0.07734 . 
## mobilityjakarta$Station   -77.33      25.29  -3.058  0.00486 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 440.5 on 28 degrees of freedom
## Multiple R-squared:  0.2504, Adjusted R-squared:  0.2236 
## F-statistic: 9.352 on 1 and 28 DF,  p-value: 0.004864
model$coefficients
##             (Intercept) mobilityjakarta$Station 
##              -1746.7915                -77.3337
model$fitted.values
##         1         2         3         4         5         6         7         8 
## 1423.8902 1501.2239 1423.8902  959.8880 1191.8891 1423.8902 1346.5565 1346.5565 
##         9        10        11        12        13        14        15        16 
## 1269.2228 1269.2228  805.2206 1037.2217 1269.2228 1655.8913 1191.8891 1269.2228 
##        17        18        19        20        21        22        23        24 
## 1191.8891  727.8869 1037.2217 1191.8891 1191.8891 1269.2228 1191.8891 1114.5554 
##        25        26        27        28        29        30 
##  495.8858  727.8869 1114.5554 1114.5554 1037.2217  882.5543
mobilityjakarta$prediksi_linear <- model$fitted.values
# plot method
plot(mobilityjakarta$Tanggal, mobilityjakarta$Dirawat, lty = 10, col = 'black', lwd =5)
#plot(mod.ss)
# add lm fit
abline(coef(lm( mobilityjakarta$Dirawat ~ mobilityjakarta$Tanggal  , data = mobilityjakarta)), lty = 10, col = 'blue', lwd =5)
rug(mobilityjakarta$Tanggal)  # add rug to plot
#points(mobilityjakarta$Tanggal, mobilityjakarta$DIRAWAT, lty = 2, col = 2, lwd = 2)
lines(mobilityjakarta$Tanggal,mod.ss$y , lty = 2, col = 'red', lwd = 4)
#points(mobilityjakarta$Tanggal,mod.smsp$y , lty = 2, col = 4, lwd = 2)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_linear, lty = 2, col = 6, lwd = 4)
legend("topright", 
       legend = c("Real", "Model SS", "Trends", "Linear Regresi"), 
       lty = 1:4, col = c("Black","Red","Blue","Pink"), lwd = 3, bty = "p")

Dari visualisasi diatas terlihat bahwa prediksi kasus dengan smooth splines adalah prediksi terbaik yang paling mendekati data real, Dari nilai trend didapatkan bahwa jumlah kasus dirawat selama bulan September 2020 terus megalami penurunan.

9. Korelasi Pearson

Korelasi Pearson adalah alat analisis statistik yang digunakan untuk melihat keeratan hubungan linier antara 2 variabel yang skala datanya adalah interval atau rasio.

Korelasi Pearson Untuk melakukan analisis korelasi Pearson cukup sederhana. Semisal kita ingin mengetahui hubungan antara Kasus Positif Covid-19 dan Mobilitas di Tempat Kerja (Workplace)

cor.test(mobilityjakarta$Dirawat,mobilityjakarta$Station)
## 
##  Pearson's product-moment correlation
## 
## data:  mobilityjakarta$Dirawat and mobilityjakarta$Station
## t = -3.058, df = 28, p-value = 0.004864
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.729188 -0.170907
## sample estimates:
##        cor 
## -0.5003673

10. Kesimpulan

  • Bahwa prediksi mobilitas di Terminal/stasiun cenderung mengalamai kenaikan selama bulan September.
  • Berdasarkan nilai korelasi diapatkan bahwa nilai mobilitas di Terminal/stasiun dengan kasus dirawat mempunyai hubungan yang lemah.