Packages yang digunakan

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 bersumber dari: https://www.kaggle.com/datasets/anuragbantu/new-delhi-air-quality

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

Penyiapan data

#Membuat peubah periode
Periode = seq(from = 1, to = 72, by = 1)
data = data[,-c(1,4,10,11,12)] #menghapus beberapa kolom yang tidak diperlukan

#Menggabungkan semua peubah ke dalam satu data frame bernama df
df = data.frame(Periode, data)

Eksplorasi Data

#Matriks korelasi antar peubah
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ readr     2.1.4
## ✔ ggplot2   3.4.2     ✔ stringr   1.5.0
## ✔ lubridate 1.9.2     ✔ tibble    3.2.1
## ✔ purrr     1.0.1     ✔ tidyr     1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
df %>%
  as_tibble() %>%
  select(-Periode) %>%
  cor() %>%
  ggcorrplot::ggcorrplot(type = "upper", lab = TRUE, lab_size = 3) +
  theme_light() +
  labs(title = "Korelasi antar peubah",
       subtitle = "Peubah respon: AQI.", 
       caption = "Source: https://www.kaggle.com/datasets/anuragbantu/new-delhi-air-quality",
       x = NULL, y = NULL)

Dari eksplorasi menggunakan plot korelasi, terlihat bahwa peubah CO, no2, o3, dan so2 memiliki nilai korelasi yang tinggi terhadap AQI. Oleh karena itu, peubah-peubah inilah yang akan digunakan sebagai peubah penjelas dalam tahapan analisis berikutnya.

Sebagai peubah respon: AQI Sebagai peubah penjelas: CO, no2, o3, dan so2

Pembentukan Penduga Model Regresi

model_aqi = lm(AQI~CO+o3+so2, data=df)
summary(model_aqi)
## 
## Call:
## lm(formula = AQI ~ CO + o3 + so2, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6297 -0.3961 -0.0265  0.1540  4.3012 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.130858   7.389425  -0.288    0.774    
## CO          -0.002784   0.036556  -0.076    0.940    
## o3           0.485279   0.024961  19.442   <2e-16 ***
## so2          3.894798   3.338980   1.166    0.248    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7034 on 68 degrees of freedom
## Multiple R-squared:  0.9492, Adjusted R-squared:  0.947 
## F-statistic: 423.6 on 3 and 68 DF,  p-value: < 2.2e-16

Model yang dihasilkan adalah \[AQI_i=-2.130858-0.002784CO_i+0.485279O3_i+3.894798SO2_i\] Dari hasil ringkasan model diperoleh beberapa informasi sebagai berikut:

  1. Uji F \(p-value = 0.000 < \alpha= 5% (Tolak H_0)\) Artinya: Minimal terdapat satu peubah penjelas yang berpengaruh signifikan terhadap nilai AQI pada taraf 5%.

  2. Uji t Hasil uji-t parsial parameter regresi koefisien regresi so3 juga menunjukkan hal yang sama, yaitu memiliki p-value < \(\alpha\) (5%) sehingga nyata dalam taraf 5%.

  3. \(R^2\) Diperoleh nilai \(R^2=0.9492\) Artinya, sebesar 94.92% keragaman nilai AQI dapat dijelaskan oleh model yang telah dibentuk. Hasil ini menunjukkan hasil yang bagus, seolah mendapatkan hasil terbaik.

Namun, pengujian asumsi tetap perlu dilakukan.

Pengujian Asumsi

Uji Multikolinearitas

car::vif(model_aqi)
##       CO       o3      so2 
## 3.909863 3.697950 3.191173

Nilai VIF pada setiap peubah penjelas < 10. Artinya, tidak terjadi multikolinearitas pada peubah penjelas yang digunakan

Sisaan menyebar normal

#Dengan eksplorasi
plot(model_aqi,2);

Hasil QQ-Plot memperlihatkan bahwa titik-titik nya cenderung mengikuti garis kenormalan. Oleh karena itu, dapat disimpulkan bahwa sisaan menyebar normal.

Nilai harapan sisaan sama dengan nol

\(H_0 : E[\varepsilon]=0\) \(H_1 : E[\varepsilon]\ne0\)

# Uji t
t.test(resid(model_aqi), mu = 0,)
## 
##  One Sample t-test
## 
## data:  resid(model_aqi)
## t = 1.5344e-15, df = 71, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.1617684  0.1617684
## sample estimates:
##    mean of x 
## 1.244882e-16

\(p-value=1 > 0.1\) (tidak tolak \(H_0\)) Artinya: Nilai harapan sisaan sama dengan nol

