Pemanggilan Packages

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.5.3

## 
## 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(lmtest)

## Warning: package 'lmtest' was built under R version 4.5.3

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.5.2

## 
## Attaching package: 'zoo'

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

library(orcutt)
library(readxl)

## Warning: package 'readxl' was built under R version 4.5.2

Input Data

Data yang digunakan adalah data Indeks Pembangunan Manusia (IPM) Provinsi Bali tahun 2010–2021. Variabel yang digunakan dalam analisis adalah Tahun sebagai variabel bebas (X) dan IPM sebagai variabel terikat (Y).

dtBali <- read_excel("C:/Users/hp/Downloads/data_ipm_bali.xlsx")

dtBali <- dtBali[, c("Tahun", "IPM")]

dtBali

## # A tibble: 12 × 2
##    Tahun   IPM
##    <dbl> <dbl>
##  1  2010  70.1
##  2  2011  70.9
##  3  2012  71.6
##  4  2013  72.1
##  5  2014  72.5
##  6  2015  73.3
##  7  2016  73.6
##  8  2017  74.3
##  9  2018  74.8
## 10  2019  75.4
## 11  2020  75.5
## 12  2021  75.7

str(dtBali)

## tibble [12 × 2] (S3: tbl_df/tbl/data.frame)
##  $ Tahun: num [1:12] 2010 2011 2012 2013 2014 ...
##  $ IPM  : num [1:12] 70.1 70.9 71.6 72.1 72.5 ...

summary(dtBali)

##      Tahun           IPM       
##  Min.   :2010   Min.   :70.10  
##  1st Qu.:2013   1st Qu.:71.97  
##  Median :2016   Median :73.46  
##  Mean   :2016   Mean   :73.31  
##  3rd Qu.:2018   3rd Qu.:74.92  
##  Max.   :2021   Max.   :75.69

Data terdiri atas 12 pengamatan, yaitu IPM Provinsi Bali dari tahun 2010 sampai dengan 2021. Dalam analisis ini, Tahun digunakan sebagai variabel independen (X), sedangkan IPM digunakan sebagai variabel dependen (Y).

Eksplorasi Data

Plot Data IPM terhadap Waktu

plot(
  dtBali$Tahun,
  dtBali$IPM,
  type = "o",
  pch = 16,
  main = "Perkembangan IPM Provinsi Bali Tahun 2010–2021",
  xlab = "Tahun",
  ylab = "IPM"
)

Berdasarkan plot, nilai IPM Provinsi Bali menunjukkan pola yang cenderung meningkat dari tahun ke tahun. Hal tersebut menunjukkan adanya hubungan positif antara waktu (Tahun) dengan nilai IPM.

Namun, karena data merupakan data runtun waktu, pola yang meningkat secara berurutan juga dapat menyebabkan nilai pada suatu periode berkorelasi dengan nilai pada periode sebelumnya. Oleh karena itu, selain melihat hubungan antara Tahun dan IPM, perlu dilakukan pemeriksaan terhadap autokorelasi pada sisaan model regresi.

Korelasi Tahun dan IPM

cor(dtBali$Tahun, dtBali$IPM)

## [1] 0.9923534

Nilai korelasi digunakan untuk melihat arah dan kekuatan hubungan linear antara Tahun dan IPM. Berdasarkan hasil korelasi, hubungan antara Tahun dan IPM bersifat positif, yang menunjukkan bahwa peningkatan tahun cenderung diikuti oleh peningkatan nilai IPM.

Namun, korelasi saja belum cukup untuk menjelaskan hubungan secara matematis. Oleh karena itu, selanjutnya dilakukan analisis regresi linear sederhana.

Membentuk Model Regresi

model <- lm(IPM ~ Tahun, data = dtBali)

summary(model)

## 
## Call:
## lm(formula = IPM ~ Tahun, data = dtBali)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4800 -0.1125  0.0800  0.1725  0.2500 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -974.75000   41.22291  -23.65 4.15e-10 ***
## Tahun          0.52000    0.02045   25.42 2.03e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2446 on 10 degrees of freedom
## Multiple R-squared:  0.9848, Adjusted R-squared:  0.9832 
## F-statistic: 646.4 on 1 and 10 DF,  p-value: 2.033e-10

Berdasarkan hasil regresi linear sederhana, diperoleh persamaan:

\[ \widehat{IPM} = -974.75 + 0.52(Tahun) \]

