Regresi dengan Penanganan Autokorelasi

Pada pertemuan kali ini, analisis data deret waktu akan lebih membahas tentang regresi. Dalam regresi terdapat asumsi tidak adanya autokorelasi yang harus terpenuhi namun kerap kali autokorelasi muncul pada data deret waktu sehingga kali ini akan dilakukan penanganan autokorelasi pada model regresi sebagai bahan bahasan utamanya.

Adapun, pengerjaan dimulai dengan panggil packages library untuk analisis data deret waktu dan regresi yang diperlukan sebagai berikut.

library(dplyr)
library(TTR)
library(forecast)
library(lmtest) #digunakan untuk uji formal pendeteksian autokorelasi
library(orcutt) #untuk membuat model regresi Cochrane-Orcutt
library(HoRM) #untuk membuat model regresi Hildreth-Lu

Preprocessing

Sebelum masuk ke analisis akan ada tahapan sebagai berikut.

Input Data

Data yang digunakan dalam kesempatan kali ini adalah data emisi karbon co2 yang berasal dari link berikut “https://www.kaggle.com/datasets/ggsri123/co2-emissions-from-fossil-fuels” dan sudah dirapihkan terlebih dahulu dengan pemilihan negara amatan yaitu di amerika serikat (USA).

#install.packages("rio") #install jika belum ada
library(rio)
data<- import("https://raw.githubusercontent.com/fax17/MPDW/main/pertemuan%202/data%20fossil.csv")
head(data) #hanya menampilkan sebagian data paling awal

##   Year                  Country  Total Solid Fuel Liquid Fuel Gas Fuel Cement
## 1 1950 UNITED STATES OF AMERICA 692124     343131      244760    87123   5341
## 2 1951 UNITED STATES OF AMERICA 712983     330676      262189   102677   5786
## 3 1952 UNITED STATES OF AMERICA 694618     293170      273190   109931   5860
## 4 1953 UNITED STATES OF AMERICA 711175     290927      286594   115541   6209
## 5 1954 UNITED STATES OF AMERICA 677689     249303      290178   121179   6398
## 6 1955 UNITED STATES OF AMERICA 742828     280082      313292   130784   7304
##   Gas Flaring Per Capita Bunker fuels (Not in Total)
## 1       11769       4.32                        8100
## 2       11654       4.39                       10253
## 3       12468       4.21                       10784
## 4       11905       4.24                       10973
## 5       10631       3.97                       10060
## 6       11366       4.28                       10842

data_fossil <- data[c(-1,-2)] #hanya mengambil peubah numeriknya
head(data_fossil)

##    Total Solid Fuel Liquid Fuel Gas Fuel Cement Gas Flaring Per Capita
## 1 692124     343131      244760    87123   5341       11769       4.32
## 2 712983     330676      262189   102677   5786       11654       4.39
## 3 694618     293170      273190   109931   5860       12468       4.21
## 4 711175     290927      286594   115541   6209       11905       4.24
## 5 677689     249303      290178   121179   6398       10631       3.97
## 6 742828     280082      313292   130784   7304       11366       4.28
##   Bunker fuels (Not in Total)
## 1                        8100
## 2                       10253
## 3                       10784
## 4                       10973
## 5                       10060
## 6                       10842

Eksplorasi Data

Selanjutnya dari data yang diperoleh akan dilakukan eksplorasi dan visualisasinya menggunakan matriks korelasi untuk melihat ada tidaknya hubungan antar peubahnya.

library(corrplot)

## corrplot 0.92 loaded

matriks <- cor(data_fossil)
corrplot.mixed(matriks, upper = 'ellipse', lower = 'number', order = "original",
               tl.col="black", tl.pos = "lt",diag = 'l',
               number.digits=2, number.cex=0.55)

