Analisa Data dengan Regresi Nonparametrik Menggunakan Pendekatan Smooting Spline Pada Data Google Mobility Index dan Covid-19 di Jakarta Juni 2020
Vivin Octavia Cahyani
4/20/2022
Dosen Pengampu : Prof. Dr. Suhartono, M.Kom
Lembaga : Universitas Islam Negeri Maulana Malik Ibrahim Malang
Jurusan : Teknik Informatika
Fakultas : Sains dan Teknologi
Pengertian Smooting 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
Untuk melihat bagaimana kinerja fungsi-fungsi ini dalam praktiknya, mari kita lihat contoh simulasinya. Secara khusus, mari simulasikan beberapa data dengan hubungan fungsional (berkala) yang memiliki beberapa gangguan. Berikut Analisa Data dengan Regresi Nonparametrik dengan Pendekatan Smooting Spline Pada Data Google Mobility Index dan Covid-19 di Jakarta Juni 2020.
Data Google Mobility Index dan Covid-19 di Jakarta Juni 2020
## Import the data
library(readxl)## Warning: package 'readxl' was built under R version 4.1.2mobilityjakarta <- read_excel(path = "dataPositifMobilityJuni2020.xlsx")
mobilityjakartaplot(mobilityjakarta$Tanggal,mobilityjakarta$POSITIF)
Data Deskriptif
library(jmv)## Warning: package 'jmv' was built under R version 4.1.3# Use the descritptives function to get the descritptive data
descriptives(mobilityjakarta, vars = vars(POSITIF, retail_and_recreation_percent_change_from_baseline), freq = TRUE)##
## DESCRIPTIVES
##
## Descriptives
## ----------------------------------------------------------------------------------------
## POSITIF retail_and_recreation_percent_change_from_baseline
## ----------------------------------------------------------------------------------------
## N 30 30
## Missing 0 0
## Mean 9393.100 -38.63333
## Median 9176.500 -38.00000
## Standard deviation 1537.402 7.141348
## Minimum 7383.000 -52.00000
## Maximum 13040.00 -27.00000
## ----------------------------------------------------------------------------------------hist(mobilityjakarta$POSITIF, main="Distribusi Positif di RS untuk Covid 19", xlab="Positif di RS Untuk Covid 19", xlim=c(3000,6000),ylim=c(0,15)) # histogram variabel wt
mean=mean(mobilityjakarta$POSITIF) # Menghitung nilai mean variabel wt
abline(v=mean, col="green",lwd=2) # Menambahkan garis vertical pada plot
# 2. Density plot
d<-density(mobilityjakarta$POSITIF) # menghitung density variabel wt
plot(d, xlab="Positif di RS untuk Covid 19", main="Distribusi Positif di RS untuk Covid 19")
Perbandingan Pendekatan
library(npreg)## Warning: package 'npreg' was built under R version 4.1.3# fit using ss
mod.ss <- ss( mobilityjakarta$Tanggal,mobilityjakarta$POSITIF, nknots = 10)
mod.ss##
## Call:
## ss(x = mobilityjakarta$Tanggal, y = mobilityjakarta$POSITIF,
## nknots = 10)
##
## Smoothing Parameter spar = 0.1905552 lambda = 1.419018e-06
## Equivalent Degrees of Freedom (Df) 9.180708
## Penalized Criterion (RSS) 11100214
## Generalized Cross-Validation (GCV) 768282.1# fit using smooth.spline
mod.smsp <- smooth.spline( mobilityjakarta$Tanggal,mobilityjakarta$POSITIF, nknots = 10)
mod.smsp## Call:
## smooth.spline(x = mobilityjakarta$Tanggal, y = mobilityjakarta$POSITIF,
## nknots = 10)
##
## Smoothing Parameter spar= 0.1819329 lambda= 5.045269e-05 (12 iterations)
## Equivalent Degrees of Freedom (Df): 9.134215
## Penalized Criterion (RSS): 11298277
## GCV: 778509.7# rmse between solutions
sqrt(mean(( mod.ss$y - mod.smsp$y )^2))## [1] 24.69569# rmse between solutions and f(x)
sqrt(mean(( mobilityjakarta$POSITIF - mod.ss$y )^2))## [1] 608.2821sqrt(mean(( mobilityjakarta$POSITIF - mod.smsp$y )^2))## [1] 613.685mobilityjakarta$prediksi_ss <- mod.ss$y
mobilityjakarta$prediksi_smsp <- mod.smsp$y
mobilityjakarta# plot results
plot(mobilityjakarta$Tanggal, mobilityjakarta$POSITIF)
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$POSITIF, lty = 10, col = 'black', 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
#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")
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 = 1.5 lambda = 4095.998
## Equivalent Degrees of Freedom (Df) 2.000001
## Penalized Criterion (RSS) 356.2382
## Generalized Cross-Validation (GCV) 13.63156# 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= 1.479131 lambda= 118781.3 (25 iterations)
## Equivalent Degrees of Freedom (Df): 2
## Penalized Criterion (RSS): 356.2382
## GCV: 13.63156mobilityjakarta$prediksi_ss1 <- mod.ss1$y
mobilityjakarta$prediksi_smsp1 <- mod.smsp1$y
mobilityjakarta# plot results
plot(mobilityjakarta$Tanggal, mobilityjakarta$retail_and_recreation_percent_change_from_baseline, lty = 10, col = 'black', lwd =5)
# add lm fit
abline(coef(lm( mobilityjakarta$retail_and_recreation_percent_change_from_baseline ~ mobilityjakarta$Tanggal , data = mobilityjakarta)), lty = 10, col = 'blue', lwd =5)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_ss1 , lty = 2, col = 'black', lwd = 4)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_smsp1, lty = 4, col = 'red', lwd = 2)
Korelasi Pearson
Untuk melakukan analisis korelasi Pearson cukup sederhana. Semisal kita ingin mengetahui hubungan antara mobilityjakartaPOSITIFdanmobilityjakartaretail_and_recreation_percent_change_from_baseline,
cor.test(mobilityjakarta$POSITIF,mobilityjakarta$retail_and_recreation_percent_change_from_baseline)##
## Pearson's product-moment correlation
##
## data: mobilityjakarta$POSITIF and mobilityjakarta$retail_and_recreation_percent_change_from_baseline
## t = 6.212, df = 28, p-value = 1.037e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5525143 0.8801389
## sample estimates:
## cor
## 0.7612531Dari output tersebut dapat kita simpulkan bahwa tidak ada hubungan yang signifikan antara mobilityjakartaPOSITIF dan mobilityjakartaretail_and_recreation_percent_change_from_baseline (p-value<0,01). Nilai p-value dalam output R dituliskan 1.037 Nilai koefisien korelasi r adalah sebesar 0.7612531 yang menunjukkan hubungan yang lemah dan terbalik.
Hasil Regresi dengan Pendekatan Smoothing Spline
model <- lm(mobilityjakarta$POSITIF ~ mobilityjakarta$retail_and_recreation_percent_change_from_baseline, data = mobilityjakarta)
summary(model)##
## Call:
## lm(formula = mobilityjakarta$POSITIF ~ mobilityjakarta$retail_and_recreation_percent_change_from_baseline,
## data = mobilityjakarta)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1348.3 -661.8 -182.7 186.4 3051.4
##
## Coefficients:
## Estimate
## (Intercept) 15724.48
## mobilityjakarta$retail_and_recreation_percent_change_from_baseline 163.88
## Std. Error
## (Intercept) 1035.92
## mobilityjakarta$retail_and_recreation_percent_change_from_baseline 26.38
## t value
## (Intercept) 15.179
## mobilityjakarta$retail_and_recreation_percent_change_from_baseline 6.212
## Pr(>|t|)
## (Intercept) 4.85e-15 ***
## mobilityjakarta$retail_and_recreation_percent_change_from_baseline 1.04e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1015 on 28 degrees of freedom
## Multiple R-squared: 0.5795, Adjusted R-squared: 0.5645
## F-statistic: 38.59 on 1 and 28 DF, p-value: 1.037e-06Rumus linear regresi menjadi y = 15724.48 + (163.88) mobilityjakarta$retail_and_recreation_percent_change_from_baseline
model$coefficients## (Intercept)
## 15724.483
## mobilityjakarta$retail_and_recreation_percent_change_from_baseline
## 163.884model$fitted.values## 1 2 3 4 5 6 7 8
## 7366.402 8185.822 8021.938 8185.822 8513.589 7858.054 7202.518 9005.241
## 9 10 11 12 13 14 15 16
## 8513.589 8841.357 9005.241 9496.893 8513.589 8021.938 10316.313 9988.545
## 17 18 19 20 21 22 23 24
## 9988.545 10152.429 10316.313 9660.777 9169.125 10644.081 10316.313 10316.313
## 25 26 27 28 29 30
## 10644.081 11135.733 10480.197 9496.893 11135.733 11299.617mobilityjakarta$prediksi_linear <- model$fitted.values
mobilityjakarta# 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$POSITIF~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$POSITIF, lty = 10, col = 'black', 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
#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")