Pendeteksian dan Penanganan Autokorelasi pada Data Indeks Pembangunan Manusia (IPM) Provinsi Riau pada Tahun 2010 - 2021
Packages & Library
Packages yang digunakan ada 9, yaitu :
dplyrreadxlforecastTTRtserieslmtestlawstatorcuttHoRM
lapply(c("dplyr","readxl","forecast","TTR","tseries",
"lmtest","orcutt","HoRM","lawstat"),
library,character.only=T)Deklarasi Data IPM Riau 2010-2021
ipmriau<- read_excel("D:/Semester 5/MPDW/Pertemuan 3/Tugas Individu MPDW.xlsx", sheet = 2)
knitr::kable(ipmriau,align = "c")| Tahun | IPM |
|---|---|
| 2010 | 68.65 |
| 2011 | 68.90 |
| 2012 | 69.15 |
| 2013 | 69.91 |
| 2014 | 70.33 |
| 2015 | 70.84 |
| 2016 | 71.20 |
| 2017 | 71.79 |
| 2018 | 72.44 |
| 2019 | 73.00 |
| 2020 | 72.71 |
| 2021 | 72.94 |
Eksplorasi Data
x <- ipmriau$Tahun
y <- ipmriau$IPM
#diagram pencar identifikasi model
plot(x,y,pch = 20, col = "#A10035", main = "Scatter Plot IPM Riau Tahun 2010-2021",
ylab = "IPM", xlab = "Tahun")
Berdasarkan plot di atas, terlihat bahwa antara tahun dengan nilai IPM memiliki hubungan yang linear positif. Artinya, semakin bertambahnya tahun, maka semakin bertambah nilai IPM.
#korelasi x dan y
cor(x,y)## [1] 0.9864531Nilai korelasi sebesar 0.9864531. Artinya, sesuai dengan plot di atas, bahwa antara kedua peubah memiliki hubungan linear positif kuat.
#model regresi
model <- lm(y~x, data = ipmriau)
summary(model)##
## Call:
## lm(formula = y ~ x, data = ipmriau)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.46795 -0.14599 0.01153 0.09640 0.47191
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -815.69072 46.62834 -17.49 7.91e-09 ***
## x 0.43993 0.02313 19.02 3.51e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2767 on 10 degrees of freedom
## Multiple R-squared: 0.9731, Adjusted R-squared: 0.9704
## F-statistic: 361.6 on 1 and 10 DF, p-value: 3.513e-09Deteksi Autokorelasi
Grafik
1. Residual Plot
#sisaan dan fitted value
resi1 <- residuals(model)
fit <- predict(model)
#Diagnostik dengan eksploratif
par(mfrow = c(2,2))
qqnorm(resi1)
qqline(resi1, col = "#A10035", lwd = 2)
plot(fit, resi1, col = "#A10035", pch = 20, xlab = "Sisaan",
ylab = "Fitted Values", main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)
hist(resi1, col = "#A10035")
plot(seq(1,12,1), resi1, col = "#A10035", pch = 20,
xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,12,1), resi1, col = "red")
abline(a = 0, b = 0, lwd = 2)
Berdasarkan Plot Sisaan vs Order membentuk pola tertentu, maka asumsi kebebasan tidak terpenuhi. Artinya, berdasarkan hasil uji eksplorasi terdeteksi adanya gejala autokolerasi.
2. ACF dan PACF Plot
#ACF dan PACF identifikasi autokorelasi
par(mfrow = c(2,1))
acf(resi1)
pacf(resi1)
Berdasarkan plot ACF dan PACF tidak ada yang signifikan. Artinya, sisaan saling bebas atau tidak ada autokorelasi. Namun, perlu dilakukan uji formal untuk analisis lebih lanjut.
Uji Statistik
1. Durbin Watson Test
H0: tidak ada autokorelasi
H1: ada autokorelasi
lmtest::dwtest(model, alternative = 'two.sided') #ada autokorelasi##
## Durbin-Watson test
##
## data: model
## DW = 1.103, p-value = 0.03271
## alternative hypothesis: true autocorrelation is not 0p-value = 0.03271 < 0.05 maka H0 ditolak. Artinya, cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf 5%.
2. Breusch-Godfrey Test
H0: tidak ada autokorelasi
H1: ada autokorelasi
lmtest::bgtest(y ~ x, data=ipmriau, order=1) #Perform Breusch-Godfrey test for first-order serial correlation##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: y ~ x
## LM test = 1.7969, df = 1, p-value = 0.1801p-value = 0.1801 > 0.05, maka tak tolak H0. Artinya, tidak cukup bukti untuk menyatakan bahwa ada autokorelasi pada taraf 5% (tidak terdapat autokorelasi).
3. Run’s Test
H0: tidak ada autokorelasi
H1: ada autokorelasi
lawstat::runs.test(resid(model), alternative = 'two.sided')##
## Runs Test - Two sided
##
## data: resid(model)
## Standardized Runs Statistic = 0.60553, p-value = 0.5448p-value = 0.5448 > 0.05, maka tak tolak H0. Artinya, tidak cukup bukti untuk menyatakan bahwa ada autokorelasi pada taraf 5% (tidak terdapat autokorelasi).
Berdasarkan Uji Statistik di atas didapatkan bahwa p-value = 0.03271 < 0.05 pada uji Durbin Watson yang artinya H0 ditolak (Terdapat Autokorelasi). Sedangkan pada Breusch-Godfrey Test dan Run’s Test p-value > 0.05 yang artinya Tak tolak H0 (Terdapat autokorelasi). Karena pada uji Durbin-Watson terdeteksi adanya Autokorelasi, maka diperlukan adanya penanganan menggunakan Cochran Orcutt dan Hildreth-lu.
Penanganan Autokorelasi
1. Cochrane-Orcutt
#Interactive method using to solve first order autocorrelation problems
modelco <- orcutt::cochrane.orcutt(model)
modelco ## Cochrane-orcutt estimation for first order autocorrelation
##
## Call:
## lm(formula = y ~ x, data = ipmriau)
##
## number of interaction: 16
## rho 0.486948
##
## Durbin-Watson statistic
## (original): 1.10303 , p-value: 1.635e-02
## (transformed): 1.42722 , p-value: 7.532e-02
##
## coefficients:
## (Intercept) x
## -784.312909 0.424346#rho optimum
rho <- modelco$rho
rho## [1] 0.4869478y## [1] 68.65 68.90 69.15 69.91 70.33 70.84 71.20 71.79 72.44 73.00 72.71 72.94y[-1]## [1] 68.90 69.15 69.91 70.33 70.84 71.20 71.79 72.44 73.00 72.71 72.94y[-12]## [1] 68.65 68.90 69.15 69.91 70.33 70.84 71.20 71.79 72.44 73.00 72.71#transformasi terhadap y dan x
(y.trans <- y[-1]-y[-12]*rho)## [1] 35.47103 35.59929 36.23756 36.28748 36.59296 36.70462 37.11931 37.48202
## [9] 37.72550 37.16281 37.53402(x.trans <- x[-1]-x[-12]*rho)## [1] 1032.235 1032.748 1033.261 1033.774 1034.287 1034.800 1035.313 1035.826
## [9] 1036.339 1036.852 1037.365Model dengan Peubah Hasil Transformasi
modelcorho <- lm(y.trans~x.trans)
summary(modelcorho)##
## Call:
## lm(formula = y.trans ~ x.trans)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42773 -0.20483 0.00321 0.17646 0.35267
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -402.39344 51.60786 -7.797 2.72e-05 ***
## x.trans 0.42435 0.04987 8.509 1.35e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2684 on 9 degrees of freedom
## Multiple R-squared: 0.8894, Adjusted R-squared: 0.8771
## F-statistic: 72.4 on 1 and 9 DF, p-value: 1.348e-05Mendeteksi Autokorelasi
lmtest::dwtest(modelcorho,alternative = 'two.sided')##
## Durbin-Watson test
##
## data: modelcorho
## DW = 1.4272, p-value = 0.1506
## alternative hypothesis: true autocorrelation is not 0bgtest(modelcorho)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: modelcorho
## LM test = 0.57296, df = 1, p-value = 0.4491runs.test(resid(modelcorho), alternative = 'two.sided')##
## Runs Test - Two sided
##
## data: resid(modelcorho)
## Standardized Runs Statistic = -0.93314, p-value = 0.3507Berdasarkan hasil uji statistik di atas, didapatkan bahwa p-value > 0.05, maka tak tolak H0. Artinya, tidak cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf 5% (Tidak terdapat autokorelasi). Oleh karena itu dapat disimpulkan bahwa hasil transformasi sudah memenuhi asumsi bebas autokorelasi.
Transformasi Balik
cat("IPM = ", coef(modelcorho)[1]/(1-0.4869478), "+", coef(model)[2]," Tahun", sep = "") #persamaan regresi setelah di transformasi ke persamaan awal## IPM = -784.3129+0.4399301 Tahun2. Hildreth-lu
hildreth.lu.func<- function(r, model){
x <- model.matrix(model)[,-1]
y <- model.response(model.frame(model))
n <- length(y)
t <- 2:n
y <- y[t]-r*y[t-1]
x <- x[t]-r*x[t-1]
return(lm(y~x))
}#mencari rho yang meminimumkan SSE (iteratif)
r <- c(seq(0.1,0.999, by= 0.001))
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
tab$rho[which.min(tab$SSE)]#rho optimal## [1] 0.487#grafik rho dan SSE
plot(tab$SSE ~ tab$rho , type = "l")
abline(v = tab[tab$SSE==min(tab$SSE),"rho"], lty = 3)
#model Hildreth-lu
modelhl <- hildreth.lu.func(0.487, model)
summary(modelhl)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42775 -0.20483 0.00321 0.17647 0.35266
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -402.34820 51.60787 -7.796 2.72e-05 ***
## x 0.42434 0.04988 8.508 1.35e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2684 on 9 degrees of freedom
## Multiple R-squared: 0.8894, Adjusted R-squared: 0.8771
## F-statistic: 72.38 on 1 and 9 DF, p-value: 1.349e-05Mendeteksi Autokorelasi
dwtest(modelhl)##
## Durbin-Watson test
##
## data: modelhl
## DW = 1.4273, p-value = 0.07533
## alternative hypothesis: true autocorrelation is greater than 0bgtest(modelhl)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: modelhl
## LM test = 0.57285, df = 1, p-value = 0.4491runs.test(resid(modelhl), alternative = 'two.sided')##
## Runs Test - Two sided
##
## data: resid(modelhl)
## Standardized Runs Statistic = -0.93314, p-value = 0.3507Berdasarkan hasil uji statistik di atas, didapatkan bahwa p-value > 0.05, maka tak tolak H0. Artinya, tidak cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf 5% (Tidak terdapat autokorelasi). Oleh karena itu dapat disimpulkan bahwa hasil transformasi sudah memenuhi asumsi bebas autokorelasi.
Transformasi Balik
cat("IPM = ", coef(modelhl)[1]/(1-0.487), "+", coef(modelhl)[2]," Tahun", sep = "") #persamaan regresi setelah di transformasi ke persamaan awal## IPM = -784.3045+0.4243417 Tahun