Hasil matriks korelasi menunjukan bahwa tiap peubah saling memiliki hubungan sehingga setalah mengetahui adanya hubungan antar peubah, maka model regresi dapat ditentukan. Pada percobaan Kali ini peneliti dengan sengaja ingin mengetahui hubungan linear berganda antar peubah bunker fuel sebagai respon (y) dan peubah yang mempengaruhinya yakni gass fuel (x1) dan gass flaring (x2). Bunker fuel sendiri merupakan data emisi karbin akibat bahan bakar bunker yang mana emisi kerap terjadi akibar pembakaran kapalnya sehingga ada co2 diudara, sementara itu gass fuel merupakan emisi karbon co2 hasil bahan bakar gas, dan gas flaring adalah emisi karbon co2 dari pembakaran gas non bahan bakar. Jadi pemilihan ini ingin mengetahui apakah emisi karbon co2 pembakaran bahan bakar bunker ini berhubungan atau bahkan dipengaruhi emisi karbon co2 akibat gas baik bahan bakar maupun pembakaran lainnya. Hal ini dicobakan karena peneliti memikirkan apakah jika emisi karbon yang terjadi di udara akibat gas meningkat berarti akibat bahan bakar bunkerjuga meningkat berhubungan akibat bahan bakar bunker ini juga sama saama akibat pembakaran yang hanya saja ini pembakaran yang biasa ada pada kapal. Lebih lanjut akan dilakukan pengujian regresi dan penanganannya.

Regresi

Pemodelan Regresi Linear Berganda

#Pembuatan Model Regresi
#model regresi
model<- lm(`Bunker fuels (Not in Total)`~`Gas Fuel`+`Gas Flaring`, data = data_fossil)
summary(model)

## 
## Call:
## lm(formula = `Bunker fuels (Not in Total)` ~ `Gas Fuel` + `Gas Flaring`, 
##     data = data_fossil)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16428.3  -5040.7    135.8   5966.4  12232.1 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -17.18989 6754.85655  -0.003    0.998    
## `Gas Fuel`       0.10182    0.01971   5.167 2.68e-06 ***
## `Gas Flaring`   -0.54377    0.42440  -1.281    0.205    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7385 on 62 degrees of freedom
## Multiple R-squared:  0.6124, Adjusted R-squared:  0.5999 
## F-statistic: 48.99 on 2 and 62 DF,  p-value: 1.731e-13

cat("y = ", coef(model)[1], " + ", coef(model)[2]," x1",coef(model)[3]," x2", sep = "")

## y = -17.18989 + 0.1018224 x1-0.5437749 x2

Model yang dihasilkan adalah \[y_i=-17.18989+0.1018224x_1-0.5437749x_2\] Berdasarkan ringkasan model dapat diketahui bahwa hasil uji F memiliki p-value < $\alpha$ (5%). Artinya, minimal terdapat satu variabel yang berpengaruh nyata terhadap model. Selanjutnya dapat dilihat juga nilai $R^2=0.6124$. Artinya, sebesar 61.24% keragaman nilai y dapat dijelaskan oleh peubah x nya. Hasil ini perlu kita uji terhadap sisaannya seperti berikut ini.

Pengujian Asumsi

#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,65,1), sisaan, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,65,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
ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan))

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

Berdasarkan uji formal kolmogorov-Smirnov didapatkan nilai p-value > $\alpha$ (5%). Artinya, cukup bukti untuk menyatakan sisaan berdistribusi normal.

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

Berdasarkan plot ACF dan PACF, terlihat seperti ada autokorelasi. 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 = 0.16641, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

Berdasarkan hasil DW Test, didapatkan nilai $DW = 0.16641$ dan p-value < $2.2e-16$. Berdasarkan tabel Durbin-Watson diperoleh nilai $DL = 1.5355$ dan $DU = 1.6621$. Nilai DW berada dibawah nilai DU. Artinya, berada di daerah autokorelasi positif. Selanjutnya, 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)

## Warning in lmtest::dwtest(regCO): exact p value cannot be computed (not in
## [0,1]), approximate p value will be used

## Warning in lmtest::dwtest(regCO): exact p value cannot be computed (not in
## [0,1]), approximate p value will be used

modelCO

## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = `Bunker fuels (Not in Total)` ~ `Gas Fuel` + `Gas Flaring`, 
##     data = data_fossil)
## 
##  number of interaction: 8
##  rho 0.935292
## 
## Durbin-Watson statistic 
## (original):    0.16641 , p-value: 4.396e-30
## (transformed): 1.88284 , p-value: 2.878e-01
##  
##  coefficients: 
##   (Intercept)    `Gas Fuel` `Gas Flaring` 
##  -2779.621296      0.095645      0.186458

