Regresi Nonparametrik dengan Pendekatan Smoothing Spline Pada Data Google Mobility Index dan Covid-19 di Jakarta Februari 2021
Shafira Halmahera
21 April 2022
Dosen Pengampu : Prof. Dr. Suhartono, M.Kom
Lembaga : Universitas Islam Negeri Maulana Malik Ibrahim Malang
Mata kuliah : Linier Algebra / C
Jurusan : Teknik Informatika
1. Pengertian Smoothing Spline
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
2. Data Google Mobility Index dan Covid-19 di Jakarta Februari 2021
## Import the data
library(readxl)
mobilityjakarta <- read_excel(path = "C:/Users/shafira halmahera/Documents/LINEAR ALGEBRA/Data_Exel/Data Google Mobility Index dan Covid-19 Februari 2021.xlsx")
mobilityjakartaplot(mobilityjakarta$Tanggal,mobilityjakarta$Dirawat)
3. Data Deskriptif
library(jmv)## Warning: package 'jmv' was built under R version 4.1.3## Warning: package 'jmv' was built under R version 4.1.3# Mendapatkan data descriptive menggunakan fungsi descritptive
descriptives(mobilityjakarta, vars = vars(Dirawat, retail_and_recreation_percent_change_from_baseline), freq = TRUE)##
## DESCRIPTIVES
##
## Descriptives
## ----------------------------------------------------------------------------------------
## Dirawat retail_and_recreation_percent_change_from_baseline
## ----------------------------------------------------------------------------------------
## N 28 28
## Missing 0 0
## Mean 5992.429 -32.39286
## Median 5709.500 -32.50000
## Standard deviation 1689.662 4.121662
## Minimum 3917.000 -39.00000
## Maximum 9888.000 -24.00000
## ----------------------------------------------------------------------------------------hist(mobilityjakarta$Dirawat, main="Distribusi Dirawat di RS untuk Covid 19", xlab="Dirawat di RS Untuk Covid 19 Bulan Februari 2021", xlim=c(3000,6000),ylim=c(0,15)) # histogram variabel wt
mean=mean(mobilityjakarta$Dirawat) # Menghitung nilai mean variabel wt
abline(v=mean, col="blue",lwd=2) # Menambahkan garis vertical pada plot
#Density plot
d<-density(mobilityjakarta$Dirawat) # menghitung density variabel wt
plot(d, xlab="Dirawat di RS untuk Covid 19", main="Distribusi Dirawat di RS untuk Covid 19")
3. Data Deskriptif
library(npreg)## Warning: package 'npreg' was built under R version 4.1.3## Warning: package 'npreg' was built under R version 4.1.3# fit using ss
mod.ss <- ss( mobilityjakarta$Tanggal,mobilityjakarta$Dirawat, nknots = 10)
mod.ss##
## Call:
## ss(x = mobilityjakarta$Tanggal, y = mobilityjakarta$Dirawat,
## nknots = 10)
##
## Smoothing Parameter spar = 0.2230408 lambda = 2.436055e-06
## Equivalent Degrees of Freedom (Df) 8.506742
## Penalized Criterion (RSS) 3682341
## Generalized Cross-Validation (GCV) 271339.6# fit using smooth.spline
mod.smsp <- smooth.spline( mobilityjakarta$Tanggal,mobilityjakarta$Dirawat, nknots = 10)
mod.smsp## Call:
## smooth.spline(x = mobilityjakarta$Tanggal, y = mobilityjakarta$Dirawat,
## nknots = 10)
##
## Smoothing Parameter spar= 0.2077412 lambda= 7.98616e-05 (14 iterations)
## Equivalent Degrees of Freedom (Df): 8.503743
## Penalized Criterion (RSS): 3748757
## GCV: 276148.6# rmse between solutions
sqrt(mean(( mod.ss$y - mod.smsp$y )^2))## [1] 11.89723# rmse between solutions and f(x)
sqrt(mean(( mobilityjakarta$DIRAWAT - mod.ss$y )^2))## Warning: Unknown or uninitialised column: `DIRAWAT`.## [1] NaNsqrt(mean(( mobilityjakarta$DIRAWAT - mod.smsp$y )^2))## Warning: Unknown or uninitialised column: `DIRAWAT`.## [1] NaNmobilityjakarta$prediksi_ss <- mod.ss$y
mobilityjakarta$prediksi_smsp <- mod.smsp$y
mobilityjakarta# plot results
plot(mobilityjakarta$Tanggal, mobilityjakarta$Dirawat)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_ss , lty = 2, col = 'black', lwd = 4)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_smsp, lty = 4, col = 'purple', lwd = 4)
# plot method
plot(mobilityjakarta$Tanggal, mobilityjakarta$Dirawat, lty = 10, col = 'black', lwd =5)
#plot(mod.ss)
# add lm fit
abline(coef(lm( mobilityjakarta$Dirawat ~ mobilityjakarta$Tanggal , data = mobilityjakarta)), lty = 10, col = 'red', 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 = 'orange', lwd = 4)
#points(mobilityjakarta$Tanggal,mod.smsp$y , lty = 2, col = 4, lwd = 2)
legend("topleft",
legend = c("Real", "Model SS", "Trends"),
lty = 1:3, col = c("Black","orange","red"), 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.4772456 lambda = 0.0001672036
## Equivalent Degrees of Freedom (Df) 4.148249
## Penalized Criterion (RSS) 276.639
## Generalized Cross-Validation (GCV) 13.61542# 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.4522355 lambda= 0.004663811 (12 iterations)
## Equivalent Degrees of Freedom (Df): 4.155284
## Penalized Criterion (RSS): 276.3759
## GCV: 13.6105mobilityjakarta$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 =4)
# add lm fit
abline(coef(lm( mobilityjakarta$retail_and_recreation_percent_change_from_baseline ~ mobilityjakarta$Tanggal , data = mobilityjakarta)), lty = 10, col = 'Brown', lwd =5)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_ss1 , lty = 2, col = 'salmon', lwd = 6)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_smsp1, lty = 4, col = 'Orchid', lwd = 4)
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$Dirawat,mobilityjakarta$retail_and_recreation_percent_change_from_baseline)##
## Pearson's product-moment correlation
##
## data: mobilityjakarta$Dirawat and mobilityjakarta$retail_and_recreation_percent_change_from_baseline
## t = -1.7876, df = 26, p-value = 0.0855
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.62657966 0.04818325
## sample estimates:
## cor
## -0.3308409Dari output tersebut dapat kita simpulkan bahwa tidak ada hubungan yang signifikan antara mobilityjakarta Dirawat danmobilityjakarta retail_and_recreation_percent_change_from_baseline (p-value<0,01). Nilai p-value dalam output R dituliskan 0.0855 Nilai koefisien korelasi r adalah sebesar -0.3308409 yang menunjukkan hubungan yang sedang dan terbalik.
6. Hasil Regresi dengan Pendekatan Smoothing Spline
model <- lm(mobilityjakarta$Dirawat ~ mobilityjakarta$retail_and_recreation_percent_change_from_baseline, data = mobilityjakarta)
summary(model)##
## Call:
## lm(formula = mobilityjakarta$Dirawat ~ mobilityjakarta$retail_and_recreation_percent_change_from_baseline,
## data = mobilityjakarta)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2042.4 -1069.2 -501.3 881.8 3677.6
##
## Coefficients:
## Estimate
## (Intercept) 1599.08
## mobilityjakarta$retail_and_recreation_percent_change_from_baseline -135.63
## Std. Error
## (Intercept) 2476.75
## mobilityjakarta$retail_and_recreation_percent_change_from_baseline 75.87
## t value
## (Intercept) 0.646
## mobilityjakarta$retail_and_recreation_percent_change_from_baseline -1.788
## Pr(>|t|)
## (Intercept) 0.5242
## mobilityjakarta$retail_and_recreation_percent_change_from_baseline 0.0855 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1625 on 26 degrees of freedom
## Multiple R-squared: 0.1095, Adjusted R-squared: 0.0752
## F-statistic: 3.196 on 1 and 26 DF, p-value: 0.0855Rumus linear regresi menjadi y = 1599.08 + (-135.63) mobilityjakarta$retail_and_recreation
mobilityjakarta$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(1, 1, 0, alpha = 0.5), ...)
# lines(density(x), col = 2, lwd = 2) # Uncomment to add density lines
}
# Equivalent with a formula
pairs(mobilityjakarta$Dirawat~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$Dirawat, lty = 10, col = 'black', lwd =4)
#plot(mod.ss)
# add lm fit
abline(coef(lm( mobilityjakarta$Dirawat ~ mobilityjakarta$Tanggal , data = mobilityjakarta)), lty = 10, col = 'brown', 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 = 'turquoise', 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 = 3)
legend("topright",
legend = c("Real", "Model SS", "Trends", "Linear Regresi"),
lty = 1:4, col = c("Black","turquoise","brown","purple"), lwd = 3, bty = "p")