Project Regresi Linear Berganda
Sandra Berliana Putri - G1401221021
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.5.3## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.5.3##
## 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.5.3## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 4.5.3library(randtests)
library(nortest)
library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectlibrary(car)## Warning: package 'car' was built under R version 4.5.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.5.3##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(olsrr)## Warning: package 'olsrr' was built under R version 4.5.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.5.3## Warning: package 'tidyr' was built under R version 4.5.3## Warning: package 'readr' was built under R version 4.5.3## Warning: package 'forcats' was built under R version 4.5.3## Warning: package 'lubridate' was built under R version 4.5.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ lubridate 1.9.5 ✔ tibble 3.3.1
## ✔ purrr 1.2.1 ✔ tidyr 1.3.2
## ✔ readr 2.2.0## ── 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\\ASUS\\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-05plot(model1)
Didapatkan 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
summary(Data_Anreg$Y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 62.93 75.61 78.78 77.95 80.91 85.62boxplot(Data_Anreg$Y, col = "lightblue")
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 signifikan dengan Y sehingga dapat dikeluarkan dari model. X8 juga dapat dikeluarkan dari model dikarenakan korelasi yang tinggi dengan peubah X5.
model2= lm(Y~X1+X2+X5+X7+X9 , data=Data_Anreg)
summary(model2)##
## Call:
## lm(formula = Y ~ X1 + X2 + X5 + X7 + X9, data = Data_Anreg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.8704 -2.4087 0.2841 1.9428 5.9382
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -19.062627 14.174112 -1.345 0.18945
## X1 0.381941 0.189058 2.020 0.05302 .
## X2 0.004925 0.002122 2.322 0.02776 *
## X5 0.540278 0.172446 3.133 0.00403 **
## X7 0.162398 0.153261 1.060 0.29837
## X9 0.209361 0.123018 1.702 0.09986 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.185 on 28 degrees of freedom
## Multiple R-squared: 0.7084, Adjusted R-squared: 0.6564
## F-statistic: 13.61 on 5 and 28 DF, p-value: 8.922e-07Pengujian Multikolinearitas
vif(model2)## X1 X2 X5 X7 X9
## 1.493739 1.142603 1.472333 1.761202 1.638213Melalui uji vif didapatkan nilai VIF dari tiap peubah yang tersisa kurang dari 10 sehingga dapat disimpulkan bahwa tidak terdapat multikolinearitas
Pendeteksian 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.4449 0.6770 0.1800
## 2 2 0.3973 0.1573 0.1355
## 3 3 -0.3988 -0.0813 0.0399
## 4 4 -2.1074 -0.8531 0.1408
## 5 5 -0.6354 -0.2402 0.1251
## 6 6 1.3550 0.4820 0.1123
## 7 7 -0.9434 -0.4404 0.1789
## 8 8 0.5977 0.3037 0.2052
## 9 9 0.4201 0.1213 0.0769
## 10 10 -0.7436 -0.2533 0.1039
## 11 11 -0.8783 -0.4853 0.2339
## 12 12 -1.8877 -1.8739 0.4963
## 13 13 0.3556 0.1783 0.2008
## 14 14 0.7467 0.4783 0.2910
## 15 15 2.2649 1.2106 0.2222
## 16 16 0.4289 0.1518 0.1114
## 17 17 0.0104 0.0050 0.1870
## 18 18 -1.7191 -0.7139 0.1471
## 19 19 1.5066 0.5934 0.1343
## 20 20 0.0383 0.0191 0.1995
## 21 21 0.1080 0.0394 0.1172
## 22 22 0.7428 0.2395 0.0942
## 23 23 0.0797 0.0327 0.1442
## 24 24 0.2066 0.0570 0.0707
## 25 25 -0.9434 -0.3438 0.1172
## 26 26 1.2276 0.3806 0.0877
## 27 27 0.8690 0.2945 0.1030
## 28 28 0.8002 0.2436 0.0848
## 29 29 -1.1423 -0.6925 0.2687
## 30 30 -0.1956 -0.0911 0.1783
## 31 31 0.4891 0.1971 0.1397
## 32 32 -2.1598 -1.8303 0.4180
## 33 33 -0.9237 -0.5970 0.2946
## 34 34 -0.5701 -0.4270 0.3594for (i in 1:dim(hasil)[1]){
absri <- abs(hasil$ri)
pencilan <- which(absri > 2)
}
pencilan## [1] 4 15 32titik_leverage <- vector("list", dim(hasil)[1])
for (i in 1:dim(hasil)[1]) {
cutoff <- 2 * 6/ 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((6/ 34))
amatan_berpengaruh[[i]] <- which(abs(DFFITSi) > cutoff)
}
berpengaruh <- unlist(amatan_berpengaruh)
amatan_berpengaruh <- sort(unique(berpengaruh))
amatan_berpengaruh## [1] 4 12 15 32ols_plot_resid_lev(model2)
Dari proses analisis sebelumnya didapatkan 3 buah pencilan yaitu amatan ke-4, ke-15, dan ke-32(Riau, Jawa Timur, Maluku Utara). Terdapat pula 3 buah leverage yaitu amatan ke-12, ke-32, dan ke-34 (Jawa Barat, Maluku Utara, dan Papua). dan 4 buah amatan berpengaruh yaittu amatan ke-4, ke-12, ke-15, dan amatan ke-32. Saat dicobakan nilai adj R-square tertinggi didapatkan dengan menghapus amatan ke-4, ke-12, dan ke-15
# model dengan mengeluarkan amatan ke-4, ke-12 dan ke-15
db1= Data_Anreg[-c(34),]
db= Data_Anreg[-c(4,15,12),]
model3 = lm(Y~X1+X2+X5+X7+X9 , data=db)
summary(model3)##
## Call:
## lm(formula = Y ~ X1 + X2 + X5 + X7 + X9, data = db)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.3766 -1.3634 -0.1341 1.7843 4.4156
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -34.549356 12.210193 -2.830 0.00906 **
## X1 0.623983 0.169406 3.683 0.00111 **
## X2 0.007434 0.002465 3.016 0.00581 **
## X5 0.579666 0.142637 4.064 0.00042 ***
## X7 0.201312 0.128735 1.564 0.13044
## X9 0.067596 0.107376 0.630 0.53472
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.607 on 25 degrees of freedom
## Multiple R-squared: 0.806, Adjusted R-squared: 0.7672
## F-statistic: 20.77 on 5 and 25 DF, p-value: 3.57e-08anova(model3)## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 397.25 397.25 58.4561 5.313e-08 ***
## X2 1 64.21 64.21 9.4487 0.005053 **
## X5 1 222.63 222.63 32.7595 5.810e-06 ***
## X7 1 18.98 18.98 2.7935 0.107122
## X9 1 2.69 2.69 0.3963 0.534717
## Residuals 25 169.89 6.80
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Uji 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(model3,1)
Terlihat sisaan menyebar disekitar 0 sehingga nilai harapan sisaan 0
Formal
t.test(model3$residuals,mu = 0,conf.level = 0.95)##
## One Sample t-test
##
## data: model3$residuals
## t = 1.6741e-17, df = 30, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.872893 0.872893
## sample estimates:
## mean of x
## 7.155297e-18Uji-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(model3,1)
Lebar pita hampir sama untuk setiap dugaan sehingga asumsi kehomogenan ragam terpenuhi
Formal
bptest(model3)##
## studentized Breusch-Pagan test
##
## data: model3
## BP = 2.4192, df = 5, p-value = 0.7886Melalui 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(db)[1],
y = model3$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(model3)##
## Durbin-Watson test
##
## data: model3
## DW = 2.3538, p-value = 0.7988
## alternative hypothesis: true autocorrelation is greater than 0Pengujian autokorelasi menggunakan Uji Durbin-Watson menghasilkan nilai p-value 0.1119 > 0.05 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(model3,2)
Melalui grafik tersebut terlihat data mengikuti garis, sehingga asumsi kenormalan sisaan terpenuhi.
Formal
ks.test(model3$residuals, "pnorm", mean=mean(model3$residuals), sd=sd(model3$residuals))##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: model3$residuals
## D = 0.079584, p-value = 0.9806
## alternative hypothesis: two-sidedPemilihan model terbaik
Forward
step(model3,direction="forward")## Start: AIC=64.74
## Y ~ X1 + X2 + X5 + X7 + X9##
## Call:
## lm(formula = Y ~ X1 + X2 + X5 + X7 + X9, data = db)
##
## Coefficients:
## (Intercept) X1 X2 X5 X7 X9
## -34.549356 0.623983 0.007434 0.579666 0.201312 0.067596Backward
step(model3,direction="backward")## Start: AIC=64.74
## Y ~ X1 + X2 + X5 + X7 + X9
##
## Df Sum of Sq RSS AIC
## - X9 1 2.693 172.59 63.224
## <none> 169.89 64.737
## - X7 1 16.618 186.51 65.630
## - X2 1 61.826 231.72 72.358
## - X1 1 92.199 262.09 76.176
## - X5 1 112.235 282.13 78.460
##
## Step: AIC=63.22
## Y ~ X1 + X2 + X5 + X7
##
## Df Sum of Sq RSS AIC
## <none> 172.59 63.224
## - X7 1 18.984 191.57 64.459
## - X2 1 73.525 246.11 72.226
## - X1 1 113.180 285.77 76.857
## - X5 1 145.273 317.86 80.156##
## Call:
## lm(formula = Y ~ X1 + X2 + X5 + X7, data = db)
##
## Coefficients:
## (Intercept) X1 X2 X5 X7
## -35.838295 0.657121 0.007834 0.612837 0.212938Stepwise
step(model3,direction="both")## Start: AIC=64.74
## Y ~ X1 + X2 + X5 + X7 + X9
##
## Df Sum of Sq RSS AIC
## - X9 1 2.693 172.59 63.224
## <none> 169.89 64.737
## - X7 1 16.618 186.51 65.630
## - X2 1 61.826 231.72 72.358
## - X1 1 92.199 262.09 76.176
## - X5 1 112.235 282.13 78.460
##
## Step: AIC=63.22
## Y ~ X1 + X2 + X5 + X7
##
## Df Sum of Sq RSS AIC
## <none> 172.59 63.224
## - X7 1 18.984 191.57 64.459
## + X9 1 2.693 169.89 64.737
## - X2 1 73.525 246.11 72.226
## - X1 1 113.180 285.77 76.857
## - X5 1 145.273 317.86 80.156##
## Call:
## lm(formula = Y ~ X1 + X2 + X5 + X7, data = db)
##
## Coefficients:
## (Intercept) X1 X2 X5 X7
## -35.838295 0.657121 0.007834 0.612837 0.212938best <- ols_step_best_subset(model3)
best## Best Subsets Regression
## -----------------------------
## Model Index Predictors
## -----------------------------
## 1 X1
## 2 X1 X5
## 3 X1 X2 X5
## 4 X1 X2 X5 X7
## 5 X1 X2 X5 X7 X9
## -----------------------------
##
## Subsets Regression Summary
## ----------------------------------------------------------------------------------------------------------------------------------
## Adj. Pred
## Model R-Square R-Square R-Square C(p) AIC SBIC SBC MSEP FPE HSP APC
## ----------------------------------------------------------------------------------------------------------------------------------
## 1 0.4537 0.4348 0.3523 43.3980 178.8049 87.9657 183.1069 511.4768 17.5611 0.5892 0.6217
## 2 0.6996 0.6782 0.6122 13.7055 162.2613 73.0133 167.9972 291.6307 10.3031 0.3479 0.3648
## 3 0.7812 0.7569 0.7066 5.1898 154.4336 67.2370 161.6035 220.5678 8.0107 0.2729 0.2836
## 4 0.8029 0.7726 0.7034 4.3963 153.1985 67.3335 161.8025 206.6590 7.7086 0.2655 0.2729
## 5 0.8060 0.7672 0.6941 6.0000 154.7110 69.5016 164.7489 211.9106 8.1111 0.2832 0.2871
## ----------------------------------------------------------------------------------------------------------------------------------
## 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 CriteriaMetode backward, stepwise dan best subtest regression menghasilkan peubah penjelas yang sama sebagai model terbaik yaitu X1, X2, X5, X7 sedangkan metode forward X9 tetap dimasukkan. namun berdasarkan hasil dari best subtest regression kombinasi prediktor dengan nilai adj. R-square tertinggi adalah model dengan 4 peubah penjelas yaitu X1, X2, X5, X7. sehingga didapatkan model \[ \hat Y = -35.84+0.657X_1+0.008X_2+0.613X_5+0.213X_7 \]
model4=lm(formula = Y ~ X1 + X2 + X5 + X7, data = db)
summary(model4)##
## Call:
## lm(formula = Y ~ X1 + X2 + X5 + X7, data = db)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.3391 -1.4916 -0.0157 1.5579 4.5025
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -35.838295 11.896733 -3.012 0.005710 **
## X1 0.657121 0.159139 4.129 0.000334 ***
## X2 0.007834 0.002354 3.328 0.002618 **
## X5 0.612837 0.131000 4.678 7.86e-05 ***
## X7 0.212938 0.125916 1.691 0.102769
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.576 on 26 degrees of freedom
## Multiple R-squared: 0.8029, Adjusted R-squared: 0.7726
## F-statistic: 26.48 on 4 and 26 DF, p-value: 7.746e-09t.test(model4$residuals,mu = 0,conf.level = 0.95)##
## One Sample t-test
##
## data: model4$residuals
## t = -2.4941e-16, df = 30, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.8797844 0.8797844
## sample estimates:
## mean of x
## -1.074409e-16bptest(model4)##
## studentized Breusch-Pagan test
##
## data: model4
## BP = 2.2134, df = 4, p-value = 0.6966dwtest(model4)##
## Durbin-Watson test
##
## data: model4
## DW = 2.324, p-value = 0.7867
## alternative hypothesis: true autocorrelation is greater than 0ks.test(model4$residuals, "pnorm", mean=mean(model4$residuals), sd=sd(model4$residuals))##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: model4$residuals
## D = 0.092738, p-value = 0.9301
## alternative hypothesis: two-sidedanova(model4)## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 397.25 397.25 59.8457 3.299e-08 ***
## X2 1 64.21 64.21 9.6733 0.004497 **
## X5 1 222.63 222.63 33.5383 4.225e-06 ***
## X7 1 18.98 18.98 2.8599 0.102769
## Residuals 26 172.59 6.64
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pengujian asumsi kembali model terbaik didapatkan bahwa keseluruhan asumsi terpenuhi dengan nilai p-value > 0.05.
Untuk pengujian kelayakan model sendiri melalui uji F-Simultan dan uji t-parsial. Uji f- simultan mengahsilkan nilai F- statistik dari model sebesar 26.48 dengan p-value $7.746^{-9} $ dan untuk uji t parsial didapatkan nilai intersep, X1, X2, X5 yang signifikan pada taraf 5% sedangkan X7 tidak signifikan.