Penanganan Autokorelasi dengan Metode Cochrane-Orcutt dan Hildreth-Lu

Dinda Khamila Nurfatimah

2023-09-11

Pemanggilan Packages

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, union

library(TTR)
library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(lmtest) #digunakan untuk uji formal pendeteksian autokorelasi

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(orcutt) #untuk membuat model regresi Cochrane-Orcutt
library(HoRM) #untuk membuat model regresi Hildreth-Lu

Input Data

Data yang digunakan dalam kesempatan kali ini adalah data New Delhi Air Quality (Kualitas Udara di New Delhi) tahun 2022.

library(rio)
data <- import("https://raw.githubusercontent.com/DindaKhamila/mpdw/main/Pertemuan%202/data2.csv")
View(data)

Eksplorasi Data

str(data)     #struktur data

## 'data.frame':    72 obs. of  7 variables:
##  $ AQI : num  30.2 28.2 26.6 25 26 26 26 24 23 24 ...
##  $ CO  : num  199 198 199 202 205 ...
##  $ no2 : num  0.0469 0.0465 0.0469 0.0482 0.0489 ...
##  $ o3  : num  55.8 54.9 54.6 55.1 55.8 ...
##  $ pm10: num  10.5 10.7 11.2 11.1 10.4 ...
##  $ pm25: num  5.64 4.62 3.52 2.23 1.98 ...
##  $ so2 : num  0.387 0.41 0.402 0.376 0.339 ...

dim(data)     #dimensi data

## [1] 72  7

summary(data$AQI)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   19.00   25.00   27.00   26.18   28.00   32.00

Melakukan eksplorasi data dengan membuat matriks korelasi dan histogram antar peubah untuk mengetahui pola sebaran data dan hubungan antar peubah

#pair plot
library(ggplot2)
library(GGally)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

ggpairs(data)

Gambar di atas menunjukkan matriks korelasi dan scatter plot yang memuat hubungan linier antara peubah respon dan peubah penjelas serta nilai koefisien korelasi. Tanda bintang (*) pada koefisien korelasi menunjukkan adanya hubungan linier dan hubungan signifikan antara peubah respon dengan peubah penjelas. Berdasarkan matriks korelasi tersebut diketahui bahwa peubah PM10 tidak memiliki hubungan linier dan signifikan terhadap peubah respon AQI. Hal tersebut dapat diartikan bahwa pada tahun 2022 kualitas udara di New Delhi tidak memiliki pengaruh pada peubah PM10. Oleh karena itu, pada penelitian ini, peubah PM10 dapat terwakili oleh peubah-peubah lain yang signifikan.

# Memasukkan pustaka ggplot2 (pastikan telah menginstalnya terlebih dahulu)
library(ggplot2)

# Membuat histogram terhadap seluruh peubah dengan ggplot2
ggplot(data, aes(x = AQI)) +
  geom_histogram(binwidth = 1, fill = "blue", color = "black", alpha=0.7) +
  labs(title = "Histogram",
       x = "AQI",
       y = "Frekuensi")

ggplot(data, aes(x = CO)) +
  geom_histogram(binwidth = 1, fill = "blue", color = "black", alpha=0.7) +
  labs(title = "Histogram",
       x = "CO",
       y = "Frekuensi")

ggplot(data, aes(x = no2)) +
  geom_histogram(binwidth = 1, fill = "blue", color = "black", alpha=0.7) +
  labs(title = "Histogram",
       x = "NO2",
       y = "Frekuensi")

ggplot(data, aes(x = o3)) +
  geom_histogram(binwidth = 1, fill = "blue", color = "black", alpha=0.7) +
  labs(title = "Histogram",
       x = "O3",
       y = "Frekuensi")

ggplot(data, aes(x = pm25)) +
  geom_histogram(binwidth = 1, fill = "blue", color = "black", alpha=0.7) +
  labs(title = "Histogram",
       x = "PM25",
       y = "Frekuensi")

ggplot(data, aes(x = so2)) +
  geom_histogram(binwidth = 1, fill = "blue", color = "black", alpha=0.7) +
  labs(title = "Histogram",
       x = "SO2",
       y = "Frekuensi")

Menamakan Peubah

