Analisa Data dengan Regresi Linear dengan tiga model
Prof Dr Suhartono M.Kom UIN Maulana Malik Ibrahim Malang
4/16/2022
## Import the data
library(readxl)
mobilityjakarta <- read_excel(path = "DataDirawatJanuari2021.xlsx")
mobilityjakarta## # A tibble: 30 x 12
## Tanggal DIRAWAT retail_and_recreation grocery_and_pharmacy parks
## <dttm> <dbl> <dbl> <dbl> <dbl>
## 1 2021-01-01 00:00:00 5789 -25 -8 18
## 2 2021-01-02 00:00:00 4599 -20 3 0
## 3 2021-01-03 00:00:00 4410 -21 -2 -6
## 4 2021-01-04 00:00:00 4499 -15 5 -5
## 5 2021-01-05 00:00:00 4479 -18 3 -10
## 6 2021-01-06 00:00:00 4254 -15 4 -8
## 7 2021-01-07 00:00:00 4601 -16 4 -10
## 8 2021-01-08 00:00:00 4236 -18 0 -16
## 9 2021-01-09 00:00:00 5237 -22 -4 -22
## 10 2021-01-10 00:00:00 4863 -20 -1 -16
## # ... with 20 more rows, and 7 more variables: transit_stations <dbl>,
## # workplaces <dbl>, residential <dbl>, Meninggal <dbl>, POSITIF <dbl>,
## # Self_isolation <dbl>, Sembuh <dbl>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.
mobilityjakarta## # A tibble: 30 x 12
## Tanggal DIRAWAT retail_and_recreation grocery_and_pharmacy parks
## <dttm> <dbl> <dbl> <dbl> <dbl>
## 1 2021-01-01 00:00:00 5789 -25 -8 18
## 2 2021-01-02 00:00:00 4599 -20 3 0
## 3 2021-01-03 00:00:00 4410 -21 -2 -6
## 4 2021-01-04 00:00:00 4499 -15 5 -5
## 5 2021-01-05 00:00:00 4479 -18 3 -10
## 6 2021-01-06 00:00:00 4254 -15 4 -8
## 7 2021-01-07 00:00:00 4601 -16 4 -10
## 8 2021-01-08 00:00:00 4236 -18 0 -16
## 9 2021-01-09 00:00:00 5237 -22 -4 -22
## 10 2021-01-10 00:00:00 4863 -20 -1 -16
## # ... with 20 more rows, and 7 more variables: transit_stations <dbl>,
## # workplaces <dbl>, residential <dbl>, Meninggal <dbl>, POSITIF <dbl>,
## # Self_isolation <dbl>, Sembuh <dbl>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(DIRAWAT, retail_and_recreation), freq = TRUE)##
## DESCRIPTIVES
##
## Descriptives
## -----------------------------------------------------------
## DIRAWAT retail_and_recreation
## -----------------------------------------------------------
## N 30 30
## Missing 0 0
## Mean 4545.500 -24.33333
## Median 4532.000 -25.00000
## Standard deviation 532.6330 5.128240
## Minimum 3056.000 -32.00000
## Maximum 5789.000 -15.00000
## -----------------------------------------------------------hist(mobilityjakarta$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(mobilityjakarta$DIRAWAT) # Menghitung nilai mean variabel wt
abline(v=mean, col="red",lwd=2) # Menambahkan garis vertical pada plot
# 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)## 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.3360082 lambda = 1.59531e-05
## Equivalent Degrees of Freedom (Df) 6.392432
## Penalized Criterion (RSS) 2923328
## Generalized Cross-Validation (GCV) 157360.7# 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.3068372 lambda= 0.0004029795 (13 iterations)
## Equivalent Degrees of Freedom (Df): 6.610864
## Penalized Criterion (RSS): 2864717
## GCV: 157099.4# rmse between solutions
sqrt(mean(( mod.ss$y - mod.smsp$y )^2))## [1] 7.993845# rmse between solutions and f(x)
sqrt(mean(( mobilityjakarta$DIRAWAT - mod.ss$y )^2))## [1] 312.1607sqrt(mean(( mobilityjakarta$DIRAWAT - mod.smsp$y )^2))## [1] 309.0155mobilityjakarta$prediksi_ss <- mod.ss$y
mobilityjakarta$prediksi_smsp <- mod.smsp$y
mobilityjakarta## # A tibble: 30 x 14
## Tanggal DIRAWAT retail_and_recreation grocery_and_pharmacy parks
## <dttm> <dbl> <dbl> <dbl> <dbl>
## 1 2021-01-01 00:00:00 5789 -25 -8 18
## 2 2021-01-02 00:00:00 4599 -20 3 0
## 3 2021-01-03 00:00:00 4410 -21 -2 -6
## 4 2021-01-04 00:00:00 4499 -15 5 -5
## 5 2021-01-05 00:00:00 4479 -18 3 -10
## 6 2021-01-06 00:00:00 4254 -15 4 -8
## 7 2021-01-07 00:00:00 4601 -16 4 -10
## 8 2021-01-08 00:00:00 4236 -18 0 -16
## 9 2021-01-09 00:00:00 5237 -22 -4 -22
## 10 2021-01-10 00:00:00 4863 -20 -1 -16
## # ... with 20 more rows, and 9 more variables: transit_stations <dbl>,
## # workplaces <dbl>, residential <dbl>, Meninggal <dbl>, POSITIF <dbl>,
## # Self_isolation <dbl>, Sembuh <dbl>, prediksi_ss <dbl>, prediksi_smsp <dbl># 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 = 'red', lwd = 2)
# 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 = '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, nknots = 10)
mod.ss1##
## Call:
## ss(x = mobilityjakarta$Tanggal, y = mobilityjakarta$retail_and_recreation,
## nknots = 10)
##
## Smoothing Parameter spar = 0.1755807 lambda = 1.106118e-06
## Equivalent Degrees of Freedom (Df) 9.421911
## Penalized Criterion (RSS) 66.91026
## Generalized Cross-Validation (GCV) 4.740279# fit using smooth.spline
mod.smsp1 <- smooth.spline( mobilityjakarta$Tanggal,mobilityjakarta$retail_and_recreation, nknots = 10)
mod.smsp1## Call:
## smooth.spline(x = mobilityjakarta$Tanggal, y = mobilityjakarta$retail_and_recreation,
## nknots = 10)
##
## Smoothing Parameter spar= 0.1710485 lambda= 4.209671e-05 (13 iterations)
## Equivalent Degrees of Freedom (Df): 9.34048
## Penalized Criterion (RSS): 68.23911
## GCV: 4.796386mobilityjakarta$prediksi_ss1 <- mod.ss1$y
mobilityjakarta$prediksi_smsp1 <- mod.smsp1$y
mobilityjakarta## # A tibble: 30 x 16
## Tanggal DIRAWAT retail_and_recreation grocery_and_pharmacy parks
## <dttm> <dbl> <dbl> <dbl> <dbl>
## 1 2021-01-01 00:00:00 5789 -25 -8 18
## 2 2021-01-02 00:00:00 4599 -20 3 0
## 3 2021-01-03 00:00:00 4410 -21 -2 -6
## 4 2021-01-04 00:00:00 4499 -15 5 -5
## 5 2021-01-05 00:00:00 4479 -18 3 -10
## 6 2021-01-06 00:00:00 4254 -15 4 -8
## 7 2021-01-07 00:00:00 4601 -16 4 -10
## 8 2021-01-08 00:00:00 4236 -18 0 -16
## 9 2021-01-09 00:00:00 5237 -22 -4 -22
## 10 2021-01-10 00:00:00 4863 -20 -1 -16
## # ... with 20 more rows, and 11 more variables: transit_stations <dbl>,
## # workplaces <dbl>, residential <dbl>, Meninggal <dbl>, POSITIF <dbl>,
## # Self_isolation <dbl>, Sembuh <dbl>, prediksi_ss <dbl>, prediksi_smsp <dbl>,
## # prediksi_ss1 <dbl>, prediksi_smsp1 <dbl># plot results
plot(mobilityjakarta$Tanggal, mobilityjakarta$retail_and_recreation, lty = 10, col = 'black', lwd =5)
# add lm fit
abline(coef(lm( mobilityjakarta$retail_and_recreation ~ 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 mobilityjakarta\(DIRAWAT dan mobilityjakarta\)retail_and_recreation,
cor.test(mobilityjakarta$DIRAWAT,mobilityjakarta$retail_and_recreation)##
## Pearson's product-moment correlation
##
## data: mobilityjakarta$DIRAWAT and mobilityjakarta$retail_and_recreation
## t = -0.84284, df = 28, p-value = 0.4065
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4898100 0.2151644
## sample estimates:
## cor
## -0.1572981Dari output tersebut dapat kita simpulkan bahwa tidak ada hubungan yang signifikan antara mobilityjakarta\(DIRAWAT dan mobilityjakarta\)retail_and_recreation (p-value<0,01). Nilai p-value dalam output R dituliskan 0,4065 Nilai koefisien korelasi r adalah sebesar -0,1572981 yang menunjukkan hubungan yang sedang dan terbalik.
model <- lm(mobilityjakarta$DIRAWAT ~ mobilityjakarta$retail_and_recreation, data = mobilityjakarta)
summary(model)##
## Call:
## lm(formula = mobilityjakarta$DIRAWAT ~ mobilityjakarta$retail_and_recreation,
## data = mobilityjakarta)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1533.07 -189.28 -5.75 198.61 1232.61
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4147.96 481.69 8.611 2.34e-09 ***
## mobilityjakarta$retail_and_recreation -16.34 19.38 -0.843 0.406
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 535.3 on 28 degrees of freedom
## Multiple R-squared: 0.02474, Adjusted R-squared: -0.01009
## F-statistic: 0.7104 on 1 and 28 DF, p-value: 0.4065Rumus linear regresi menjadi y = 4147.96 + (-16.34) mobilityjakarta$retail_and_recreation
model$coefficients## (Intercept) mobilityjakarta$retail_and_recreation
## 4147.95629 -16.33741model$fitted.values## 1 2 3 4 5 6 7 8
## 4556.392 4474.705 4491.042 4393.017 4442.030 4393.017 4409.355 4442.030
## 9 10 11 12 13 14 15 16
## 4507.379 4474.705 4442.030 4507.379 4523.717 4556.392 4638.079 4670.753
## 17 18 19 20 21 22 23 24
## 4670.753 4638.079 4621.741 4605.404 4589.066 4621.741 4589.066 4670.753
## 25 26 27 28 29 30
## 4556.392 4589.066 4572.729 4540.054 4589.066 4589.066mobilityjakarta$prediksi_linear <- model$fitted.values
mobilityjakarta## # A tibble: 30 x 17
## Tanggal DIRAWAT retail_and_recreation grocery_and_pharmacy parks
## <dttm> <dbl> <dbl> <dbl> <dbl>
## 1 2021-01-01 00:00:00 5789 -25 -8 18
## 2 2021-01-02 00:00:00 4599 -20 3 0
## 3 2021-01-03 00:00:00 4410 -21 -2 -6
## 4 2021-01-04 00:00:00 4499 -15 5 -5
## 5 2021-01-05 00:00:00 4479 -18 3 -10
## 6 2021-01-06 00:00:00 4254 -15 4 -8
## 7 2021-01-07 00:00:00 4601 -16 4 -10
## 8 2021-01-08 00:00:00 4236 -18 0 -16
## 9 2021-01-09 00:00:00 5237 -22 -4 -22
## 10 2021-01-10 00:00:00 4863 -20 -1 -16
## # ... with 20 more rows, and 12 more variables: transit_stations <dbl>,
## # workplaces <dbl>, residential <dbl>, Meninggal <dbl>, POSITIF <dbl>,
## # Self_isolation <dbl>, Sembuh <dbl>, prediksi_ss <dbl>, prediksi_smsp <dbl>,
## # prediksi_ss1 <dbl>, prediksi_smsp1 <dbl>, prediksi_linear <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, 1, 1, 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, 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 = '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")