Regresi Linear 3 model (Smoothing Spline Model)
~ untuk Kasus Isolasi Mandiri Covid-19 DKI Jakarta dibandingkan dengan Mobilitas di Park/Ruang Terbuka Hijau ~
| Oleh | Ihsan Bagus Fahad Arafat |
| Dosen Pengampu | Prof. Dr. Suhartono, M.Kom |
| Program | Magister Informatika |
| Tanggal | 6 Juni 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.3mobilityjakarta <- read_excel(path = "Datagab.xlsx")
mobilityjakarta3. Aktivasi Library
library(npreg)## Warning: package 'npreg' was built under R version 4.1.3library(jmv)## Warning: package 'jmv' was built under R version 4.1.3library(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.3library(stats)
library(npreg)
descriptives(mobilityjakarta, vars = vars(Isoman, Park), freq = TRUE)##
## DESCRIPTIVES
##
## Descriptives
## -----------------------------------------------
## Isoman Park
## -----------------------------------------------
## N 31 31
## Missing 0 0
## Mean 983.0645 -30.77419
## Median 978.0000 -32.00000
## Standard deviation 209.1432 6.651865
## Minimum 618.0000 -46.00000
## Maximum 1285.000 -15.00000
## -----------------------------------------------4. Histrogram plot
hist(mobilityjakarta$Isoman, main="Distribusi Isolasi Mandiri untuk Covid 19", xlab="Isolasi Mandiri Untuk Covid 19", xlim=c(0,2500),ylim=c(0,15)) # histogram variabel wt
mean=mean(mobilityjakarta$Isoman) # Menghitung nilai mean variabel wt
abline(v=mean, col="red",lwd=2) # Menambahkan garis vertical pada plot
5. Density plot
d<-density(mobilityjakarta$Isoman) # menghitung density variabel wt
plot(d, xlab="Isolasi Mandiri untuk Covid 19", main="Distribusi Isolasi Mandiri untuk Covid 19")
6. Pemodelan Mobilitas
Pada pemodelan ini dilakukan dengan linier dan smooth spline
6.1 Linier Model Mobilitas di Park/Ruang Terbuka Hijau
# linear / tren
model.linear <- lm(mobilityjakarta$Park ~ mobilityjakarta$Tanggal, data = mobilityjakarta)
model.linear##
## Call:
## lm(formula = mobilityjakarta$Park ~ mobilityjakarta$Tanggal,
## data = mobilityjakarta)
##
## Coefficients:
## (Intercept) mobilityjakarta$Tanggal
## -7.094e+03 4.322e-06VieW Data Prediksi linier model Mobilitas di Park/Ruang Terbuka Hijau
mobilityjakarta$prediksi_model.linear <- model.linear$fitted.values
mobilityjakarta$prediksi_model.linear## 1 2 3 4 5 6 7 8
## -36.37500 -36.00161 -35.62823 -35.25484 -34.88145 -34.50806 -34.13468 -33.76129
## 9 10 11 12 13 14 15 16
## -33.38790 -33.01452 -32.64113 -32.26774 -31.89435 -31.52097 -31.14758 -30.77419
## 17 18 19 20 21 22 23 24
## -30.40081 -30.02742 -29.65403 -29.28065 -28.90726 -28.53387 -28.16048 -27.78710
## 25 26 27 28 29 30 31
## -27.41371 -27.04032 -26.66694 -26.29355 -25.92016 -25.54677 -25.173396.2 Smooth Spline Model Mobilitas di Park
# fit using smooth.spline
model.sspline <- smooth.spline(mobilityjakarta$Tanggal, mobilityjakarta$Park, nknots = 10)
model.sspline## Call:
## smooth.spline(x = mobilityjakarta$Tanggal, y = mobilityjakarta$Park,
## nknots = 10)
##
## Smoothing Parameter spar= 0.08753934 lambda= 1.121128e-05 (15 iterations)
## Equivalent Degrees of Freedom (Df): 10.63323
## Penalized Criterion (RSS): 391.4619
## GCV: 29.25547VieW Data Prediksi Smooth Spline Model Kasus Isolasi Mandiri Covid-19
mobilityjakarta$prediksi_model.sspline <- model.sspline$y
mobilityjakarta$prediksi_model.sspline## [1] -28.58376 -33.22636 -38.15322 -40.37724 -37.97096 -33.24555 -29.57185
## [8] -29.42120 -31.66703 -34.28328 -35.24390 -33.31304 -30.41566 -29.26694
## [15] -31.63952 -35.53596 -38.01629 -36.93151 -33.29642 -28.91680 -25.59842
## [22] -24.69844 -25.77959 -27.95598 -30.27661 -31.52996 -30.43939 -26.42179
## [29] -21.66823 -19.06335 -21.491757. Plot Regresi Linear Mobility Index DKI Jakarta di Park/Ruang Terbuka Hijau
Mobility Index di Park/Ruang Terbuka Hijau 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 park/taman antara -40 dan -20 dengan kurva naik yang lebih
stabil. Prediksi ini dapat menggambarkan bahwa mobilitas di Park/Ruang
Terbuka Hijau mempunyai tren kenaikan selama bulan Oktober 2021
8. Hasil Regresi antara Kasus Isolasi Mandiri Covid-19 dengan Mobility Index Park
# fit using ss
mod.ss <- ss( mobilityjakarta$Tanggal,mobilityjakarta$Isoman, nknots = 10)
mod.ss##
## Call:
## ss(x = mobilityjakarta$Tanggal, y = mobilityjakarta$Isoman, nknots = 10)
##
## Smoothing Parameter spar = 0.08869552 lambda = 2.606636e-07
## Equivalent Degrees of Freedom (Df) 10.35887
## Penalized Criterion (RSS) 39727.2
## Generalized Cross-Validation (GCV) 2890.563# fit using smooth.spline
mod.smsp <- smooth.spline( mobilityjakarta$Tanggal,mobilityjakarta$Isoman, nknots = 10)
mod.smsp## Call:
## smooth.spline(x = mobilityjakarta$Tanggal, y = mobilityjakarta$Isoman,
## nknots = 10)
##
## Smoothing Parameter spar= 0.2279273 lambda= 0.000115856 (14 iterations)
## Equivalent Degrees of Freedom (Df): 8.161813
## Penalized Criterion (RSS): 51536.66
## GCV: 3063.054# rmse between solutions
sqrt(mean(( mod.ss$y - mod.smsp$y )^2))## [1] 14.55239# rmse between solutions and f(x)
sqrt(mean(( mobilityjakarta$Isoman - mod.ss$y )^2))## [1] 35.79836sqrt(mean(( mobilityjakarta$Isoman - mod.smsp$y )^2))## [1] 40.77343mobilityjakarta$prediksi_ss <- mod.ss$y
mobilityjakarta$prediksi_smsp <- mod.smsp$y
mobilityjakarta# plot results
plot(mobilityjakarta$Tanggal, mobilityjakarta$Isoman)
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$Isoman, lty = 10, col = 'black', lwd =5)
#plot(mod.ss)
# add lm fit
abline(coef(lm( mobilityjakarta$Isoman ~ mobilityjakarta$Tanggal , data = mobilityjakarta)), lty = 10, col = 'blue', lwd =5)
rug(mobilityjakarta$Tanggal) # add rug to plot
#points(mobilityjakarta$Tanggal, mobilityjakarta$ISOMAN, 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$Isoman~mobilityjakarta$Park, data = mobilityjakarta,
upper.panel = NULL, # Disabling the upper panel
diag.panel = panel.hist) # Adding the histograms
model <- lm(mobilityjakarta$Isoman ~ mobilityjakarta$Park, data = mobilityjakarta)
summary(model)##
## Call:
## lm(formula = mobilityjakarta$Isoman ~ mobilityjakarta$Park, data = mobilityjakarta)
##
## Residuals:
## Min 1Q Median 3Q Max
## -397.02 -92.77 -22.66 146.84 305.83
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 541.247 163.424 3.312 0.00249 **
## mobilityjakarta$Park -14.357 5.194 -2.764 0.00982 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 189.2 on 29 degrees of freedom
## Multiple R-squared: 0.2085, Adjusted R-squared: 0.1812
## F-statistic: 7.639 on 1 and 29 DF, p-value: 0.009819model$coefficients## (Intercept) mobilityjakarta$Park
## 541.24656 -14.35677model$fitted.values## 1 2 3 4 5 6 7 8
## 957.5928 1000.6631 1101.1605 1101.1605 1101.1605 1043.7334 1000.6631 900.1658
## 9 10 11 12 13 14 15 16
## 971.9496 1043.7334 1043.7334 1043.7334 1015.0199 1000.6631 900.1658 1000.6631
## 17 18 19 20 21 22 23 24
## 1086.8037 1201.6579 1043.7334 814.0252 986.3064 885.8090 885.8090 928.8793
## 25 26 27 28 29 30 31
## 1000.6631 1000.6631 1015.0199 928.8793 756.5981 857.0955 857.0955mobilityjakarta$prediksi_linear <- model$fitted.values# plot method
plot(mobilityjakarta$Tanggal, mobilityjakarta$Isoman, lty = 10, col = 'black', lwd =5)
#plot(mod.ss)
# add lm fit
abline(coef(lm( mobilityjakarta$Isoman ~ mobilityjakarta$Tanggal , data = mobilityjakarta)), lty = 10, col = 'blue', lwd =5)
rug(mobilityjakarta$Tanggal) # add rug to plot
#points(mobilityjakarta$Tanggal, mobilityjakarta$ISOMAN, 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 isolasi mandiri selama bulan
Oktober 2021 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$Isoman,mobilityjakarta$Park)##
## Pearson's product-moment correlation
##
## data: mobilityjakarta$Isoman and mobilityjakarta$Park
## t = -2.764, df = 29, p-value = 0.009819
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6980222 -0.1220251
## sample estimates:
## cor
## -0.456621610. Kesimpulan
- Bahwa prediksi mobilitas di Park/Ruang Terbuka Hijau cenderung mengalamai kenaikan selama bulan Oktober.
- Berdasarkan nilai korelasi diapatkan bahwa nilai mobilitas di Park dengan kasus isolasi mandiri mempunyai hubungan yang lemah.
11. Referensi
https://github.com/juba/rmdformats : Format Visual Rpubs
https://medium.com/machine-learning-id/simple-linear-regression-r-ea42aed17c52
http://users.stat.umn.edu/~helwig/notes/smooth-spline-notes.html
https://www.rumusstatistik.com/2019/06/korelasi-pearson.html
https://rpubs.com/suhartono-uinmaliki/891461