y <- data[,1] #y untuk peubah AQI
x1 <- data[,2] #x1 untuk peubah CO
x2 <- data[,3] #x2 untuk peubah NO2
x3 <- data[,4] #x3 untuk peubah O3
x4 <- data[,6] #x4 untuk peubah PM25
x5 <- data[,7] #x5 untuk peubah SO2

data <- cbind(y, x1, x2, x3, x4,x5)
data <- as.data.frame(data)

Regresi

#Pembuatan Model Regresi
#model regresi
model<- lm(y~x1+x2+x3+x4+x5, data = data)
summary(model)

## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.00865 -0.19031  0.05126  0.24060  0.48683 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.98852    3.47797  -1.722   0.0898 .  
## x1           0.02813    0.01721   1.634   0.1070    
## x2           9.37637    5.24174   1.789   0.0782 .  
## x3           0.45991    0.01284  35.820  < 2e-16 ***
## x4           1.25523    0.08068  15.558  < 2e-16 ***
## x5          -7.57517    1.71003  -4.430 3.63e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3258 on 66 degrees of freedom
## Multiple R-squared:  0.9894, Adjusted R-squared:  0.9886 
## F-statistic:  1235 on 5 and 66 DF,  p-value: < 2.2e-16

Model yang dihasilkan adalah \[y=-5.98852+0.02813x1+9.37637x2+0.45991x3+1.25523x4-7.57517x5\] Berdasarkan ringkasan model dapat diketahui bahwa hasil uji F memiliki p-value(< 2.2e-16) < \(\alpha\) (5%). Artinya, minimal terdapat satu variabel yang berpengaruh nyata terhadap model. Hasil uji-t parsial kedua parameter regresi, yaitu intersep dan koefisien regresi juga menunjukkan hal yang sama, yaitu memiliki p-value < \(\alpha\) (5%) sehingga nyata dalam taraf 5%, kecuali peubah x1 yang memiliki p-value (0.1070) > \(\alpha\) (5%). Selanjutnya dapat dilihat juga nilai \(R^2=0.9894\). Artinya, sebesar 98.94% keragaman nilai AQI dapat dijelaskan oleh semua peubah penjelas. Hasil ini menunjukkan hasil yang bagus, seolah mendapatkan hasil terbaik. Namun, kita perlu melakukan uji terhadap sisaannya seperti berikut ini.

#sisaan dan fitted value
 sisaan<- residuals(model)
fitValue<- predict(model)

#Diagnostik dengan eksploratif
par(mfrow = c(2,2))
qqnorm(sisaan)
qqline(sisaan, col = "steelblue", lwd = 2)
plot(fitValue, sisaan, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Fitted Values", main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)
hist(sisaan, col = "steelblue")
plot(seq(1,72,1), sisaan, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,72,1), sisaan, col = "red")
abline(a = 0, b = 0, lwd = 2)

Dua plot di samping kiri digunakan untuk melihat apakah sisaan menyebar normal. Normal Q-Q Plot di atas menunjukkan bahwa sisaan cenderung menyebar normal, tetapi histogram dari sisaan tidak menunjukkan demikian. Selanjutnya, dua plot di samping kanan digunakan untuk melihat autokorelasi. Plot Sisaan vs Fitted Value dan Plot Sisaan vs Order menunjukkan adanya pola pada sisaan. Untuk lebih lanjut akan digunakan uji formal melihat normalitas sisaan dan plot ACF dan PACF untuk melihat apakah ada autokorelasi atau tidak.

#Melihat Sisaan Menyebar Normal/Tidak
#H0: sisaan mengikuti sebaran normal
#H1: sisaan tidak mengikuti sebaran normal
shapiro.test(sisaan)

## 
##  Shapiro-Wilk normality test
## 
## data:  sisaan
## W = 0.96301, p-value = 0.03273

ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan))

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.092866, p-value = 0.5331
## alternative hypothesis: two-sided

Berdasarkan uji formal Saphiro-Wilk didapatkan nilai p-value(0.03273)< \(\alpha\) (5%) dan Kolmogorov-Smirnov didapatkan nilai p-value(0.5331) > \(\alpha\) (5%).

#ACF dan PACF identifikasi autokorelasi
par(mfrow = c(1,2))
acf(sisaan)
pacf(sisaan)

Berdasarkan plot ACF dan PACF, terlihat semua dalam rentang batas dan tidak ada yang signifikan. Namun, untuk lebih memastikan akan dilakukan uji formal dengan uji Durbin Watson.

