Packages

Beberapa packages yang digunakan adalah lmtest, orcutt, HoRM, dan readxl

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(orcutt)
## Warning: package 'orcutt' was built under R version 4.1.2
library(HoRM)
## Warning: package 'HoRM' was built under R version 4.1.2
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(readxl)

Data yang digunakan adalah data IPM Provinsi Riau Tahun 2010-2021.

IPM adalah indikator penting untuk mengukur keberhasilan pembangunan kualitas hidup manusia.

Dimensi dasar IPM:

  1. Umur panjang dan sehat

  2. Pengetahuan

  3. Standar hidup layak

sumber : https://www.bps.go.id/subject/26/indeks-pembangunan-manusia.html

#Menentukan Working Direktori
setwd("E:/Semester 6/MPDW/Pertemuan 3/Kuliah")

#Memasukkan File Data
dataregresi<-read_excel("IPM Riau.xlsx")
attach(dataregresi)

Diagram Pencar IPM Provinsi Riau 2010-2021

plot(x,y, xlab="Tahun", ylab="IPM", pch = 20, col = "blue", main = "Scatter Plot Tahun vs IPM Provinsi Riau")

Berdasarkan diagram pencar tersebut, terlihat bahwa nilai IPM Provinsi Riau mengalami peningkatan setiap tahunnya. Namun pada tahun 2020 terjadi penurunan dibandingkan tahun sebelumnya dan pada tahun 2021 kembali mengalami kenaikan walaupun nilainya tidak setinggi pada tahun 2019. Hubungan antara tahun dan nilai IPM cenderung linear.

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

Korelasi antara tahun dengan IPM bernilai 0.9864531. Artinya hubungan antara tahun dan IPM adalah hubungan positif. Hal tesebut juga ditunjukkan pada diagram pencar sebelumnya yang menggambarkan kenaikan nilai IPM setiap tahunnya.

Model Regresi

model<- lm(y~x, data = dataregresi)
summary(model)
## 
## Call:
## lm(formula = y ~ x, data = dataregresi)
## 
## 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-09

Berdasarkan hasil regresi tersebut, didapatkan persamaan regresi yaitu y = -815.69072 + 0.43993x. Parameter regresi yang dihasilkan juga signifikan.

Diagnostik Residual secara Ekplorasi

#Sisaan dan Fitted Value
resi1<- residuals(model)
fit<- predict(model)

