Analisa Data dengan Regresi Nonparametrik Menggunakan Pendekatan Smoothing Spline Pada Data Google Mobility Index dan Covid-19 di Jakarta Agustus 2020
Adila Qurrota A’yun
4/19/2022
Dosen Pengampu : Prof. Dr. Suhartono, M.Kom
Mata Kuliah : Linear Algebra
Prodi : Teknik Informatika
Lembaga : Universitas Islam Negeri Maulana Malik Ibrahim Malang
1 Pengertian Smoothing Spline
Smoothing spline merupakan pendekatan yang ampuh untuk memperkirakan hubungan fungsional antara prediktor X dan respons Y. Smoothing spline dapat ditampung baik menggunakan fungsi smooth.spline (dalam paket stats) atau fungsi ss (dalam paket npreg). Dokumen ini memberikan latar belakang teoritis mengenai smoothing spline serta contoh-contoh yang menggambarkan bagaimana menggunakan fungsi smooth.spline dan ss.
Kedua fungsi tersebut 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 serta memiliki lebih banyak jenis spline (linier, kubik, quintic) opsi untuk batasan periodisitas opsi dan nilai simpul yang ditentukan pengguna ringkasan yang mana sesuai dengan metode plot. Berikut analisa data dengan regresi nonparametrik menggunakan pendekatan smoothing spline pada data Google Mobility Index dan Covid-19 di Jakarta Agustus 2020.
2 Data Google Mobility Index dan Covid-19 di Jakarta Agustus 2020
## Import the data
library (readxl)mobilityJakarta <- read_excel(path = "Data Google Mobility Index dan Covid-19 DKI Agustus 2020.xlsx")
mobilityJakartaplot(mobilityJakarta$Tanggal,mobilityJakarta$`Self Isolation`)
3 Data Deskriptif
library (jmv)## Warning: package 'jmv' was built under R version 4.1.3# 1. Use the descritptives function to get the descritptive data
descriptives(mobilityJakarta, vars = vars('Self Isolation', retail_and_recreation_percent_change_from_baseline), freq = TRUE)##
## DESCRIPTIVES
##
## Descriptives
## ----------------------------------------------------------------------------------------------
## Self Isolation retail_and_recreation_percent_change_from_baseline
## ----------------------------------------------------------------------------------------------
## N 31 31
## Missing 0 0
## Mean 5670.290 -28.32258
## Median 6040.000 -27.00000
## Standard deviation 746.4055 4.061093
## Minimum 4128.000 -37.00000
## Maximum 6602.000 -21.00000
## ----------------------------------------------------------------------------------------------hist(mobilityJakarta$`Self Isolation`, main="Distribusi Self Isolation di RS untuk Covid-19", xlab="Self Isolation di RS Untuk Covid-19", xlim=c(3000,6000),ylim=c(0,15)) # histogram variabel wt
mean=mean(mobilityJakarta$`Self Isolation`) # Menghitung nilai mean variabel wt
abline(v=mean, col="blue",lwd=2) # Menambahkan garis vertical pada plot
# 2. Density plot
d<-density(mobilityJakarta$`Self Isolation`) # menghitung density variabel wt
plot(d, xlab="Self Isolation di RS untuk Covid-19", main="Distribusi Self Isolation di RS untuk 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$`Self Isolation`, nknots = 10)
mod.ss##
## Call:
## ss(x = mobilityJakarta$Tanggal, y = mobilityJakarta$`Self Isolation`,
## nknots = 10)
##
## Smoothing Parameter spar = 0.2052099 lambda = 1.810773e-06
## Equivalent Degrees of Freedom (Df) 8.879359
## Penalized Criterion (RSS) 640096.1
## Generalized Cross-Validation (GCV) 40551.92# fit using smooth.spline
mod.smsp <- smooth.spline(mobilityJakarta$Tanggal,mobilityJakarta$`Self Isolation`, nknots = 10)
mod.smsp## Call:
## smooth.spline(x = mobilityJakarta$Tanggal, y = mobilityJakarta$`Self Isolation`,
## nknots = 10)
##
## Smoothing Parameter spar= 0.1640644 lambda= 3.999867e-05 (15 iterations)
## Equivalent Degrees of Freedom (Df): 9.407103
## Penalized Criterion (RSS): 586573.5
## GCV: 38999.79# rmse between solutions
sqrt(mean(( mod.ss$y - mod.smsp$y )^2))## [1] 15.91299# rmse between solutions and f(x)
sqrt(mean((mobilityJakarta$`Self Isolation` - mod.ss$y )^2))## [1] 143.695sqrt(mean((mobilityJakarta$`Self Isolation` - mod.smsp$y )^2))## [1] 137.5563mobilityJakarta$prediksi_ss <- mod.ss$y
mobilityJakarta$prediksi_smsp <- mod.smsp$y
mobilityJakarta# plot results
plot(mobilityJakarta$Tanggal, mobilityJakarta$`Self Isolation`)
lines(mobilityJakarta$Tanggal, mobilityJakarta$prediksi_ss , lty = 2, col = 'red', lwd = 4)
lines(mobilityJakarta$Tanggal, mobilityJakarta$prediksi_smsp, lty = 4, col = 'blue', lwd = 2)
# plot method
plot(mobilityJakarta$Tanggal, mobilityJakarta$`Self Isolation`, lty = 10, col = 'black', lwd =5)
# plot(mod.ss)
# add lm fit
abline(coef(lm(mobilityJakarta$`Self Isolation` ~ mobilityJakarta$Tanggal , data = mobilityJakarta)), lty = 10, col = 'blue', lwd =5)
rug(mobilityJakarta$Tanggal) # add rug to plot
#points(mobilityJakarta$Tanggal, mobilityJakarta$'Self Isolation', 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 = 0.3686358 lambda = 2.74518e-05
## Equivalent Degrees of Freedom (Df) 5.788914
## Penalized Criterion (RSS) 249.0038
## Generalized Cross-Validation (GCV) 12.14464# 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.3479396 lambda= 0.0008521042 (10 iterations)
## Equivalent Degrees of Freedom (Df): 5.797622
## Penalized Criterion (RSS): 248.853
## GCV: 12.14567mobilityJakarta$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)
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.
Untuk melakukan analisis korelasi pearson cukup sederhana. Semisal kita ingin mengetahui hubungan antara mobilityJakarta Self Isolation dan mobilityJakarta retail_and_recreation_percent_change_from_baseline yakni sebagai berikut :
cor.test(mobilityJakarta$`Self Isolation`,mobilityJakarta$retail_and_recreation_percent_change_from_baseline)##
## Pearson's product-moment correlation
##
## data: mobilityJakarta$`Self Isolation` and mobilityJakarta$retail_and_recreation_percent_change_from_baseline
## t = -1.3474, df = 29, p-value = 0.1883
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5497841 0.1221125
## sample estimates:
## cor
## -0.2427305Dari output tersebut dapat disimpulkan bahwa tidak ada hubungan yang signifikan antara mobilityJakarta SelfIsolation dan mobilityJakarta retail_and_recreation (p-value<0,01). Nilai p-value dalam output R dituliskan 0.1883 dan nilai koefisien korelasi r adalah sebesar -0.2427305 yang menunjukkan hubungan sedang dan terbalik.
6 Hasil Regresi dengan Pendekatan Smoothing Spline
model <- lm(mobilityJakarta$`Self Isolation` ~ mobilityJakarta$retail_and_recreation_percent_change_from_baseline, data = mobilityJakarta)
summary(model)##
## Call:
## lm(formula = mobilityJakarta$`Self Isolation` ~ mobilityJakarta$retail_and_recreation_percent_change_from_baseline,
## data = mobilityJakarta)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1374.4 -524.3 256.6 610.7 921.3
##
## Coefficients:
## Estimate
## (Intercept) 4406.75
## mobilityJakarta$retail_and_recreation_percent_change_from_baseline -44.61
## Std. Error
## (Intercept) 947.02
## mobilityJakarta$retail_and_recreation_percent_change_from_baseline 33.11
## t value
## (Intercept) 4.653
## mobilityJakarta$retail_and_recreation_percent_change_from_baseline -1.347
## Pr(>|t|)
## (Intercept) 6.64e-05 ***
## mobilityJakarta$retail_and_recreation_percent_change_from_baseline 0.188
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 736.5 on 29 degrees of freedom
## Multiple R-squared: 0.05892, Adjusted R-squared: 0.02647
## F-statistic: 1.816 on 1 and 29 DF, p-value: 0.1883rumus linear regresi menjadi y = 4406.75 + (-44.61) mobilityJakarta$retail_and_recreation_percent_change_from_baseline
model$coefficients## (Intercept)
## 4406.75016
## mobilityJakarta$retail_and_recreation_percent_change_from_baseline
## -44.61247model$fitted.values## 1 2 3 4 5 6 7 8
## 6057.411 6057.411 5611.287 5655.899 5611.287 5611.287 5432.837 5611.287
## 9 10 11 12 13 14 15 16
## 5834.349 5522.062 5655.899 5566.674 5834.349 5432.837 5655.899 5878.962
## 17 18 19 20 21 22 23 24
## 5834.349 5789.737 5566.674 5834.349 5700.512 5834.349 5968.186 5611.287
## 25 26 27 28 29 30 31
## 5566.674 5566.674 5522.062 5343.612 5522.062 5655.899 5432.837mobilityJakarta$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$`Self Isolation`~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$`Self Isolation`, lty = 10, col = 'black', lwd =5)
# plot(mod.ss)
# add lm fit
abline(coef(lm(mobilityJakarta$`Self Isolation` ~ mobilityJakarta$Tanggal , data = mobilityJakarta)), lty = 10, col = 'blue', lwd =5)
rug(mobilityJakarta$Tanggal) # add rug to plot
#points(mobilityJakarta$Tanggal, mobilityJakarta$'Self Isolation', 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","Purple"), lwd = 3, bty = "p")