#Deteksi autokorelasi dengan uji-Durbin Watson
#H0: tidak ada autokorelasi
#H1: ada autokorelasi
dwtest(model)

## 
##  Durbin-Watson test
## 
## data:  model
## DW = 1.6803, p-value = 0.02539
## alternative hypothesis: true autocorrelation is greater than 0

Berdasarkan hasil DW Test, didapatkan nilai \(DW = 1.6803\) dan p-value = \(0.02539\). Berdasarkan tabel Durbin-Watson diperoleh nilai \(DL = 1.4732\) dan \(DU = 1.7688\). Nilai DW masih berada di antara nilai DL dan DU. Artinya, berada di daerah inkonklusif, tidak dapat dikatakan berada di daerah autokorelasi positif maupun bebas dari autokorelasi. Namun, dengan nilai p-value < 0.05 dapat disimpulkan bahwa tolak H0, cukup bukti mengatakan adanya autokorelasi. Oleh karena itu, diperlukan penangan autokorelasi. Penanganan yang akan digunakan menggunakan dua metode, yaitu Cochrane-Orcutt dan Hildret-Lu.

Penanganan Autokorelasi

Metode Cochrane-Orcutt

Penanganan metode Cochrane-Orcutt dapat dilakukan dengan bantuan packages Orcutt pada aplikasi R maupun secara manual. Berikut ini ditampilkan cara menggunakan bantuan library packages Orcutt.

#Penanganan Autokorelasi Cochrane-Orcutt
modelCO<-cochrane.orcutt(model)
modelCO

## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5, data = data)
## 
##  number of interaction: 11
##  rho 0.147689
## 
## Durbin-Watson statistic 
## (original):    1.68029 , p-value: 2.539e-02
## (transformed): 2.13850 , p-value: 5.042e-01
##  
##  coefficients: 
## (Intercept)          x1          x2          x3          x4          x5 
##   -4.224160    0.018503    8.227037    0.462981    1.024578   -6.070250

Hasil keluaran model setelah dilakukan penanganan adalah sebagai berikut. \[y=-4.224160+0.018503x1+8.227037x2+0.462981x3+1.024578x4-6.070250x5\] Hasil juga menunjukkan bahwa nilai DW dan p-value meningkat menjadi \(2.13850\) dan \(0.5042\). Nilai DW sudah berada pada rentang DU < DW < 4-DU atau \(1.7688 < DW < 2.2312\). Hal tersebut juga didukung dengan nilai p-value(0.5042) > 0.05, artinya belum cukup bukti menyatakan bahwa sisaan terdapat autokorelasi pada taraf nyata 5%. Untuk nilai \(ρ ̂\) optimum yang digunakan adalah \(0.1476895\). Nilai tersebut dapat diketahui dengan syntax berikut.

#Rho optimum
rho<- modelCO$rho
rho

## [1] 0.1476895

Selanjutnya akan dilakukan transformasi secara manual dengan syntax berikut ini.

#Transformasi Manual
y.trans <- y[-1]-y[-72]*rho
x1.trans <- x1[-1]-x1[-72]*rho
x2.trans <- x2[-1]-x2[-72]*rho
x3.trans <- x3[-1]-x3[-72]*rho
x4.trans <- x4[-1]-x4[-72]*rho
x5.trans <- x5[-1]-x5[-72]*rho
modelCOmanual <- lm(y.trans ~ x1.trans + x2.trans + x3.trans + x4.trans + x5.trans)
summary(modelCOmanual)

## 
## Call:
## lm(formula = y.trans ~ x1.trans + x2.trans + x3.trans + x4.trans + 
##     x5.trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.83902 -0.25457  0.03448  0.25319  0.54531 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.60030    3.34238  -1.077   0.2854    
## x1.trans     0.01850    0.01939   0.954   0.3436    
## x2.trans     8.22704    5.88076   1.399   0.1666    
## x3.trans     0.46298    0.01434  32.277  < 2e-16 ***
## x4.trans     1.02458    0.15577   6.578 9.61e-09 ***
## x5.trans    -6.07025    2.12271  -2.860   0.0057 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3171 on 65 degrees of freedom
## Multiple R-squared:  0.9863, Adjusted R-squared:  0.9853 
## F-statistic: 937.7 on 5 and 65 DF,  p-value: < 2.2e-16

