Input Data
#Data Peta
peta <- readOGR("D:/MY COLLEGE/SEMESTER 6/DC/DATA/JABAR_KAB.shp")
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## OGR data source with driver: ESRI Shapefile
## Source: "D:\MY COLLEGE\SEMESTER 6\DC\DATA\JABAR_KAB.shp", layer: "JABAR_KAB"
## with 27 features
## It has 6 fields
#Data Excel
df = read.xlsx("D:/MY COLLEGE/SEMESTER 6/DC/DATA/Jabar Data.xlsx")
View(df)
head(df)
## PROVNO KABKOTNO KODE2010 PROVINSI KABKOT IDSP2010 Y X1 X2 X3 X4
## 1 32 1 3201 JAWA BARAT BOGOR 3201 4 71 29.50 77938 20
## 2 32 2 3202 JAWA BARAT SUKABUMI 3202 2 54 33.82 222497 19
## 3 32 3 3203 JAWA BARAT CIANJUR 3203 1 88 37.45 785294 23
## 4 32 4 3204 JAWA BARAT BANDUNG 3204 3 70 32.59 722933 21
## 5 32 5 3205 JAWA BARAT GARUT 3205 1 50 36.94 861324 30
## 6 32 6 3206 JAWA BARAT TASIKMALAYA 3206 2 68 36.53 638521 22
## X5 X6 X7
## 1 13809 0.779 45339
## 2 17793 0.735 25616
## 3 15906 0.439 20002
## 4 15314 0.477 35590
## 5 18059 0.345 23356
## 6 17152 0.321 20855str(df)## 'data.frame': 27 obs. of 14 variables:
## $ PROVNO : num 32 32 32 32 32 32 32 32 32 32 ...
## $ KABKOTNO: num 1 2 3 4 5 6 7 8 9 10 ...
## $ KODE2010: num 3201 3202 3203 3204 3205 ...
## $ PROVINSI: chr "JAWA BARAT" "JAWA BARAT" "JAWA BARAT" "JAWA BARAT" ...
## $ KABKOT : chr "BOGOR" "SUKABUMI" "CIANJUR" "BANDUNG" ...
## $ IDSP2010: num 3201 3202 3203 3204 3205 ...
## $ Y : num 4 2 1 3 1 2 2 4 4 3 ...
## $ X1 : num 71 54 88 70 50 68 137 48 130 72 ...
## $ X2 : num 29.5 33.8 37.5 32.6 36.9 ...
## $ X3 : num 77938 222497 785294 722933 861324 ...
## $ X4 : num 20 19 23 21 30 22 30 18 30 24 ...
## $ X5 : num 13809 17793 15906 15314 18059 ...
## $ X6 : num 0.779 0.735 0.439 0.477 0.345 0.321 0.554 0.838 0.617 0.692 ...
## $ X7 : num 45339 25616 20002 35590 23356 ...dt.new <- df %>% dplyr::select(Y, X1, X2, X3, X4, X5, X6, X7)Data
Tabel 1 Peubah-Peubah yang Digunakan
Prosedur Analisis
Gambar 1 Prosedur Penelitian
1. Eksplorasi Data
Boxplot
par(mfrow=c(2,2))
boxplot(df$Y, xlab = "Y")
boxplot(df$X1, xlab = "X1")
boxplot(df$X2, xlab = "X2")
boxplot(df$X3, xlab = "X3")
boxplot(df$X4, xlab = "X4")
boxplot(df$X5, xlab = "X5")
boxplot(df$X6, xlab = "X6")
boxplot(df$X7, xlab = "X7")
Berdasarkan boxplot diatas, sebaran Y, X1, X3, X5 cenderung menjulur ke kanan. Artinya, terdapat sebagin besar kabupaten/kota yang memiliki nilai prevalensi pengidap DM, tingkat konsumsi minuman, jumlah pengidap hipertensi, dan rataan pembelian beras yang cenderung kecil. Sementara peubah X2, X6 dan X7 cenderung simetris. Hal ini berarti nilai persentase perokok, rataan konsumsi gula pasir, dan PDRB perkapita tiap kabupaten/kota cenderung sama. Peubah X4 cenderung menjulur ke kiri. Artinya, sebagian besar kabupaten/kota yang memiliki skor PPH cenderung besar. Selain itu terdeteksi adanya pencilan pada peubah Y dan X1, artinya terdapat nilai kedua peubah tersebut sangat jauh dari nilai rata-rata sebagian kabupaten/kota lainnya.
Histogram
pdc <- df[,c(-1,-2,-3,-4,-5,-6)]
plot_histogram(data = pdc, nrow = 2, ncol = 2, geom_histogram_args = list (fill="#D27685"))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Scatter Plot
plot(df$X1, df$Y,
xlab="Tingkat Konsumsi Minuman",
ylab="Prevalensi Pengidap DM",
pch=20, col="orange", cex=2)
reg.ols = lm(Y~X1, data = df)
lines.lm(reg.ols, col=2, add=T)
plot(df$X2, df$Y,
xlab="Persentase Penduduk Perokok",
ylab= "Prevalensi Pengidap DM",
pch=20, col="orange", cex=2)
reg.ols = lm(Y~X2, data = df)
lines.lm(reg.ols, col=2, add=T)
plot(df$X3, df$Y,
xlab="Jumlah Pengidap Hipertensi",
ylab="Prevalensi Pengidap DM",
pch=20, col="orange", cex=2)
reg.ols = lm(Y~X3, data = df)
lines.lm(reg.ols, col=2, add=T)
plot(df$X4, df$Y,
xlab="Skor PPH Konsumsi Sayur dan Buah",
ylab="Prevalensi Pengidap DM",
pch=20, col="orange", cex=2)
reg.ols = lm(Y~X4, data = df)
lines.lm(reg.ols, col=2, add=T)
plot(df$X5, df$Y,
xlab="Pengeluaran Pembelian Beras",
ylab="Prevalensi Pengidap DM",
pch=20, col="orange", cex=2)
reg.ols = lm(Y~X5, data = df)
lines.lm(reg.ols, col=2, add=T)
plot(df$X6, df$Y,
xlab="Konsumsi Gula Pasir",
ylab="Prevalensi Pengidap DM",
pch=20, col="orange", cex=2)
reg.ols = lm(Y~X6, data = df)
lines.lm(reg.ols, col=2, add=T)
plot(df$X7, df$Y,
xlab="PDRB per Kapita",
ylab="Prevalensi Pengidap DM",
pch=20, col="orange", cex=2)
reg.ols = lm(Y~X7, data = df)
lines.lm(reg.ols, col=2, add=T)
Matriks Korelasi
corrplot::corrplot(cor(df[,c(7:14)]), method = "number", type="upper", sig.level = 0.05)
kor <- cor(pdc)
ggcorrplot(kor, type="lower", lab = TRUE)
correlation heatmap dapat menunjukkan kekuatan hubungan antarpeubah. Dapat terlihat bahwa tingkat konsumsi minuman (X1), konsumsi gula pasir (X6),dan PDRB perkapita atas dasar harga berlaku (X7) memiliki korelasi positif terhadap prevalensi diabetes melitus (Y). Sementara itu, peubah lainnya yaitu persentase penduduk merokok (X2), jumlah orang yang mengidap hipertensi (X3), tingkat konsumsi sayur dan buah (X4), dan pengeluaran pembelian beras (X5) memiliki korelasi negatif terhadap prevalensi diabetes melitus (Y). Kemudian, terlihat bahwa persentase perokok dan rataan konsumsi gula pasir memiliki korelasi negatif terkuat terhadap peubah respon.
Peta Geospasial
Diabetes Melitus Tipe 1 dan 2
k=16
colfunc <- colorRampPalette(c("green", "yellow","red"))
color <- colfunc(k)
peta$y <- df$Y
spplot(peta, "y", col.regions=color, main="Prevalensi Orang yang Mengidap Penyakit\nDiabetes Melitus di Jawa Barat Tahun 2021")
Peta geospasial di samping merupakan peta geospasial prevalensi kasus diabetes melitus tipe 1 dan 2 di Provinsi Jawa Barat pada tahun 2021. Dapat dilihat bahwa mayoritas kabupaten/kota di Provinsi Jawa Barat memiliki angka kasus diabetes melitus relatif rendah. Umumnya kabupaten/kota dengan kasus diabetes melitus yang relatif tinggi berada pada Kota Depok, Kota Bekasi, Kota Bogor, Kota Bandung, dan Kota Cirebon.
Penduduk Perokok
k=16
colfunc <- colorRampPalette(c("green", "yellow","red"))
color <- colfunc(k)
peta$x2 <- df$X2
spplot(peta, "x2", col.regions=color, main="Persentase Perokok\nUsia 15-24 & >=55 Tahun dalam Sebulan Terakhir\ndi Jawa Barat Tahun 2021")
Konsumsi Sayur dan Buah-Buahan
k=16
colfunc <- colorRampPalette(c("green", "yellow","red"))
color <- colfunc(k)
peta$x3 <- df$X3
spplot(peta, "x3", col.regions=color, main="Jumlah Penderita Penyakit Hipertensi\ndi Jawa Barat Tahun 2021")
Rataan Konsumsi Gula Pasir
k=16
colfunc <- colorRampPalette(c("green", "yellow","red"))
color <- colfunc(k)
peta$x6 <- df$X6
spplot(peta, "x6", col.regions=color, main="Rata-Rata Pengeluaran Pembelian Beras Perkapita\nSeminggu di Jawa Barat Tahun 2021")
2. Pemodelan Penuh
RLB OLS Penuh
model.awal <- lm(Y~., data = dt.new)
summary(model.awal)##
## Call:
## lm(formula = Y ~ ., data = dt.new)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.93392 -0.91283 -0.00431 0.55453 2.86382
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.117e+01 4.862e+00 4.355 0.000341 ***
## X1 1.159e-02 1.651e-02 0.702 0.491081
## X2 -2.967e-01 1.097e-01 -2.704 0.014070 *
## X3 -2.199e-06 1.252e-06 -1.757 0.095085 .
## X4 -8.111e-03 9.197e-02 -0.088 0.930646
## X5 -4.447e-04 2.763e-04 -1.609 0.124005
## X6 -4.184e-01 2.252e+00 -0.186 0.854545
## X7 1.084e-05 1.332e-05 0.814 0.425737
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.534 on 19 degrees of freedom
## Multiple R-squared: 0.7196, Adjusted R-squared: 0.6164
## F-statistic: 6.967 on 7 and 19 DF, p-value: 0.0003475Berdasarkan persamaan regresi di atas, dapat diinterpretasikan bahwa ketika nilai peubah lain sama dengan nol maka prevalensi orang yang mengidap penyakit diabetes melitus di Jawa Barat akan bernilai sebesar 21.17 jiwa. Lalu, untuk setiap pertambahan tingkat konsumsi makanan sebesar 1 gram/kap/hari, prevalensi orang yang mengidap penyakit diabetes melitus akan meningkat sebesar 0.0116 jiwa dengan asumsi peubah lain tetap. Peubah lainnya berhubungan negatif sehingga peningkatan pada masing-masing persentase penduduk usia 15-24 &>= 55 tahun yang merokok dalam sebulan terakhir, prevalensi orang yang mengidap penyakit hipertensi, skor PPH tingkat konsumsi sayur dan buah-buahan, rata-rata pengeluaran pembelian beras perkapita seminggu, rata-rata konsumsi gula pasir perkapita seminggu, dan PDRB perkapita atas dasar harga berlaku akan menurunkan prevalensi orang yang mengidap penyakit diabetes melitus berurut sebesar 0.2967 jiwa, 0.0000022 jiwa, 0.0081 jiwa, 0.00045 jiwa, 0.4184 jiwa,dan 0.000011 jiwa dengan asumsi peubah lainnya tetap.
anova(model.awal)## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 0.316 0.316 0.1343 0.71802
## X2 1 91.034 91.034 38.7021 5.633e-06 ***
## X3 1 14.506 14.506 6.1670 0.02252 *
## X4 1 0.097 0.097 0.0412 0.84132
## X5 1 7.149 7.149 3.0395 0.09742 .
## X6 1 0.055 0.055 0.0235 0.87980
## X7 1 1.558 1.558 0.6626 0.42574
## Residuals 19 44.691 2.352
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13. Identifikasi Multikolinieritas
ols_vif_tol(model.awal)## Variables Tolerance VIF
## 1 X1 0.8142157 1.228176
## 2 X2 0.3906496 2.559838
## 3 X3 0.8071965 1.238856
## 4 X4 0.6792024 1.472315
## 5 X5 0.4739803 2.109792
## 6 X6 0.8142741 1.228088
## 7 X7 0.6651699 1.503375Berdasarkan model penuh, seluruh peubah yang digunakan tidak terdeteksi multikolinieritas karena memiliki nilai VIF < 10 dan tolerance > 0.1.
4. Uji Asumsi Klasik
Pada RLB OLS Penuh
uji_asumsi <- function(modelreg, autocol.test = c("runs", "dw", "bgodfrey"),
homos.test = c("bp", "glejser", "ncv"),
normal.test = c("ks", "shapiro.w", "jb"), taraf.nyata=0.05){
#Hipoteisis: H0: Asumsi terpenuhi vs H1: Asumsi tidak terpenuhi
sisaan <- modelreg$residuals
uji_t <- c()
autokol <- c()
homoskedastisitas <- c()
ks <- c()
p_value <- matrix(NA, ncol = 1, nrow = 4)
a <- t.test(sisaan,
mu = 0,
conf.level = 0.95)
p_value[1,1] <- round(a$p.value,3)
if(a$p.value<taraf.nyata){hasil = "Tolak H0"}else{hasil="Tak Tolak H0"}
uji_t <- hasil
if (autocol.test == "runs"){
b <- randtests::runs.test(sisaan)
p_value[2,1] <- round(b$p.value,3)
if(b$p.value<taraf.nyata){hasil = "Tolak H0"}else{hasil="Tak Tolak H0"}
autokol <- hasil
}else if (autocol.test == "dw") {
b <- lmtest::dwtest(modelreg)
p_value[2,1] <- round(b$p.value,3)
if(b$p.value<taraf.nyata){hasil = "Tolak H0"}else{hasil="Tak Tolak H0"}
autokol <- hasil
}else if (autocol.test == "bgodfrey") {
b <- lmtest::bgtest(modelreg)
p_value[2,1] <- round(b$p.value,3)
if(b$p.value<taraf.nyata){hasil = "Tolak H0"}else{hasil="Tak Tolak H0"}
autokol <- hasil
}
if (homos.test == "bp") {
c <- lmtest::bptest(modelreg)
p_value[3,1] <- round(c$p.value,3)
if(c$p.value<taraf.nyata){hasil = "Tolak H0"}else{hasil="Tak Tolak H0"}
homoskedastisitas <- hasil
} else if (homos.test == "glejser") {
c <- skedastic::glejser(modelreg)
p_value[3,1] <- round(c$p.value,3)
if(c$p.value<taraf.nyata){hasil = "Tolak H0"}else{hasil="Tak Tolak H0"}
homoskedastisitas <- hasil
} else if (homos.test == "ncv") {
c <- car::ncvTest(modelreg)
p_value[3,1] <- round(c$p,3)
if(c$p<taraf.nyata){hasil = "Tolak H0"}else{hasil="Tak Tolak H0"}
homoskedastisitas <- hasil
}
if (normal.test == "ks"){
d <- ks.test(sisaan, "pnorm")
p_value[4,1] <- round(d$p.value,3)
if(d$p.value<taraf.nyata){hasil = "Tolak H0"}else{hasil="Tak Tolak H0"}
ks <- hasil
} else if (normal.test == "shapiro.w") {
d <- shapiro.test(sisaan)
p_value[4,1] <- round(d$p.value,3)
if(d$p.value<taraf.nyata){hasil = "Tolak H0"}else{hasil="Tak Tolak H0"}
ks <- hasil
} else if (normal.test == "jb") {
d <- tseries::jarque.bera.test(sisaan)
p_value[4,1] <- round(d$p.value,3)
if(d$p.value<taraf.nyata){hasil = "Tolak H0"}else{hasil="Tak Tolak H0"}
ks <- hasil
}
keputusan <- rbind(uji_t, autokol, homoskedastisitas, ks)
tabel <- data.frame(p_value, keputusan)
colnames(tabel) <- c("P-value", "Keputusan")
rownames(tabel) <- c("E(sisaan) = 0", "Non-autokorelasi", "Homoskedastisitas",
"Normality")
duga <- plot(modelreg,1)
kuantil <- plot(modelreg,2)
histo <- hist(sisaan)
urutan <- plot(x = 1:length(sisaan),
y = sisaan,
type = 'b',
ylab = "Residuals",
xlab = "Observation")
print(tabel)
}
par(mfrow=c(2,2))
uji_asumsi(model.awal, "runs", "bp", "ks")
## P-value Keputusan
## E(sisaan) = 0 1.000 Tak Tolak H0
## Non-autokorelasi 0.423 Tak Tolak H0
## Homoskedastisitas 0.202 Tak Tolak H0
## Normality 0.858 Tak Tolak H0Berdasarkan uji formal dan secara eksploratif, dapat disimpulkan bahwa semua asumsi klasik dalam model regresi terpenuhi pada taraf nyata 5%. Hal ini dapat dilihat dari semua nilai p-value dalam uji formal yang lebih dari alpha (tak tolak H0). Lalu, plot sisaan vs Y duga, QQ-plot, histogram sisaan, dan plot sisaan vs urutan juga menunjukkan bahwa semua asumsi klasik terpenuhi.
5. Pemilihan Peubah Penjelas
Menggunakan stepwise elimination.
stepw <- step(lm(Y~1, data=dt.new), direction="both", scope=formula(dt.new), trace=1)## Start: AIC=49.94
## Y ~ 1
##
## Df Sum of Sq RSS AIC
## + X2 1 90.364 69.043 29.350
## + X5 1 83.353 76.055 31.962
## + X7 1 40.095 119.312 44.119
## + X6 1 18.042 141.365 48.699
## + X4 1 13.229 146.179 49.603
## <none> 159.407 49.942
## + X3 1 9.603 149.805 50.264
## + X1 1 0.316 159.091 51.888
##
## Step: AIC=29.35
## Y ~ X2
##
## Df Sum of Sq RSS AIC
## + X5 1 15.860 53.183 24.304
## + X3 1 12.715 56.329 25.855
## <none> 69.043 29.350
## + X7 1 2.042 67.002 30.540
## + X6 1 1.078 67.965 30.925
## + X1 1 0.986 68.057 30.962
## + X4 1 0.086 68.958 31.317
## - X2 1 90.364 159.407 49.942
##
## Step: AIC=24.3
## Y ~ X2 + X5
##
## Df Sum of Sq RSS AIC
## + X3 1 5.7070 47.476 23.239
## <none> 53.183 24.304
## + X7 1 0.8497 52.334 25.869
## + X1 1 0.2257 52.958 26.189
## + X4 1 0.0023 53.181 26.302
## + X6 1 0.0017 53.182 26.303
## - X5 1 15.8600 69.043 29.350
## - X2 1 22.8712 76.055 31.962
##
## Step: AIC=23.24
## Y ~ X2 + X5 + X3
##
## Df Sum of Sq RSS AIC
## <none> 47.476 23.239
## - X3 1 5.7070 53.183 24.304
## + X7 1 1.4934 45.983 24.376
## + X1 1 1.1635 46.313 24.569
## + X4 1 0.2004 47.276 25.124
## + X6 1 0.0480 47.428 25.211
## - X5 1 8.8522 56.329 25.855
## - X2 1 27.3690 74.845 33.529stepAIC(model.awal, direction = "both")## Start: AIC=29.61
## Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7
##
## Df Sum of Sq RSS AIC
## - X4 1 0.0183 44.710 27.618
## - X6 1 0.0812 44.773 27.655
## - X1 1 1.1598 45.851 28.298
## - X7 1 1.5585 46.250 28.532
## <none> 44.691 29.606
## - X5 1 6.0930 50.784 31.057
## - X3 1 7.2583 51.950 31.670
## - X2 1 17.1972 61.889 36.396
##
## Step: AIC=27.62
## Y ~ X1 + X2 + X3 + X5 + X6 + X7
##
## Df Sum of Sq RSS AIC
## - X6 1 0.0789 44.789 25.665
## - X1 1 1.2048 45.914 26.335
## - X7 1 1.5469 46.257 26.536
## <none> 44.710 27.618
## - X5 1 6.0749 50.785 29.057
## + X4 1 0.0183 44.691 29.606
## - X3 1 7.4132 52.123 29.760
## - X2 1 21.0768 65.786 36.046
##
## Step: AIC=25.67
## Y ~ X1 + X2 + X3 + X5 + X7
##
## Df Sum of Sq RSS AIC
## - X1 1 1.1945 45.983 24.376
## - X7 1 1.5243 46.313 24.569
## <none> 44.789 25.665
## - X5 1 6.0331 50.822 27.077
## + X6 1 0.0789 44.710 27.618
## + X4 1 0.0159 44.773 27.655
## - X3 1 7.3350 52.123 27.760
## - X2 1 21.0966 65.885 34.086
##
## Step: AIC=24.38
## Y ~ X2 + X3 + X5 + X7
##
## Df Sum of Sq RSS AIC
## - X7 1 1.4934 47.476 23.239
## <none> 45.983 24.376
## + X1 1 1.1945 44.789 25.665
## - X3 1 6.3508 52.334 25.869
## + X6 1 0.0686 45.914 26.335
## + X4 1 0.0665 45.917 26.337
## - X5 1 7.4941 53.477 26.452
## - X2 1 20.0212 66.004 32.135
##
## Step: AIC=23.24
## Y ~ X2 + X3 + X5
##
## Df Sum of Sq RSS AIC
## <none> 47.476 23.239
## - X3 1 5.7070 53.183 24.304
## + X7 1 1.4934 45.983 24.376
## + X1 1 1.1635 46.313 24.569
## + X4 1 0.2004 47.276 25.124
## + X6 1 0.0480 47.428 25.211
## - X5 1 8.8522 56.329 25.855
## - X2 1 27.3690 74.845 33.529##
## Call:
## lm(formula = Y ~ X2 + X3 + X5, data = dt.new)
##
## Coefficients:
## (Intercept) X2 X3 X5
## 2.348e+01 -3.139e-01 -1.857e-06 -5.073e-04
stepAIC(model.awal, direction = "backward")## Start: AIC=29.61
## Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7
##
## Df Sum of Sq RSS AIC
## - X4 1 0.0183 44.710 27.618
## - X6 1 0.0812 44.773 27.655
## - X1 1 1.1598 45.851 28.298
## - X7 1 1.5585 46.250 28.532
## <none> 44.691 29.606
## - X5 1 6.0930 50.784 31.057
## - X3 1 7.2583 51.950 31.670
## - X2 1 17.1972 61.889 36.396
##
## Step: AIC=27.62
## Y ~ X1 + X2 + X3 + X5 + X6 + X7
##
## Df Sum of Sq RSS AIC
## - X6 1 0.0789 44.789 25.665
## - X1 1 1.2048 45.914 26.335
## - X7 1 1.5469 46.257 26.536
## <none> 44.710 27.618
## - X5 1 6.0749 50.785 29.057
## - X3 1 7.4132 52.123 29.760
## - X2 1 21.0768 65.786 36.046
##
## Step: AIC=25.67
## Y ~ X1 + X2 + X3 + X5 + X7
##
## Df Sum of Sq RSS AIC
## - X1 1 1.1945 45.983 24.376
## - X7 1 1.5243 46.313 24.569
## <none> 44.789 25.665
## - X5 1 6.0331 50.822 27.077
## - X3 1 7.3350 52.123 27.760
## - X2 1 21.0966 65.885 34.086
##
## Step: AIC=24.38
## Y ~ X2 + X3 + X5 + X7
##
## Df Sum of Sq RSS AIC
## - X7 1 1.4934 47.476 23.239
## <none> 45.983 24.376
## - X3 1 6.3508 52.334 25.869
## - X5 1 7.4941 53.477 26.452
## - X2 1 20.0212 66.004 32.135
##
## Step: AIC=23.24
## Y ~ X2 + X3 + X5
##
## Df Sum of Sq RSS AIC
## <none> 47.476 23.239
## - X3 1 5.7070 53.183 24.304
## - X5 1 8.8522 56.329 25.855
## - X2 1 27.3690 74.845 33.529##
## Call:
## lm(formula = Y ~ X2 + X3 + X5, data = dt.new)
##
## Coefficients:
## (Intercept) X2 X3 X5
## 2.348e+01 -3.139e-01 -1.857e-06 -5.073e-04stepAIC(model.awal, direction = "forward")## Start: AIC=29.61
## Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7, data = dt.new)
##
## Coefficients:
## (Intercept) X1 X2 X3 X4 X5
## 2.117e+01 1.159e-02 -2.967e-01 -2.199e-06 -8.111e-03 -4.447e-04
## X6 X7
## -4.184e-01 1.084e-05model.bestsubsets <- regsubsets(Y~., dt.new)
summary(model.bestsubsets)## Subset selection object
## Call: regsubsets.formula(Y ~ ., dt.new)
## 7 Variables (and intercept)
## Forced in Forced out
## X1 FALSE FALSE
## X2 FALSE FALSE
## X3 FALSE FALSE
## X4 FALSE FALSE
## X5 FALSE FALSE
## X6 FALSE FALSE
## X7 FALSE FALSE
## 1 subsets of each size up to 7
## Selection Algorithm: exhaustive
## X1 X2 X3 X4 X5 X6 X7
## 1 ( 1 ) " " "*" " " " " " " " " " "
## 2 ( 1 ) " " "*" " " " " "*" " " " "
## 3 ( 1 ) " " "*" "*" " " "*" " " " "
## 4 ( 1 ) " " "*" "*" " " "*" " " "*"
## 5 ( 1 ) "*" "*" "*" " " "*" " " "*"
## 6 ( 1 ) "*" "*" "*" " " "*" "*" "*"
## 7 ( 1 ) "*" "*" "*" "*" "*" "*" "*"6. Pendeteksian Pencilan, Titik Leverage, dam Amatan Berpengaruh
Pencilan
#Mendeteksi Pencilan
which(abs(rstandard(model.awal))>2)## 1 20 22
## 1 20 22df[which(abs(rstandard(model.awal))>2),]## PROVNO KABKOTNO KODE2010 PROVINSI KABKOT IDSP2010 Y X1 X2 X3 X4
## 1 32 1 3201 JAWA BARAT BOGOR 3201 4 71 29.5 77938 20
## 20 32 72 3272 JAWA BARAT SUKABUMI 3272 3 66 32.8 77257 30
## 22 32 74 3274 JAWA BARAT CIREBON 3274 10 63 31.0 77723 21
## X5 X6 X7
## 1 13809 0.779 45339
## 20 13829 0.558 37209
## 22 12874 0.696 72749Diketahui bahwa terdapat tiga amatan pencilan, yaitu pada data ke-1 (Kab Bogor), data ke-20 (Kota Sukabumi), dan data ke-22 (Kota Cirebon).
Titik Leverage
#Mendeteksi Leverage
which(abs(hatvalues(model.awal))>0.52)## named integer(0)Diketahui bahwa tidak ada amatan data yang merupakan titik leverage.
df[which(abs(hatvalues(model.awal))>0.52),] #2p/n=2x7/27## [1] PROVNO KABKOTNO KODE2010 PROVINSI KABKOT IDSP2010 Y X1
## [9] X2 X3 X4 X5 X6 X7
## <0 rows> (or 0-length row.names)#Batas pencilan dan leverage
ri = rstandard(model.awal)
hii = hatvalues(model.awal)
Obs = c(1:27)
summ <- cbind.data.frame(Obs, hii, ri)
View(summ)Amatan Berpengaruh
#Eksploratif amatan berpengaruh (jarak cook)
plot(model.awal, 5)
Dengan metode jarak Cook (Cook’s distance) secara eksploratif, terlihat bahwa tidak ada amatan berpengaruh pada taraf nyata 5% karena tidak ada titik plot yang berada di luar garis putus-putus skala 0.5.
#mendeteksi amatan berpengaruh
influence.measures(model.awal)## Influence measures of
## lm(formula = Y ~ ., data = dt.new) :
##
## dfb.1_ dfb.X1 dfb.X2 dfb.X3 dfb.X4 dfb.X5 dfb.X6
## 1 -0.364948 -1.38e-01 0.214111 0.495444 0.383129 4.67e-02 -0.34748
## 2 0.239548 3.24e-02 0.193821 0.235025 0.259874 -5.41e-01 -0.32561
## 3 -0.143860 -1.25e-01 -0.239652 -0.326736 0.276653 2.21e-01 0.16664
## 4 -0.180672 2.11e-02 -0.025253 -0.191729 0.153679 1.02e-01 0.12923
## 5 0.009162 -1.12e-01 -0.046793 0.076385 0.093606 3.13e-02 -0.06033
## 6 -0.000574 2.51e-05 -0.000121 -0.000486 0.000772 -4.35e-05 0.00113
## 7 0.027757 -5.24e-02 -0.018339 0.010015 0.009046 -5.87e-03 -0.00214
## 8 0.027582 -1.07e-01 0.252428 0.067749 -0.228408 -1.56e-01 0.20448
## 9 0.001106 -7.12e-02 0.008754 -0.019473 -0.014488 1.64e-02 -0.00474
## 10 0.007177 -2.55e-02 0.054903 0.106990 -0.030320 -5.28e-02 0.05299
## 11 -0.577348 6.61e-02 0.063003 -0.363538 0.296961 4.11e-01 0.19040
## 12 0.138191 -1.51e-02 0.009919 -0.006613 -0.052163 -8.97e-02 -0.13969
## 13 -0.007113 -7.36e-02 -0.002469 0.031736 0.078545 -1.70e-02 0.03531
## 14 -0.096717 -4.79e-02 -0.202603 0.048574 -0.015437 2.86e-01 0.08261
## 15 0.103004 -1.33e-02 -0.118310 -0.082675 0.026714 1.26e-02 0.01766
## 16 0.224628 1.57e-01 0.234184 -0.108852 -0.151484 -2.99e-01 -0.21504
## 17 0.005594 -7.37e-02 0.022956 0.014206 0.057640 -2.29e-02 -0.01880
## 18 0.038442 5.22e-03 -0.032774 0.036002 0.010594 -2.63e-02 0.00195
## 19 0.249947 8.60e-02 -0.274703 -0.368639 -0.048504 1.00e-01 -0.23980
## 20 -0.711479 6.16e-01 0.447174 0.625793 -1.146645 5.40e-01 0.43554
## 21 -0.071253 1.26e-02 0.123112 0.165418 -0.063988 -3.41e-02 -0.13670
## 22 0.510933 -1.71e-01 0.450613 -0.353028 -0.419559 -6.63e-01 -0.04750
## 23 0.282595 3.57e-02 -0.384167 0.285455 -0.081453 8.91e-03 0.23396
## 24 0.231700 9.60e-02 -0.457758 -0.006833 0.209516 5.80e-02 -0.02030
## 25 0.064178 4.50e-02 -0.088088 -0.086677 0.000852 1.24e-02 -0.01625
## 26 -0.001488 -1.07e-02 0.013697 0.022874 0.012481 -2.55e-02 0.02679
## 27 0.045008 -8.38e-02 0.043302 -0.057224 0.097860 -9.02e-02 -0.03097
## dfb.X7 dffit cov.r cook.d hat inf
## 1 0.206537 -1.17380 0.215 1.38e-01 0.196
## 2 0.058975 -0.78173 1.660 7.65e-02 0.388
## 3 0.093726 -0.54221 1.555 3.73e-02 0.283
## 4 0.067213 -0.32595 1.573 1.37e-02 0.193
## 5 -0.037669 0.19823 2.668 5.17e-03 0.435 *
## 6 0.000232 -0.00182 2.055 4.39e-07 0.250
## 7 -0.007019 -0.06315 3.097 5.26e-04 0.503 *
## 8 -0.040204 0.43428 2.706 2.46e-02 0.480 *
## 9 0.036453 -0.10462 2.442 1.44e-03 0.374 *
## 10 -0.038982 0.13768 1.926 2.49e-03 0.223
## 11 0.142923 0.83635 0.600 7.98e-02 0.199
## 12 -0.014311 -0.17609 2.640 4.08e-03 0.427 *
## 13 -0.038963 0.11804 1.914 1.83e-03 0.213
## 14 -0.153777 -0.42614 1.570 2.32e-02 0.238
## 15 -0.314340 -0.37246 1.678 1.79e-02 0.242
## 16 -0.453500 -0.81174 1.426 8.15e-02 0.355
## 17 -0.026100 0.13237 1.623 2.29e-03 0.109
## 18 -0.024098 -0.08288 1.998 9.05e-04 0.236
## 19 -0.051138 0.54459 1.238 3.69e-02 0.216
## 20 0.465643 -1.68130 0.207 2.77e-01 0.312 *
## 21 0.560497 0.68569 1.946 5.98e-02 0.410
## 22 0.452469 1.33728 0.223 1.79e-01 0.238
## 23 -0.374972 0.77572 1.131 7.32e-02 0.285
## 24 -0.316362 0.57100 1.739 4.15e-02 0.333
## 25 -0.005855 0.16744 1.723 3.67e-03 0.162
## 26 -0.004548 -0.04713 2.602 2.93e-04 0.409 *
## 27 -0.039870 0.19820 2.084 5.16e-03 0.291summary(influence.measures(model.awal))## Potentially influential observations of
## lm(formula = Y ~ ., data = dt.new) :
##
## dfb.1_ dfb.X1 dfb.X2 dfb.X3 dfb.X4 dfb.X5 dfb.X6 dfb.X7 dffit cov.r
## 5 0.01 -0.11 -0.05 0.08 0.09 0.03 -0.06 -0.04 0.20 2.67_*
## 7 0.03 -0.05 -0.02 0.01 0.01 -0.01 0.00 -0.01 -0.06 3.10_*
## 8 0.03 -0.11 0.25 0.07 -0.23 -0.16 0.20 -0.04 0.43 2.71_*
## 9 0.00 -0.07 0.01 -0.02 -0.01 0.02 0.00 0.04 -0.10 2.44_*
## 12 0.14 -0.02 0.01 -0.01 -0.05 -0.09 -0.14 -0.01 -0.18 2.64_*
## 20 -0.71 0.62 0.45 0.63 -1.15_* 0.54 0.44 0.47 -1.68 0.21
## 26 0.00 -0.01 0.01 0.02 0.01 -0.03 0.03 0.00 -0.05 2.60_*
## cook.d hat
## 5 0.01 0.44
## 7 0.00 0.50
## 8 0.02 0.48
## 9 0.00 0.37
## 12 0.00 0.43
## 20 0.28 0.31
## 26 0.00 0.41Berdasarkan ringkasan di atas dengan pendekatan Jarak Cook, tidak ada amatan berpengaruh.
dt.new[c(1,20,22),1:8]## Y X1 X2 X3 X4 X5 X6 X7
## 1 4 71 29.5 77938 20 13809 0.779 45339
## 20 3 66 32.8 77257 30 13829 0.558 37209
## 22 10 63 31.0 77723 21 12874 0.696 72749D1 <- cooks.distance(model.awal)
r <- stdres(model.awal)
a <- cbind(dt.new, D1, r)
a[D1 > 4/51, ] # D1 = jarak pada GLS yang diaproksimasikan kok bisa dapet 4/51## Y X1 X2 X3 X4 X5 X6 X7 D1 r
## 1 4 71 29.50 77938 20 13809 0.779 45339 0.13826500 -2.132915
## 11 4 81 39.82 248174 30 17566 0.593 32130 0.07978193 1.604910
## 16 5 66 27.18 658978 25 15143 0.750 107788 0.08151170 -1.089310
## 20 3 66 32.80 77257 30 13829 0.558 37209 0.27701774 -2.210954
## 22 10 63 31.00 77723 21 12874 0.696 72749 0.17910375 2.139637r_absolut <- abs(r)
a <- cbind(dt.new, D1, r, r_absolut)
sort_a <- a[order(-r_absolut), ]
sort_a[1:27, ]## Y X1 X2 X3 X4 X5 X6 X7 D1 r
## 20 3 66 32.80 77257 30 13829 0.558 37209 2.770177e-01 -2.210953862
## 22 10 63 31.00 77723 21 12874 0.696 72749 1.791038e-01 2.139636839
## 1 4 71 29.50 77938 20 13809 0.779 45339 1.382650e-01 -2.132915048
## 11 4 81 39.82 248174 30 17566 0.593 32130 7.978193e-02 1.604909625
## 23 8 78 26.32 667811 21 14186 0.770 39527 7.324726e-02 1.213061647
## 16 5 66 27.18 658978 25 15143 0.750 107788 8.151170e-02 -1.089310291
## 19 8 73 29.81 55386 22 14741 0.526 45921 3.693005e-02 1.034064940
## 2 2 54 33.82 222497 19 17793 0.735 25616 7.653676e-02 -0.982332987
## 3 1 88 37.45 785294 23 15906 0.439 20002 3.725054e-02 -0.868554790
## 21 7 79 30.00 722933 25 14507 0.558 121126 5.979113e-02 0.829292029
## 24 8 92 26.08 513142 26 14157 0.653 35659 4.151516e-02 0.814950337
## 14 5 87 34.09 231691 27 13066 0.610 69976 2.321194e-02 -0.770234240
## 4 3 70 32.59 722933 21 15314 0.477 35590 1.368113e-02 -0.675726456
## 15 4 80 32.84 623205 26 15064 0.601 98726 1.787164e-02 -0.669991367
## 8 4 48 37.08 222360 18 14617 0.838 22805 2.460501e-02 0.462024993
## 25 8 82 28.03 175156 23 13838 0.666 59906 3.669519e-03 0.390387398
## 17 3 63 37.62 423891 28 16226 0.522 26879 2.293570e-03 0.387379668
## 27 4 62 38.40 61015 30 14728 0.546 22892 5.155837e-03 0.316976030
## 10 3 72 35.75 768968 24 15432 0.692 25930 2.491795e-03 0.263534375
## 13 3 60 36.34 512669 30 15838 0.655 26292 1.833268e-03 0.232774231
## 5 1 50 36.94 861324 30 18059 0.345 23356 5.170323e-03 0.231637474
## 12 2 88 36.60 550316 30 17417 0.883 44072 4.081832e-03 -0.209342201
## 18 1 65 41.65 147446 26 17817 0.508 28366 9.051440e-04 -0.153062808
## 9 4 130 34.41 648010 30 14892 0.617 22833 1.442642e-03 -0.139016931
## 7 2 137 41.07 398281 30 16752 0.554 27218 5.260376e-04 -0.064460511
## 26 3 73 35.39 215761 22 18021 0.299 31556 2.930294e-04 -0.058259114
## 6 2 68 36.53 638521 22 17152 0.321 20855 4.394392e-07 -0.003246992
## r_absolut
## 20 2.210953862
## 22 2.139636839
## 1 2.132915048
## 11 1.604909625
## 23 1.213061647
## 16 1.089310291
## 19 1.034064940
## 2 0.982332987
## 3 0.868554790
## 21 0.829292029
## 24 0.814950337
## 14 0.770234240
## 4 0.675726456
## 15 0.669991367
## 8 0.462024993
## 25 0.390387398
## 17 0.387379668
## 27 0.316976030
## 10 0.263534375
## 13 0.232774231
## 5 0.231637474
## 12 0.209342201
## 18 0.153062808
## 9 0.139016931
## 7 0.064460511
## 26 0.058259114
## 6 0.0032469927. Kelayakan Model
Uji F Simultan
Model di atas merupakan model OLS setelah dilakukan metode stepwise elimination. Terlihat bahwa model tersebut memiliki p-value sebesar 3.004e-06 mendekati 0 < 0.05 maka tolak H0. Artinya, minimal ada satu peubah penjelas yang berpengaruh linier terhadap peubah respon sehingga model layak untuk digunakan.
Uji T Parsial
Terlihat bahwa peubah X2 dan X5 signifikan terhadap taraf nyata 5%. Artinya, peubah X2 dan X5 secara individu berpengaruh nyata terhadap tingkat prevalensi diabetes melitus di Provinsi Jawa Barat. Sementara itu, peubah X3 memiliki p-value > 0.05 sehingga dikatakan bahwa peubah X3 tidak signifikan pada taraf nyata 5% dan selanjutnya dikeluarkan dari model.
8. Pemodelan Regresi Robust Estimasi-M
Pemodelan Reduksi
OLS Berdasarkan Stepwise
model.ols.stepw <- lm(Y~X2+X3+X5, data = dt.new)
summary(model.ols.stepw)##
## Call:
## lm(formula = Y ~ X2 + X3 + X5, data = dt.new)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.06930 -0.54561 -0.01149 0.65769 2.92688
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.348e+01 2.832e+00 8.292 2.31e-08 ***
## X2 -3.139e-01 8.622e-02 -3.641 0.00136 **
## X3 -1.857e-06 1.117e-06 -1.663 0.10993
## X5 -5.073e-04 2.450e-04 -2.071 0.04978 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.437 on 23 degrees of freedom
## Multiple R-squared: 0.7022, Adjusted R-squared: 0.6633
## F-statistic: 18.08 on 3 and 23 DF, p-value: 3.004e-06Berdasarkan Stepwise & Signifikansi Parameter
model.ols.fix <- lm(Y~X2+X5, data = dt.new)
summary(model.ols.fix)##
## Call:
## lm(formula = Y ~ X2 + X5, data = dt.new)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5701 -0.7340 -0.2142 0.8059 3.3177
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 23.5447298 2.9339223 8.025 2.99e-08 ***
## X2 -0.2776334 0.0864192 -3.213 0.00372 **
## X5 -0.0006413 0.0002397 -2.675 0.01323 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.489 on 24 degrees of freedom
## Multiple R-squared: 0.6664, Adjusted R-squared: 0.6386
## F-statistic: 23.97 on 2 and 24 DF, p-value: 1.902e-06Regresi Robust: Pembobot Huber
model.huber.fix <- rlm(Y~X2+X5, data=dt.new, psi=psi.huber)
summary(model.huber.fix)##
## Call: rlm(formula = Y ~ X2 + X5, data = dt.new, psi = psi.huber)
## Residuals:
## Min 1Q Median 3Q Max
## -2.5236 -0.6820 -0.1490 0.8004 3.3519
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 23.7399 2.8219 8.4127
## X2 -0.2870 0.0831 -3.4526
## X5 -0.0006 0.0002 -2.7611
##
## Residual standard error: 1.076 on 24 degrees of freedomRegresi Robust: Pembobot Tukey Bisquare
model.tukbis.fix <- rlm(Y~X2+X5, data=dt.new, psi = psi.bisquare)
summary(model.tukbis.fix)##
## Call: rlm(formula = Y ~ X2 + X5, data = dt.new, psi = psi.bisquare)
## Residuals:
## Min 1Q Median 3Q Max
## -2.4403 -0.6917 -0.1811 0.8178 3.4632
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 23.2737 3.1227 7.4530
## X2 -0.2904 0.0920 -3.1572
## X5 -0.0006 0.0003 -2.3548
##
## Residual standard error: 1.13 on 24 degrees of freedomTerlihat bahwa nilai intersep dan kemiringan pada ketiga model regresi tersebut tidak berbeda signifikan dan memiliki tanda masing-masing parameter yang sama.
Robust: Pembobot Huber
Tidak usah dimasukin PPT
model.huber.awal <- rlm(Y~., data=dt.new, psi=psi.huber)
summary(model.huber.awal)##
## Call: rlm(formula = Y ~ ., data = dt.new, psi = psi.huber)
## Residuals:
## Min 1Q Median 3Q Max
## -3.5079 -0.5906 0.0360 0.4532 2.7635
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 23.9518 4.6254 5.1783
## X1 0.0076 0.0157 0.4821
## X2 -0.3670 0.1044 -3.5154
## X3 0.0000 0.0000 -2.3366
## X4 0.0345 0.0875 0.3946
## X5 -0.0005 0.0003 -1.7910
## X6 -0.4053 2.1421 -0.1892
## X7 0.0000 0.0000 0.0201
##
## Residual standard error: 0.701 on 19 degrees of freedomRobust: Pembobot Tukey Bisquare
Tidak usah dimasukin PPT
model.tukbis.awal <- rlm(Y~., data=dt.new, psi = psi.bisquare)## Warning in rlm.default(x, y, weights, method = method, wt.method = wt.method, :
## 'rlm' failed to converge in 20 stepssummary(model.tukbis.awal)##
## Call: rlm(formula = Y ~ ., data = dt.new, psi = psi.bisquare)
## Residuals:
## Min 1Q Median 3Q Max
## -3.92766 -0.49885 -0.04701 0.39539 2.92536
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 25.4379 3.3529 7.5868
## X1 0.0058 0.0114 0.5080
## X2 -0.4201 0.0757 -5.5511
## X3 0.0000 0.0000 -3.8198
## X4 0.0618 0.0634 0.9742
## X5 -0.0004 0.0002 -2.3613
## X6 -0.0124 1.5528 -0.0080
## X7 0.0000 0.0000 -1.4158
##
## Residual standard error: 0.6634 on 19 degrees of freedomBerdasarkan Matriks Korelasi
Tidak usah dimasukin PPT
model4 <- lm(Y~X2+X3+X4+X5+X6+X7, data = dt.new)
summary(model4)##
## Call:
## lm(formula = Y ~ X2 + X3 + X4 + X5 + X6 + X7, data = dt.new)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.03325 -0.91865 0.08759 0.61461 2.79651
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.191e+01 4.686e+00 4.676 0.000145 ***
## X2 -2.965e-01 1.083e-01 -2.737 0.012701 *
## X3 -2.033e-06 1.214e-06 -1.676 0.109399
## X4 1.416e-02 8.523e-02 0.166 0.869698
## X5 -4.765e-04 2.691e-04 -1.771 0.091797 .
## X6 -3.754e-01 2.222e+00 -0.169 0.867521
## X7 1.018e-05 1.311e-05 0.776 0.446853
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.514 on 20 degrees of freedom
## Multiple R-squared: 0.7124, Adjusted R-squared: 0.6261
## F-statistic: 8.255 on 6 and 20 DF, p-value: 0.0001397Korelasi >= 0.5
Tidak usah dimasukin PPT
model5 <- lm(Y~X2+X5+X7, data = dt.new)
summary(model4)##
## Call:
## lm(formula = Y ~ X2 + X3 + X4 + X5 + X6 + X7, data = dt.new)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.03325 -0.91865 0.08759 0.61461 2.79651
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.191e+01 4.686e+00 4.676 0.000145 ***
## X2 -2.965e-01 1.083e-01 -2.737 0.012701 *
## X3 -2.033e-06 1.214e-06 -1.676 0.109399
## X4 1.416e-02 8.523e-02 0.166 0.869698
## X5 -4.765e-04 2.691e-04 -1.771 0.091797 .
## X6 -3.754e-01 2.222e+00 -0.169 0.867521
## X7 1.018e-05 1.311e-05 0.776 0.446853
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.514 on 20 degrees of freedom
## Multiple R-squared: 0.7124, Adjusted R-squared: 0.6261
## F-statistic: 8.255 on 6 and 20 DF, p-value: 0.00013979. Ukuran Kebaikan Model
Berdasarkan nilai Galat Baku dan AIC.
Nilai AIC
AIC(model.ols.fix)## [1] 102.9263AIC(model.huber.fix)## [1] 102.9786AIC(model.tukbis.fix)## [1] 103.0437Galat Baku
Model OLS stepwise : 1.489
Model Regresi Robust: Pembobot Huber Stepwise : 1.076
Model Regresi Robust: Pembobot Tukey Bisquare Stepwise : 1.13
Dalam penentuan kebaikan model, dilakukan perbandingan galat baku sisaan dan nilai AIC setiap model yang diterapkan. Semakin kecil nilai galat baku sisaan dan AIC, maka model semakin baik sehingga nilai prediksi Y akan sedekat mungkin dengan aktualnya. Berdasarkan tabel di atas, didapatkan bahwa model regresi robust dengan pembobot Huber memiliki galat baku sisaan terkecil yaitu sebesar 1.076 dan nilai AIC yang relatif kecil. Dengan demikian, model regresi linier terbaik yang terpilih adalah model regresi robust Estimasi-M dengan pembobot Huber.
R-Squared
preds.huber <- model.huber.fix$pred
preds.tukbis <- model.tukbis.fix$pred
actual <- dt.new$Yrss.huber <- sum((preds.huber - actual) ^ 2)
tss.huber <- sum((actual - mean(actual)) ^ 2)
rsq.huber <- 1 - rss.huber/tss.huber
rsq.huber## [1] 1rss.huber <- sum((preds.tukbis - actual) ^ 2)
tss.huber <- sum((actual - mean(actual)) ^ 2)
rsq.huber <- 1 - rss.huber/tss.huber
rsq.huber## [1] 1Model OLS stepwise : 0.6664
Model Regresi Robust: Pembobot Huber Stepwise : 1
Model Regresi Robust: Pembobot Tukey Bisquare Stepwise : 1
10. Model Terbaik
Model regresi linier terbaik yang terpilih adalah Model Regresi Robust Estimasi-M dengan pembobot Huber.
11. Interpretasi Model
Persamaan Model Regresi Robust Estimasi-M dengan pembobot Huber dinyatakan sengai berikut.
\[\widehat{Y} = 23.7399 - 0.2870X2 - 0.0006X5\]
Interpetasi :
Intersep
Dugaan rata-rata prevalensi Orang yang mengidap penyakit diabetes melitus tipe 1 dan 2 di Jawa Barat tahun 2021 jika tidak dapat dapat dijelaskan oleh persentase penduduk usia 15-24 dan >= 55 tahun yang merokok dalam sebulan terakhir dan rata-rata pengeluaran pembelian beras perkapita seminggu sebesar 23.7399 jiwa/1000 orang.
X2
Setiap kenaikan satu persen penduduk usia 15-24 dan >= 55 tahun yang merokok dalam sebulan terakhir (\(X2\)), maka akan menurunkan prevalensi Orang yang mengidap penyakit diabetes melitus tipe 1 dan 2 di Jawa Barat tahun 2021 sebesar 0.2870 apabila rata-rata pengeluaran pembelian beras perkapita seminggu (\(X5\)) tetap.
X5
Setiap kenaikan satu rupiah/kapita/minggu rata-rata pengeluaran pembelian beras perkapita seminggu (\(X2\)), maka akan menurunkan prevalensi Orang yang mengidap penyakit diabetes melitus tipe 1 dan 2 di Jawa Barat tahun 2021 sebesar 0.0006 apabila persentase penduduk usia 15-24 dan >= 55 tahun yang merokok dalam sebulan terakhir (\(X5\)) tetap.