Penerapan Ridge Regression dan Regresi Least Absolute Shrinkage and Selection Operator (LASSO) Determinan Kejadian Stunting di Provinsi Papua Barat tahun 2022
Sapta Hastho Ponco
2024-01-10
Percepatan penurunan stunting pada anak balita merupakan salah satu agenda utama Pemerintah Indonesia. Pemerintah memandang bahwa sebagai upaya mewujudkan sumber daya manusia yang sehat, cerdas, dan produktif, serta pencapaian Sustainable Development Goals (SDGs) khususnya pada Tujuan 2, yaitu menghilangkan kelaparan, mencapai ketahanan pangan dan gizi yang baik, serta meningkatkan pertanian berkelanjutan, perlu dilakukan percepatan penurunan stunting.
Kondisi stunting di Provinsi Papua Barat hingga tahun 2022 masih tergolong tinggi. Oleh karena itu ingin dilihat determinannya dari indikator-indikator sebagai berikut:
X1: Imunisasi Dasar Lengkap
X2: Penolong persalinan di fasilitas kesehatan oleh tenaga kesehatan
X3: Penggunaan alat KB Modern
X4: ASI Eksklusif**
X5: Complementary Feeding
X6: Akses Air Minum Layak
X7: Akses Sanitasi Layak
X8: Pendidikan Anak Usia Dini (PAUD)
X9: Kepemilikan JKN/jamkesda
X10: Penerima KPS/KKS atau bantuan pangan
Pada kesempatan ini akan dilakukan perbandingan performa model dari regresi OLS, regresi ridge dan regresi LASSO. Dari perbandingan ini ingin diketahui manakah dari ketiga model tersebut yang lebih baik dalam mengidentifikasi determinan stunting di Provinsi Papua Barat. Perbandingan dilakukan dengan melihat nilai \(R^{2}\) pada setiap model.
Load Packages
pacman::p_load(car, corrplot, glmnet, hrbrthemes, factoextra, readxl, dplyr, lmtest)Setting Working Directory
setwd("D:\\Share\\")Load Dataset
data_penelitian <- read_excel("Stunt2022_9100.xlsx")
data = data_penelitian
View(data)A. PERSIAPAN DATA
Checking Distribution & Missing Data
str(data)## tibble [13 × 12] (S3: tbl_df/tbl/data.frame)
## $ kabupaten_kota: chr [1:13] "9101. FAKFAK" "9102. KAIMANA" "9103. TELUK WONDAMA" "9104. TELUK BINTUNI" ...
## $ y : num [1:13] 29 29.2 26.1 22.8 36.6 36.7 23.8 31.1 39.1 27.3 ...
## $ x1 : num [1:13] 63.8 69.8 39.8 61.4 56.6 ...
## $ x2 : num [1:13] 76.7 51.4 45.5 71.3 88.7 ...
## $ x3 : num [1:13] 36 22.3 22.3 34.2 39.4 ...
## $ x4 : num [1:13] 63 49.6 50.7 80.3 61.7 ...
## $ x5 : num [1:13] 64.6 52.3 51.9 60.6 61.4 ...
## $ x6 : num [1:13] 98.5 85.9 24 87.5 81.1 ...
## $ x7 : num [1:13] 69.2 76.5 53 78.6 82.6 ...
## $ x8 : num [1:13] 31.8 21.5 32.2 23.2 26.8 ...
## $ x9 : num [1:13] 74 62.4 82.1 99.9 60 ...
## $ x10 : num [1:13] 33.4 23.8 30.5 19.5 12.7 ...summary(data)## kabupaten_kota y x1 x2
## Length:13 Min. :22.80 Min. : 0.00 Min. :26.25
## Class :character 1st Qu.:27.20 1st Qu.:45.16 1st Qu.:45.50
## Mode :character Median :29.00 Median :55.59 Median :71.28
## Mean :31.35 Mean :51.40 Mean :62.57
## 3rd Qu.:36.60 3rd Qu.:63.76 3rd Qu.:79.04
## Max. :51.50 Max. :70.64 Max. :88.67
## x3 x4 x5 x6
## Min. : 3.034 Min. :24.83 Min. :33.06 Min. :23.97
## 1st Qu.:20.840 1st Qu.:50.68 1st Qu.:51.91 1st Qu.:69.92
## Median :30.594 Median :61.66 Median :61.45 Median :80.91
## Mean :26.312 Mean :59.25 Mean :60.90 Mean :73.86
## 3rd Qu.:34.206 3rd Qu.:71.69 3rd Qu.:70.97 3rd Qu.:87.54
## Max. :39.447 Max. :80.32 Max. :84.64 Max. :98.49
## x7 x8 x9 x10
## Min. :33.74 Min. : 2.231 Min. :44.01 Min. : 2.737
## 1st Qu.:56.57 1st Qu.:21.496 1st Qu.:62.40 1st Qu.:19.487
## Median :76.48 Median :24.092 Median :73.39 Median :21.216
## Mean :68.52 Mean :23.631 Mean :72.42 Mean :22.292
## 3rd Qu.:82.58 3rd Qu.:31.753 3rd Qu.:78.74 3rd Qu.:29.919
## Max. :88.42 Max. :35.507 Max. :99.85 Max. :41.600Pemeriksaan missing data:
colSums(is.na(data))## kabupaten_kota y x1 x2 x3
## 0 0 0 0 0
## x4 x5 x6 x7 x8
## 0 0 0 0 0
## x9 x10
## 0 0Dari output di atas tampak tidak ditemukan missing data untuk seluruh variabel. Jika ditemukan missing data, maka perlu dilakukan imputasi.
Plot Data
Plot masing-masing variabel setelah dilakukan imputasi data
par(mfrow=c(3,4))
ploty = plot(density(data$y))
plotx1 = plot(density(data$x1))
plotx2 = plot(density(data$x2))
plotx3 = plot(density(data$x3))
plotx4 = plot(density(data$x4))
plotx5 = plot(density(data$x5))
plotx6 = plot(density(data$x6))
plotx7 = plot(density(data$x7))
plotx8 = plot(density(data$x8))
plotx9 = plot(density(data$x9))
plotx10 = plot(density(data$x10)) 
Plot variabel menurut Kabupaten/kota
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
library(grid)
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols))
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
library(ggplot2)
ploty <- ggplot(data = data, aes(x = factor(kabupaten_kota), y = y)) + geom_point(colour="red") +
theme(axis.text.x = element_text(angle = 45, size = 6, hjust = 1)) + ggtitle("Stunting")+ labs(x="",y="")
plotx1 <- ggplot(data = data, aes(x = factor(kabupaten_kota), y = x1)) + geom_point(colour="red") +
theme(axis.text.x = element_text(angle = 45, size = 6, hjust = 1)) + ggtitle("X1 Imunisasi")+ labs(x="",y="")
plotx2 <- ggplot(data = data, aes(x = factor(kabupaten_kota), y = x2)) + geom_point(colour="red") +
theme(axis.text.x = element_text(angle = 45, size = 6, hjust = 1)) + ggtitle("X2 Persalinan")+ labs(x="",y="")
plotx3 <- ggplot(data = data, aes(x = factor(kabupaten_kota), y = x3)) + geom_point(colour="red") +
theme(axis.text.x = element_text(angle = 45, size = 6, hjust = 1)) + ggtitle("X3 KB Modern")+ labs(x="",y="")
plotx4 <- ggplot(data = data, aes(x = factor(kabupaten_kota), y = x4)) + geom_point(colour="red") +
theme(axis.text.x = element_text(angle = 45, size = 6, hjust = 1)) + ggtitle("X4 ASI Eks")+ labs(x="",y="")
plotx5 <- ggplot(data = data, aes(x = factor(kabupaten_kota), y = x5)) + geom_point(colour="red") +
theme(axis.text.x = element_text(angle = 45, size = 6, hjust = 1)) + ggtitle("X5 MP-ASI")+ labs(x="",y="")
multiplot(ploty,plotx1,plotx2,plotx3,plotx4,plotx5,cols=3)
plotx6 <- ggplot(data = data, aes(x = factor(kabupaten_kota), y = x6)) + geom_point(colour="red") +
theme(axis.text.x = element_text(angle = 45, size = 6, hjust = 1)) + ggtitle("X6 Air Minum Layak")+ labs(x="",y="")
plotx7 <- ggplot(data = data, aes(x = factor(kabupaten_kota), y = x7)) + geom_point(colour="red") +
theme(axis.text.x = element_text(angle = 45, size = 6, hjust = 1)) + ggtitle("X7 Sanitasi Layak")+ labs(x="",y="")
plotx8 <- ggplot(data = data, aes(x = factor(kabupaten_kota), y = x8)) + geom_point(colour="red") +
theme(axis.text.x = element_text(angle = 45, size = 6, hjust = 1)) + ggtitle("X8 APK PAUD")+ labs(x="",y="")
plotx9 <- ggplot(data = data, aes(x = factor(kabupaten_kota), y = x9)) + geom_point(colour="red") +
theme(axis.text.x = element_text(angle = 45, size = 6, hjust = 1)) + ggtitle("X9 Kepemilikan JKN/Jamkesda")+ labs(x="",y="")
plotx10 <- ggplot(data = data, aes(x = factor(kabupaten_kota), y = x10)) + geom_point(colour="red") +
theme(axis.text.x = element_text(angle = 45, size = 6, hjust = 1)) + ggtitle("X10 Penerima KPS/KKS atau Bantuan Pangan")+ labs(x="",y="")
multiplot(plotx6,plotx7,plotx8,plotx9,plotx10,cols=3)
Dari visualisasi plot di atas, tampak data bersifat menyebar dan cenderung tidak menunjukkan hubungan yang linier.
library(hrbrthemes)
data %>%
filter(!is.na(data$y)) %>%
arrange(data$y) %>%
tail(50) %>%
mutate(Nama.Daerah=factor(kabupaten_kota, kabupaten_kota)) %>%
ggplot( aes(x=kabupaten_kota, y=y) ) +
geom_segment( aes(x=kabupaten_kota,xend=kabupaten_kota, y=0, yend=y), color="#3B1877") +
geom_point(size=3, color="#E69A8D") +
coord_flip() +
theme_ipsum() +
theme(
panel.grid.minor.y = element_blank(),
panel.grid.major.y = element_blank(),
legend.position="none"
) +
xlab("") +
ylab("Stunting tahun 2022")
Kabupaten Pegunungan Arfak dan Tambrauw merupakan dua kabupaten/kota yang memiliki tingkat stunting tertinggi dibanding kabupaten/kota lain di Provinsi Papua Barat.
Hubungan antar variabel
Eksplorasi dan visualisasi dapat menggunakan matriks korelasi untuk melihat ada tidaknya hubungan antar variabel.
library(corrplot)
cordata = data[,2:12]
cordata = cor(cordata)
corrplot(cordata,method = 'number')
Berdasarkan matriks korelasi di atas terlihat bahwa Y memiliki korelasi negatif yang cukup erat dengan variabel X8, X1, X9, X3 dan X7. Sementara hubungan keeratan antara variabel Y dengan X4 dan X6 sangat rendah.
B.1 ORDINARY LEAST SQUARES (OLS) REGRESSION
Analisis regresi merupakan metode yang digunakan untuk menentukan hubungan antara satu variabel respon dengan satu atau lebih variabel bebas. Salah satu metode yang digunakan untuk mengestimasi parameter model regresi yaitu metode kuadrat terkecil atau Ordinary Least Square (OLS). Tujuan utama dari metode ini adalah mengestimasi koefisien regresi untuk meminimumkan jumlah kuadrat error. Penggunaan metode OLS memerlukan beberapa asumsi yang harus dipenuhi. Apabila asumsi-asumsi tersebut tidak terpenuhi, maka estimator yang dihasilkan bersifat bias dan interpretasi hasil yang diberikan pun menjadi tidak valid.
Model_klasik = lm(formula = y~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10, data = data)
sum = summary(Model_klasik)
print(sum)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 +
## x10, data = data)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## -0.37258 -0.30965 -0.90576 -0.35698 -0.72126 0.66174 0.40937 -0.14932
## 9 10 11 12 13
## 0.03711 0.62571 1.82504 0.32919 -1.07260
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 91.76738 8.94636 10.258 0.00937 **
## x1 -0.01222 0.10571 -0.116 0.91853
## x2 0.12764 0.06423 1.987 0.18524
## x3 0.25624 0.10412 2.461 0.13295
## x4 -0.27621 0.08825 -3.130 0.08871 .
## x5 -0.07935 0.08329 -0.953 0.44131
## x6 0.12811 0.05642 2.271 0.15117
## x7 -0.56719 0.09991 -5.677 0.02965 *
## x8 -0.07613 0.11048 -0.689 0.56194
## x9 -0.11717 0.07020 -1.669 0.23703
## x10 -0.61139 0.10684 -5.723 0.02920 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.912 on 2 degrees of freedom
## Multiple R-squared: 0.9901, Adjusted R-squared: 0.9408
## F-statistic: 20.06 on 10 and 2 DF, p-value: 0.04838cat("AIC = ",AIC(Model_klasik))## AIC = 53.41361cat("Adj R2 = ",sum$adj.r.squared)## Adj R2 = 0.9407846Berdasarkan hasil uji signifikansi di atas, tidak seluruh variabel signifikan pada taraf nyata 5%. Variabel yang memiliki p-value < α = 0.05 adalah X7 dan X10 sehingga kedua variabel tersebut signifikan terhadap model yang digunakan, yaitu model penuh dengan seluruh variabel bebas.
Adjusted R-Squared yang dihasilkan sebesar 94.08%, yang berarti variabel-variabel yang dilibatkan mampu menjelaskan 94% varians stunting di Papua Barat sedangkan sisanya berasal dari faktor lain.
\(H_{0}\): Tidak ada pengaruh variabel bebas terhadap variabel Y
\(H_{1}\): Paling sedikit ada satu variabel bebas mempengaruhi variabel Y
Dari uji signifikansi F, diperoleh p-value sebesar 0.04838 (p-value < α) sehingga memberikan informasi bahwa cukup bukti untuk menolak \(H_{0}\).
Pemeriksaan Asumsi Multikolinearitas
Multikolinearitas adalah suatu kondisi adanya hubungan yang linear diantara variabel-variabel bebas dalam model regresi. Multikolinearitas dapat dideteksi menggunakan Variance Inflation Factors (VIF). Nilai VIF dapat dihitung dengan rumus: \(VIF_{j}=\frac{1}{1-R_{j}^{2}}\) dengan \(R_{j}\) merupakan koefisien determinasi ke-j, j=1,2,…,k.
# VIF
vif(Model_klasik)## x1 x2 x3 x4 x5 x6 x7 x8
## 13.049053 5.817837 4.093461 5.686733 5.175505 4.801289 11.353699 3.815464
## x9 x10
## 3.185579 4.404467# Tolerance
1/vif(Model_klasik)## x1 x2 x3 x4 x5 x6 x7
## 0.07663391 0.17188517 0.24429204 0.17584790 0.19321787 0.20827742 0.08807702
## x8 x9 x10
## 0.26209133 0.31391468 0.22704224# rata-rata VIF
mean(vif(Model_klasik))## [1] 6.138309Dari hasil di atas terlihat ditemukan adanya multikolinearitas yang ditandai dengan nilai VIF > 10, yaitu pada variabel X1 dan X7 yang memiliki nilai VIF sebesar 13.0490 dan 11.3537.
B.2 REGRESI OLS TANPA X1 dan X7
modelnew <- lm(formula = y~x2+x3+x4+x5+x6+x8+x9+x10, data = data)
sum = summary(modelnew)
print(sum)##
## Call:
## lm(formula = y ~ x2 + x3 + x4 + x5 + x6 + x8 + x9 + x10, data = data)
##
## Residuals:
## 1 2 3 4 5 6 7 8 9 10 11
## 6.623 -4.371 -2.604 -3.351 1.917 2.763 -5.658 1.162 1.958 2.160 1.003
## 12 13
## 1.386 -2.988
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 73.42660 26.13444 2.810 0.0483 *
## x2 -0.10666 0.15488 -0.689 0.5289
## x3 0.13268 0.32211 0.412 0.7015
## x4 0.09049 0.14492 0.624 0.5662
## x5 -0.06259 0.20849 -0.300 0.7790
## x6 -0.08750 0.10904 -0.802 0.4673
## x8 -0.37617 0.31207 -1.205 0.2945
## x9 -0.25629 0.17189 -1.491 0.2102
## x10 -0.29280 0.27748 -1.055 0.3508
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.042 on 4 degrees of freedom
## Multiple R-squared: 0.8029, Adjusted R-squared: 0.4088
## F-statistic: 2.037 on 8 and 4 DF, p-value: 0.2567cat("AIC = ",AIC(modelnew))## AIC = 88.337cat("Adj R2 = ",sum$adj.r.squared)## Adj R2 = 0.4088072Jika variabel X1 dan X7 dihilangkan dari model (karena mengandung multikolinearitas), tampak nilai \(R^{2}\) menurun. Oleh karena itu, perlu kehati-hatian dalam menghilangkan suatu variabel.
B.3 SELEKSI VARIABEL REGRESI
Salah satu cara untuk mengatasai masalah multikolinearitas dalam model adalah dengan menghilangkan salah satu atau beberapa variabel yang kolinier, terutama yang memiliki hubungan koliner kuat dengan variabel lain. Pengeluaran variabel bebas ini harus dilakukan secara hati-hati karena tidak memungkinkan variabel yang dihilangkan justru variabel yang penting (spesification bias).
backward <- step(Model_klasik, direction="backward")## Start: AIC=14.52
## y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10
##
## Df Sum of Sq RSS AIC
## - x1 1 0.049 7.362 12.608
## <none> 7.313 14.521
## - x8 1 1.737 9.050 15.291
## - x5 1 3.318 10.632 17.385
## - x9 1 10.187 17.500 23.864
## - x2 1 14.441 21.754 26.693
## - x6 1 18.852 26.165 29.093
## - x3 1 22.148 29.461 30.635
## - x4 1 35.817 43.130 35.591
## - x7 1 117.845 125.159 49.440
## - x10 1 119.748 127.061 49.636
##
## Step: AIC=12.61
## y ~ x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10
##
## Df Sum of Sq RSS AIC
## <none> 7.362 12.608
## - x8 1 1.802 9.164 13.454
## - x5 1 4.531 11.893 16.843
## - x2 1 15.087 22.449 25.102
## - x9 1 15.737 23.099 25.473
## - x3 1 22.109 29.471 28.640
## - x6 1 23.943 31.305 29.425
## - x4 1 48.340 55.702 36.916
## - x10 1 132.346 139.708 48.870
## - x7 1 138.661 146.023 49.445summary(backward)##
## Call:
## lm(formula = y ~ x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10,
## data = data)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## -0.38980 -0.42572 -0.84984 -0.36890 -0.74082 0.68343 0.38447 -0.15565
## 9 10 11 12 13
## 0.06662 0.52409 1.91392 0.32499 -0.96678
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 91.56897 7.19290 12.730 0.00105 **
## x2 0.12564 0.05067 2.480 0.08931 .
## x3 0.25542 0.08509 3.002 0.05760 .
## x4 -0.27074 0.06100 -4.438 0.02126 *
## x5 -0.07348 0.05407 -1.359 0.26735
## x6 0.12483 0.03996 3.124 0.05232 .
## x7 -0.57147 0.07602 -7.517 0.00488 **
## x8 -0.07726 0.09015 -0.857 0.45447
## x9 -0.12163 0.04803 -2.532 0.08525 .
## x10 -0.61498 0.08374 -7.344 0.00522 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.567 on 3 degrees of freedom
## Multiple R-squared: 0.9901, Adjusted R-squared: 0.9603
## F-statistic: 33.22 on 9 and 3 DF, p-value: 0.007526forward <- step(Model_klasik, direction="forward")## Start: AIC=14.52
## y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10summary(forward)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 +
## x10, data = data)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## -0.37258 -0.30965 -0.90576 -0.35698 -0.72126 0.66174 0.40937 -0.14932
## 9 10 11 12 13
## 0.03711 0.62571 1.82504 0.32919 -1.07260
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 91.76738 8.94636 10.258 0.00937 **
## x1 -0.01222 0.10571 -0.116 0.91853
## x2 0.12764 0.06423 1.987 0.18524
## x3 0.25624 0.10412 2.461 0.13295
## x4 -0.27621 0.08825 -3.130 0.08871 .
## x5 -0.07935 0.08329 -0.953 0.44131
## x6 0.12811 0.05642 2.271 0.15117
## x7 -0.56719 0.09991 -5.677 0.02965 *
## x8 -0.07613 0.11048 -0.689 0.56194
## x9 -0.11717 0.07020 -1.669 0.23703
## x10 -0.61139 0.10684 -5.723 0.02920 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.912 on 2 degrees of freedom
## Multiple R-squared: 0.9901, Adjusted R-squared: 0.9408
## F-statistic: 20.06 on 10 and 2 DF, p-value: 0.04838stepwise <- step(Model_klasik, direction="both")## Start: AIC=14.52
## y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10
##
## Df Sum of Sq RSS AIC
## - x1 1 0.049 7.362 12.608
## <none> 7.313 14.521
## - x8 1 1.737 9.050 15.291
## - x5 1 3.318 10.632 17.385
## - x9 1 10.187 17.500 23.864
## - x2 1 14.441 21.754 26.693
## - x6 1 18.852 26.165 29.093
## - x3 1 22.148 29.461 30.635
## - x4 1 35.817 43.130 35.591
## - x7 1 117.845 125.159 49.440
## - x10 1 119.748 127.061 49.636
##
## Step: AIC=12.61
## y ~ x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10
##
## Df Sum of Sq RSS AIC
## <none> 7.362 12.608
## - x8 1 1.802 9.164 13.454
## + x1 1 0.049 7.313 14.521
## - x5 1 4.531 11.893 16.843
## - x2 1 15.087 22.449 25.102
## - x9 1 15.737 23.099 25.473
## - x3 1 22.109 29.471 28.640
## - x6 1 23.943 31.305 29.425
## - x4 1 48.340 55.702 36.916
## - x10 1 132.346 139.708 48.870
## - x7 1 138.661 146.023 49.445summary(stepwise)##
## Call:
## lm(formula = y ~ x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10,
## data = data)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## -0.38980 -0.42572 -0.84984 -0.36890 -0.74082 0.68343 0.38447 -0.15565
## 9 10 11 12 13
## 0.06662 0.52409 1.91392 0.32499 -0.96678
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 91.56897 7.19290 12.730 0.00105 **
## x2 0.12564 0.05067 2.480 0.08931 .
## x3 0.25542 0.08509 3.002 0.05760 .
## x4 -0.27074 0.06100 -4.438 0.02126 *
## x5 -0.07348 0.05407 -1.359 0.26735
## x6 0.12483 0.03996 3.124 0.05232 .
## x7 -0.57147 0.07602 -7.517 0.00488 **
## x8 -0.07726 0.09015 -0.857 0.45447
## x9 -0.12163 0.04803 -2.532 0.08525 .
## x10 -0.61498 0.08374 -7.344 0.00522 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.567 on 3 degrees of freedom
## Multiple R-squared: 0.9901, Adjusted R-squared: 0.9603
## F-statistic: 33.22 on 9 and 3 DF, p-value: 0.007526Seleksi variabel menggunakan metode backward, forward, dan stepwise hasilnya hampir sama, yaitu menggunakan variabel X2, X3, X4, X5, X6, X7, X8, X9, dan X10.
\(R^{2}\) =96.03%
Oleh karena itu, untuk selanjutnya akan digunakan model dengan menggunakan kesembilan variabel bebas ini, tanpa melibatkan variabel X1.
modelnew1 <- lm(formula = y~x2+x3+x4+x5+x6+x7+x8+x9+x10, data = data)
sum = summary(modelnew1)
print(sum)##
## Call:
## lm(formula = y ~ x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10,
## data = data)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## -0.38980 -0.42572 -0.84984 -0.36890 -0.74082 0.68343 0.38447 -0.15565
## 9 10 11 12 13
## 0.06662 0.52409 1.91392 0.32499 -0.96678
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 91.56897 7.19290 12.730 0.00105 **
## x2 0.12564 0.05067 2.480 0.08931 .
## x3 0.25542 0.08509 3.002 0.05760 .
## x4 -0.27074 0.06100 -4.438 0.02126 *
## x5 -0.07348 0.05407 -1.359 0.26735
## x6 0.12483 0.03996 3.124 0.05232 .
## x7 -0.57147 0.07602 -7.517 0.00488 **
## x8 -0.07726 0.09015 -0.857 0.45447
## x9 -0.12163 0.04803 -2.532 0.08525 .
## x10 -0.61498 0.08374 -7.344 0.00522 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.567 on 3 degrees of freedom
## Multiple R-squared: 0.9901, Adjusted R-squared: 0.9603
## F-statistic: 33.22 on 9 and 3 DF, p-value: 0.007526cat("AIC = ",AIC(modelnew1))## AIC = 51.50018cat("Adj R2 = ",sum$adj.r.squared)## Adj R2 = 0.96025931. Pemeriksaan Asumsi Kenormalan
layout(matrix(c(1,2,3,4),2,2))
plot(modelnew1) #Normal Q-Q
Hipotesis
\(H_{0}\): residual menyebar normal
\(H_{1}\): residual tidak menyebar normal
library(nortest)
lillie.test(modelnew1[['residuals']])##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: modelnew1[["residuals"]]
## D = 0.14271, p-value = 0.6649ks.test(modelnew$residuals, "pnorm", mean=mean(modelnew$residuals), sd=sd(modelnew$residuals))##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: modelnew$residuals
## D = 0.22852, p-value = 0.4407
## alternative hypothesis: two-sidedshapiro.test(modelnew1$residuals)##
## Shapiro-Wilk normality test
##
## data: modelnew1$residuals
## W = 0.91745, p-value = 0.2317ad.test(modelnew1$residuals)##
## Anderson-Darling normality test
##
## data: modelnew1$residuals
## A = 0.36277, p-value = 0.3855Berdasarkan Uji Liliefors, Kolmogorov-Smirnov, Shapiro-Wilk dan Anderson-Darling, residual data menyebar normal dengan p-value > 5%.
2. Pemeriksaan Asumsi Autokorelasi
Autokorelasi dapat didefinisikan terjadinya korelasi diantara data pengamatan sebelumnya, atau dengan kata lain bahwa munculnya suatu data dipengaruhi oleh data sebelumnya.
Gujarati (2006) menyebut keberadaan autokorelasi pada metode OLS memiliki konsekuensi antara lain: estimasi masih linier dan tidak bias, namun estimator-estimator tersebut tidak lagi efisien (memiliki varian terkecil). Oleh karena itu, interval estimasi maupun uji hipotesis yang didasarkan pada distribusi t maupun F tidak tidak dapat digunakan untuk evaluasi hasil regresi.
Hipotesis
\(H_{0}\): \(Cov[\)\(ε_{i}, ε_{j}\)] = 0 (tidak terdapat autokorelasi pada residual)
\(H_{1}\): \(Cov[\)\(ε_{i}, ε_{j}\)] ≠ 0 (terdapat autokorelasi pada residual)
car::dwt(modelnew1)## lag Autocorrelation D-W Statistic p-value
## 1 0.2915785 1.269244 0.572
## Alternative hypothesis: rho != 0randtests::runs.test(modelnew1$residuals)##
## Runs Test
##
## data: modelnew1$residuals
## statistic = -2.4221, runs = 3, n1 = 6, n2 = 6, n = 12, p-value =
## 0.01543
## alternative hypothesis: nonrandomness#plot(x = 1:dim(data)[1],y = modelnew1$residuals,type = 'b', ylab = "Residuals", xlab = "Observation", pch=20)3. Pemeriksaan Asumsi Multikolinearitas
car::vif(modelnew1)## x2 x3 x4 x5 x6 x7 x8 x9
## 5.395591 4.074306 4.048140 3.250397 3.589544 9.795426 3.785949 2.222182
## x10
## 4.032143Nilai VIF pada seluruh variabel bebas < 10 sehingga model sudah terbebas dari multikolinearitas
4. Pemeriksaan Asumsi Heteroskedastisitas
Heteroskedastisitas terjadi jika varians dari error/residual suatu pengamatan ke pengamatan lain terjadi ketidaksamaan (tidak konstan). Bila dalam suatu pengamatan terdapat kasus heteroskedastisitas sedangkan asumsi lain dipenuhi, maka parameter regresi yang diduga dengan metode OLS tidak lagi memiliki sifat efisien (varians minimum) (Draper and Smith, 1998).
Hipotesis
\(H_{0}\): \(Var[ε]\) = \(σ^{2}I\) (varians residual bersifat homogen)
\(H_{1}\): \(Var[ε]\) ≠ \(σ^{2}I\) (varians residual tidak homogen)
Uji statistik yang dapat digunakan untuk menguji apakah varian dari error bersifat homoskedastik atau tidak adalah uji Breusch Pagan Godfrey. Pada prinsipnya, uji ini mencoba mengukur varians \(ε_{i}^{2}\) akibat perubahan nilai-nilai variabel bebasnya.
lmtest::bptest(modelnew1)##
## studentized Breusch-Pagan test
##
## data: modelnew1
## BP = 4.5954, df = 9, p-value = 0.8681p−value > α=0.05 (tidak cukup bukti tolak \(H_{0}\)), artinya varians residual homogen
#plot(modelnew1, 1, pch=20)C. RIDGE REGRESSION
Menurut Montgomery dkk (2006), salah satu cara untuk mengatasi masalah multikolinearitas dalam model adalah menggunakan estimasi regresi Ridge. Regresi ridge diperkenalkan pertama kali oleh Hoer dan R.W. Kennard pada tahun 1962. Regresi ridge adalah satu diantara metode yang digunakan untuk mengatasi masalah multikolinearitas yang merupakan modifikasi dari metode kuadrat terkecil. Modifikasi tersebut dilakukan dengan menambahkan konstanta bias c pada diagonal matriks (\(X^{t}X\)) yang mempengaruhi besarnya koefisien penduga ridge dan penduga yang dihasilkan adalah penduga yang bias (Hoerl dan Kennard, 1970).
Sifat dari estimator ridge pada umumnya memiliki varians yang minimum sehingga nilai VIF pada regresi ridge adalah diagonal utama dari matriks berikut: \(((Z^{t}Z+cI)^{-1}Z^{t}Z((Z^{t}Z+cI)^{-1})\)
library(glmnet)
x<-cbind(data$x1,data$x2,data$x3,data$x4,data$x5,data$x6,data$x7,data$x8,data$x9,data$x10)
y<-data$y
modelRidge <- glmnet::cv.glmnet(x=x,y=y,alpha=0)
summary(modelRidge)## Length Class Mode
## lambda 100 -none- numeric
## cvm 100 -none- numeric
## cvsd 100 -none- numeric
## cvup 100 -none- numeric
## cvlo 100 -none- numeric
## nzero 100 -none- numeric
## call 4 -none- call
## name 1 -none- character
## glmnet.fit 12 elnet list
## lambda.min 1 -none- numeric
## lambda.1se 1 -none- numeric
## index 2 -none- numericprint(modelRidge)##
## Call: glmnet::cv.glmnet(x = x, y = y, alpha = 0)
##
## Measure: Mean-Squared Error
##
## Lambda Index Measure SE Nonzero
## min 5.48 76 40.18 13.10 10
## 1se 51.11 52 52.78 27.29 10bestlambda <- modelRidge$lambda.min
cat("Lambda Terbaik:", bestlambda, "\n")## Lambda Terbaik: 5.480201coef(modelRidge, s = bestlambda)## 11 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 55.632327167
## V1 -0.066191046
## V2 -0.009023555
## V3 -0.064041001
## V4 -0.016042586
## V5 0.029230076
## V6 0.002151203
## V7 -0.077009815
## V8 -0.149501888
## V9 -0.105967169
## V10 -0.140464985Output di atas merupakan nilai intersep dan koefisien dari setiap variabel bebas terhadap variabel respon (stunting). Dapat terlihat bahwa variabel X8 merupakan variabel dengan koefisien negatif tertinggi yang berpengaruh terhadap model. Dengan demikian dapat disimpulkan bahwa variabel ini merupakan variabel yang paling berpengaruh terhadap variabel respon jika dibandingkan variabel bebas lainnya.
Y.est <- predict(modelRidge, s = bestlambda, newx = x)
r_squared_ridge <- 1 - sum((y - Y.est)^2) / sum((y - mean(y))^2)
r_squared_ridge## [1] 0.7994744D. REGRESI LASSO
Metode Least Absolute Shrinkage and Selection Operator (LASSO) diperkenalkan pertama kali oleh Tibshirani pada tahun 1996. Sesuai namanya regresi LASSO merupakan metode regresi berganda yang digunakan untuk shrinkage yaitu menyusutkan koefisien taksiran mendekati angka nol dan selection operator yaitu menyeleksi variabel-variabel bebas sehingga menghasilkan model dengan variabel terbaik. LASSO menyusutkan koefisien regresi dari variabel bebas yang memiliki korelasi tinggi dengan residual menjadi tepat pada nol atau mendekati nol (Tibshirani, 1996). Selain itu, regresi LASSO juga digunakan untuk data yang kontinu dan memerlukan variabel independen yang berdistribusi normal baku.
modelLasso <- glmnet::cv.glmnet(x=x,y=y,alpha=1)
print(modelLasso)##
## Call: glmnet::cv.glmnet(x = x, y = y, alpha = 1)
##
## Measure: Mean-Squared Error
##
## Lambda Index Measure SE Nonzero
## min 0.027 59 58.67 20.51 9
## 1se 5.876 1 72.15 37.96 0bestlambdalasso <- modelLasso$lambda.min
cat("Lambda terbaik:", bestlambdalasso, "\n")## Lambda terbaik: 0.02664804coef(modelLasso, s = bestlambdalasso)## 11 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 88.85195942
## V1 .
## V2 0.11213812
## V3 0.22996875
## V4 -0.24433642
## V5 -0.05912117
## V6 0.11103131
## V7 -0.53536945
## V8 -0.07810096
## V9 -0.12537495
## V10 -0.58672960Output di atas merupakan nilai intersep dan koefisien regresi dari masing-masing variabel bebas yang selanjutnya digunakan untuk membuat persamaan model regresi LASSO.
Y.estlasso <- predict(modelLasso, s = bestlambdalasso, newx = x)
r_squared_lasso <- 1 - sum((y - Y.estlasso)^2) / sum((y - mean(y))^2)
r_squared_lasso## [1] 0.9889791E. EVALUASI MODEL
Didasarkan pada perbandingan nilai R-Square, diperoleh hasil:
cat("Performa Model Regresi Klasik:", sum$adj.r.squared, "\n")## Performa Model Regresi Klasik: 0.9602593cat("Performa Model Regresi Ridge:", r_squared_ridge, "\n")## Performa Model Regresi Ridge: 0.7994744cat("Performa Model Regresi LASSO:", r_squared_lasso, "\n")## Performa Model Regresi LASSO: 0.9889791Berdasarkan hasil analisis yang dilakukan melalui variable selection, ridge regression dan LASSO regression, model terbaik yang dihasilkan adalah LASSO regression.
References
- Badan Pusat Statistik (BPS). 2021. Laporan Indeks Khusus Penanganan Stunting. Jakarta: Badan Pusat Stastitik.
- Draper, N.R., & Smith, H. 1998. Applied Regression Analysis. New York: John Wiley and sons.
- Gujarati, Damodar N. 2006. Ekonometrika Dasar. Jakarta: Penerbit Erlangga.
- Hoerl, A.E. & Kennard, R.W. 1970. Ridge Regression: Biased Estimation for Non-Orthogonal Problems. Technometrics, 12, 55-67.
- Montgomery, D. C., Peck, E. A., & Vining, G. G. 2006. Introduction to Linear Regression Analysis. 4th Ed. Canada: John Wiley & Sons.
- Tibshirani, R. 1996. Regression shrinkage and selection via the lasso. Journal of The Royal Statistical Society Series B Methodological, Vol.58, Nomor 1, p: 267-288.
Direktorat Statistik Kesejahteraan Rakyat, BPS, saptahas@bps.go.id