Penanganan Autokorelasi Pada Model Regresi Linear Sederhana

Data IPM Provinsi Jawa Barat Tahun 2010-2021

Library

library(rmarkdown)
library(dplyr)
library(readxl)
library(forecast)
library(TTR)
library(tseries)
library(lmtest) #uji-Durbin Watson
library(orcutt) #Cochrane-Orcutt
library(HoRM) #Hildreth Lu
library(lawstat)
library(ggplot2)
library(knitr)

Importing Data

Dataset yang digunakan adalah IPM (Indeks Pembangunan Manusia) dari Badan Pusat Statistika tahun 2010 sampai 2021 di Provinsi Jawa Barat.

Menurut Badan Pusat Statistik (BPS), 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.

data <- readxl::read_xlsx("C:/Users/user/Downloads/Tugas MPDW.xlsx", sheet ="Jawa Barat")
data
## # A tibble: 12 x 2
##    Tahun   IPM
##    <dbl> <dbl>
##  1  2010  66.2
##  2  2011  66.7
##  3  2012  67.3
##  4  2013  68.2
##  5  2014  68.8
##  6  2015  69.5
##  7  2016  70.0
##  8  2017  70.7
##  9  2018  71.3
## 10  2019  72.0
## 11  2020  72.1
## 12  2021  72.4

Eksplorasi Data

Time Series Plot

x <- data$Tahun
y <- data$IPM
p <- ggplot(data, aes(x=Tahun, y=IPM)) + geom_line(lwd=1.5,col="blue") + geom_point(size=2,col="black")
p +labs(x="Tahun",y = "Indeks Pembangunan Manusia",
         title="Time Series Plot Indeks Pembangunan Manusia Di Provinsi Jawa Barat", subtitle="Tahun 2010-2021") +
  theme_bw()+
  theme(
    plot.title = element_text(size = 14L,
                              face = "bold",
                              hjust = 0.5))

Indeks Pembangunan Manusia di Jawa Barat mengalami kenaikan secara konstan dari tahun ke tahun. Namun pada tahun 2019 ke 2020 mengalami perlambatan kenaikan Indeks Pembangunan Manusia (IPM) jika dibandingkan dengan tahun-tahun sebelumnya.

Korelasi Tahun (x) dan IPM (y)

cor(x,y)
## [1] 0.9939809

Diperoleh nilai Korelasi sangat tinggi yaitu 0.9939809 menunjukkan adanya hubungan linear kuat positif antara tahun dengan IPM.

Model Regresi

model <- lm(y~x, data = data)
summary(model)
## 
## Call:
## lm(formula = y ~ x, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4760 -0.1888  0.1183  0.1785  0.3104 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.146e+03  4.237e+01  -27.05 1.10e-10 ***
## x            6.032e-01  2.102e-02   28.69 6.16e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2514 on 10 degrees of freedom
## Multiple R-squared:  0.988,  Adjusted R-squared:  0.9868 
## F-statistic: 823.2 on 1 and 10 DF,  p-value: 6.159e-11

Diperoleh persamaan regresi awal yaitu IPM = -1146 + 0.6032 Tahun. Pengaruh peubah x (Tahun) terhadap Y (IPM) bernilai positif sebesar 0.603, hal ini menunjukkan bahwa jika waktunya bertambah satu tahun maka dugaan nilai IPM bertambah sebesar 0.603.

Diperoleh Nilai R-Squared yaitu 98.8%, artinya keragaman Y yang mampu dijelaskan oleh peubah penjelas X adalah sebesar 98.8%.

Deteksi Autokorelasi

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 = "red", lwd = 2)

plot(fit, resi1, col = "red", pch = 20, xlab = "Sisaan", 
     ylab = "Fitted Values", main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)

hist(resi1, col = "red",main = "Histogram Residual")