Ragam sisaan homogen

\(H_0: Var[\varepsilon]=\sigma^2I\) \(H_1:Var[\varepsilon]\ne \sigma^2I\)

# Uji Breusch-Pagan
lmtest::bptest(model_aqi)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_aqi
## BP = 1.76, df = 3, p-value = 0.6237

\(p-value=0.6237> 0.1\) (tidak tolak \(H_0\)) Artinya: Ragam sisaan homogen

Antar sisaan tidak saling berkorelasi (tidak terjadi autokorelasi)

\(H_0 : Cov[\varepsilon_i,\varepsilon_j]=0\) (tidak terjadi autokorelasi pada sisaan) \(H_1 : Cov[\varepsilon_i,\varepsilon_j]\neq0\) (terjadi autokorelasi pada sisaan)

# Uji Durbin Watson
library(lmtest)
dwtest(model_aqi)
## 
##  Durbin-Watson test
## 
## data:  model_aqi
## DW = 0.55676, p-value = 1.13e-14
## alternative hypothesis: true autocorrelation is greater than 0
#ACF dan PACF identifikasi autokorelasi
sisaan = model_aqi$residuals
par(mfrow = c(1,2))
acf(sisaan)
pacf(sisaan)

\(p-value=0.00<0.1\) (tolak \(H_0\)) Artinya: Terjadi autkorelasi pada sisaan pada taraf 5%

Selain dari itu, plot ACF juga memperlihatkan adanya autokorelasi. Oleh karena itu, diperlukan penanganan untuk hal tersebut. Terdapat 2 metode yang akan dicobakan: 1. Cochrane-Orcutt 2. Hildreth-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 dengan Packages Orcutt
library(orcutt)
modelCO<-cochrane.orcutt(model_aqi)
modelCO
## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = AQI ~ CO + o3 + so2, data = df)
## 
##  number of interaction: 9
##  rho 0.46121
## 
## Durbin-Watson statistic 
## (original):    0.55676 , p-value: 1.13e-14
## (transformed): 2.71731 , p-value: 9.976e-01
##  
##  coefficients: 
## (Intercept)          CO          o3         so2 
##    4.989103   -0.031450    0.483940   -0.027116

Hasil keluaran model setelah dilakukan penanganan adalah sebagai berikut. \[y_i=4.989103 -0.031450CO_i+0.483940O3_i-0.027116SO2_i\] Hasil juga menunjukkan bahwanilai DW dan p-value meningkat menjadi \(2.71731\) dan \(9.976e-01\). 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 terdapat autokorelasi pada sisaan pada taraf nyata 5%. Untuk nilai \(ρ ̂\) optimum yang digunakan adalah \(0.46121\).

Selanjutnya akan dilakukan penghitungan secara manual

#Rho optimum
rho<- modelCO$rho
rho
## [1] 0.4612096

Selanjutnya akan dilakukan transformasi secara manual dengan syntax berikut ini.

#Transformasi Manual
#tanpa data pertama - tanpa data terakhir
aqi.trans <- df$AQI[-1]-df$AQI[-72]*rho 
CO.trans <- df$CO[-1]-df$CO[-72]*rho
o3.trans <- df$o3[-1]-df$o3[-72]*rho
so2.trans <- df$so2[-1]-df$so2[-72]*rho

#Membentuk model dengan peubah yang sudah ditransformasi
modelCOmanual<- lm(aqi.trans~CO.trans+o3.trans+so2.trans)
summary(modelCOmanual)
## 
## Call:
## lm(formula = aqi.trans ~ CO.trans + o3.trans + so2.trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.93684 -0.23232  0.04205  0.18381  0.79208 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.68808    3.13685   0.857    0.395    
## CO.trans    -0.03145    0.02819  -1.115    0.269    
## o3.trans     0.48394    0.01940  24.943   <2e-16 ***
## so2.trans   -0.02712    2.69876  -0.010    0.992    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3435 on 67 degrees of freedom
## Multiple R-squared:  0.9635, Adjusted R-squared:  0.9619 
## F-statistic: 589.6 on 3 and 67 DF,  p-value: < 2.2e-16

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

#Mencari Penduga Koefisien Regresi setelah Transformasi ke Persamaan Awal
#b0 bintang = menghapus koefisien b1,b2,b3 dari model
b0bintang <- modelCOmanual$coefficients[-c(2,3,4)]
b0 <- b0bintang/(1-rho)

#b1 = menghapus koefisien b0,b2,b3 dari model
b1 <- modelCOmanual$coefficients[-c(1,3,4)]

#b2 = menghapus koefisien b0,b1,b3 dari model
b2 <- modelCOmanual$coefficients[-c(1,2,4)]

#b3 = menghapus koefisien b0,b1,b2, b3 dari model
b3 <- modelCOmanual$coefficients[-c(1,2,3)]

b0;b1;b2;b3
## (Intercept) 
##    4.989103
##   CO.trans 
## -0.0314504
##  o3.trans 
## 0.4839403
##   so2.trans 
## -0.02711614

Hasil perhitungan secara manual maupun dengan package menghasilkan kesimpulan yang sama.

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]
  
  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]
  
  
  return(lm(y~x1+x2+x3))
}

#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_aqi))}))
round(tab, 4)
##   rho     SSE
## 1 0.1 12.0739
## 2 0.2 10.0700
## 3 0.3  8.7231
## 4 0.4  8.0226
## 5 0.5  7.9523
## 6 0.6  8.4915
## 7 0.7  9.6241
## 8 0.8 11.3621
## 9 0.9 13.7733

Pertama-tama akan dicari di mana kira-kira \(ρ\) yang menghasilkan SSE minimum. Pada hasil di atas terlihat \(ρ\) minimum ketika 0.4. 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.4 sampai dengan 0.5.

#Rho optimal di sekitar 0.4
rOpt <- seq(0.4,0.5, by= 0.001)
tabOpt <- data.frame("rho" = rOpt, "SSE" = sapply(rOpt, function(i){deviance(hildreth.lu.func(i, model_aqi))}))
head(tabOpt[order(tabOpt$SSE),])
##      rho      SSE
## 62 0.461 7.905986
## 63 0.462 7.906004
## 61 0.460 7.906030
## 64 0.463 7.906084
## 60 0.459 7.906136
## 65 0.464 7.906225
#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.461, y=7.905986, labels = "rho=0.461", cex = 0.8)

Perhitungan yang dilakukan aplikasi R menunjukkan bahwa nilai \(ρ\) optimum, yaitu saat SSE terkecil terdapat pada nilai \(ρ=0.461\). 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.461, model_aqi)
summary(modelHL)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.93664 -0.23231  0.04208  0.18384  0.79297 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.68575    3.13731   0.856    0.395    
## x1          -0.03142    0.02819  -1.115    0.269    
## x2           0.48393    0.01940  24.948   <2e-16 ***
## x3          -0.02506    2.69806  -0.009    0.993    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3435 on 67 degrees of freedom
## Multiple R-squared:  0.9635, Adjusted R-squared:  0.9619 
## F-statistic:   590 on 3 and 67 DF,  p-value: < 2.2e-16
#Transformasi Balik
cat("y = ", coef(modelHL)[1]/(1-0.461), "+", coef(modelHL)[2],"x", sep = "")
## y = 4.98283+-0.03142007x

Setelah dilakukan tranformasi balik, didapatkan model dengan metode Hildreth-Lu sebagai berikut. \[y_i=2.68575-0.03142x1_i+0.48393x2_i-0.02506x3_i\]

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

Hasil uji Durbin-Watson juga menunjukkan bawah nilai DW sebesar \(0.9976\) berada pada selang daerah tidak ada autokorelasi, yaitu pada rentang DU < DW < 4-DU. Hal tersebut juga didukung oleh p-value sebesar \(0.9976\), 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_aqi)$`Sum Sq`[-c(1,2,3)]
sseModelCO <- anova(modelCOmanual)$`Sum Sq`[-c(1,2,3)]
sseModelHL <- anova(modelHL)$`Sum Sq`[-c(1,2,3)]
mseModelawal <- sseModelawal/length(df$AQI)
mseModelCO <- sseModelCO/length(df$AQI)
mseModelHL <- sseModelHL/length(df$AQI)
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 33.6474846             7.9059848         7.9059861
## MSE  0.4673262             0.1098053         0.1098054

Hasil penanganan autokorelasi dengan metode Cochrane-Orcutt dan Hildreth-Lu memiliki nilai SSE sebesar \(7.9059861\). Jauh lebih rendah dibandingkan model awal dengan SSE sebesar $33.6474846 $. Hal ini menunjukkan bahwa model setelah penanganan lebih baik dibandingkan model awal ketika autokorelasi masih terjadi.

Kesimpulan: Pada data yang digunakan, metode Cochrane-Orcutt dan Hildreth-Lu terbukti efektif dalam menangani autokorelasi.