Koefisien Tahun sebesar 0,52 menunjukkan bahwa setiap kenaikan satu tahun berkaitan dengan peningkatan rata-rata nilai IPM sebesar sekitar 0,52 poin.

Koefisien Tahun memiliki p-value sebesar \(2.03 \times 10^{-10}\), yang lebih kecil dari taraf signifikansi 5%. Dengan demikian, terdapat bukti yang cukup bahwa terdapat hubungan linear yang signifikan antara waktu (Tahun) dan IPM.

Nilai \(R^2\) sebesar 0,9848 menunjukkan bahwa sekitar 98,48% variasi IPM dapat dijelaskan oleh komponen tren waktu dalam model regresi linear sederhana ini. Namun, nilai \(R^2\) yang tinggi tidak serta-merta menunjukkan bahwa model bebas dari masalah autokorelasi, sehingga pemeriksaan terhadap sisaan tetap diperlukan.

Plot Diagnostik Sisaan

sisaan <- residuals(model)
fitValue <- fitted(model)

par(mfrow = c(2,2))

qqnorm(
  sisaan,
  main = "Normal Q-Q Plot"
)
qqline(
  sisaan,
  col = "steelblue",
  lwd = 2
)

plot(
  fitValue,
  sisaan,
  pch = 16,
  col = "steelblue",
  xlab = "Fitted Values",
  ylab = "Sisaan",
  main = "Sisaan vs Fitted Values"
)
abline(h = 0, lwd = 2)

hist(
  sisaan,
  main = "Histogram Sisaan",
  xlab = "Sisaan"
)

plot(
  seq_along(sisaan),
  sisaan,
  pch = 16,
  col = "steelblue",
  xlab = "Urutan Pengamatan",
  ylab = "Sisaan",
  main = "Sisaan vs Urutan"
)
lines(
  seq_along(sisaan),
  sisaan
)
abline(h = 0, lwd = 2)

par(mfrow = c(1,1))

Plot diagnostik digunakan untuk melihat karakteristik sisaan model. Plot Q-Q digunakan untuk melihat kesesuaian distribusi sisaan terhadap distribusi normal, sedangkan plot sisaan terhadap fitted values digunakan untuk melihat pola tertentu pada sisaan.

Khusus untuk tugas ini, perhatian utama diberikan pada pola sisaan berdasarkan urutan waktu, karena data yang digunakan merupakan data runtun waktu. Jika sisaan pada suatu periode memiliki hubungan dengan sisaan pada periode sebelumnya, maka terdapat indikasi autokorelasi.

Uji Shapiro-Wilk

shapiro.test(sisaan)

## 
##  Shapiro-Wilk normality test
## 
## data:  sisaan
## W = 0.89467, p-value = 0.1354

Uji Shapiro-Wilk digunakan untuk menguji apakah sisaan mengikuti distribusi normal.

Hipotesis:

\[ H_0:\text{sisaan berdistribusi normal} \] \[ H_1:\text{sisaan tidak berdistribusi normal} \]

Dengan taraf signifikansi 5%, diperoleh p-value sebesar 0,1354. Karena p-value > 0,05, maka tidak terdapat cukup bukti untuk menolak \(H_0\).

Dengan demikian, berdasarkan uji Shapiro-Wilk, tidak terdapat cukup bukti bahwa sisaan menyimpang dari distribusi normal.

Pemeriksaan Autokorelasi (ACF dan PACF)

par(mfrow = c(1,2))

acf(
  sisaan,
  main = "ACF Sisaan"
)

pacf(
  sisaan,
  main = "PACF Sisaan"
)

par(mfrow = c(1,1))

Plot ACF digunakan untuk melihat korelasi antara sisaan pada suatu periode dengan sisaan pada beberapa periode sebelumnya. Adanya korelasi yang signifikan pada lag tertentu dapat menjadi indikasi adanya autokorelasi.

Untuk memperoleh keputusan formal, selanjutnya digunakan uji Durbin-Watson.

Uji Durbin-Watson

dw_awal <- dwtest(model)

dw_awal

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

Uji Durbin-Watson digunakan untuk menguji keberadaan autokorelasi positif pada sisaan model.

Hipotesis yang digunakan adalah:

\[ H_0: \rho = 0 \]

\[ H_1: \rho > 0 \]

Dengan taraf signifikansi 5%, diperoleh nilai Durbin-Watson sebesar 0,87546 dan p-value sebesar 0,003527.

