Analisa Data dengan Regresi Linear Tiga Model Menggunakan Pendekatan Smooting Spline Pada Data Covid-19 dan Google Mobility Index di Jakarta Januari 2022
Enggarani Wahyu Ekaputri, Prof. Dr. Suhartono, M.Kom
6 Juni 2022
Lembaga : Universitas Islam Negeri Maulana Malik Ibrahim Malang
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
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.
2 Data Covid-19 dan Google Mobility Index di Jakarta Januari 2022
library(readxl)covidmobilityjakarta <- read_excel(path = "DataCovid19 DKI Jakarta.xlsx")
covidmobilityjakartaplot(covidmobilityjakarta$Tanggal,covidmobilityjakarta$Dirawat)
3 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(covidmobilityjakarta, vars = vars(Dirawat, retail_and_recreation_percent_change_from_baseline), freq = TRUE)##
## DESCRIPTIVES
##
## Descriptives
## ----------------------------------------------------------------------------------------
## Dirawat retail_and_recreation_percent_change_from_baseline
## ----------------------------------------------------------------------------------------
## N 31 31
## Missing 0 0
## Mean 1634.839 7.451613
## Median 781.0000 7.000000
## Standard deviation 1930.584 3.085868
## Minimum 116.0000 1.000000
## Maximum 6809.000 14.00000
## ----------------------------------------------------------------------------------------hist(covidmobilityjakarta$Dirawat, main="Distribusi Dirawat di RS untuk Covid 19", xlab="Dirawat di RS Untuk Covid 19", xlim=c(3000,6000),ylim=c(0,15)) # histogram variabel wt
mean=mean(covidmobilityjakarta$Dirawat) # Menghitung nilai mean variabel wt
abline(v=mean, col="red",lwd=2) # Menambahkan garis vertical pada plot
# 2. Density plot
d<-density(covidmobilityjakarta$Dirawat) # menghitung density variabel wt
plot(d, xlab="Dirawat di RS untuk Covid 19", main="Distribusi Dirawat 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( covidmobilityjakarta$Tanggal,covidmobilityjakarta$Dirawat, nknots = 10)
mod.ss##
## Call:
## ss(x = covidmobilityjakarta$Tanggal, y = covidmobilityjakarta$Dirawat,
## nknots = 10)
##
## Smoothing Parameter spar = 0.01776397 lambda = 8.009766e-08
## Equivalent Degrees of Freedom (Df) 10.76019
## Penalized Criterion (RSS) 55189.81
## Generalized Cross-Validation (GCV) 4176.454# fit using smooth.spline
mod.smsp <- smooth.spline( covidmobilityjakarta$Tanggal,covidmobilityjakarta$Dirawat, nknots = 10)
mod.smsp## Call:
## smooth.spline(x = covidmobilityjakarta$Tanggal, y = covidmobilityjakarta$Dirawat,
## nknots = 10)
##
## Smoothing Parameter spar= -0.05618678 lambda= 1.0263e-06 (16 iterations)
## Equivalent Degrees of Freedom (Df): 11.7885
## Penalized Criterion (RSS): 20919.79
## GCV: 1757.101# rmse between solutions
sqrt(mean(( mod.ss$y - mod.smsp$y )^2))## [1] 32.00289# rmse between solutions and f(x)
sqrt(mean(( covidmobilityjakarta$Dirawat - mod.ss$y )^2))## [1] 42.1938sqrt(mean(( covidmobilityjakarta$Dirawat - mod.smsp$y )^2))## [1] 25.97753covidmobilityjakarta$prediksi_ss <- mod.ss$y
covidmobilityjakarta$prediksi_smsp <- mod.smsp$y
covidmobilityjakarta# plot results
plot(covidmobilityjakarta$Tanggal, covidmobilityjakarta$Dirawat)
lines(covidmobilityjakarta$Tanggal, covidmobilityjakarta$prediksi_ss , lty = 2, col = 'black', lwd = 4)
lines(covidmobilityjakarta$Tanggal, covidmobilityjakarta$prediksi_smsp, lty = 4, col = 'red', lwd = 2)
# plot method
plot(covidmobilityjakarta$Tanggal, covidmobilityjakarta$Dirawat, lty = 10, col = 'black', lwd =5)
#plot(mod.ss)
# add lm fit
abline(coef(lm( covidmobilityjakarta$Dirawat ~ covidmobilityjakarta$Tanggal , data = covidmobilityjakarta)), lty = 10, col = 'blue', lwd =5)
rug(covidmobilityjakarta$Tanggal) # add rug to plot
#points(covidmobilityjakarta$Tanggal, covidmobilityjakarta$Dirawat, lty = 2, col = 2, lwd = 2)
points(covidmobilityjakarta$Tanggal,mod.ss$y , lty = 2, col = 'red', lwd = 4)
#points(covidmobilityjakarta$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( covidmobilityjakarta$Tanggal,covidmobilityjakarta$retail_and_recreation_percent_change_from_baseline, nknots = 10)
mod.ss1##
## Call:
## ss(x = covidmobilityjakarta$Tanggal, y = covidmobilityjakarta$retail_and_recreation_percent_change_from_baseline,
## nknots = 10)
##
## Smoothing Parameter spar = 0.2344568 lambda = 2.945537e-06
## Equivalent Degrees of Freedom (Df) 8.355361
## Penalized Criterion (RSS) 96.48915
## Generalized Cross-Validation (GCV) 5.833233# fit using smooth.spline
mod.smsp1 <- smooth.spline( covidmobilityjakarta$Tanggal,covidmobilityjakarta$retail_and_recreation_percent_change_from_baseline, nknots = 10)
mod.smsp1## Call:
## smooth.spline(x = covidmobilityjakarta$Tanggal, y = covidmobilityjakarta$retail_and_recreation_percent_change_from_baseline,
## nknots = 10)
##
## Smoothing Parameter spar= -0.004828107 lambda= 2.409043e-06 (15 iterations)
## Equivalent Degrees of Freedom (Df): 11.55074
## Penalized Criterion (RSS): 55.17046
## GCV: 4.52129covidmobilityjakarta$prediksi_ss1 <- mod.ss1$y
covidmobilityjakarta$prediksi_smsp1 <- mod.smsp1$y
covidmobilityjakarta# plot results
plot(covidmobilityjakarta$Tanggal, covidmobilityjakarta$retail_and_recreation_percent_change_from_baseline, lty = 10, col = 'black', lwd =5)
# add lm fit
abline(coef(lm( covidmobilityjakarta$retail_and_recreation_percent_change_from_baseline ~ covidmobilityjakarta$Tanggal , data = covidmobilityjakarta)), lty = 10, col = 'blue', lwd =5)
lines(covidmobilityjakarta$Tanggal, covidmobilityjakarta$prediksi_ss1 , lty = 2, col = 'black', lwd = 4)
lines(covidmobilityjakarta$Tanggal, covidmobilityjakarta$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. Korelasi Pearson Untuk melakukan analisis korelasi Pearson cukup sederhana. Semisal kita ingin mengetahui hubungan antara mobilityjakartaMeninggaldanmobilityjakartaretail_and_recreation_percent_change_from_baseline,
cor.test(covidmobilityjakarta$Dirawat,covidmobilityjakarta$retail_and_recreation_percent_change_from_baseline)##
## Pearson's product-moment correlation
##
## data: covidmobilityjakarta$Dirawat and covidmobilityjakarta$retail_and_recreation_percent_change_from_baseline
## t = 1.6015, df = 29, p-value = 0.1201
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.07707823 0.58073124
## sample estimates:
## cor
## 0.2850471Dari 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 0.2033 Nilai koefisien korelasi r adalah sebesar 0.2479615 yang menunjukkan hubungan yang lemah dan terbalik.
6 Hasil Regresi dengan Pendekatan Smoothing Spline
model <- lm(covidmobilityjakarta$Dirawat ~ covidmobilityjakarta$retail_and_recreation_percent_change_from_baseline, data = covidmobilityjakarta)
summary(model)##
## Call:
## lm(formula = covidmobilityjakarta$Dirawat ~ covidmobilityjakarta$retail_and_recreation_percent_change_from_baseline,
## data = covidmobilityjakarta)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2662.6 -1193.6 -388.6 503.0 4615.4
##
## Coefficients:
## Estimate
## (Intercept) 306.0
## covidmobilityjakarta$retail_and_recreation_percent_change_from_baseline 178.3
## Std. Error
## (Intercept) 896.0
## covidmobilityjakarta$retail_and_recreation_percent_change_from_baseline 111.4
## t value
## (Intercept) 0.342
## covidmobilityjakarta$retail_and_recreation_percent_change_from_baseline 1.601
## Pr(>|t|)
## (Intercept) 0.735
## covidmobilityjakarta$retail_and_recreation_percent_change_from_baseline 0.120
##
## Residual standard error: 1882 on 29 degrees of freedom
## Multiple R-squared: 0.08125, Adjusted R-squared: 0.04957
## F-statistic: 2.565 on 1 and 29 DF, p-value: 0.1201Rumus linear regresi menjadi y = 5577.47 + 20.77 mobilityjakarta$retail_and_recreation
model$coefficients## (Intercept)
## 305.9820
## covidmobilityjakarta$retail_and_recreation_percent_change_from_baseline
## 178.3314model$fitted.values## 1 2 3 4 5 6 7 8
## 840.9763 2267.6276 2802.6218 2267.6276 2445.9590 2089.2962 1375.9705 1732.6334
## 9 10 11 12 13 14 15 16
## 1554.3019 1554.3019 484.3135 1197.6391 1732.6334 1375.9705 1197.6391 1375.9705
## 17 18 19 20 21 22 23 24
## 1197.6391 484.3135 1375.9705 1375.9705 1554.3019 1375.9705 1554.3019 1732.6334
## 25 26 27 28 29 30 31
## 1732.6334 1910.9648 2089.2962 1554.3019 1910.9648 1732.6334 2802.6218covidmobilityjakarta$prediksi_linear <- model$fitted.values
covidmobilityjakarta# 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(covidmobilityjakarta$Dirawat~covidmobilityjakarta$retail_and_recreation_percent_change_from_baseline, data = covidmobilityjakarta,
upper.panel = NULL, # Disabling the upper panel
diag.panel = panel.hist) # Adding the histograms
# plot method
plot(covidmobilityjakarta$Tanggal, covidmobilityjakarta$Dirawat, lty = 10, col = 'black', lwd =5)
#plot(mod.ss)
# add lm fit
abline(coef(lm( covidmobilityjakarta$Dirawat ~ covidmobilityjakarta$Tanggal , data = covidmobilityjakarta)), lty = 10, col = 'blue', lwd =5)
rug(covidmobilityjakarta$Tanggal) # add rug to plot
#points(covidmobilityjakarta$Tanggal, covidmobilityjakarta$Dirawat, lty = 2, col = 2, lwd = 2)
lines(covidmobilityjakarta$Tanggal,mod.ss$y , lty = 2, col = 'red', lwd = 4)
#points(covidmobilityjakarta$Tanggal,mod.smsp$y , lty = 2, col = 4, lwd = 2)
lines(covidmobilityjakarta$Tanggal, covidmobilityjakarta$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")