Dosen Pengampu : Prof. Dr. Suhartono, M.Kom

Lembaga : Universitas Islam Negeri Maulana Malik Ibrahim Malang

Jurusan : Teknik Informatika

Fakultas : Sains dan Teknologi

1. Pengertian smoothing spline

Smoothing splines adalah pendekatan yang ampuh untuk memperkirakan hubungan fungsional antara prediktor X dan respons Y. Smoothing splines dapat di-fit baik menggunakan fungsi smooth.spline (dalam paket stats) atau fungsi ss (dalam paket npreg). Dokumen ini memberikan latar belakang teoritis tentang smoothing splines, serta contoh-contoh yang menggambarkan bagaimana menggunakan fungsi smooth.spline dan ss. Seperti yang saya tunjukkan dalam tutorial ini, kedua fungsi memiliki sintaks yang sangat mirip, tetapi fungsi ss menawarkan beberapa opsi tambahan.

Regresi spline merupakan smoothing untuk memplot data dengan mempertimbangkan kemulusan kurva. Spline adalah model polinomial yang tersegmentasi atau terbagi, dan sifat segmen ini memberikan fleksibilitas yang lebih besar daripada model polinomial biasa. Properti ini memungkinkan model regresi spline untuk secara efektif disesuaikan dengan properti lokal data. Penggunaan splines menitikberatkan pada adanya perilaku atau pola data yang memiliki sifat yang berbeda pada suatu area tertentu dengan pada area lainnya. Berikut regresi nonparametrik dengan pendekatan smoothing spline pada data Google Mobility Index dan Covid-19 di Jakarta Januari 2021.

2. Data Google Mobility Index dan Covid-19 di Jakarta Januari 2021

library(readxl)
## Warning: package 'readxl' was built under R version 4.1.2
mobilityjakarta <- read_excel(path = "~/linear algebra/DataSembuhJanuari2021.xlsx")
mobilityjakarta
## # A tibble: 30 x 13
##    Tanggal             Nama_provinsi POSITIF Dirawat Sembuh Meninggal
##    <dttm>              <chr>           <dbl>   <dbl>  <dbl>     <dbl>
##  1 2021-01-01 00:00:00 DKI JAKARTA    185691    5789 166512      3308
##  2 2021-01-02 00:00:00 DKI JAKARTA    187586    4599 168781      3334
##  3 2021-01-03 00:00:00 DKI JAKARTA    189243    4410 170510      3345
##  4 2021-01-04 00:00:00 DKI JAKARTA    191075    4299 173036      3369
##  5 2021-01-05 00:00:00 DKI JAKARTA    192899    4479 174131      3392
##  6 2021-01-06 00:00:00 DKI JAKARTA    195301    4254 175441      3410
##  7 2021-01-07 00:00:00 DKI JAKARTA    197699    4601 176882      3435
##  8 2021-01-08 00:00:00 DKI JAKARTA    200658    4236 179562      3463
##  9 2021-01-09 00:00:00 DKI JAKARTA    203411    5237 181613      3485
## 10 2021-01-10 00:00:00 DKI JAKARTA    206122    4863 184576      3517
## # ... with 20 more rows, and 7 more variables: `Self Isolation` <dbl>,
## #   retail_and_recreation_percent_change_from_baseline <dbl>,
## #   grocery_and_pharmacy_percent_change_from_baseline <dbl>,
## #   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>
plot(mobilityjakarta$Tanggal,mobilityjakarta$POSITIF)

3 Data Deskriptif

library(jmv)
## Warning: package 'jmv' was built under R version 4.1.3
# Mendapatkan data descriptive menggunakan fungsi descritptive
descriptives(mobilityjakarta, vars = vars(POSITIF, grocery_and_pharmacy_percent_change_from_baseline), freq = TRUE)
## 
##  DESCRIPTIVES
## 
##  Descriptives                                                                            
##  --------------------------------------------------------------------------------------- 
##                          POSITIF     grocery_and_pharmacy_percent_change_from_baseline   
##  --------------------------------------------------------------------------------------- 
##    N                           30                                                   30   
##    Missing                      0                                                    0   
##    Mean                  223397.5                                            -8.733333   
##    Median                222202.0                                            -10.50000   
##    Standard deviation    25404.08                                             8.016936   
##    Minimum               185691.0                                            -20.00000   
##    Maximum               266244.0                                             5.000000   
##  ---------------------------------------------------------------------------------------