Karena p-value < 0,05, maka \(H_0\) ditolak. Dengan demikian, terdapat bukti yang cukup bahwa sisaan model regresi awal mengalami autokorelasi positif.

Nilai Durbin-Watson yang berada cukup jauh di bawah 2 juga konsisten dengan adanya autokorelasi positif. Oleh karena itu, diperlukan penanganan terhadap autokorelasi sebelum model digunakan sebagai model akhir.

Penanganan Autokorelasi

Metode Cochrane-Orcutt

modelCO <- cochrane.orcutt(model)

## Warning in cochrane.orcutt(model): Did not converge

modelCO

## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = IPM ~ Tahun, data = dtBali)
## 
##  number of interaction: 100
##  rho 0.734645
## 
## Durbin-Watson statistic 
## (original):    0.87546 , p-value: 3.527e-03
## (transformed): NA , p-value: NA
##  
##  coefficients: 
## [1] NA

Metode Cochrane-Orcutt digunakan untuk menangani autokorelasi dengan melakukan transformasi model berdasarkan estimasi parameter autokorelasi.

Pada data IPM Provinsi Bali, metode Cochrane-Orcutt memberikan peringatan “Did not converge” setelah mencapai 100 iterasi. Selain itu, koefisien model hasil transformasi tidak dapat diperoleh karena menghasilkan nilai NA.

Dengan demikian, metode Cochrane-Orcutt tidak mencapai konvergensi pada data yang digunakan. Oleh karena itu, hasil Cochrane-Orcutt tidak digunakan sebagai model akhir dan penanganan autokorelasi dilanjutkan menggunakan metode Hildreth-Lu.

Metode Hildreth-Lu

Metode Hildreth-Lu dilakukan dengan mencari nilai \(\rho\) yang menghasilkan SSE terkecil.

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))
}

Fungsi tersebut melakukan transformasi:

\[ Y_t^*=Y_t-\rho Y_{t-1} \]

dan

\[ X_t^*=X_t-\rho X_{t-1} \]

untuk memperoleh model yang telah ditransformasi.

Pencarian Nilai Rho

Pencarian Awal

r <- 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.3873
## 2 0.2 0.3555
## 3 0.3 0.3291
## 4 0.4 0.3082
## 5 0.5 0.2927
## 6 0.6 0.2827
## 7 0.7 0.2780
## 8 0.8 0.2788
## 9 0.9 0.2851

Pencarian awal dilakukan pada nilai \(\rho\) dari 0,1 sampai 0,9 dengan interval 0,1.

Berdasarkan hasil pencarian, SSE terkecil diperoleh pada:

\[ \rho=0.7 \]

dengan SSE sebesar:

\[ SSE=0.2780 \]

Nilai tersebut kemudian digunakan sebagai pusat pencarian yang lebih teliti.

Pencarian Lebih Teliti

Karena nilai awal terbaik berada di sekitar 0,7, dilakukan pencarian pada interval 0,6 sampai 0,8 dengan interval 0,001.

rOpt <- seq(0.6, 0.8, by = 0.001)

tabOpt <- data.frame(
  rho = rOpt,
  SSE = sapply(
    rOpt,
    function(i){
      deviance(hildreth.lu.func(i, model))
    }
  )
)

head(
  tabOpt[order(tabOpt$SSE), ],
  10
)

##       rho       SSE
## 136 0.735 0.2776919
## 137 0.736 0.2776920
## 135 0.734 0.2776924
## 138 0.737 0.2776926
## 134 0.733 0.2776934
## 139 0.738 0.2776938
## 133 0.732 0.2776949
## 140 0.739 0.2776955
## 132 0.731 0.2776970
## 141 0.740 0.2776978

rho_opt <- tabOpt$rho[which.min(tabOpt$SSE)]

rho_opt

## [1] 0.735

tabOpt[which.min(tabOpt$SSE), ]

##       rho       SSE
## 136 0.735 0.2776919

Metode Hildreth-Lu dilakukan dengan mencari nilai \(\rho\) yang menghasilkan Sum of Squared Errors (SSE) terkecil setelah transformasi data.

Pencarian awal dilakukan pada interval 0,1 sampai 0,9 dengan interval 0,1. Hasil pencarian menunjukkan bahwa SSE terkecil berada di sekitar \(\rho = 0,7\). Selanjutnya dilakukan pencarian yang lebih teliti pada interval 0,6 sampai 0,8 dengan interval 0,001.

Berdasarkan pencarian tersebut, diperoleh nilai optimum:

\[ \rho = 0.735 \]