plot(seq(1,12,1), resi1, col = "red", 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 Normal Q-Q Plot, sisaan berada disekitar garis normal dan histogram menunjukkan sisaan menjulur ke kiri artinya banyak data yang bernilai besar

Berdasarkan Plot Sisaan vs vitted values dan sisaan vs order, terlihat adanya pola khusus sehingga sisaan tidak saling bebas (adanya autokorelasi).

2. ACF dan PACF Plot

par(mfrow = c(1,2))
acf(resi1)
pacf(resi1)

Berdasarkan plot ACF di atas, dapat dilihat bahwa hanya lag 0 yang melewati garis sedangkan lag ke-1 hingga ke-12 tidak ada yang melewati garis. Begitu juga dengan plot PACF di atas, dapat dilihat bahwa hanya lag ke-1 hingga ke-12 tidak ada yang melewati garis. Dapat disimpulkan bahwa menurut Plot ACF dan PACF, sisaan Saling bebas atau tidak ada autokorelasi.

Namun kita perlu melakukan uji formal untuk memastikan apakah terdapat autokorelasi pada data IPM ini.

Uji Formal

  1. Durbin Watson Test

Hipotesis :

H0 = Sisaan Saling Bebas (Tidak Autokorelasi)

H1 = Sisaan Tidak Saling Bebas (Autokorelasi)

lmtest::dwtest(model, alternative = 'two.sided')
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 0.79521, p-value = 0.003536
## alternative hypothesis: true autocorrelation is not 0

Berdasarkan Durbin-Watson test didapatkan p-value < 0.05 yang artinya Tolak H0, sehingga cukup bukti untuk menyatakan bahwa terdapat autokorelasi (positif) atau sisaan tidak saling bebas pada taraf 5%.

  1. Run’s Test

Hipotesis

H0 = Sisaan Saling Bebas (Tidak Autokorelasi)

H1 = Sisaan Tidak Saling Bebas (Autokorelasi)

lawstat::runs.test(resid(model), alternative = 'two.sided')
## 
##  Runs Test - Two sided
## 
## data:  resid(model)
## Standardized Runs Statistic = -1.2111, p-value = 0.2259

Berdasarkan uji Run’s Test, p-value > 0.05 yang artinya Tak Tolak H0. Cukup bukti untuk menyatakan bahwa sisaan saling bebas pada taraf 5% (Tidak adanya autokorelasi)

  1. Breusch-Godfrey Test

Hipotesis :

H0 = Sisaan Saling Bebas (Tidak Autokorelasi)

H1 = Sisaan Tidak Saling Bebas (Autokorelasi)

bgtest(y ~ x, data=data, order=1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  y ~ x
## LM test = 3.9403, df = 1, p-value = 0.04714

Berdasarkan uji Breusch-Godfrey, p-value < 0.05 yang artinya Tolak H0. Cukup bukti untuk menyatakan bahwa sisaan tidak saling bebas pada taraf 5% (terdapat autokorelasi), maka diperlukan adanya penanganan menggunakan Cochran Orcutt dan Hildreth-lu.

Penanganan Autokorelasi

Cochran Orcutt

modelco <- orcutt::cochrane.orcutt(model,convergence=5,max.iter = 1000)
modelco 
## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = y ~ x, data = data)
## 
##  number of interaction: 412
##  rho 0.878545
## 
## Durbin-Watson statistic 
## (original):    0.79521 , p-value: 1.768e-03
## (transformed): 2.07953 , p-value: 4.014e-01
##  
##  coefficients: 
## (Intercept)           x 
## -638.590309    0.352236
rho <- modelco$rho
rho
## [1] 0.8785453
(y.trans <- y[-1]-y[-12]*rho)
##  [1] 8.554229 8.747386 9.106331 8.839284 9.056084 8.991103 9.147903 9.195634
##  [9] 9.389721 8.808383 9.115670
(x.trans <- x[-1]-x[-12]*rho)
##  [1] 245.1240 245.2454 245.3669 245.4883 245.6098 245.7312 245.8527 245.9742
##  [9] 246.0956 246.2171 246.3385

Model dengan peubah transformasi

modelcorho <- lm(y.trans~x.trans)
summary(modelcorho) 
## 
## Call:
## lm(formula = y.trans ~ x.trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.35835 -0.08547 -0.00451  0.11199  0.26577 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -77.5598    38.7676  -2.001   0.0765 .
## x.trans       0.3522     0.1578   2.233   0.0525 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.201 on 9 degrees of freedom
## Multiple R-squared:  0.3564, Adjusted R-squared:  0.2849 
## F-statistic: 4.985 on 1 and 9 DF,  p-value: 0.05247

Diperoleh model Cochran :

IPM = -77.5598 + 0.3522 Tahun

Pengecekan Autokorelasi

  1. Durbin-Watson test

Hipotesis :

H0 = Sisaan Saling Bebas (Tidak Autokorelasi)

H1 = Sisaan Tidak Saling Bebas (Autokorelasi)

lmtest::dwtest(modelcorho,alternative = 'two.sided') 
## 
##  Durbin-Watson test
## 
## data:  modelcorho
## DW = 2.0795, p-value = 0.8028
## alternative hypothesis: true autocorrelation is not 0

Didapatkan P-value > 5%. Tak Tolak H0, sehingga belum cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf 5%.

  1. Run’s Test

Hipotesis :

H0 = Sisaan Saling Bebas (Tidak Autokorelasi)

H1 = Sisaan Tidak Saling Bebas (Autokorelasi)

lawstat::runs.test(resid(modelcorho), alternative = 'two.sided')
## 
##  Runs Test - Two sided
## 
## data:  resid(modelcorho)
## Standardized Runs Statistic = -0.93314, p-value = 0.3507

Didapatkan P-value > 5%. Tak Tolak H0, sehingga belum cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf 5%.

  1. Breusch-Godfrey test

Hipotesis :

H0 = Sisaan Saling Bebas (Tidak Autokorelasi)

H1 = Sisaan Tidak Saling Bebas (Autokorelasi)

bgtest(modelcorho)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  modelcorho
## LM test = 0.17252, df = 1, p-value = 0.6779

P-value > 0.05. Tak Tolak H0, sehingga belum cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf 5%.

Berdasarkan hasil di atas dapat disimpulkan hasil transformasi Cochran Orcutt memenuhi asumsi bebas autokorelasi.

Transformasi Balik

Persamaan regresi setelah di transformasi ke persamaan awal

cat("IPM = ", coef(modelcorho)[1]/(1-rho), "+", coef(modelcorho)[2]," Tahun", sep = "")
## IPM = -638.5903+0.3522361 Tahun

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.9, 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))}))
optrho <- which.min(round(tab, 4)[,2])
round(tab, 4)
##     rho    SSE
## 1  0.10 0.5525
## 2  0.20 0.5071
## 3  0.30 0.4680
## 4  0.40 0.4350
## 5  0.50 0.4083
## 6  0.60 0.3878
## 7  0.70 0.3735
## 8  0.80 0.3655
## 9  0.90 0.3636
## 10 0.90 0.3636
## 11 0.91 0.3638
## 12 0.92 0.3640
## 13 0.93 0.3643
## 14 0.94 0.3646
## 15 0.95 0.3650
## 16 0.96 0.3655
## 17 0.97 0.3660
## 18 0.98 0.3666
## 19 0.99 0.3672
round(tab, 4)[optrho,]
##   rho    SSE
## 9 0.9 0.3636