4. Perbandingan Pendekatan

library(npreg)
## Warning: package 'npreg' was built under R version 4.1.3
mod.ss <- ss(mobilityjakarta$Tanggal,mobilityjakarta$POSITIF, nknots = 10)
mobilityjakarta$prediksi_ss <- mod.ss$y
mobilityjakarta
## # A tibble: 30 x 14
##    Tanggal             Nama_provinsi POSITIF Dirawat Sembuh Meninggal
##    <dttm>              <chr>           <dbl>   <dbl>  <dbl>     <dbl>
##  1 2021-01-01 00:00:00 DKI JAKARTA    185691    5789 166512      3308
##  2 2021-01-02 00:00:00 DKI JAKARTA    187586    4599 168781      3334
##  3 2021-01-03 00:00:00 DKI JAKARTA    189243    4410 170510      3345
##  4 2021-01-04 00:00:00 DKI JAKARTA    191075    4299 173036      3369
##  5 2021-01-05 00:00:00 DKI JAKARTA    192899    4479 174131      3392
##  6 2021-01-06 00:00:00 DKI JAKARTA    195301    4254 175441      3410
##  7 2021-01-07 00:00:00 DKI JAKARTA    197699    4601 176882      3435
##  8 2021-01-08 00:00:00 DKI JAKARTA    200658    4236 179562      3463
##  9 2021-01-09 00:00:00 DKI JAKARTA    203411    5237 181613      3485
## 10 2021-01-10 00:00:00 DKI JAKARTA    206122    4863 184576      3517
## # ... with 20 more rows, and 8 more variables: `Self Isolation` <dbl>,
## #   retail_and_recreation_percent_change_from_baseline <dbl>,
## #   grocery_and_pharmacy_percent_change_from_baseline <dbl>,
## #   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>, prediksi_ss <dbl>

5. Perbandingan Pendekatan

