Analisis Data dengan Relasi Tiga Model Smoothing Spline Model Kasus COVID-19 Jakarta Bulan November 2021
Date: 2022-04-25 Author : Setiyaris
Mata Kuliah : Pengantar Konsep Pemrograman Algoritma dan Struktur Data Magister Teknik Informatika
Dosen Pengampu : Prof Dr Suhartono M.Kom UIN Maulana Malik Ibrahim Malang 1. Penjelasan Smoothing Spline Model Smoothing merupakan salah satu metode yang digunakan dalam analisis data non parametrik. Tujuan dari smoothing adalah untuk meminimalkan keragaman karakteristik data dari data yang tidak memiliki pengaruh sehingga ciri-ciri dari data akan tampak lebih jelas.
Smooth Spline Model adalah Salah satu model regresi dengan pendekatan non parametrik yang dapat digunakan untuk menduga kurva regresi adalah regresi spline.
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.
Tujuan dari Analisa Data saat ini adalah membandingkan Kasus Positif Covid-19 dengan pergerakan orang di tempat Bekerja. Masing-masing unsur menggunakan Smooth Spline sehingga akan didapat Regresi Linear diantaranya dan memprediksi adakah hubungan keterikatan antara Kasus Positif Covid-19 dan Lokasi Tempat Bekerja.
Sumber Data diambil dari Data Covid-19 DKI Jakarta dari website Pemprov DKI Jakarta dan Google mobility index pada bulan November 2021.
library(readxl)
mobilityjakarta <- read_excel(path = "covidjakartanov.xlsx")## Warning in strptime(x, format, tz = tz): unable to identify current timezone 'C':
## please set environment variable 'TZ'mobilityjakarta## # A tibble: 30 x 13
## Tanggal Nama_Kota POSITIF Dirawat Sembuh Meninggal
## <dttm> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 2021-11-01 00:00:00 DKI Jakarta 861623 273 847066 13562
## 2 2021-11-02 00:00:00 DKI Jakarta 861700 266 847201 13562
## 3 2021-11-03 00:00:00 DKI Jakarta 861832 258 847353 13563
## 4 2021-11-04 00:00:00 DKI Jakarta 861943 258 847403 13564
## 5 2021-11-05 00:00:00 DKI Jakarta 862061 256 847438 13565
## 6 2021-11-06 00:00:00 DKI Jakarta 862130 231 847685 13566
## 7 2021-11-07 00:00:00 DKI Jakarta 862247 220 847808 13566
## 8 2021-11-08 00:00:00 DKI Jakarta 862276 198 847885 13566
## 9 2021-11-09 00:00:00 DKI Jakarta 862370 193 848034 13566
## 10 2021-11-10 00:00:00 DKI Jakarta 862465 198 848093 13566
## # ... with 20 more rows, and 7 more variables: `Self Isolation` <dbl>,
## # retail_and_recreation_percent_change_from_baseline <dbl>,
## # grocery_and_pharmacy_percent_change_from_baseline <dbl>,
## # parks_percent_change_from_baseline <dbl>,
## # transit_stations_percent_change_from_baseline <dbl>,
## # workplaces_percent_change_from_baseline <dbl>,
## # residential_percent_change_from_baseline <dbl>library(jmv)
descriptives(mobilityjakarta, vars = vars(Dirawat, parks_percent_change_from_baseline), freq = TRUE)##
## DESCRIPTIVES
##
## Descriptives
## ------------------------------------------------------------------------
## Dirawat parks_percent_change_from_baseline
## ------------------------------------------------------------------------
## N 30 30
## Missing 0 0
## Mean 204.5000 -22.80000
## Median 198.0000 -22.00000
## Standard deviation 30.06631 4.999310
## Minimum 168.0000 -35.00000
## Maximum 273.0000 -13.00000
## ------------------------------------------------------------------------# 2. 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")
library(npreg)
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.07486715 lambda = 2.070972e-07
## Equivalent Degrees of Freedom (Df) 10.52988
## Penalized Criterion (RSS) 1220.487
## Generalized Cross-Validation (GCV) 96.58671mod.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.05018013 lambda= 5.636477e-06 (14 iterations)
## Equivalent Degrees of Freedom (Df): 11.14365
## Penalized Criterion (RSS): 1196.037
## GCV: 100.9138sqrt(mean(( mod.ss$y - mod.smsp$y )^2))## [1] 1.507532sqrt(mean(( mobilityjakarta$Dirawat - mod.ss$y )^2))## [1] 6.378315sqrt(mean(( mobilityjakarta$Dirawat - mod.smsp$y )^2))## [1] 6.314102mobilityjakarta$prediksi_ss <- mod.ss$y
mobilityjakarta$prediksi_smsp <- mod.smsp$y
mobilityjakarta## # A tibble: 30 x 15
## Tanggal Nama_Kota POSITIF Dirawat Sembuh Meninggal
## <dttm> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 2021-11-01 00:00:00 DKI Jakarta 861623 273 847066 13562
## 2 2021-11-02 00:00:00 DKI Jakarta 861700 266 847201 13562
## 3 2021-11-03 00:00:00 DKI Jakarta 861832 258 847353 13563
## 4 2021-11-04 00:00:00 DKI Jakarta 861943 258 847403 13564
## 5 2021-11-05 00:00:00 DKI Jakarta 862061 256 847438 13565
## 6 2021-11-06 00:00:00 DKI Jakarta 862130 231 847685 13566
## 7 2021-11-07 00:00:00 DKI Jakarta 862247 220 847808 13566
## 8 2021-11-08 00:00:00 DKI Jakarta 862276 198 847885 13566
## 9 2021-11-09 00:00:00 DKI Jakarta 862370 193 848034 13566
## 10 2021-11-10 00:00:00 DKI Jakarta 862465 198 848093 13566
## # ... with 20 more rows, and 9 more variables: `Self Isolation` <dbl>,
## # retail_and_recreation_percent_change_from_baseline <dbl>,
## # grocery_and_pharmacy_percent_change_from_baseline <dbl>,
## # parks_percent_change_from_baseline <dbl>,
## # transit_stations_percent_change_from_baseline <dbl>,
## # workplaces_percent_change_from_baseline <dbl>,
## # residential_percent_change_from_baseline <dbl>, prediksi_ss <dbl>, ...# plot results
plot(mobilityjakarta$Tanggal, mobilityjakarta$Dirawat)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_ss , lty = 2, col = 'red', lwd = 4)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_smsp, lty = 4, col = 'green', lwd = 2)
# plot method
plot(mobilityjakarta$Tanggal, mobilityjakarta$Dirawat, lty = 10, col = 'Green', lwd =5)
#plot(mod.ss)
# add lm fit
abline(coef(lm( mobilityjakarta$Dirawat ~ mobilityjakarta$Tanggal , data = mobilityjakarta)), lty = 10, col = 'Yellow', 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("Green","Red","Yellow"), lwd = 3, bty = "p")
library(npreg)
# fit using ss
mod.ss1 <- ss( mobilityjakarta$Tanggal,mobilityjakarta$parks_percent_change_from_baseline, nknots = 10)
mod.ss1##
## Call:
## ss(x = mobilityjakarta$Tanggal, y = mobilityjakarta$parks_percent_change_from_baseline,
## nknots = 10)
##
## Smoothing Parameter spar = 0.7161693 lambda = 0.008900255
## Equivalent Degrees of Freedom (Df) 2.243905
## Penalized Criterion (RSS) 627.9689
## Generalized Cross-Validation (GCV) 24.45359mod.smsp1 <- smooth.spline( mobilityjakarta$Tanggal,mobilityjakarta$parks_percent_change_from_baseline, nknots = 10)
mod.smsp1## Call:
## smooth.spline(x = mobilityjakarta$Tanggal, y = mobilityjakarta$parks_percent_change_from_baseline,
## nknots = 10)
##
## Smoothing Parameter spar= 0.697281 lambda= 0.2670617 (15 iterations)
## Equivalent Degrees of Freedom (Df): 2.243935
## Penalized Criterion (RSS): 627.966
## GCV: 24.45353mobilityjakarta$prediksi_ss1 <- mod.ss1$y
mobilityjakarta$prediksi_smsp1 <- mod.smsp1$y
mobilityjakarta## # A tibble: 30 x 17
## Tanggal Nama_Kota POSITIF Dirawat Sembuh Meninggal
## <dttm> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 2021-11-01 00:00:00 DKI Jakarta 861623 273 847066 13562
## 2 2021-11-02 00:00:00 DKI Jakarta 861700 266 847201 13562
## 3 2021-11-03 00:00:00 DKI Jakarta 861832 258 847353 13563
## 4 2021-11-04 00:00:00 DKI Jakarta 861943 258 847403 13564
## 5 2021-11-05 00:00:00 DKI Jakarta 862061 256 847438 13565
## 6 2021-11-06 00:00:00 DKI Jakarta 862130 231 847685 13566
## 7 2021-11-07 00:00:00 DKI Jakarta 862247 220 847808 13566
## 8 2021-11-08 00:00:00 DKI Jakarta 862276 198 847885 13566
## 9 2021-11-09 00:00:00 DKI Jakarta 862370 193 848034 13566
## 10 2021-11-10 00:00:00 DKI Jakarta 862465 198 848093 13566
## # ... with 20 more rows, and 11 more variables: `Self Isolation` <dbl>,
## # retail_and_recreation_percent_change_from_baseline <dbl>,
## # grocery_and_pharmacy_percent_change_from_baseline <dbl>,
## # parks_percent_change_from_baseline <dbl>,
## # transit_stations_percent_change_from_baseline <dbl>,
## # workplaces_percent_change_from_baseline <dbl>,
## # residential_percent_change_from_baseline <dbl>, prediksi_ss <dbl>, ...# plot results
plot(mobilityjakarta$Tanggal, mobilityjakarta$parks_percent_change_from_baseline, lty = 10, col = 'yellow', lwd =5)
# add lm fit
abline(coef(lm( mobilityjakarta$parks_percent_change_from_baseline ~ mobilityjakarta$Tanggal , data = mobilityjakarta)), lty = 10, col = 'green', lwd =5)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_ss1 , lty = 2, col = 'yellow', lwd = 4)
lines(mobilityjakarta$Tanggal, mobilityjakarta$prediksi_smsp1, lty = 4, col = 'red', lwd = 2)
2.Korelasi Pearson Korelasi
Pearson digunakan untuk melakukan analisis korelasi Pearson cukup
sederhana. Semisal kita ingin mengetahui hubungan antara
mobilityjakartaDirawat dan
mobilityjakartaparks_percent_change_from_baseline
cor.test(mobilityjakarta$Dirawat,mobilityjakarta$parks_percent_change_from_baseline)##
## Pearson's product-moment correlation
##
## data: mobilityjakarta$Dirawat and mobilityjakarta$parks_percent_change_from_baseline
## t = -1.2987, df = 28, p-value = 0.2046
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5512856 0.1333649
## sample estimates:
## cor
## -0.2383567model <- lm(mobilityjakarta$Dirawat ~ mobilityjakarta$parks_percent_change_from_baseline, data = mobilityjakarta)
summary(model)##
## Call:
## lm(formula = mobilityjakarta$Dirawat ~ mobilityjakarta$parks_percent_change_from_baseline,
## data = mobilityjakarta)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.821 -23.061 -5.636 8.290 62.681
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 171.816 25.745 6.674
## mobilityjakarta$parks_percent_change_from_baseline -1.433 1.104 -1.299
## Pr(>|t|)
## (Intercept) 3.05e-07 ***
## mobilityjakarta$parks_percent_change_from_baseline 0.205
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 29.72 on 28 degrees of freedom
## Multiple R-squared: 0.05681, Adjusted R-squared: 0.02313
## F-statistic: 1.687 on 1 and 28 DF, p-value: 0.2046model$coefficients## (Intercept)
## 171.816225
## mobilityjakarta$parks_percent_change_from_baseline
## -1.433499model$fitted.values## 1 2 3 4 5 6 7 8
## 219.1217 213.3877 209.0872 207.6537 193.3187 194.7522 221.9887 203.3532
## 9 10 11 12 13 14 15 16
## 209.0872 207.6537 200.4862 203.3532 214.8212 200.4862 206.2202 204.7867
## 17 18 19 20 21 22 23 24
## 201.9197 197.6192 190.4517 206.2202 201.9197 201.9197 207.6537 201.9197
## 25 26 27 28 29 30
## 197.6192 203.3532 204.7867 211.9542 196.1857 201.9197mobilityjakarta$prediksi_linear <- model$fitted.values
mobilityjakarta## # A tibble: 30 x 18
## Tanggal Nama_Kota POSITIF Dirawat Sembuh Meninggal
## <dttm> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 2021-11-01 00:00:00 DKI Jakarta 861623 273 847066 13562
## 2 2021-11-02 00:00:00 DKI Jakarta 861700 266 847201 13562
## 3 2021-11-03 00:00:00 DKI Jakarta 861832 258 847353 13563
## 4 2021-11-04 00:00:00 DKI Jakarta 861943 258 847403 13564
## 5 2021-11-05 00:00:00 DKI Jakarta 862061 256 847438 13565
## 6 2021-11-06 00:00:00 DKI Jakarta 862130 231 847685 13566
## 7 2021-11-07 00:00:00 DKI Jakarta 862247 220 847808 13566
## 8 2021-11-08 00:00:00 DKI Jakarta 862276 198 847885 13566
## 9 2021-11-09 00:00:00 DKI Jakarta 862370 193 848034 13566
## 10 2021-11-10 00:00:00 DKI Jakarta 862465 198 848093 13566
## # ... with 20 more rows, and 12 more variables: `Self Isolation` <dbl>,
## # retail_and_recreation_percent_change_from_baseline <dbl>,
## # grocery_and_pharmacy_percent_change_from_baseline <dbl>,
## # parks_percent_change_from_baseline <dbl>,
## # transit_stations_percent_change_from_baseline <dbl>,
## # workplaces_percent_change_from_baseline <dbl>,
## # residential_percent_change_from_baseline <dbl>, prediksi_ss <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, 0, 1, alpha = 0.5), ...)
# lines(density(x), col = 2, lwd = 2) # Uncomment to add density lines
}
# Equivalent with a formula
pairs(mobilityjakarta$Dirawat~mobilityjakarta$parks_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 =5)
#plot(mod.ss)
# add lm fit
abline(coef(lm( mobilityjakarta$Dirawat ~ mobilityjakarta$Tanggal , data = mobilityjakarta)), lty = 10, col = 'green', 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 = 'Orange', 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","Orange","green","Pink"), lwd = 3, bty = "p")
3. Referensi https://github.com/juba/rmdformats
https://medium.com/machine-learning-id/simple-linear-regression-r-ea42aed17c52