Hasil keluaran model setelah dilakukan penanganan adalah sebagai berikut. \[y_i=-2779.621296+0.0956455x_1+0.186458x_2\] Hasil juga menunjukkan bahwa nilai DW dan p-value meningkat menjadi $1.88284$ dan $0.2878$. Nilai DW sudah berada pada rentang DU < DW < 4-DU. Hal tersebut juga didukung dengan nilai p-value > 0.05, artinya belum cukup bukti menyatakan bahwa sisaan terdapat autokorelasi pada taraf nyata 5%. Untuk nilai $ρ ̂$ optimum yang digunakan adalah $0.935292$. Nilai tersebut dapat diketahui dengan syntax berikut.

#Rho optimum
rho<- modelCO$rho
rho

## [1] 0.9352923

Selanjutnya akan dilakukan transformasi atau penanganan autokorelasi secara manual dengan syntax berikut ini.

#Transformasi Manual
Bunker.trans<- data_fossil$`Bunker fuels (Not in Total)`[-1]-data_fossil$`Bunker fuels (Not in Total)`[-65]*rho
Gasfuel.trans<- data_fossil$`Gas Fuel`[-1]-data_fossil$`Gas Fuel`[-65]*rho
Gasflaring.trans <- data_fossil$`Gas Flaring`[-1]-data_fossil$`Gas Flaring`[-65]*rho
modelCOmanual<- lm(Bunker.trans~Gasfuel.trans+Gasflaring.trans)
summary(modelCOmanual)

## 
## Call:
## lm(formula = Bunker.trans ~ Gasfuel.trans + Gasflaring.trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11773.7  -1196.1   -221.5   1232.2   9006.3 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)  
## (Intercept)      -179.86291  853.14990  -0.211   0.8337  
## Gasfuel.trans       0.09565    0.03657   2.615   0.0112 *
## Gasflaring.trans    0.18646    0.49238   0.379   0.7062  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2920 on 61 degrees of freedom
## Multiple R-squared:  0.1117, Adjusted R-squared:  0.08261 
## F-statistic: 3.836 on 2 and 61 DF,  p-value: 0.02695

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[-c(2,3)]
b0 <- b0bintang/(1-rho)
b1 <- modelCOmanual$coefficients[-c(1,3)]
b2 <- modelCOmanual$coefficients[-c(1,2)]

cat("y = ", b0, " + ", b1," x1", " + ",b2," x2", sep = "")

## y = -2779.621 + 0.0956452 x1 + 0.1864576 x2

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.multi<- 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]
  
  # Modifikasi untuk beberapa variabel bebas
  x1 <- x[,1][t]-r*x[,1][t-1]
  x2 <- x[,2][t]-r*x[,2][t-1]
  
  return(lm(y~x1+x2))
}

# Pencarian rho yang meminimumkan SSE
rho <- c(seq(0.1,1, by= 0.1))
tab <- data.frame("rho" = rho, "SSE" = sapply(rho, function(i){deviance(hildreth.lu.func.multi(i, model))}))
round(tab,1)

##    rho        SSE
## 1  0.1 2767452693
## 2  0.2 2255583591
## 3  0.3 1810398485
## 4  0.4 1431849108
## 5  0.5 1119805136
## 6  0.6  873970971
## 7  0.7  693747371
## 8  0.8  577976521
## 9  0.9  524285535
## 10 1.0  536178032

Pertama-tama akan dicari di mana kira-kira $ρ$ yang menghasilkan SSE minimum. Pada hasil di atas terlihat $ρ$ minimum ketika 0.9. 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.00001 dan dilakukan pada selang 0.9 sampai dengan 1.

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

##          rho       SSE
## 3530 0.93529 520050663
## 3531 0.93530 520050663
## 3529 0.93528 520050663
## 3532 0.93531 520050664
## 3528 0.93527 520050664
## 3533 0.93532 520050665

#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.93529, y=520050663, labels = "rho=0.93529", cex = 0.99)