library(npreg)
mod.ss <- ss(mobilityjakarta$Tanggal,mobilityjakarta$POSITIF, nknots = 10)
mobilityjakarta$prediksi_ss <- mod.ss$y
mobilityjakarta
## # A tibble: 30 x 14
##    Tanggal             Nama_provinsi POSITIF Dirawat Sembuh Meninggal
##    <dttm>              <chr>           <dbl>   <dbl>  <dbl>     <dbl>
##  1 2021-01-01 00:00:00 DKI JAKARTA    185691    5789 166512      3308
##  2 2021-01-02 00:00:00 DKI JAKARTA    187586    4599 168781      3334
##  3 2021-01-03 00:00:00 DKI JAKARTA    189243    4410 170510      3345
##  4 2021-01-04 00:00:00 DKI JAKARTA    191075    4299 173036      3369
##  5 2021-01-05 00:00:00 DKI JAKARTA    192899    4479 174131      3392
##  6 2021-01-06 00:00:00 DKI JAKARTA    195301    4254 175441      3410
##  7 2021-01-07 00:00:00 DKI JAKARTA    197699    4601 176882      3435
##  8 2021-01-08 00:00:00 DKI JAKARTA    200658    4236 179562      3463
##  9 2021-01-09 00:00:00 DKI JAKARTA    203411    5237 181613      3485
## 10 2021-01-10 00:00:00 DKI JAKARTA    206122    4863 184576      3517
## # ... with 20 more rows, and 8 more variables: `Self Isolation` <dbl>,
## #   retail_and_recreation_percent_change_from_baseline <dbl>,
## #   grocery_and_pharmacy_percent_change_from_baseline <dbl>,
## #   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>, prediksi_ss <dbl>
mod.smsp <- smooth.spline( mobilityjakarta$Tanggal,mobilityjakarta$POSITIF, nknots = 10)
mobilityjakarta$prediksi_smsp <- mod.smsp$y
mobilityjakarta
## # A tibble: 30 x 15
##    Tanggal             Nama_provinsi POSITIF Dirawat Sembuh Meninggal
##    <dttm>              <chr>           <dbl>   <dbl>  <dbl>     <dbl>
##  1 2021-01-01 00:00:00 DKI JAKARTA    185691    5789 166512      3308
##  2 2021-01-02 00:00:00 DKI JAKARTA    187586    4599 168781      3334
##  3 2021-01-03 00:00:00 DKI JAKARTA    189243    4410 170510      3345
##  4 2021-01-04 00:00:00 DKI JAKARTA    191075    4299 173036      3369
##  5 2021-01-05 00:00:00 DKI JAKARTA    192899    4479 174131      3392
##  6 2021-01-06 00:00:00 DKI JAKARTA    195301    4254 175441      3410
##  7 2021-01-07 00:00:00 DKI JAKARTA    197699    4601 176882      3435
##  8 2021-01-08 00:00:00 DKI JAKARTA    200658    4236 179562      3463
##  9 2021-01-09 00:00:00 DKI JAKARTA    203411    5237 181613      3485
## 10 2021-01-10 00:00:00 DKI JAKARTA    206122    4863 184576      3517
## # ... with 20 more rows, and 9 more variables: `Self Isolation` <dbl>,
## #   retail_and_recreation_percent_change_from_baseline <dbl>,
## #   grocery_and_pharmacy_percent_change_from_baseline <dbl>,
## #   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>, prediksi_ss <dbl>, ...
# plot method
plot(mobilityjakarta$Tanggal, mobilityjakarta$POSITIF, lty = 10, col = 'red', lwd =10)
# plot(mod.ss)
# add lm fit
abline(coef(lm( mobilityjakarta$POSITIF ~ mobilityjakarta$Tanggal  , data = mobilityjakarta)), lty = 2.5, col = 'blue', lwd =7)
rug(mobilityjakarta$Tanggal)  # add rug to plot

points(mobilityjakarta$Tanggal,mod.ss$y , lty = 2, col = 'green', lwd = 5)

legend("bottomright", 
       legend = c("Real", "Model SS", "Trends"), 
       lty = 1:3, col = c("red","green","Blue"), lwd = 3, bty = "p")

library(npreg)
mod.ss1 <- ss(mobilityjakarta$Tanggal,mobilityjakarta$grocery_and_pharmacy_percent_change_from_baselineF, nknots = 10)
## Warning: Unknown or uninitialised column:
## `grocery_and_pharmacy_percent_change_from_baselineF`.
mobilityjakarta$prediksi_ss1 <- mod.ss1$y
mobilityjakarta
## # A tibble: 30 x 16
##    Tanggal             Nama_provinsi POSITIF Dirawat Sembuh Meninggal
##    <dttm>              <chr>           <dbl>   <dbl>  <dbl>     <dbl>
##  1 2021-01-01 00:00:00 DKI JAKARTA    185691    5789 166512      3308
##  2 2021-01-02 00:00:00 DKI JAKARTA    187586    4599 168781      3334
##  3 2021-01-03 00:00:00 DKI JAKARTA    189243    4410 170510      3345
##  4 2021-01-04 00:00:00 DKI JAKARTA    191075    4299 173036      3369
##  5 2021-01-05 00:00:00 DKI JAKARTA    192899    4479 174131      3392
##  6 2021-01-06 00:00:00 DKI JAKARTA    195301    4254 175441      3410
##  7 2021-01-07 00:00:00 DKI JAKARTA    197699    4601 176882      3435
##  8 2021-01-08 00:00:00 DKI JAKARTA    200658    4236 179562      3463
##  9 2021-01-09 00:00:00 DKI JAKARTA    203411    5237 181613      3485
## 10 2021-01-10 00:00:00 DKI JAKARTA    206122    4863 184576      3517
## # ... with 20 more rows, and 10 more variables: `Self Isolation` <dbl>,
## #   retail_and_recreation_percent_change_from_baseline <dbl>,
## #   grocery_and_pharmacy_percent_change_from_baseline <dbl>,
## #   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>, prediksi_ss <dbl>, ...
mod.smsp1 <- smooth.spline( mobilityjakarta$Tanggal,mobilityjakarta$grocery_and_pharmacy_percent_change_from_baseline, nknots = 10)
mobilityjakarta$prediksi_smsp1 <- mod.smsp1$y
mobilityjakarta
## # A tibble: 30 x 17
##    Tanggal             Nama_provinsi POSITIF Dirawat Sembuh Meninggal
##    <dttm>              <chr>           <dbl>   <dbl>  <dbl>     <dbl>
##  1 2021-01-01 00:00:00 DKI JAKARTA    185691    5789 166512      3308
##  2 2021-01-02 00:00:00 DKI JAKARTA    187586    4599 168781      3334
##  3 2021-01-03 00:00:00 DKI JAKARTA    189243    4410 170510      3345
##  4 2021-01-04 00:00:00 DKI JAKARTA    191075    4299 173036      3369
##  5 2021-01-05 00:00:00 DKI JAKARTA    192899    4479 174131      3392
##  6 2021-01-06 00:00:00 DKI JAKARTA    195301    4254 175441      3410
##  7 2021-01-07 00:00:00 DKI JAKARTA    197699    4601 176882      3435
##  8 2021-01-08 00:00:00 DKI JAKARTA    200658    4236 179562      3463
##  9 2021-01-09 00:00:00 DKI JAKARTA    203411    5237 181613      3485
## 10 2021-01-10 00:00:00 DKI JAKARTA    206122    4863 184576      3517
## # ... with 20 more rows, and 11 more variables: `Self Isolation` <dbl>,
## #   retail_and_recreation_percent_change_from_baseline <dbl>,
## #   grocery_and_pharmacy_percent_change_from_baseline <dbl>,
## #   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>, prediksi_ss <dbl>, ...
# Hasil plot

plot(mobilityjakarta$Tanggal, mobilityjakarta$grocery_and_pharmacy_percent_change_from_baseline, lty = 10, col = 'Red', lwd =5)

# add lm fit
abline(coef(lm( mobilityjakarta$grocery_and_pharmacy_percent_change_from_baseline ~ mobilityjakarta$Tanggal  , data = mobilityjakarta)), lty = 10, col = 'Blue', lwd =5)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_ss1 , lty = 2, col = 'yellow', lwd = 4)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_smsp1, lty = 4, col = 'green', lwd = 2)

6. 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.

cor.test(mobilityjakarta$POSITIF,mobilityjakarta$grocery_and_pharmacy_percent_change_from_baseline)
## 
##  Pearson's product-moment correlation
## 
## data:  mobilityjakarta$POSITIF and mobilityjakarta$grocery_and_pharmacy_percent_change_from_baseline
## t = -5.0335, df = 28, p-value = 2.531e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8407334 -0.4376168
## sample estimates:
##        cor 
## -0.6892232

Dari output tersebut dapat kita simpulkan bahwa adanya hubungan yang signifikan antara mobilityjakarta POSITIF dan mobilityjakarta grocery_and_pharmacy_percent_change_from_baseline (p-value<0,01). Nilai p-value dalam output R dituliskan 0.002834. Nilai koefisien korelasi r adalah sebesar 0.5180438 yang menunjukkan hubungan yang cukup kuat dan positif (berbanding lurus) antara variabel POSITIF dan grocery_and_pharmacy_percent_change_from_baseline.

7. Hasil Regresi dengan Pendekatan Smoothing Spline