Grafik rho dan SSE

plot(tab$SSE ~ tab$rho , type = "l")
abline(v = tab[tab$SSE==min(tab$SSE),"rho"], lty = 3)

r <- seq(0.1,1, by= 0.01)
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))

# grafik rho dan SSE
plot(tab$SSE ~ tab$rho , type = "l")
abline(v = tab[tab$SSE==min(tab$SSE),"rho"], lty = 3)

Model dengan Peubah Transformasi Hidreth-lu dengan rho optimum 0.9

modelhl <- hildreth.lu.func(0.9, model)
summary(modelhl)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.36185 -0.07965 -0.00773  0.11035  0.26443 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -51.7880    38.8082  -1.334    0.215
## x             0.2928     0.1916   1.528    0.161
## 
## Residual standard error: 0.201 on 9 degrees of freedom
## Multiple R-squared:  0.206,  Adjusted R-squared:  0.1177 
## F-statistic: 2.335 on 1 and 9 DF,  p-value: 0.1609

Model hasil transformasi : IPM = - 51.7880 + 0.2928 Tahun

Pengecekan Asumsi Autokorelasi

  1. Durbin-Watson test

Hipotesis :

H0 = Sisaan Saling Bebas (Tidak Autokorelasi)

H1 = Sisaan Tidak Saling Bebas (Autokorelasi)

lmtest::dwtest(modelhl,alternative = "two.sided")
## 
##  Durbin-Watson test
## 
## data:  modelhl
## DW = 2.1207, p-value = 0.8594
## alternative hypothesis: true autocorrelation is not 0

Didapat P-value > 5%. Tak Tolak H0, sehingga belum cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf 5%.

  1. Run’s Test

Hipotesis :

H0 = Sisaan Saling Bebas (Tidak Autokorelasi)

H1 = Sisaan Tidak Saling Bebas (Autokorelasi)

lawstat::runs.test(resid(modelhl), alternative = 'two.sided')
## 
##  Runs Test - Two sided
## 
## data:  resid(modelhl)
## Standardized Runs Statistic = -0.93314, p-value = 0.3507

Berdasarkan uji Run’s Test, p-value > 0.05 yang artinya Tak Tolak H0. belum cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf 5%.

  1. Breusch-Godfrey test

Hipotesis :

H0 = Sisaan Saling Bebas (Tidak Autokorelasi)

H1 = Sisaan Tidak Saling Bebas (Autokorelasi)

bgtest(modelhl)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  modelhl
## LM test = 0.22421, df = 1, p-value = 0.6358

P-value > 5%. Tak Tolak H0, sehingga belum cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf 5%.

Berdasarkan hasil di atas dapat disimpulkan hasil transformasi Hildreth-Lu memenuhi asumsi bebas autokorelasi.

Transformasi Balik

Persamaan regresi setelah di transformasi ke persamaan awal

cat("IPM = ", coef(modelhl)[1]/(1-0.735), "+", coef(modelhl)[2]," Tahun", sep = "") 
## IPM = -195.4262+0.2928182 Tahun

Kesimpulan

Terdapat autokorelasi pada deret waktu IPM Provinsi Jawa Barat tahun 2010-2021. Hal ini dapat dideteksi setelah melakukan pengujian dengan metode grafik dan uji formal. Sehingga dilakukan penanganan dengan metode Cochrane-Orcutt dan Hildreth-Lu. Kedua metode tersebut dapat menangani autokorelasi pada data IPM Provinsi Jawa Barat tahun 2010-2021.