dengan SSE minimum sebesar 0,2776919.

Nilai \(\rho = 0,735\) kemudian digunakan sebagai parameter transformasi dalam metode Hildreth-Lu.

Grafik Hubungan Rho dan SSE

plot(
  tabOpt$rho,
  tabOpt$SSE,
  type = "l",
  xlab = expression(rho),
  ylab = "SSE",
  main = "Hubungan Rho dengan SSE"
)

abline(
  v = rho_opt,
  lty = 2,
  col = "red",
  lwd = 2
)

points(
  rho_opt,
  min(tabOpt$SSE),
  pch = 19
)

text(
  rho_opt,
  min(tabOpt$SSE),
  labels = paste0(
    "rho = ",
    round(rho_opt, 3)
  ),
  pos = 4
)

Grafik menunjukkan hubungan antara nilai \(\rho\) dengan SSE. Titik minimum pada grafik menunjukkan nilai \(\rho\) yang menghasilkan SSE paling kecil.

Berdasarkan grafik dan hasil pencarian numerik, SSE minimum diperoleh pada:

\[ \boxed{\rho=0.735} \]

Oleh karena itu, nilai tersebut digunakan dalam pembentukan model Hildreth-Lu.

Pembentukan Model Hildreth-Lu

modelHL <- hildreth.lu.func(
  rho_opt,
  model
)

summary(modelHL)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.18521 -0.15406 -0.03417  0.13618  0.24766 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -180.9202    33.8105  -5.351 0.000462 ***
## x              0.3753     0.0632   5.939 0.000218 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1757 on 9 degrees of freedom
## Multiple R-squared:  0.7967, Adjusted R-squared:  0.7741 
## F-statistic: 35.27 on 1 and 9 DF,  p-value: 0.0002183

Dengan menggunakan nilai \(\rho = 0,735\), diperoleh model pada skala transformasi:

\[ Y_t^* = -180.9202 + 0.3753X_t^* + \varepsilon_t \]

Koefisien variabel Tahun pada model transformasi bernilai positif sebesar 0,3753. Hal ini menunjukkan bahwa setelah dilakukan transformasi untuk menangani autokorelasi, hubungan antara waktu dan IPM tetap menunjukkan arah positif.

Koefisien tersebut memiliki p-value sebesar 0,000218, sehingga lebih kecil dari 0,05. Dengan demikian, terdapat bukti yang cukup bahwa hubungan linear antara komponen waktu dan IPM tetap signifikan setelah transformasi.

Transformasi Balik Model

Karena model Hildreth-Lu menggunakan data yang telah ditransformasi, intercept perlu dikembalikan ke skala model awal.

beta0_HL <- coef(modelHL)[1] / (1 - rho_opt)
beta1_HL <- coef(modelHL)[2]

cat(
  "IPM = ",
  round(beta0_HL, 4),
  " + ",
  round(beta1_HL, 7),
  " Tahun",
  sep = ""
)

## IPM = -682.7179 + 0.3753413 Tahun

Setelah dilakukan transformasi balik terhadap model Hildreth-Lu, diperoleh persamaan:

\[ \widehat{IPM} = -682.7179 + 0.3753413(Tahun) \]

Koefisien Tahun sebesar 0,3753413 menunjukkan bahwa setiap peningkatan waktu sebesar satu tahun berkaitan dengan peningkatan IPM sekitar 0,3753 poin, dengan struktur autokorelasi yang telah diperhitungkan dalam model.

Intercept sebesar -682,7179 tidak memiliki interpretasi praktis secara langsung karena Tahun = 0 berada jauh di luar rentang pengamatan, yaitu 2010–2021.

Evaluasi Autokorelasi Setelah Hildreth-Lu

Setelah model ditransformasi menggunakan \(\rho=0.735\), dilakukan kembali uji Durbin-Watson.

dw_HL <- dwtest(modelHL)

dw_HL

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

Setelah dilakukan transformasi menggunakan metode Hildreth-Lu, dilakukan kembali uji Durbin-Watson untuk mengevaluasi autokorelasi pada sisaan model.

Hipotesis yang digunakan adalah:

\[ H_0: \rho = 0 \]

\[ H_1: \rho > 0 \]

Hasil pengujian menunjukkan nilai Durbin-Watson sebesar 2,1296 dengan p-value sebesar 0,4359.

Karena p-value > 0,05, maka \(H_0\) tidak ditolak. Dengan demikian, tidak terdapat cukup bukti adanya autokorelasi positif pada sisaan model setelah transformasi Hildreth-Lu.