Hasil model transformasi bukan merupakan model sesungguhnya. Koefisien regresi masih perlu dicari kembali mengikuti \(β_0^*=β_0+ρ ̂β_0\) dan \(β_1^*=β_1\).

#Mencari Penduga Koefisien Regresi setelah Transformasi ke Persamaan Awal
b0bintang <- modelCOmanual$coefficients[-2]
b0 <- b0bintang/(1-rho)
b1 <- modelCOmanual$coefficients[-1]
b0

## (Intercept)    x2.trans    x3.trans    x4.trans    x5.trans 
##  -4.2241601   9.6526281   0.5432069   1.2021180  -7.1221102

b1

##    x1.trans    x2.trans    x3.trans    x4.trans    x5.trans 
##  0.01850297  8.22703665  0.46298100  1.02457782 -6.07024960

Hasil perhitungan koefisien regresi tersebut akan menghasilkan hasil yang sama dengan model yang dihasilkan menggunakan packages.

Metode Hildreth-Lu

Penanganan kedua adalah menggunakan metode Hildreth-Lu. Metode ini akan mencari nilai SSE terkecil dan dapat dicari secara manual maupun menggunakan packages. Jika menggunakan packages, gunakan library packages HORM.

#Penanganan Autokorelasi Hildreth lu
# Hildreth-Lu
hildreth.lu.func<- function(r, model){
  x1 <- model.matrix(model)[,2]
  x2 <- model.matrix(model)[,3]
  x3 <- model.matrix(model)[,4]
  x4 <- model.matrix(model)[,5]
  x5 <- model.matrix(model)[,6]
  
  y <- model.response(model.frame(model))
  n <- length(y)
  t <- 2:n
  y <- y[t]-r*y[t-1]
  x1 <- x1[t]-r*x1[t-1]
  x2 <- x2[t]-r*x2[t-1]
  x3 <- x3[t]-r*x3[t-1]
  x4 <- x4[t]-r*x4[t-1]
  x5 <- x5[t]-r*x5[t-1]
  
  return(lm(y~x1+x2+x3+x4+x5))
}

#Pencariab 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 6.5505
## 2 0.2 6.5532
## 3 0.3 6.6750
## 4 0.4 6.9165
## 5 0.5 7.2891
## 6 0.6 7.7976
## 7 0.7 8.3710
## 8 0.8 8.9575
## 9 0.9 9.6240

Pertama-tama akan dicari di mana kira-kira \(ρ\) yang menghasilkan SSE minimum. Pada hasil di atas terlihat \(ρ\) minimum ketika 0.1. Namun, hasil tersebut masih kurang teliti sehingga akan dicari kembali \(ρ\) yang lebih optimum dengan ketelitian yang lebih. Jika sebelumnya jarak antar \(ρ\) yang dicari adalah 0.1, kali ini jarak antar \(ρ\) adalah 0.001 dan dilakukan pada selang 0.2 sampai dengan 0.9.

#Rho optimal di sekitar 0.1
rOpt <- seq(0.1,0.9, by= 0.001)
tabOpt <- data.frame("rho" = rOpt, "SSE" = sapply(rOpt, function(i){deviance(hildreth.lu.func(i, model))}))
head(tabOpt[order(tabOpt$SSE),])

##      rho      SSE
## 49 0.148 6.536743
## 48 0.147 6.536745
## 50 0.149 6.536752
## 47 0.146 6.536759
## 51 0.150 6.536774
## 46 0.145 6.536786

#Grafik SSE optimum
par(mfrow = c(1,1))
plot(tab$SSE ~ tab$rho , type = "l", xlab = "Rho", ylab = "SSE")
abline(v = tabOpt[tabOpt$SSE==min(tabOpt$SSE),"rho"], lty = 2, col="red",lwd=2)
text(x=0.1476895, y=6.5505, labels = "rho=0.1476895", cex = 0.8)

Perhitungan yang dilakukan aplikasi R menunjukkan bahwa nilai \(ρ\) optimum, yaitu saat SSE terkecil terdapat pada nilai \(ρ=0.1476895\). Hal tersebut juga ditunjukkan pada plot. Selanjutnya, model dapat didapatkan dengan mengevaluasi nilai \(ρ\) ke dalam fungsi hildreth.lu.func, serta dilanjutkan dengan pengujian autokorelasi dengan uji Durbin-Watson. Namun, setelah pengecekan tersebut tidak lupa koefisien regresi tersebut digunakan untuk transformasi balik. Persamaan hasil transformasi itulah yang menjadi persamaan sesungguhnya.