model <- lm(mobilityjakarta$POSITIF ~ mobilityjakarta$grocery_and_pharmacy_percent_change_from_baseline, data = mobilityjakarta)
summary(model)
## 
## Call:
## lm(formula = mobilityjakarta$POSITIF ~ mobilityjakarta$grocery_and_pharmacy_percent_change_from_baseline, 
##     data = mobilityjakarta)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -36105 -10052  -1669   2134  44448 
## 
## Coefficients:
##                                                                   Estimate
## (Intercept)                                                       204323.8
## mobilityjakarta$grocery_and_pharmacy_percent_change_from_baseline  -2184.0
##                                                                   Std. Error
## (Intercept)                                                           5104.5
## mobilityjakarta$grocery_and_pharmacy_percent_change_from_baseline      433.9
##                                                                   t value
## (Intercept)                                                        40.028
## mobilityjakarta$grocery_and_pharmacy_percent_change_from_baseline  -5.034
##                                                                   Pr(>|t|)    
## (Intercept)                                                        < 2e-16 ***
## mobilityjakarta$grocery_and_pharmacy_percent_change_from_baseline 2.53e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18730 on 28 degrees of freedom
## Multiple R-squared:  0.475,  Adjusted R-squared:  0.4563 
## F-statistic: 25.34 on 1 and 28 DF,  p-value: 2.531e-05
mobilityjakarta$prediksi_model <- model$fitted.values
mobilityjakarta
## # A tibble: 30 x 18
##    Tanggal             Nama_provinsi POSITIF Dirawat Sembuh Meninggal
##    <dttm>              <chr>           <dbl>   <dbl>  <dbl>     <dbl>
##  1 2021-01-01 00:00:00 DKI JAKARTA    185691    5789 166512      3308
##  2 2021-01-02 00:00:00 DKI JAKARTA    187586    4599 168781      3334
##  3 2021-01-03 00:00:00 DKI JAKARTA    189243    4410 170510      3345
##  4 2021-01-04 00:00:00 DKI JAKARTA    191075    4299 173036      3369
##  5 2021-01-05 00:00:00 DKI JAKARTA    192899    4479 174131      3392
##  6 2021-01-06 00:00:00 DKI JAKARTA    195301    4254 175441      3410
##  7 2021-01-07 00:00:00 DKI JAKARTA    197699    4601 176882      3435
##  8 2021-01-08 00:00:00 DKI JAKARTA    200658    4236 179562      3463
##  9 2021-01-09 00:00:00 DKI JAKARTA    203411    5237 181613      3485
## 10 2021-01-10 00:00:00 DKI JAKARTA    206122    4863 184576      3517
## # ... with 20 more rows, and 12 more variables: `Self Isolation` <dbl>,
## #   retail_and_recreation_percent_change_from_baseline <dbl>,
## #   grocery_and_pharmacy_percent_change_from_baseline <dbl>,
## #   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>, prediksi_ss <dbl>, ...
# Menambahkan 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, 0, alpha = 0.5), ...)
    # lines(density(x), col = 2, lwd = 2) # Uncomment to add density lines
}
# Menyetarakan berdasarkan formula
pairs(mobilityjakarta$POSITIF~mobilityjakarta$grocery_and_pharmacy_percent_change_from_baseline, data = mobilityjakarta,
      upper.panel = NULL,         # Disabling the upper panel
      diag.panel = panel.hist)    # Adding the histograms

# plot method
plot(mobilityjakarta$Tanggal, mobilityjakarta$POSITIF, lty = 10, col = 'Red', lwd =5)

#plot(mod.ss)
# add lm fit
abline(coef(lm( mobilityjakarta$POSITIF ~ mobilityjakarta$Tanggal  , data = mobilityjakarta)), lty = 10, col = 'Blue', lwd =5)
rug(mobilityjakarta$Tanggal)  # add rug to plot

lines(mobilityjakarta$Tanggal,mod.ss$y , lty = 2, col = 'Orange', lwd = 4)

#plot(mod.smsp)

lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_model, lty = 2, col = 'Green', lwd = 4)
legend("topleft", 
       legend = c("Real", "Model SS", "Trends", "Regresi Nonparametrik"), 
       lty = 1:4, col = c("Red","Orange","Blue","Green"), lwd = 3, bty = "p")

Daftar Pustaka

https://rpubs.com/suhartono-uinmaliki/891461

http://users.stat.umn.edu/~helwig/notes/smooth-spline-notes.html