Perhitungan yang dilakukan aplikasi R menunjukkan bahwa nilai $ρ$ optimum, yaitu saat SSE terkecil terdapat pada nilai $ρ=0.93529$. 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 terbaaaik
modelHL <- hildreth.lu.func.multi(0.93529, model)
summary(modelHL)

## 
## Call:
## lm(formula = y ~ x1 + x2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11773.7  -1196.2   -221.5   1232.2   9006.3 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -179.87246  853.16575  -0.211   0.8337  
## x1             0.09565    0.03657   2.615   0.0112 *
## x2             0.18645    0.49237   0.379   0.7062  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2920 on 61 degrees of freedom
## Multiple R-squared:  0.1117, Adjusted R-squared:  0.08261 
## F-statistic: 3.836 on 2 and 61 DF,  p-value: 0.02695

#Transformasi Balik
cat("y = ", coef(modelHL)[1]/(1-0.93529), " + ", coef(modelHL)[2]," x1", " + ",coef(modelHL)[3]," x2", sep = "")

## y = -2779.67 + 0.09564552 x1 + 0.1864518 x2

Setelah dilakukan tranformasi balik, didapatkan model dengan metode Hildreth-Lu sebagai berikut. \[y_i=-2619.329 + 0.09564552x_1 +0.1864518x_2\]

#Deteksi autokorelasi
lmtest::dwtest(modelHL, alternative = 'two.sided')

## Warning in lmtest::dwtest(modelHL, alternative = "two.sided"): exact p value
## cannot be computed (not in [0,1]), approximate p value will be used

## 
##  Durbin-Watson test
## 
## data:  modelHL
## DW = 1.8828, p-value = 0.5756
## alternative hypothesis: true autocorrelation is not 0

Hasil uji Durbin-Watson juga menunjukkan bawah nilai DW sebesar $1.8828$ berada pada selang daerah tidak ada autokorelasi, yaitu pada rentang DU < DW < 4-DU. Hal tersebut juga didukung oleh p-value sebesar $0.5756$, di mana p-value > $\alpha$=5%. Artinya tak tolak $H_0$ atau belum cukup bukti menyatakan bahwa ada autokorelasi dalam data nilai IPM 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)]
sseModelCO <- anova(modelCOmanual)$`Sum Sq`[-c(1,2)]
sseModelHL <- anova(modelHL)$`Sum Sq`[-c(1,2)]
mseModelawal <- sseModelawal/length(Bunker.trans)
mseModelCO <- sseModelCO/length(Bunker.trans)
mseModelHL <- sseModelHL/length(Bunker.trans)

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 3381549788             520050663         520050663
## MSE   52836715               8125792           8125792

Berdasarkan hasil tersebut dapat diketahui bahwa hasil penanganan autokorelasi dengan metode Cochrane-Orcutt dan Hildreth-Lu memiliki SSE yang sama, sebesar $520050663$ dan lebih baik dibandingkan model awal ketika autokorelasi masih terjadi, yaitu sebesar $3381549788 (ingat dalam SSE maupun MSE, semakin kecil nilainya semakin baik).

Simpulan

Autokorelasi yang terdapat pada data yang dianalisis terjadi akibat adanya korelasi di antara unsur penyusunnya. Indikator bunker fuel yang erat hubungannya dengan gas (fuel & flaring) sangat rawan menjadi penyebab adanya autokorelasi. 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. Kemudian, autokorelasi tersebut dapat ditangani dengan metode Cochrane-Orcutt dan Hildreth-Lu. Kedua metode menghasilkan nilai SSE yang sama dan lebih kecil dari model awalnya, artinya keduanya baik untuk digunakan.

Daftar Pustaka

Aprianto A, Debataraja NN, Imro’ah N. 2020. Metode cochrane-orcutt untuk mengatasi autokorelasi pada estimasi parameter ordinary least squares. Bimaster : Buletin Ilmiah Matematika, Statistika dan Terapannya. 9(1):95–102. doi:10.26418/bbimst.v9i1.38590.

BPS. 2021a. Indeks Pembangunan Manusia 2020. Jakarta (ID): Badan Pusat Statistik.

BPS. 2021b. Indeks Pembangunan Manusia (IPM) 2021. Berita Resmi Statistik., siap terbit.