#Model terbaik
modelHL <- hildreth.lu.func(0.1476895, model)
summary(modelHL)

## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.83902 -0.25457  0.03448  0.25319  0.54531 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.60030    3.34238  -1.077   0.2854    
## x1           0.01850    0.01939   0.954   0.3436    
## x2           8.22704    5.88076   1.399   0.1666    
## x3           0.46298    0.01434  32.277  < 2e-16 ***
## x4           1.02458    0.15577   6.578 9.61e-09 ***
## x5          -6.07025    2.12271  -2.860   0.0057 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3171 on 65 degrees of freedom
## Multiple R-squared:  0.9863, Adjusted R-squared:  0.9853 
## F-statistic: 937.7 on 5 and 65 DF,  p-value: < 2.2e-16

Setelah dilakukan tranformasi balik, didapatkan model dengan metode Hildreth-Lu sebagai berikut. \[y=-3.60030+0.01850x1+8.22704x2+0.46298x3+1.02458x4-6.07025x5\]

#Deteksi autokorelasi
dwtest(modelHL)

## 
##  Durbin-Watson test
## 
## data:  modelHL
## DW = 2.1385, p-value = 0.5042
## alternative hypothesis: true autocorrelation is greater than 0

Hasil uji Durbin-Watson juga menunjukkan bawah nilai DW sebesar \(2.1385\) berada pada selang daerah tidak ada autokorelasi, yaitu pada rentang DU < DW < 4-DU atau \(1.7688 < DW < 2.2312\). Hal tersebut juga didukung oleh p-value sebesar \(0.5042\), di mana p-value > \(\alpha\)=5%. Artinya tak tolak \(H_0\) atau belum cukup bukti menyatakan bahwa ada autokorelasi dalam data nilai AQI dengan metode Hildreth-Lu pada taraf nyata 5%.

Terakhir, akan dibandingkan nilai SSE dari ketiga metode (metode awal, metode Cochrane-Orcutt, dan Hildreth-Lu).

#Perbandingan
sseModelawal <- anova(model)$`Sum Sq`[-c(1,2,3,4,5)]
sseModelCO <- anova(modelCOmanual)$`Sum Sq`[-c(1,2,3,4,5)]
sseModelHL <- anova(modelHL)$`Sum Sq`[-c(1,2,3,4,5)]
mseModelawal <- sseModelawal/length(y)
mseModelCO <- sseModelCO/length(y)
mseModelHL <- sseModelHL/length(y)
akurasi <- matrix(c(sseModelawal,sseModelCO,sseModelHL,
                    mseModelawal,mseModelCO,mseModelHL),nrow=2,ncol=3,byrow = T)
colnames(akurasi) <- c("Model Awal", "Model Cochrane-Orcutt", "Model Hildreth-Lu")
row.names(akurasi) <- c("SSE","MSE")
akurasi

##     Model Awal Model Cochrane-Orcutt Model Hildreth-Lu
## SSE 7.00358086            6.53674213        6.53674213
## MSE 0.09727196            0.09078809        0.09078809

Berdasarkan hasil tersebut dapat diketahui bahwa hasil penanganan autokorelasi dengan metode Cochrane-Orcutt dan Hildreth-Lu memiliki SSE yang sama, sebesar \(6.53674213\) dan lebih baik dibandingkan model awal ketika autokorelasi masih terjadi, yaitu sebesar \(7.00358086\).

Simpulan

Autokorelasi yang terdapat pada data AQI terjadi akibat adanya korelasi di antara unsur penyusunnya. Adanya autokorelasi menyebabkan model regresi kurang baik karena akan meingkatkan galatnya. Autokorelasi dapat dideteksi secara eksploratif melalui plot sisaan, ACF, dan PACF, serta dengan uji formal Durbin-Watson. Namun, autokorelasi tersebut dapat ditangani dengan metode Cochrane-Orcutt dan Hildreth-Lu. Kedua metode menghasilkan nilai SSE yang sama, artinya keduanya baik untuk digunakan.