#Diagnostik secara Eksploratif
par(mfrow = c(2,2))
qqnorm(resi1)
qqline(resi1, col = "steelblue", lwd = 2)
plot(fit, resi1, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Fitted Values", main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)
hist(resi1, col = "steelblue")
plot(seq(1,12,1), resi1, col = "steelblue", 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)

Pengecekan Diagnostik Residual secara Ekplorasi:

  1. Berdasarkan normal Q-Q Plot dan Histogram, terlihat bahwa sisaan cenderung menyebar normal.

  2. Berdasarkan plot Sisaan vs Fitted Values, terlihat bahwa sisaan tidak acak artinya sisaan belum homogen. Namun hal tersebut perlu pengecekan melalui uji formal.

  3. Berdasarkan Plot sisaan vs Order, terlihat bahwa sisaan tidak menyebar acak atau terdapat pola naik turun. Hal tersebut mengindikasikan adanya autokorelasi.

Identifikasi autokorelasi dengan plot ACF dan PACF

#Mengidentifikasi Autokorelasi dengan ACF dan PACF
par(mfrow = c(1,2))
acf(resi1)
pacf(resi1)

Berdasarkan plot ACF dan PACF diatas, terlihat bahwa selain pada lag 0, tidak terdapat garis vertikal yang melewati batas signifikan. Oleh karena itu dapat disimpulkan bahwa tidak terdapat autokorelasi.

Terdapat perbedaan kesimpulan dari hasil eksplorasi pada sisaan vs order dengan plot ACF dan PACF. Oleh karena itu diperlukan uji formal untuk memastikan ada tidaknya autokorelasi.

Uji Durbin Warson

dwtest(model)
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 1.103, p-value = 0.01635
## alternative hypothesis: true autocorrelation is greater than 0

Hipotesis uji Durbin Watson

H0 : tidak terdapat autokorelasi

H1 : ada autokorelasi positif

Berdasarkan uji Durbin Watson, diperoleh nilai p-value sebesar 0.01635 < alpha 0.05 sehingga tolak H0, artinya cukup bukti untuk menyatakan bahwa data tersebut memiliki autokorelasi positif pada taraf nyata 5%.

Uji Breusch-Godfrey

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

Hipotesis uji Breusch-Godfrey

H0 : ada autokorelasi

H1 : tidak ada autokorelasi

Berdasarkan uji Breusch-Godfrey, diperoleh nilai p-value sebesar 0.1801 > alpha 0.05 sehingga tak tolak H0, artinya tidak cukup bukti untuk menyatakan bahwa data tersebut tidak ada autokorelasi pada taraf nyata 5%.

Hasil uji formal baik pada uji Durbin Watson maupun uji Breusch-Godfrey keduanya menghasilkan kesimpulan bahwa terdapat autokorelasi. Oleh karena itu diperlukan penanganan autokorelasi pada model tersebut.

Penanganan Autokorelasi dengan Cochrane-Orcutt

Penentuan rho terbaik pada metode Cochrane-Orcutt yaitu dengan cara iteratif.

modelco<-cochrane.orcutt(model)
modelco
## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = y ~ x, data = dataregresi)
## 
##  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
summary(modelco)
## Call:
## lm(formula = y ~ x, data = dataregresi)
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept) -784.312909  100.589892  -7.797 2.716e-05 ***
## x              0.424346    0.049872   8.509 1.348e-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-05
## 
## Durbin-Watson statistic 
## (original):    1.10303 , p-value: 1.635e-02
## (transformed): 1.42722 , p-value: 7.532e-02

Berdasarkan hasil diatas, diperoleh nilai rho sebesar 0.486948 dan persamaan yang diperoleh yaitu y = -784.312909 + 0.424346x. Persamaan tersebut adalah persamaan yang sudah ditransformasi balik. Setelah transformasi diperoleh nilai Durbin-Watson statistic sebesar 1.42722 dengan p-value: 7.532e-02. Artinya model hasil penanganan dengan Cochrane-Orcutt tidak ada autokorelasi.

#Cara Manual
#Menentukan rho optimum
rho<- modelco$rho

#Mentransformasi
y.trans<- y[-1]-y[-12]*rho
x.trans<- x[-1]-x[-12]*rho
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-05

Berdasarkan hasil diatas diperoleh persamaan y = -402.39344 + 0.42435x. Namun hasil tersebut belum ditransformasi balik sehingga harus ditransformasi balik terlebih dahulu.

#Transformasi balik
cat("y = ", coef(modelcorho)[1]/(1-rho), "+", coef(modelcorho)[2],"x", sep = "")
## y = -784.3129+0.4243458x

Setelah ditansformasi balik diperoleh persamaan y = -784.3129 + 0.4243458x. Hasil tersebut sama seperti hasil Cochrane-Orcutt sebelumnya yang menggunakan fungsi langsung.

Penanganan Autokorelasi dengan Hildreth lu

Penentuan rho terbaik pada metode Hildreth lu yaitu dengan memilih nilai rho yang menghasilkan SSE terkecil.

# 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
r <- c(seq(0.1,0.9, by= 0.1))
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4)
##   rho    SSE
## 1 0.1 0.7162
## 2 0.2 0.6856
## 3 0.3 0.6640
## 4 0.4 0.6516
## 5 0.5 0.6482
## 6 0.6 0.6540
## 7 0.7 0.6688
## 8 0.8 0.6927
## 9 0.9 0.7257
#Membuat Grafik rho dan SSE
plot(tab$SSE ~ tab$rho , type = "l")
abline(v = tab[tab$SSE==min(tab$SSE),"rho"], lty = 3)

