Projek Regresi
Sintong M.N Purba - G1401221077
2024-04-21
library(readxl)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(lmtest)## Warning: package 'lmtest' was built under R version 4.3.2## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.3.2##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(GGally)## Warning: package 'GGally' was built under R version 4.3.3## Loading required package: ggplot2## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2library(randtests)
library(nortest)
library(MASS)## Warning: package 'MASS' was built under R version 4.3.3##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectlibrary(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(olsrr)## Warning: package 'olsrr' was built under R version 4.3.3##
## Attaching package: 'olsrr'## The following object is masked from 'package:MASS':
##
## cement## The following object is masked from 'package:datasets':
##
## riverslibrary(tidyverse)## Warning: package 'tidyverse' was built under R version 4.3.2## Warning: package 'lubridate' was built under R version 4.3.2## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.0
## ✔ readr 2.1.4## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ car::recode() masks dplyr::recode()
## ✖ MASS::select() masks dplyr::select()
## ✖ purrr::some() masks car::some()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsData
Menggunakan data dengan peubah sebagai berikut :
\(Y\) = Indeks Demokrasi
\(X_1\)= Indeks Kemerdekaan Pers
\(X_2\) = Jumlah Ormas/LSM
\(X_3\) = Jumlah Aturan Yang Membatasi
Kebebasan
\(X_4\) = Persentase
Perempuan di DPRD
\(X_5\) = IPM
\(X_6\) = PDRB (Ribu Rupiah)
\(X_7\) = Indeks Keterbukaan Informasi
Publik
\(X_8\) = Indeks
Pembangunan Gender
\(X_9\) = Nilai
Integritas Pemerintah
Data_Anreg <- read_excel("C:/Users/Poerba/Downloads/Data_Anreg_Kel 9.xlsx")
head(Data_Anreg)## # A tibble: 6 × 10
## Y X1 X2 X3 X4 X5 X6 X7 X8 X9
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 80.8 76.4 238 2 11.1 72.8 38900 79.1 92.2 63.4
## 2 79.5 75.9 676 3 13 72.7 62922 73.4 91.1 66.2
## 3 77.4 78.7 204 16 6.15 73.3 50264 75.4 94.7 70.6
## 4 73.6 82.0 334 3 18.5 73.5 151259 76.7 88.7 64.2
## 5 77.2 83.7 139 4 14.6 72.1 76164 74.0 89.0 69.4
## 6 80.6 81.4 200 2 21.3 70.9 68237 71.0 93.0 65.6Model Awal
model1<- lm(Y~. , data= Data_Anreg)
summary(model1)##
## Call:
## lm(formula = Y ~ ., data = Data_Anreg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.7352 -1.8433 0.2983 2.0116 5.7584
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.220e+01 1.964e+01 -2.149 0.0419 *
## X1 4.323e-01 2.055e-01 2.103 0.0461 *
## X2 5.024e-03 2.213e-03 2.270 0.0325 *
## X3 -4.704e-02 1.811e-01 -0.260 0.7973
## X4 1.330e-02 8.483e-02 0.157 0.8768
## X5 4.199e-01 2.960e-01 1.418 0.1690
## X6 -4.197e-06 1.460e-05 -0.288 0.7762
## X7 1.188e-01 1.681e-01 0.707 0.4866
## X8 3.528e-01 2.776e-01 1.271 0.2159
## X9 1.984e-01 1.275e-01 1.556 0.1329
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.236 on 24 degrees of freedom
## Multiple R-squared: 0.7419, Adjusted R-squared: 0.6451
## F-statistic: 7.666 on 9 and 24 DF, p-value: 3.146e-05Didapatkan model regresi sebagai berikut : \[ \hat Y = -4.220+0.432X_1+0.005X_2-0.047X_3+0.0133X_4-0.42X_5-0.000004X_6+0.119X_7+0.353X_8+0.198X_9 \] Namun, model tersebut belum dapat digunakan karena belum dilakukan pengujian multikolinearitas, asumsi dan kriteria model terbaik.
Eksplorasi Data
ggpairs(Data_Anreg,
upper = list(continuous = wrap('cor', size = 2)),
title = "Matriks Scatterplot dan Korelasi")
Pada matriks diatas tanda * artinya nilai p-value < 0.05, yang artinya terdapat nilai korelasi antara 2 peubah. Pada matriks tersebut, peubah-peubah yang memiliki korelasi terhadap respon Y adalah X1, X2, X5, X7, X8, X9. Sedangkan peubah X3, X4, X6 tidak berkorelasi dengan Y sehingga dapat dikeluarkan dari model.
model2= lm(Y~X1+X2+X5+X7+X8+X9 , data=Data_Anreg)
summary(model2)##
## Call:
## lm(formula = Y ~ X1 + X2 + X5 + X7 + X8 + X9, data = Data_Anreg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.4679 -2.2469 0.5111 1.9693 5.7142
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -41.659321 18.477099 -2.255 0.0325 *
## X1 0.417464 0.182853 2.283 0.0305 *
## X2 0.005055 0.002041 2.476 0.0198 *
## X5 0.363267 0.192478 1.887 0.0699 .
## X7 0.124633 0.148843 0.837 0.4098
## X8 0.390766 0.215735 1.811 0.0812 .
## X9 0.206939 0.118302 1.749 0.0916 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.062 on 27 degrees of freedom
## Multiple R-squared: 0.74, Adjusted R-squared: 0.6823
## F-statistic: 12.81 on 6 and 27 DF, p-value: 8.158e-07Pengujian Multikolinearitas
vif(model2)## X1 X2 X5 X7 X8 X9
## 1.511123 1.143998 1.983678 1.796453 1.618779 1.638422Melalui uji vif didapatkan nilai VIF dari tiap peubah yang tersisa kurang dari 10 sehingga dapat disimpulkan bahwa tidak terdapat multikolinearitas
Uji asumsi
Pengujian asumsi dilakukan dengan metode non-formal menggunakan
visualisasi dan metode formal. Pada analisis regresi terdapat beberapa
asumsi yang harus dipenuhi yaitu :
Nilai Harapan sisaan 0 (E[e]=0)
H0: \(E[\epsilon]=0\)
H1: \(E[\epsilon]\not=0\)
Non-Formal
plot(model2,1)
Terlihat sisaan menyebar disekitar 0 sehingga nilai harapan sisaan 0
Formal
t.test(model2$residuals,mu = 0,conf.level = 0.95)##
## One Sample t-test
##
## data: model2$residuals
## t = -4.2279e-16, df = 33, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.9664885 0.9664885
## sample estimates:
## mean of x
## -2.008445e-16Uji-T menghasilkan p-value 1 > dibanding 0.05 sehingga pada taraf 5% H0 tidak ditolak. Asumsi nilai harapan sisaan 0 terpenuhi.
Ragam sisaan homogen (Homoscedasticity)
H0: \(Var[\epsilon]=\sigma^2 I\)
H1: \(Var[\epsilon]\not=\sigma^2 I\)
Non-Formal
plot(model2,1)
Lebar pita hampir sama untuk setiap dugaan sehingga asumsi kehomogenan ragam terpenuhi
Formal
bptest(model2)##
## studentized Breusch-Pagan test
##
## data: model2
## BP = 3.869, df = 6, p-value = 0.6944Melalui uji Breusch-Pagan didapatkan p-value = 0.6944 > 0.05 sehingga pada taraf 5% H0 tidak ditolak. Asumsi Homoscedasticity terpenuhi.
Sisaan saling bebas/ tidak ada autokorelasi
H0: \(E[\epsilon_i , \epsilon_j]=0\)
H1: \(E[\epsilon_i , \epsilon_j]\not=0\)
Non-Formal
plot(x = 1:dim(Data_Anreg)[1],
y = model2$residuals,
type = 'b',
ylab = "Residuals",
xlab = "Observation")
Melalui plot sisaan dan urutan dapat dilihat bahwa data menyebar acak dan tidak berpola. Sehingga tidak terdapat autokorelasi pada data
Formal
dwtest(model2)##
## Durbin-Watson test
##
## data: model2
## DW = 1.7, p-value = 0.1119
## alternative hypothesis: true autocorrelation is greater than 0Pengujian autokorelasi menggunakan Uji Durbin-Watson menghasilkan nilai p-value 0.1119 > 0.005 sehingga pada taraf 5% tak tolak H0. Asumsi sisaan saling bebas terpenuhi.
Sisaan menyebar normal
H0: Sisaan menyebar normal
H1: Sisaan tidak menyebar normal
Non-Formal
plot(model2,2)
Melalui grafik tersebut terlihat data mengikuti garis, sehingga asumsi kenormalan sisaan terpenuhi.
Formal
shapiro.test(model2$residuals)##
## Shapiro-Wilk normality test
##
## data: model2$residuals
## W = 0.96272, p-value = 0.2914Pendeteksian pencilan, leverage, dan amatan berpengaruh
index <- c(1:34)
ri <- studres(model2)
DFFITSi <- dffits(model2)
hii <- hatvalues(model2)
hasil <- data.frame(index,ri,DFFITSi,hii); round(hasil,4)## index ri DFFITSi hii
## 1 1 1.4441 0.6786 0.1809
## 2 2 0.4249 0.1682 0.1355
## 3 3 -0.8573 -0.2710 0.0909
## 4 4 -1.8093 -0.8361 0.1760
## 5 5 -0.4621 -0.1833 0.1359
## 6 6 0.9828 0.4341 0.1632
## 7 7 -0.9662 -0.4511 0.1790
## 8 8 0.4653 0.2409 0.2113
## 9 9 0.6595 0.2078 0.0903
## 10 10 -0.8943 -0.3106 0.1076
## 11 11 -0.8428 -0.4672 0.2350
## 12 12 -1.6503 -1.6845 0.5103
## 13 13 0.0939 0.0499 0.2201
## 14 14 0.7606 0.4874 0.2910
## 15 15 2.2761 1.2224 0.2239
## 16 16 0.4061 0.1441 0.1118
## 17 17 -0.1144 -0.0556 0.1911
## 18 18 -1.9547 -0.8253 0.1513
## 19 19 0.9298 0.5384 0.2511
## 20 20 0.2548 0.1317 0.2110
## 21 21 0.3048 0.1164 0.1273
## 22 22 0.9098 0.3015 0.0989
## 23 23 1.1551 0.8655 0.3595
## 24 24 0.5833 0.2024 0.1075
## 25 25 -1.4845 -0.6638 0.1666
## 26 26 0.9411 0.3469 0.1196
## 27 27 0.6016 0.2314 0.1289
## 28 28 0.8459 0.2576 0.0849
## 29 29 -0.8435 -0.5468 0.2959
## 30 30 -0.4553 -0.2234 0.1941
## 31 31 0.1123 0.0531 0.1825
## 32 32 -2.5985 -2.2443 0.4272
## 33 33 -0.3089 -0.2478 0.3915
## 34 34 0.0789 0.0711 0.4483for (i in 1:dim(hasil)[1]){
absri <- abs(hasil$ri)
pencilan <- which(absri > 2)
}
pencilan## [1] 15 32titik_leverage <- vector("list", dim(hasil)[1])
for (i in 1:dim(hasil)[1]) {
cutoff <- 2 * 7 / 34
titik_leverage[[i]] <- which(hii > cutoff)
}
leverages <- unlist(titik_leverage)
titik_leverage <- sort(unique(leverages))
titik_leverage## [1] 12 32 34amatan_berpengaruh <- vector("list", dim(hasil)[1])
for (i in 1:dim(hasil)[1]) {
cutoff <- 2 * sqrt((7/ 34))
amatan_berpengaruh[[i]] <- which(abs(DFFITSi) > cutoff)
}
berpengaruh <- unlist(amatan_berpengaruh)
amatan_berpengaruh <- sort(unique(berpengaruh))
amatan_berpengaruh## [1] 12 15 32db= Data_Anreg[-c(12,15),]ols_plot_resid_lev(model2)
Terdapat 2 amatan yang merupakan pencilan yaitu amatan ke-15 dan ke-32, dan 3 amatan yang merupakan leverage yaitu amatan ke-12, ke-32 dan ke-34. Sedangkan amatan berpengaruh terdiri dari amatan ke-12, ke-15, dan ke-32. Namun saat seluruh amatan berpengaruh dimasukkan maupun dikeluarkan dari model terjadi penurunan adjusted R-Squared. setelah dicobakan mengkombinasikan amatan yang dimasukkan dan dikeluarkan didapatkan nilai adjusted R-Squared tertinggi saat amatan ke-12 dan ke-15 dikeluarkan.
# model dengan mengeluarkan amatan ke-12 dan ke-15
db= Data_Anreg[-c(12,15),]
model3 = lm(Y~X1+X2+X5+X7+X8+X9 , data=db)
summary(model3)##
## Call:
## lm(formula = Y ~ X1 + X2 + X5 + X7 + X8 + X9, data = db)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.6775 -1.9039 0.3456 2.1890 4.0467
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -46.547628 17.261624 -2.697 0.0124 *
## X1 0.549607 0.180422 3.046 0.0054 **
## X2 0.006140 0.002659 2.309 0.0295 *
## X5 0.375846 0.178264 2.108 0.0452 *
## X7 0.146967 0.141765 1.037 0.3098
## X8 0.332473 0.202305 1.643 0.1128
## X9 0.161275 0.110841 1.455 0.1581
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.832 on 25 degrees of freedom
## Multiple R-squared: 0.7752, Adjusted R-squared: 0.7212
## F-statistic: 14.36 on 6 and 25 DF, p-value: 4.86e-07Pemilihan model terbaik
Forward
step(model3,direction="forward")## Start: AIC=72.72
## Y ~ X1 + X2 + X5 + X7 + X8 + X9##
## Call:
## lm(formula = Y ~ X1 + X2 + X5 + X7 + X8 + X9, data = db)
##
## Coefficients:
## (Intercept) X1 X2 X5 X7 X8
## -46.54763 0.54961 0.00614 0.37585 0.14697 0.33247
## X9
## 0.16127Backward
step(model3,direction="backward")## Start: AIC=72.72
## Y ~ X1 + X2 + X5 + X7 + X8 + X9
##
## Df Sum of Sq RSS AIC
## - X7 1 8.618 209.08 72.063
## <none> 200.46 72.716
## - X9 1 16.976 217.44 73.318
## - X8 1 21.657 222.12 73.999
## - X5 1 35.644 236.11 75.953
## - X2 1 42.752 243.21 76.903
## - X1 1 74.408 274.87 80.818
##
## Step: AIC=72.06
## Y ~ X1 + X2 + X5 + X8 + X9
##
## Df Sum of Sq RSS AIC
## <none> 209.08 72.063
## - X9 1 20.822 229.90 73.101
## - X8 1 27.317 236.40 73.993
## - X2 1 40.267 249.35 75.700
## - X5 1 43.637 252.72 76.129
## - X1 1 120.300 329.38 84.607##
## Call:
## lm(formula = Y ~ X1 + X2 + X5 + X8 + X9, data = db)
##
## Coefficients:
## (Intercept) X1 X2 X5 X8 X9
## -48.651646 0.630356 0.005944 0.409078 0.368002 0.176945Stepwise
step(model3,direction="both")## Start: AIC=72.72
## Y ~ X1 + X2 + X5 + X7 + X8 + X9
##
## Df Sum of Sq RSS AIC
## - X7 1 8.618 209.08 72.063
## <none> 200.46 72.716
## - X9 1 16.976 217.44 73.318
## - X8 1 21.657 222.12 73.999
## - X5 1 35.644 236.11 75.953
## - X2 1 42.752 243.21 76.903
## - X1 1 74.408 274.87 80.818
##
## Step: AIC=72.06
## Y ~ X1 + X2 + X5 + X8 + X9
##
## Df Sum of Sq RSS AIC
## <none> 209.08 72.063
## + X7 1 8.618 200.46 72.716
## - X9 1 20.822 229.90 73.101
## - X8 1 27.317 236.40 73.993
## - X2 1 40.267 249.35 75.700
## - X5 1 43.637 252.72 76.129
## - X1 1 120.300 329.38 84.607##
## Call:
## lm(formula = Y ~ X1 + X2 + X5 + X8 + X9, data = db)
##
## Coefficients:
## (Intercept) X1 X2 X5 X8 X9
## -48.651646 0.630356 0.005944 0.409078 0.368002 0.176945best <- ols_step_best_subset(model3)
best## Best Subsets Regression
## --------------------------------
## Model Index Predictors
## --------------------------------
## 1 X5
## 2 X1 X5
## 3 X1 X2 X5
## 4 X1 X2 X5 X8
## 5 X1 X2 X5 X8 X9
## 6 X1 X2 X5 X7 X8 X9
## --------------------------------
##
## Subsets Regression Summary
## ----------------------------------------------------------------------------------------------------------------------------------
## Adj. Pred
## Model R-Square R-Square R-Square C(p) AIC SBIC SBC MSEP FPE HSP APC
## ----------------------------------------------------------------------------------------------------------------------------------
## 1 0.4244 0.4052 0.3442 36.0017 185.6096 92.2975 190.0068 547.4837 18.1757 0.5899 0.6524
## 2 0.6398 0.6150 0.5449 14.0482 172.6072 80.5086 178.4702 354.8160 12.1114 0.3955 0.4347
## 3 0.7105 0.6795 0.6237 8.1874 167.6149 76.7562 174.9436 295.7331 10.3699 0.3414 0.3722
## 4 0.7421 0.7039 0.6411 6.6715 165.9134 76.2353 174.7078 273.5613 9.8453 0.3275 0.3534
## 5 0.7655 0.7204 0.6583 6.0747 164.8755 76.6870 175.1356 258.7365 9.5493 0.3217 0.3428
## 6 0.7752 0.7212 0.6164 7.0000 165.5286 78.4797 177.2544 258.4084 9.7725 0.3341 0.3508
## ----------------------------------------------------------------------------------------------------------------------------------
## AIC: Akaike Information Criteria
## SBIC: Sawa's Bayesian Information Criteria
## SBC: Schwarz Bayesian Criteria
## MSEP: Estimated error of prediction, assuming multivariate normality
## FPE: Final Prediction Error
## HSP: Hocking's Sp
## APC: Amemiya Prediction Criteria