Penanganan Autokorelasi Data IPM Provinsi Sumatera Utara Tahun 2010 - 2021
Diftania Putri Imani
2/16/2022
Pendahuluan
IPM (Indeks Pembangunan Manusia) merupakan indikator penting untuk mengukur keberhasilan dalam upaya membangun kualitas hidup manusia (masyarakat/penduduk). IPM menjelaskan bagaimana penduduk dapat mengakses hasil pembangunan dalam memperoleh pendapatan, kesehatan, pendidikan, dan sebagainya. IPM merupakan data strategis bagi Indonesia karena IPM digunakan sebagai salah satu alokato penentuan Dana Alokasi Umum (DAU).
Nilai IPB berkisar antara 0 sampai dengan 100. Semakin tinggi nilai IPM suatu negara/daerah, menunjukkan pencapaian pembangunan manusia di kawasan tersebut semakin baik. Capaian IPM di suatu wilayah dapat dikelompokkan menjadi empat kategori yaitu Rendah (IPM < 60), Sedang (60 ≤ IPM < 70), Tinggi (70 ≤ IPM < 80), dan Sangat Tinggi (IPM ≥ 80) (www.bps.go.id)
Identifikasi Autokorelasi
Pertama-tama, set beberapa library seperti di bawah ini.
library(readxl)## Warning: package 'readxl' was built under R version 4.0.5library(orcutt)## Warning: package 'orcutt' was built under R version 4.0.5## Loading required package: lmtest## Warning: package 'lmtest' was built under R version 4.0.5## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.0.5##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(lmtest)
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.0.5## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.4 v dplyr 1.0.7
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 2.0.1 v forcats 0.5.1## Warning: package 'ggplot2' was built under R version 4.0.5## Warning: package 'tibble' was built under R version 4.0.5## Warning: package 'tidyr' was built under R version 4.0.5## Warning: package 'readr' was built under R version 4.0.5## Warning: package 'purrr' was built under R version 4.0.5## Warning: package 'dplyr' was built under R version 4.0.5## Warning: package 'stringr' was built under R version 4.0.5## Warning: package 'forcats' was built under R version 4.0.5## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(forecast)## Warning: package 'forecast' was built under R version 4.0.5## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zooKemudian, tentukan directory dari file data yang akan digunakan. Data yang digunakan pada kasus ini adalah data IPM Provinsi Sumatera Utara pada tahun 2010 - 2021.
setwd("C:/Users/LENOVO/Downloads/Project-20220127")
sumut <- read_excel("sumut.xlsx")
attach(sumut)Selanjutnya, data IPM Sumatera Utara diubah menjadi data time series dan didapatkan peramalan 10 periode ke depan sebagai berikut.
sumut.ts<-ts(sumut$ipm)
desregresi.opt<- HoltWinters(sumut.ts, gamma = FALSE)
desregresi.opt## Holt-Winters exponential smoothing with trend and without seasonal component.
##
## Call:
## HoltWinters(x = sumut.ts, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 1
## beta : 0.7131977
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 72.000000
## b 0.217053plot(desregresi.opt)
ramalandesopt<- forecast(desregresi.opt, h=10)
ramalandesopt## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 13 72.21705 71.94012 72.49399 71.79352 72.64059
## 14 72.43411 71.88475 72.98346 71.59394 73.27427
## 15 72.65116 71.78322 73.51909 71.32377 73.97855
## 16 72.86821 71.63969 74.09674 70.98934 74.74708
## 17 73.08526 71.45809 74.71244 70.59671 75.57382
## 18 73.30232 71.24159 75.36305 70.15070 76.45394
## 19 73.51937 70.99271 76.04604 69.65517 77.38357
## 20 73.73642 70.71352 76.75933 69.11328 78.35956
## 21 73.95348 70.40574 77.50122 68.52768 79.37928
## 22 74.17053 70.07083 78.27023 67.90059 80.44047plot(ramalandesopt)
Data nilai IPM Sumatera Utara Tahun 2010 - 2021 merupakan data trend sehingga nilai peramalan 10 periode ke depan tidak konstan.
Untuk mengidentifikasi sebaran data, dibuat diagram pencar seperti di bawah ini.
plot(sumut$ipm~sumut$tahun, main = "Scatter Plot IPM Sumatera Utara vs Tahun", xlab = "Tahun", ylab = "IPM",
pch = 19, col="#ab6881")
cor(tahun,ipm)## [1] 0.9932939Berdasarkan diagram pencar di atas, terlihat bahwa nilai IPM cenderung meningkat setiap tahunnya dan terdapat hubungan positif antara nilai IPM dengan tahun. Diketahui pula nilai korelasi yang cukup tinggi antara peubah IPM dan Tahun yaitu sebesar 0.9932939.
Apabla data diregresikan, maka model regresi sederhana yang dihasilkan adalah sebagai berikut
model <- lm(ipm~tahun, data=sumut)
summary(model)##
## Call:
## lm(formula = ipm ~ tahun, data = sumut)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.40333 -0.11958 0.00167 0.13542 0.32667
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -927.99167 36.72254 -25.27 2.16e-10 ***
## tahun 0.49500 0.01822 27.17 1.06e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2179 on 10 degrees of freedom
## Multiple R-squared: 0.9866, Adjusted R-squared: 0.9853
## F-statistic: 738.1 on 1 and 10 DF, p-value: 1.056e-10Berdasarkan hasil regresi di atas, model yang terbentuk adalah y = -927.99167 + 0.49500x. Hasil uji F dengan p-value = 1.056e-10 < 0.05 menunjukkan bahwa terdapat minimal satu variabel yang berpengaruh nyata pada model pada taraf nyata 5%. Hasil uji parsial juga menunjukkan baik intersep maupun peubah tahun signifikan pada taraf nyata 5%. Nilai R-squared = 0.9853 menunjukkan bahwa sebesar 98.53% keragaman IPM Sumatera Utara dapat dijelaskan oleh peubah tahun.
Sisaan model kemudian diplotkan untuk mengidentifikasi autokorelasi secara eksploratif.
residu <- residuals(model)
fit <- predict(model)
par(mfrow = c(2,2))
qqnorm(residu)
qqline(residu, col = "#ab6881", lwd = 2)
plot(fit, residu, col = "#ab6881", pch = 20, xlab = "Sisaan", ylab = "Fitted Values", main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)
hist(residu, col = "#ab6881")
plot(seq(1,12,1), residu, col = "#ab6881", pch = 20, xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,12,1), residu, col = "red")
abline(a = 0, b = 0, lwd = 2)
Berdasarkan grafik di ats, Normal Q-Q Plot menunjukkan sisaan dari model regresi cenderung menyebar normal. Plot Sisaan vs Fitted Values juga menunjukkan bahwa ragam sisaan homogen karena tidak menunjukkan adanya pola. Autokorelasi dapat diidentifikasi dari plot Sisaan vs Order yang memiliki pola naik-turun.
par(mfrow = c(1,2))
acf(residu)
pacf(residu)
Berdasarkan hasil plot di atas, data tersebut tidak memiliki autokorelasi karena data masih terlihat berada di antara batas autokorelasi. Untuk memastikan hal tersebut, dilakukan uji lain dalam mengidentifikasi autokorelasi.
dwtest(model)##
## Durbin-Watson test
##
## data: model
## DW = 0.87521, p-value = 0.00352
## alternative hypothesis: true autocorrelation is greater than 0H0: Tidak ada autokorelasi vs H1: Ada autokorelasi. Berdasarkan hasil uji Durbin-Watson, didapatkan p-value = 0.00352 < 0.05 yang artinya Tolak H0. Cukup bukti untuk menyatakan bahwa data IPM Sumatera Utara tahun 2010-2021 memiliki autokorelasi pada taraf nyata 5%.
bgtest(ipm~tahun, data=sumut, order=1)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: ipm ~ tahun
## LM test = 3.1615, df = 1, p-value = 0.07539H0: Ada autokorelasi vs H1: Tidak ada autokorelasi. Berdasarkan hasil uji Breusch-Godgrey, didapatkan p-value = 0.07539 > 0.05 yang artinya tak tolak H0. Tidak cukup bukti untuk menyatakan bahwa data IPM Sumatera Utara tahun 2010-2021 tidak memiliki autokorelasi pada taraf nyata 5%.
Penanganan Autokorelasi dengan Cochrane-Orcutt
modelco <- cochrane.orcutt(model)
modelco## Cochrane-orcutt estimation for first order autocorrelation
##
## Call:
## lm(formula = ipm ~ tahun, data = sumut)
##
## number of interaction: 28
## rho 0.673825
##
## Durbin-Watson statistic
## (original): 0.87521 , p-value: 3.52e-03
## (transformed): 1.01601 , p-value: 1.053e-02
##
## coefficients:
## (Intercept) tahun
## -880.82056 0.47157summary(modelco)## Call:
## lm(formula = ipm ~ tahun, data = sumut)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -880.820555 114.547419 -7.690 3.032e-05 ***
## tahun 0.471570 0.056761 8.308 1.635e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1942 on 9 degrees of freedom
## Multiple R-squared: 0.8846 , Adjusted R-squared: 0.8718
## F-statistic: 69 on 1 and 9 DF, p-value: < 1.635e-05
##
## Durbin-Watson statistic
## (original): 0.87521 , p-value: 3.52e-03
## (transformed): 1.01601 , p-value: 1.053e-02rho <- modelco$rho
rho## [1] 0.6738248y.trans<- ipm[-1]-ipm[-12]*rho
x.trans<- tahun[-1]-tahun[-12]*rho
modelcorho<- lm(y.trans~x.trans)
summary(modelcorho)##
## Call:
## lm(formula = y.trans ~ x.trans)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.29120 -0.16472 0.05668 0.14722 0.21480
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -287.30179 37.36252 -7.690 3.03e-05 ***
## x.trans 0.47157 0.05676 8.308 1.63e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1942 on 9 degrees of freedom
## Multiple R-squared: 0.8846, Adjusted R-squared: 0.8718
## F-statistic: 69.02 on 1 and 9 DF, p-value: 1.635e-05Berdasarkan hasil tersebut, didapatkan persamaan hasil yang belum ditransformasi balik sebesar y = -287.30179 + 0.47157.
cat("y = ", coef(modelcorho)[1]/(1-rho), "+", coef(modelcorho)[2], "x", sep = "")## y = -880.8206+0.4715699xBerdasarkan hasil tersebut, didapatkan persamaan hasil yang sudah ditransformasi balik sebesar y = -880.82056 + 0.47157x dengan rho optimum sebesar 0.6738248. Didapatkan pula nilai p-value Durbin-Watson yang baru hasil transformasi sebesar 1.053e-02. Sehingga dapat disimpulkan bahwa cukup bukti untuk menyatakan bahwa data memiliki autokorelasi pada taraf nyata 5%. Oleh karena itu, Cochrane-Orcutt belum mampu menangani masalah autokorelasi pada kasus ini.
Penanganan Autokorelasi dengan 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))
}
r <- c(seq(0.1,0.8, by= 0.1), seq(0.9,0.99, by= 0.01))
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4)## rho SSE
## 1 0.10 0.4197
## 2 0.20 0.3941
## 3 0.30 0.3734
## 4 0.40 0.3576
## 5 0.50 0.3467
## 6 0.60 0.3407
## 7 0.70 0.3395
## 8 0.80 0.3432
## 9 0.90 0.3518
## 10 0.91 0.3530
## 11 0.92 0.3541
## 12 0.93 0.3554
## 13 0.94 0.3566
## 14 0.95 0.3580
## 15 0.96 0.3593
## 16 0.97 0.3607
## 17 0.98 0.3622
## 18 0.99 0.3637Berdasarkan hasil di atas, rho optimum ada di kisaran 0.7 karena memiliki SSE yang paling kecil yaitu sebesar 0.3395. Untuk melihat secara spesifik rho yang optimum, maka selang di sekitar 0.7 dapat diperkecil.
r <- seq(0.65,0.70, by= 0.01)
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4)## rho SSE
## 1 0.65 0.3395
## 2 0.66 0.3394
## 3 0.67 0.3393
## 4 0.68 0.3393
## 5 0.69 0.3394
## 6 0.70 0.3395Pada tabel di atas terlihat bahwa nilai rho optimum sebesar 0.67. Plot antara rho dan SSE dapat dilihat di bawah ini.
plot(tab$SSE ~ tab$rho , type = "l", xlab = "Nilai Rho", ylab = "SSE",
main = "Plot SSE vs Rho Hildreth-Lu", lwd = 1.5)
abline(v = 0.67, lty = 2, col = "#ab6881", lwd = 2)
text(x = 0.67, y = 0.35, labels = "rho = 0.67", cex = 0.72)
modelhlu <- hildreth.lu.func(0.67, model)
summary(modelhlu)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.29043 -0.16445 0.05684 0.14709 0.21505
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -291.0211 37.3621 -7.789 2.74e-05 ***
## x 0.4721 0.0561 8.415 1.47e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1942 on 9 degrees of freedom
## Multiple R-squared: 0.8872, Adjusted R-squared: 0.8747
## F-statistic: 70.81 on 1 and 9 DF, p-value: 1.475e-05Dengan memasukkan nilai rho optimum 0.67, didapatkan persamaan yang belum ditransformasi balik yaitu y = -291.0211 + 0.4721x
# Transformasi Balik
cat("y = ", coef(modelhlu)[1]/(1-0.67), "+", coef(modelhlu)[2],"x", sep = "")## y = -881.8823+0.4720964xDengan memasukkan nilai rho optimum 0.67, didapatkan persamaan yang sudah ditransformasi balik yaitu y = -881.8823 + 0.4720964x
dwtest(modelhlu)##
## Durbin-Watson test
##
## data: modelhlu
## DW = 1.0135, p-value = 0.01037
## alternative hypothesis: true autocorrelation is greater than 0Setelah dilakukan uji Durbin Watson pada data yang telah ditangani, didapatkan p-value sebesar 0.01037 < 0.05 yang artinya tolak H0. Cukup bukti untuk menyatakan bahwa data tersebut memiiki autokorelasi pada taraf nyata 5%. Oleh karena itu, metode Hildreth Lu belum mampu menangani masalah autokorelasi pada kasus ini.
Untuk membandingkan antara model OLS, model hasil penanganan metode Cochrane-Orcutt, dan model hasil penanganan metode Hildreth-Lu, dicari nilai MSE dari masing-masing model.
modelolser <- mean(model$residual^2)
modelolser## [1] 0.03955972modelcoerr <- mean(modelco$residual^2)
modelcoerr## [1] 0.0488397modelhluer <- mean(modelhlu$residual^2)
modelhluer## [1] 0.03084947Berdasarkan hasil di atas, nilai MSE model hasil penanganan metode Hildreth-Lu lebih kecil dibandingkan model lain sehingga model ini dianggap sebagai model yang paling baik.
Simpulan
Baik penanganan Cochrane-Orcutt maupun Hildreth-Lu belum dapat menghilangkan autokorelasi pada kasus ini. Namun, jika dihitung nilai MSE-nya, model hasil penanganan Hildreth-Lu memiliki nilai yang lebih kecil dibandingkan model lain sehingga model ini dianggap sebagai model yang paling baik.