Berdasarkan hasil diatas nilai SSE terkecil terdapat pada rho 0.5. Selanjutnya untuk lebih rinci akan dilakukan pemilihan rho terbaik disekitar 0.4 dan 0.5 dengan jarak 0.01.

#Rho Optimal di sekitar 0.4 dan 0.5
r <- seq(0.45,0.55, 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.45 0.6488
## 2  0.46 0.6485
## 3  0.47 0.6483
## 4  0.48 0.6482
## 5  0.49 0.6482
## 6  0.50 0.6482
## 7  0.51 0.6484
## 8  0.52 0.6486
## 9  0.53 0.6490
## 10 0.54 0.6494
## 11 0.55 0.6500

Berdasarkan hasil diatas, terlihat bahwa rho terbaik adalah 0.48 karena memiliki nilai SSE terkecil.

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

Nilai rho tersebut masih dapat diperinci lagi dengan memperkecil increment yang digunakan.

# Model Terbaik
modelhl <- hildreth.lu.func(0.48, model)
summary(modelhl)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42545 -0.20442  0.00342  0.17527  0.35429 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -408.41831   51.60807  -7.914 2.41e-05 ***
## x              0.42490    0.04921   8.635 1.20e-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.8923, Adjusted R-squared:  0.8803 
## F-statistic: 74.56 on 1 and 9 DF,  p-value: 1.197e-05

Persamaan yang dihasilkan dengan menggunakan rho 0.48 adalah y = -408.41831 + 0.42490x. Hasil tersebut belum ditransformasi balik.

#Transformasi Balik
cat("y = ", coef(modelhl)[1]/(1-0.48), "+", coef(modelhl)[2],"x", sep = "")
## y = -785.4198+0.4248951x

Setelah ditransformasi balik diperoleh persamaan y = -785.4198+0.4248951x.

# Deteksi autokorelasi
dwtest(modelhl)
## 
##  Durbin-Watson test
## 
## data:  modelhl
## DW = 1.4204, p-value = 0.07346
## alternative hypothesis: true autocorrelation is greater than 0

Setelah dilakukan penanganan dengan Hildreth lu, diperoleh p-value untuk uji Durbin Watson sebesar 0.07346 > alpha 0.05 sehingga tak tolak H0. Artinya tidak cukup bukti untuk menyatakan bahwa ada autokorelasi pada taraf nyata 5%

Selain dengan cara fungsi yang telah dibentuk, penanganan autokorelasi dengan Hildreth lu juga dapat menggunakan fungsi “hildreth.lu” pada package HoRM.

#mengidentifikasi vektor x dan y
y<-dataregresi$y
x<-dataregresi$x

#fungsi hildreth lu dengan rho tertentu
modelhl2<-hildreth.lu(y,x,0.48)
deviance(modelhl2)
## [1] 0.6481727
summary(modelhl2)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42545 -0.20442  0.00342  0.17527  0.35429 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -408.41831   51.60807  -7.914 2.41e-05 ***
## x              0.42490    0.04921   8.635 1.20e-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.8923, Adjusted R-squared:  0.8803 
## F-statistic: 74.56 on 1 and 9 DF,  p-value: 1.197e-05
#Transformasi Balik
cat("y = ", coef(modelhl2)[1]/(1-0.48), "+", coef(modelhl2)[2],"x", sep = "")
## y = -785.4198+0.4248951x

Hasil di atas sama seperti hasil Hildreth lu dengan fungsi yang dibuat sendiri.

Simpulan

Hasil regresi IPM Provinsi Riau tahun 2010-2021 menunjukkan adanya autokorelasi sehingga diperlukan penanganan. Penanganan autokorelasi dilakukan dengan metode Cochare-orcutt dan Hildreth lu. Kedua metode tersebut memiliki prinsip yang sama yaitu menentukan rho optimum terlebih dahulu. Namun cara memperoleh nilai rho-nya yang berbeda. Hasil dari penanganan dengan kedua metode tersebut relatif sama. Persamaan yang dihasilkan setelah dilakukan penanganan adalah IPM = -784.312909 + 0.424346(tahun). Model hasil penanganan juga telah terbebas dari masalah autokorelasi.