Nilai Durbin-Watson juga meningkat dari 0,87546 pada model awal menjadi 2,1296 setelah transformasi, sehingga hasil tersebut menunjukkan bahwa masalah autokorelasi positif yang terdeteksi pada model awal telah tertangani.

# Sisaan setelah transformasi Hildreth-Lu
sisaan_HL <- residuals(modelHL)

# ACF dan PACF Sisaan Setelah Transformasi Hildreth-Lu
par(mfrow = c(1, 2))

acf(
  sisaan_HL,
  main = "ACF Sisaan Setelah Hildreth-Lu"
)

pacf(
  sisaan_HL,
  main = "PACF Sisaan Setelah Hildreth-Lu"
)

par(mfrow = c(1, 1))

Perbandingan Model Sebelum dan Sesudah Penanganan

SSE_awal <- deviance(model)
MSE_awal <- mean(residuals(model)^2)

SSE_awal

## [1] 0.5982

MSE_awal

## [1] 0.04985

SSE_HL <- deviance(modelHL)
MSE_HL <- mean(residuals(modelHL)^2)

SSE_HL

## [1] 0.2776919

MSE_HL

## [1] 0.02524472

perbandingan <- data.frame(
  Model = c(
    "Regresi Awal",
    "Hildreth-Lu"
  ),
  SSE = c(
    SSE_awal,
    SSE_HL
  ),
  MSE = c(
    MSE_awal,
    MSE_HL
  ),
  DW = c(
    dw_awal$statistic,
    dw_HL$statistic
  ),
  p_value_DW = c(
    dw_awal$p.value,
    dw_HL$p.value
  )
)

perbandingan

##          Model       SSE        MSE        DW  p_value_DW
## 1 Regresi Awal 0.5982000 0.04985000 0.8754597 0.003526813
## 2  Hildreth-Lu 0.2776919 0.02524472 2.1295833 0.435883328

Berdasarkan hasil perbandingan, model Hildreth-Lu menghasilkan SSE dan MSE yang lebih kecil dibandingkan model regresi awal. Namun, nilai SSE dan MSE kedua model tidak dapat dijadikan satu-satunya dasar perbandingan karena model Hildreth-Lu menggunakan data yang telah ditransformasi dan kehilangan satu pengamatan akibat proses transformasi.

Oleh karena itu, keberhasilan penanganan autokorelasi terutama dievaluasi berdasarkan perubahan hasil uji Durbin-Watson.

Pada model regresi awal diperoleh DW sebesar 0,87546 dengan p-value 0,003527, yang menunjukkan adanya bukti autokorelasi positif. Setelah dilakukan transformasi Hildreth-Lu, nilai DW meningkat menjadi 2,1296 dengan p-value 0,4359. Karena p-value setelah transformasi lebih besar dari 0,05, tidak terdapat cukup bukti adanya autokorelasi positif pada sisaan.

Dengan demikian, berdasarkan uji Durbin-Watson, metode Hildreth-Lu berhasil menangani masalah autokorelasi positif pada model regresi IPM Provinsi Bali.

Kesimpulan

Berdasarkan seluruh tahapan analisis, diperoleh beberapa kesimpulan sebagai berikut.

Model regresi linear sederhana antara Tahun dan IPM Provinsi Bali menunjukkan hubungan linear positif yang signifikan. Koefisien Tahun sebesar 0,52 menunjukkan bahwa peningkatan satu tahun berkaitan dengan peningkatan IPM sekitar 0,52 poin.
Pemeriksaan terhadap sisaan menunjukkan adanya autokorelasi positif. Hal ini ditunjukkan oleh hasil uji Durbin-Watson dengan nilai DW sebesar 0,87546 dan p-value sebesar 0,003527. Karena p-value < 0,05, terdapat bukti yang cukup adanya autokorelasi positif pada model awal.
Penanganan autokorelasi menggunakan metode Cochrane-Orcutt tidak mencapai konvergensi setelah 100 iterasi. Oleh karena itu, hasil Cochrane-Orcutt tidak digunakan sebagai model akhir.
Penanganan selanjutnya dilakukan menggunakan metode Hildreth-Lu. Berdasarkan pencarian nilai \(\rho\), diperoleh nilai optimum sebesar 0,735 yang menghasilkan SSE minimum sebesar 0,2776919.
Setelah dilakukan transformasi dan transformasi balik, diperoleh model: