Penerapan Analisa Data dengan Regresi Linier Menggunakan Pendekatan smoothing spline pada Data Google Mobility Index Sembuh Covid-19 di DKI Jakarta Bulan Juli 2021
Nurul Izzah, Prof. Dr. Suhartono, M.Kom
7 Juni 2022
Lembaga : Universitas Islam Negeri Maulana Malik Ibrahim Malang
Jurusan : Teknik Informatika
1. Smoothing splines
Smoothing splines adalah pendekatan yang ampuh untuk memperkirakan hubungan fungsional antara prediktor X dan respons Y. Smoothing splines dapat ditampung 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. Dibandingkan dengan fungsi smooth.spline, fungsi ss memiliki lebih banyak metode pemilihan parameter penghalusan lebih banyak jenis spline (linier, kubik, quintic) opsi untuk batasan periodisitas opsi untuk nilai simpul yang ditentukan pengguna ringkasan yang sesuai dan metode plot.
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 Sembuh Covid-19 di DKI Jakarta bulan Juli 2021.
2. Data Google Mobility Index Sembuh Covid-19 di DKI Jakarta pada Bulan Juli 2021
library(readxl)## Warning: package 'readxl' was built under R version 4.1.2mobilityjakarta <- read_excel(path = "data sembuh covid-19 DKI Jakarta Juli 2021.xlsx")## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'tibble'## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'pillar'mobilityjakarta## # A tibble: 31 x 8
## Tanggal Sembuh retail_and_recr~ grocery_and_pha~ parks_percent_c~
## <dttm> <dbl> <dbl> <dbl> <dbl>
## 1 2021-07-01 00:00:00 468461 -42 -5 -67
## 2 2021-07-02 00:00:00 473467 -32 2 -51
## 3 2021-07-03 00:00:00 479150 -34 -1 -53
## 4 2021-07-04 00:00:00 484949 -34 -2 -52
## 5 2021-07-05 00:00:00 491556 -37 -6 -55
## 6 2021-07-06 00:00:00 497492 -31 2 -45
## 7 2021-07-07 00:00:00 501199 -36 2 -56
## 8 2021-07-08 00:00:00 512085 -39 -5 -63
## 9 2021-07-09 00:00:00 526941 -32 -2 -50
## 10 2021-07-10 00:00:00 543867 -35 -6 -55
## # ... with 21 more rows, and 3 more variables:
## # transit_stations_percent_change_from_baseline <dbl>,
## # workplaces_percent_change_from_baseline <dbl>,
## # residential_percent_change_from_baseline <dbl>plot(mobilityjakarta$Tanggal,mobilityjakarta$Sembuh)
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(Sembuh, retail_and_recreation_percent_change_from_baseline), freq = TRUE)##
## DESCRIPTIVES
##
## Descriptives
## ----------------------------------------------------------------------------------------
## Sembuh retail_and_recreation_percent_change_from_baseline
## ----------------------------------------------------------------------------------------
## N 31 31
## Missing 0 0
## Mean 614735.5 -31.16129
## Median 604060.0 -31.00000
## Standard deviation 101497.5 4.831127
## Minimum 468461.0 -42.00000
## Maximum 784668.0 -22.00000
## ----------------------------------------------------------------------------------------hist(mobilityjakarta$Sembuh, main="Distribusi Sembuh Covid 19", xlab="Sembuh Covid 19", xlim=c(3000,6000),ylim=c(0,15)) # histogram variabel wt
mean=mean(mobilityjakarta$Sembuh) # Menghitung nilai mean variabel wt
abline(v=mean, col="red",lwd=2) # Menambahkan garis vertical pada plot
# 2. Density plot
d<-density(mobilityjakarta$Sembuh) # menghitung density variabel wt
plot(d, xlab="Sembuh Covid 19", main="Distribusi Sembuh Covid 19")
4. Perbandingan Pendekatan
library(npreg)## Warning: package 'npreg' was built under R version 4.1.3# fit using ss
mod.ss <- ss( mobilityjakarta$Tanggal,mobilityjakarta$Sembuh, nknots = 10)
mod.ss##
## Call:
## ss(x = mobilityjakarta$Tanggal, y = mobilityjakarta$Sembuh, nknots = 10)
##
## Smoothing Parameter spar = 0.1317937 lambda = 5.338915e-07
## Equivalent Degrees of Freedom (Df) 9.937134
## Penalized Criterion (RSS) 195331500
## Generalized Cross-Validation (GCV) 13648945# fit using smooth.spline
mod.smsp <- smooth.spline( mobilityjakarta$Tanggal,mobilityjakarta$Sembuh, nknots = 10)
mod.smsp## Call:
## smooth.spline(x = mobilityjakarta$Tanggal, y = mobilityjakarta$Sembuh,
## nknots = 10)
##
## Smoothing Parameter spar= 0.09125248 lambda= 1.191242e-05 (16 iterations)
## Equivalent Degrees of Freedom (Df): 10.58322
## Penalized Criterion (RSS): 179058131
## GCV: 13316234# rmse between solutions
sqrt(mean(( mod.ss$y - mod.smsp$y )^2))## [1] 386.8594# rmse between solutions and f(x)
sqrt(mean(( mobilityjakarta$Sembuh - mod.ss$y )^2))## [1] 2510.182sqrt(mean(( mobilityjakarta$Sembuh - mod.smsp$y )^2))## [1] 2403.345mobilityjakarta$prediksi_ss <- mod.ss$y
mobilityjakarta$prediksi_smsp <- mod.smsp$y
mobilityjakarta## # A tibble: 31 x 10
## Tanggal Sembuh retail_and_recr~ grocery_and_pha~ parks_percent_c~
## <dttm> <dbl> <dbl> <dbl> <dbl>
## 1 2021-07-01 00:00:00 468461 -42 -5 -67
## 2 2021-07-02 00:00:00 473467 -32 2 -51
## 3 2021-07-03 00:00:00 479150 -34 -1 -53
## 4 2021-07-04 00:00:00 484949 -34 -2 -52
## 5 2021-07-05 00:00:00 491556 -37 -6 -55
## 6 2021-07-06 00:00:00 497492 -31 2 -45
## 7 2021-07-07 00:00:00 501199 -36 2 -56
## 8 2021-07-08 00:00:00 512085 -39 -5 -63
## 9 2021-07-09 00:00:00 526941 -32 -2 -50
## 10 2021-07-10 00:00:00 543867 -35 -6 -55
## # ... with 21 more rows, and 5 more variables:
## # transit_stations_percent_change_from_baseline <dbl>,
## # workplaces_percent_change_from_baseline <dbl>,
## # residential_percent_change_from_baseline <dbl>, prediksi_ss <dbl>,
## # prediksi_smsp <dbl># plot results
plot(mobilityjakarta$Tanggal, mobilityjakarta$Sembuh)
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$Sembuh, lty = 10, col = 'black', lwd =5)
#plot(mod.ss)
# add lm fit
abline(coef(lm( mobilityjakarta$Sembuh ~ mobilityjakarta$Tanggal , data = mobilityjakarta)), lty = 10, col = 'blue', lwd =5)
rug(mobilityjakarta$Tanggal) # add rug to plot
#points(mobilityjakarta$Tanggal, mobilityjakarta$Sembuh, 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("purple","pink","light blue"), lwd = 3, bty = "p")
library(npreg)
# fit using ss
mod.ss1 <- ss( mobilityjakarta$Tanggal,mobilityjakarta$retail_and_recreation_percent_change_from_baseline, nknots = 10)
mod.ss1##
## Call:
## ss(x = mobilityjakarta$Tanggal, y = mobilityjakarta$retail_and_recreation_percent_change_from_baseline,
## nknots = 10)
##
## Smoothing Parameter spar = 0.6931016 lambda = 0.006063805
## Equivalent Degrees of Freedom (Df) 2.331711
## Penalized Criterion (RSS) 282.6364
## Generalized Cross-Validation (GCV) 10.66071# fit using smooth.spline
mod.smsp1 <- smooth.spline( mobilityjakarta$Tanggal,mobilityjakarta$retail_and_recreation_percent_change_from_baseline, nknots = 10)
mod.smsp1## Call:
## smooth.spline(x = mobilityjakarta$Tanggal, y = mobilityjakarta$retail_and_recreation_percent_change_from_baseline,
## nknots = 10)
##
## Smoothing Parameter spar= 0.6722594 lambda= 0.1879722 (14 iterations)
## Equivalent Degrees of Freedom (Df): 2.331817
## Penalized Criterion (RSS): 282.6331
## GCV: 10.66067mobilityjakarta$prediksi_ss1 <- mod.ss1$y
mobilityjakarta$prediksi_smsp1 <- mod.smsp1$y
mobilityjakarta## # A tibble: 31 x 12
## Tanggal Sembuh retail_and_recr~ grocery_and_pha~ parks_percent_c~
## <dttm> <dbl> <dbl> <dbl> <dbl>
## 1 2021-07-01 00:00:00 468461 -42 -5 -67
## 2 2021-07-02 00:00:00 473467 -32 2 -51
## 3 2021-07-03 00:00:00 479150 -34 -1 -53
## 4 2021-07-04 00:00:00 484949 -34 -2 -52
## 5 2021-07-05 00:00:00 491556 -37 -6 -55
## 6 2021-07-06 00:00:00 497492 -31 2 -45
## 7 2021-07-07 00:00:00 501199 -36 2 -56
## 8 2021-07-08 00:00:00 512085 -39 -5 -63
## 9 2021-07-09 00:00:00 526941 -32 -2 -50
## 10 2021-07-10 00:00:00 543867 -35 -6 -55
## # ... with 21 more rows, and 7 more variables:
## # transit_stations_percent_change_from_baseline <dbl>,
## # workplaces_percent_change_from_baseline <dbl>,
## # residential_percent_change_from_baseline <dbl>, prediksi_ss <dbl>,
## # prediksi_smsp <dbl>, prediksi_ss1 <dbl>, prediksi_smsp1 <dbl># plot results
plot(mobilityjakarta$Tanggal, mobilityjakarta$retail_and_recreation_percent_change_from_baseline, lty = 10, col = 'blue', lwd =5)
# add lm fit
abline(coef(lm( mobilityjakarta$retail_and_recreation_percent_change_from_baseline ~ mobilityjakarta$Tanggal , data = mobilityjakarta)), lty = 10, col = 'purple', lwd =5)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_ss1 , lty = 2, col = 'pink', lwd = 4)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_smsp1, lty = 4, col = 'orange', lwd = 2)
5. 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$Sembuh, mobilityjakarta$retail_and_recreation_percent_change_from_baseline)##
## Pearson's product-moment correlation
##
## data: mobilityjakarta$Sembuh and mobilityjakarta$retail_and_recreation_percent_change_from_baseline
## t = 6.7568, df = 29, p-value = 2.045e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5916026 0.8897917
## sample estimates:
## cor
## 0.7820109Dari output tersebut dapat kita simpulkan bahwa adanya hubungan yang signifikan antara mobilityjakarta Sembuh Covid-19 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 Sembuh dan retail_and_recreation_percent_change_from_baseline.
6. Hasil Regresi dengan Pendekatan Smoothing Splines
model <- lm(mobilityjakarta$Sembuh ~ mobilityjakarta$retail_and_recreation_percent_change_from_baseline, data = mobilityjakarta)
summary(model)##
## Call:
## lm(formula = mobilityjakarta$Sembuh ~ mobilityjakarta$retail_and_recreation_percent_change_from_baseline,
## data = mobilityjakarta)
##
## Residuals:
## Min 1Q Median 3Q Max
## -127489 -41471 15543 36479 118318
##
## Coefficients:
## Estimate
## (Intercept) 1126694
## mobilityjakarta$retail_and_recreation_percent_change_from_baseline 16429
## Std. Error
## (Intercept) 76646
## mobilityjakarta$retail_and_recreation_percent_change_from_baseline 2432
## t value
## (Intercept) 14.700
## mobilityjakarta$retail_and_recreation_percent_change_from_baseline 6.757
## Pr(>|t|)
## (Intercept) 5.65e-15 ***
## mobilityjakarta$retail_and_recreation_percent_change_from_baseline 2.05e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 64340 on 29 degrees of freedom
## Multiple R-squared: 0.6115, Adjusted R-squared: 0.5981
## F-statistic: 45.65 on 1 and 29 DF, p-value: 2.045e-07model$coefficients## (Intercept)
## 1126694.32
## mobilityjakarta$retail_and_recreation_percent_change_from_baseline
## 16429.32model$fitted.values## 1 2 3 4 5 6 7 8
## 436662.8 600956.1 568097.4 568097.4 518809.4 617385.4 535238.8 485950.8
## 9 10 11 12 13 14 15 16
## 600956.1 551668.1 568097.4 633814.7 650244.0 584526.7 535238.8 683102.7
## 17 18 19 20 21 22 23 24
## 518809.4 600956.1 617385.4 666673.3 617385.4 568097.4 666673.3 633814.7
## 25 26 27 28 29 30 31
## 699532.0 715961.3 765249.3 715961.3 650244.0 748819.9 732390.6mobilityjakarta$prediksi_linear <- model$fitted.values
mobilityjakarta## # A tibble: 31 x 13
## Tanggal Sembuh retail_and_recr~ grocery_and_pha~ parks_percent_c~
## <dttm> <dbl> <dbl> <dbl> <dbl>
## 1 2021-07-01 00:00:00 468461 -42 -5 -67
## 2 2021-07-02 00:00:00 473467 -32 2 -51
## 3 2021-07-03 00:00:00 479150 -34 -1 -53
## 4 2021-07-04 00:00:00 484949 -34 -2 -52
## 5 2021-07-05 00:00:00 491556 -37 -6 -55
## 6 2021-07-06 00:00:00 497492 -31 2 -45
## 7 2021-07-07 00:00:00 501199 -36 2 -56
## 8 2021-07-08 00:00:00 512085 -39 -5 -63
## 9 2021-07-09 00:00:00 526941 -32 -2 -50
## 10 2021-07-10 00:00:00 543867 -35 -6 -55
## # ... with 21 more rows, and 8 more variables:
## # transit_stations_percent_change_from_baseline <dbl>,
## # workplaces_percent_change_from_baseline <dbl>,
## # residential_percent_change_from_baseline <dbl>, prediksi_ss <dbl>,
## # prediksi_smsp <dbl>, prediksi_ss1 <dbl>, prediksi_smsp1 <dbl>,
## # prediksi_linear <dbl># 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$Sembuh~mobilityjakarta$retail_and_recreation_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$Sembuh, lty = 10, col = 'black', lwd =5)
#plot(mod.ss)
# add lm fit
abline(coef(lm( mobilityjakarta$Sembuh ~ 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")
Sumber:
https://rpubs.com/suhartono-uinmaliki/891461
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00