Pertemuan 2 - Regresi
Kheni Hikmah Lestari
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, unionlibrary(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(TTR)
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(lmtest) #digunakan untuk uji formal pendeteksian autokorelasi
library(orcutt) #untuk membuat model regresi Cochrane-Orcutt
library(HoRM) #untuk membuat model regresi Hildreth-LuInput Data
Data yang digunakan dalam kesempatan kali ini adalah data IPM Provinsi Gorontalo periode tahun 2010-2021.
library(rio)
data <- import("https://raw.githubusercontent.com/khenihikmah/praktikummpdw/main/Pertemuan2/Salesmonthly.csv")
data## Y X1 X2
## 1 878.030 354.0 50
## 2 1001.900 347.0 31
## 3 779.275 232.0 20
## 4 698.500 209.0 18
## 5 628.780 270.0 23
## 6 548.225 323.0 23
## 7 491.900 348.0 21
## 8 583.850 420.0 29
## 9 887.820 399.0 14
## 10 1856.815 472.0 30
## 11 723.800 489.0 19
## 12 1015.660 492.0 25
## 13 1044.240 463.0 24
## 14 953.252 243.0 9
## 15 1084.850 208.0 13
## 16 940.170 192.0 5
## 17 765.900 194.0 10
## 18 746.788 217.0 12
## 19 708.828 203.0 6
## 20 790.788 265.5 15
## 21 852.125 243.5 11
## 22 1574.335 222.0 8
## 23 1277.725 228.0 18
## 24 1258.349 286.0 28
## 25 1476.324 248.0 24
## 26 1224.862 239.0 20
## 27 1150.700 250.0 13
## 28 998.337 318.0 18
## 29 997.150 275.0 18
## 30 760.050 311.0 20
## 31 652.362 240.0 8
## 32 753.050 275.5 12
## 33 1118.699 307.0 18
## 34 1617.275 312.0 11
## 35 1062.686 246.0 27
## 36 1624.335 257.0 18
## 37 0.000 1.0 0
## 38 526.350 144.0 7
## 39 612.500 165.0 9
## 40 540.200 132.0 9
## 41 547.940 148.0 23
## 42 496.100 163.0 8
## 43 479.350 219.0 15
## 44 549.300 239.0 12
## 45 863.750 223.0 23
## 46 1184.350 226.0 15
## 47 867.899 192.0 15
## 48 1007.180 226.0 6
## 49 1134.325 229.0 11
## 50 1255.374 268.0 12
## 51 999.123 381.0 42
## 52 836.037 289.0 21
## 53 644.648 259.0 13
## 54 584.343 248.0 18
## 55 679.350 283.0 19
## 56 733.838 253.0 20
## 57 1058.262 263.0 12
## 58 1129.275 287.0 25
## 59 995.150 252.2 22
## 60 1213.950 254.0 27
## 61 1660.612 295.2 23
## 62 1001.212 249.4 12
## 63 941.050 301.4 19
## 64 647.650 299.4 22
## 65 703.562 265.8 26
## 66 610.000 193.0 25
## 67 649.800 250.6 20
## 68 518.100 237.0 26
## 69 984.480 227.8 16
## 70 295.150 86.0 7Eksplorasi Data
Sebelum melakukan regresi, akan diperlihatkan plot time-series dari IPM Provinsi Gorontalo Periode 2010-2021
#Membentuk objek time series
data.ts<-ts(data$Y)
data.ts## Time Series:
## Start = 1
## End = 70
## Frequency = 1
## [1] 878.030 1001.900 779.275 698.500 628.780 548.225 491.900 583.850
## [9] 887.820 1856.815 723.800 1015.660 1044.240 953.252 1084.850 940.170
## [17] 765.900 746.788 708.828 790.788 852.125 1574.335 1277.725 1258.349
## [25] 1476.324 1224.862 1150.700 998.337 997.150 760.050 652.362 753.050
## [33] 1118.699 1617.275 1062.686 1624.335 0.000 526.350 612.500 540.200
## [41] 547.940 496.100 479.350 549.300 863.750 1184.350 867.899 1007.180
## [49] 1134.325 1255.374 999.123 836.037 644.648 584.343 679.350 733.838
## [57] 1058.262 1129.275 995.150 1213.950 1660.612 1001.212 941.050 647.650
## [65] 703.562 610.000 649.800 518.100 984.480 295.150#Membuat plot time series
ts.plot(data.ts, xlab="Time Period ", ylab="Y", main= "Time Series Plot of Y")
points(data.ts)
Selanjutnya akan dilakukan ramalan dan pemulusan dengan metode DMA dan DES karena terlihat pada plot di atas menunjukkan adanya trend.
dt.sma <- SMA(data.ts, n=3)
dma <- SMA(dt.sma, n = 3)
At <- 2*dt.sma - dma
Bt <- 2/(4-1)*(dt.sma - dma)
dt.dma<- At+Bt
dt.ramal<- c(NA, dt.sma)
t = 1:14
f = c()
for (i in t) {
f[i] = At[length(At)] + Bt[length(Bt)]*(i)
}dt.gab <- cbind(aktual = c(data.ts,rep(NA,14)),
pemulusan1 = c(dt.sma,rep(NA,14)),
pemulusan2 = c(dt.dma, rep(NA,14)),
At = c(At, rep(NA,14)),
Bt = c(Bt,rep(NA,14)),
ramalan = c(dt.ramal, f[-1]))
as.data.frame(dt.gab)## aktual pemulusan1 pemulusan2 At Bt ramalan
## 1 878.030 NA NA NA NA NA
## 2 1001.900 NA NA NA NA NA
## 3 779.275 886.4017 NA NA NA NA
## 4 698.500 826.5583 NA NA NA 886.4017
## 5 628.780 702.1850 530.74611 599.32167 -68.575556 826.5583
## 6 548.225 625.1683 470.49796 532.36611 -61.868148 702.1850
## 7 491.900 556.3017 436.99611 484.71833 -47.722222 625.1683
## 8 583.850 541.3250 486.42500 508.38500 -21.960000 556.3017
## 9 887.820 654.5233 771.97889 724.99667 46.982222 541.3250
## 10 1856.815 1109.4950 1677.90704 1450.54222 227.364815 654.5233
## 11 723.800 1156.1450 1460.74037 1338.90222 121.838148 1109.4950
## 12 1015.660 1198.7583 1272.02315 1242.71722 29.305926 1156.1450
## 13 1044.240 927.9000 650.62037 761.53222 -110.911852 1198.7583
## 14 953.252 1004.3840 938.88937 965.08722 -26.197852 927.9000
## 15 1084.850 1027.4473 1095.56437 1068.31756 27.246815 1004.3840
## 16 940.170 992.7573 967.02585 977.31844 -10.292593 1027.4473
## 17 765.900 930.3067 841.64481 877.10956 -35.464741 992.7573
## 18 746.788 817.6193 657.71637 721.67756 -63.961185 930.3067
## 19 708.828 740.5053 592.21904 651.53356 -59.314519 817.6193
## 20 790.788 748.8013 715.17800 728.62733 -13.449333 740.5053
## 21 852.125 783.9137 827.53626 810.08722 17.449037 748.8013
## 22 1574.335 1072.4160 1412.48100 1276.45500 136.026000 783.9137
## 23 1277.725 1234.7283 1575.35444 1439.10400 136.250444 1072.4160
## 24 1258.349 1370.1363 1610.76319 1514.51244 96.250741 1234.7283
## 25 1476.324 1337.4660 1376.39230 1360.82178 15.570519 1370.1363
## 26 1224.862 1319.8450 1282.11593 1297.20756 -15.091630 1337.4660
## 27 1150.700 1283.9620 1234.30256 1254.16633 -19.863778 1319.8450
## 28 998.337 1124.6330 927.66578 1006.45267 -78.786889 1283.9620
## 29 997.150 1048.7290 875.87511 945.01667 -69.141556 1124.6330
## 30 760.050 918.5123 731.65826 806.39989 -74.741630 1048.7290
## 31 652.362 803.1873 602.70585 682.89844 -80.192593 918.5123
## 32 753.050 721.8207 567.34381 629.13456 -61.790741 803.1873
## 33 1118.699 841.3703 928.99959 893.94789 35.051704 721.8207
## 34 1617.275 1163.0080 1586.79967 1417.28300 169.516667 841.3703
## 35 1062.686 1266.2200 1559.58759 1442.24056 117.347037 1163.0080
## 36 1624.335 1434.7653 1679.37793 1581.53289 97.845037 1266.2200
## 37 0.000 895.6737 390.31922 592.46100 -202.141778 1434.7653
## 38 526.350 716.8950 218.75667 418.01200 -199.255333 895.6737
## 39 612.500 379.6167 -94.45852 95.17156 -189.630074 716.8950
## 40 540.200 559.6833 572.38056 567.30167 5.078889 379.6167
## 41 547.940 566.8800 674.91333 631.70000 43.213333 559.6833
## 42 496.100 528.0800 488.96704 504.61222 -15.645185 566.8800
## 43 479.350 507.7967 463.70407 481.34111 -17.637037 528.0800
## 44 549.300 508.2500 497.48519 501.79111 -4.305926 507.7967
## 45 863.750 630.8000 767.21852 712.65111 54.567407 508.2500
## 46 1184.350 865.8000 1194.99444 1063.31667 131.677778 630.8000
## 47 867.899 971.9997 1220.55485 1121.13278 99.422074 865.8000
## 48 1007.180 1019.8097 1131.93170 1087.08289 44.848815 971.9997
## 49 1134.325 1003.1347 1011.16800 1007.95467 3.213333 1019.8097
## 50 1255.374 1132.2930 1266.53837 1212.84022 53.698148 1003.1347
## 51 999.123 1129.6073 1198.37789 1170.86967 27.508222 1132.2930
## 52 836.037 1030.1780 918.20893 962.99656 -44.787630 1129.6073
## 53 644.648 826.6027 545.16933 657.74267 -112.573333 1030.1780
## 54 584.343 688.3427 421.62304 528.31089 -106.687852 826.6027
## 55 679.350 636.1137 501.27033 555.20767 -53.937333 688.3427
## 56 733.838 665.8437 669.86089 668.25400 1.606889 636.1137
## 57 1058.262 823.8167 1015.85889 939.04200 76.816889 665.8437
## 58 1129.275 973.7917 1228.19333 1126.43267 101.760667 823.8167
## 59 995.150 1060.8957 1240.99733 1168.95667 72.040667 973.7917
## 60 1213.950 1112.7917 1218.84500 1176.42367 42.421333 1060.8957
## 61 1660.612 1289.9040 1515.52659 1425.27756 90.249037 1112.7917
## 62 1001.212 1291.9247 1392.56559 1352.30922 40.256370 1289.9040
## 63 941.050 1200.9580 1101.00652 1140.98711 -39.980593 1291.9247
## 64 647.650 863.3040 437.59585 607.87911 -170.283259 1200.9580
## 65 703.562 764.0873 466.26104 585.39156 -119.130519 863.3040
## 66 610.000 653.7373 476.00585 547.09844 -71.092593 764.0873
## 67 649.800 654.4540 593.94474 618.14844 -24.203704 653.7373
## 68 518.100 592.6333 524.34185 551.65844 -27.316593 654.4540
## 69 984.480 717.4600 821.81148 780.07089 41.740593 592.6333
## 70 295.150 599.2433 537.23963 562.04111 -24.801481 717.4600
## 71 NA NA NA NA NA 599.2433
## 72 NA NA NA NA NA 512.4381
## 73 NA NA NA NA NA 487.6367
## 74 NA NA NA NA NA 462.8352
## 75 NA NA NA NA NA 438.0337
## 76 NA NA NA NA NA 413.2322
## 77 NA NA NA NA NA 388.4307
## 78 NA NA NA NA NA 363.6293
## 79 NA NA NA NA NA 338.8278
## 80 NA NA NA NA NA 314.0263
## 81 NA NA NA NA NA 289.2248
## 82 NA NA NA NA NA 264.4233
## 83 NA NA NA NA NA 239.6219
## 84 NA NA NA NA NA 214.8204#Plot time series
ts.plot(dt.gab[,1], xlab="Time Period ", ylab="IPM",
main= "DMA N=6 Data IPM")
points(dt.gab[,1])
points(dt.gab[,3])
points(dt.gab[,6])
lines(dt.gab[,3],col="green",lwd=2)
lines(dt.gab[,6],col="red",lwd=2)
legend("topleft",c("data aktual","data pemulusan","data peramalan"),
lty=8, col=c("black","green","red"), cex=0.8)
Selanjutnya akan dilihat keakuratan dari metode DMA
Menghitung nilai keakuratan
error.dma = data.ts-dt.ramal[1:length(data.ts)]
SSE.dma = sum(error.dma[7:length(data.ts)]^2)
MSE.dma = mean(error.dma[7:length(data.ts)]^2)
MAPE.dma = mean(abs((error.dma[7:length(data.ts)]/data.ts[7:length(data.ts)])*100))
akurasi.dma <- matrix(c(SSE.dma, MSE.dma, MAPE.dma))
row.names(akurasi.dma)<- c("SSE", "MSE", "MAPE")
colnames(akurasi.dma) <- c("Akurasi m = 6")
akurasi.dma## Akurasi m = 6
## SSE 8329966.6
## MSE 130155.7
## MAPE InfSelanjutnya akan digunakan metode Double Exponential Smoothing dengan cara sebagai berikut.
Pertama akan data akan dibagi menjadi data training dan data testing.
#membagi training dan testing
training<-data[1:56,1]
testing<-data[57:70,1]
#data time series
training.ts<-ts(training)
testing.ts<-ts(testing,start=57)
#eksplorasi data
plot(data.ts, col="red",main="Plot semua data")
points(data.ts)
plot(training.ts, col="blue",main="Plot semua data")
points(training.ts)
Selanjutnya akan dilakukan pemulusan dengan DES, kali ini langsung dicari lambda dan gamma optimum sebagai berikut. Nilai lambda dan gamma optimum dapat dilihat pada smoothing parameters alpha untuk nilai lambda dan beta untuk nilai gamma.
#Lamda dan gamma optimum
des.opt<- HoltWinters(training.ts, gamma = FALSE)
des.opt## Holt-Winters exponential smoothing with trend and without seasonal component.
##
## Call:
## HoltWinters(x = training.ts, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0.59612
## beta : 0.06118864
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 702.699765
## b -6.551854plot(des.opt)
legend("topleft", c("Data Aktual", "Peramalan", "Data Aktual", "Peramalan"), col = c("black", "red"),
lty = c(1,1))
#ramalan
ramalandesopt<- forecast(des.opt, h=14)
ramalandesopt## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 57 696.1479 234.50633 1157.789 -9.87196 1402.168
## 58 689.5961 143.33976 1235.852 -145.83086 1525.023
## 59 683.0442 55.51142 1310.577 -276.68438 1642.773
## 60 676.4924 -30.53269 1383.517 -404.80913 1757.794
## 61 669.9405 -115.69557 1455.577 -531.58619 1871.467
## 62 663.3886 -200.54618 1527.323 -657.88566 1984.663
## 63 656.8368 -285.46334 1599.137 -784.28690 2097.960
## 64 650.2849 -370.70947 1671.279 -911.19127 2211.761
## 65 643.7331 -456.47172 1743.938 -1038.88496 2326.351
## 66 637.1812 -542.88643 1817.249 -1167.57652 2441.939
## 67 630.6294 -630.05453 1891.313 -1297.42028 2558.679
## 68 624.0775 -718.05152 1966.207 -1428.53172 2676.687
## 69 617.5257 -806.93429 2041.986 -1560.99784 2796.049
## 70 610.9738 -896.74577 2118.693 -1694.88431 2916.832Selanjutnya akan dicari akurasi dari metode DES.
ssedes.train<-des.opt$SSE
msedes.train<-ssedes.train/length(training.ts)
sisaandes<-ramalandesopt$residuals
head(sisaandes)## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## [1] NA NA -346.4950 -331.9487 -302.9107 -290.9689mapedes.train <- sum(abs(sisaandes[3:length(training.ts)]/training.ts[3:length(training.ts)])*100)/length(training.ts)
akurasides.opt <- matrix(c(ssedes.train,msedes.train,mapedes.train))
row.names(akurasides.opt)<- c("SSE", "MSE", "MAPE")
colnames(akurasides.opt) <- c("Akurasi lamda dan gamma optimum")
akurasides.opt## Akurasi lamda dan gamma optimum
## SSE 7113984.1
## MSE 127035.4
## MAPE InfAkurasi data testing
selisihdesopt<-ramalandesopt$mean-testing.ts
selisihdesopt## Time Series:
## Start = 57
## End = 70
## Frequency = 1
## [1] -362.114089 -439.678942 -312.105796 -537.457649 -990.671503 -337.823356
## [7] -284.213210 2.634937 -59.828917 27.181230 -19.170624 105.977523
## [13] -366.954331 315.823815SSEtestingdesopt<-sum(selisihdesopt^2)
SSEtestingdesopt<-SSEtestingdesopt/length(testing.ts)
MAPEtestingdesopt<-sum(abs(selisihdesopt/testing.ts)*100)/length(testing.ts)
akurasiDesTesting <- matrix(c(SSEtestingdesopt,SSEtestingdesopt,MAPEtestingdesopt))
row.names(akurasiDesTesting)<- c("SSE", "MSE", "MAPE")
colnames(akurasiDesTesting) <- c("Akurasi lamda dan gamma optimum")
akurasiDesTesting## Akurasi lamda dan gamma optimum
## SSE 152669.3504
## MSE 152669.3504
## MAPE 32.3885Setelah didapatkan nilai akurasi untuk metode DMA dan DES, selanjutnya akan dibandingkan keakuratan antar metode keduanya.
cbind(akurasi.dma, akurasiDesTesting)## Akurasi m = 6 Akurasi lamda dan gamma optimum
## SSE 8329966.6 152669.3504
## MSE 130155.7 152669.3504
## MAPE Inf 32.3885Berdasarkan perbandingan akurasi tersebut, terlihat nilai SSE, MSE, dan MAPE metode DES lebih kecil dibandingkan dengan metode DMA. Oleh karena itu, metode peramalan dan pemulusan yang terbaik antara keduanya adalah dengan metode DES.
Setelah melakukan peramalan, data yang telah dimasukkan kemudian dieksplorasi. Eksplorasi pertama yang dilakukan adalah dengan menggunakan scatter plot.
GGally::ggpairs(data[sapply(data,is.numeric)])## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
Berdasarkan scatter plot di atas, terlihat adanya hubungan / korelasi
positif antara peubah tahun dengan nilai IPM, terlihat titik-titik pada
plot yang naik ke arah kanan atas. Hal tersebut juga diperkuat dengan
hasil perhitungan aplikasi R di mana didapatkan nilai
korelasi sebesar \(0.9966848\).
Setalah mengetahui adanya hubungan antar dua peubah, maka model regresi dapat ditentukan.
Regresi
#Pembuatan Model Regresi
#model regresi
data## Y X1 X2
## 1 878.030 354.0 50
## 2 1001.900 347.0 31
## 3 779.275 232.0 20
## 4 698.500 209.0 18
## 5 628.780 270.0 23
## 6 548.225 323.0 23
## 7 491.900 348.0 21
## 8 583.850 420.0 29
## 9 887.820 399.0 14
## 10 1856.815 472.0 30
## 11 723.800 489.0 19
## 12 1015.660 492.0 25
## 13 1044.240 463.0 24
## 14 953.252 243.0 9
## 15 1084.850 208.0 13
## 16 940.170 192.0 5
## 17 765.900 194.0 10
## 18 746.788 217.0 12
## 19 708.828 203.0 6
## 20 790.788 265.5 15
## 21 852.125 243.5 11
## 22 1574.335 222.0 8
## 23 1277.725 228.0 18
## 24 1258.349 286.0 28
## 25 1476.324 248.0 24
## 26 1224.862 239.0 20
## 27 1150.700 250.0 13
## 28 998.337 318.0 18
## 29 997.150 275.0 18
## 30 760.050 311.0 20
## 31 652.362 240.0 8
## 32 753.050 275.5 12
## 33 1118.699 307.0 18
## 34 1617.275 312.0 11
## 35 1062.686 246.0 27
## 36 1624.335 257.0 18
## 37 0.000 1.0 0
## 38 526.350 144.0 7
## 39 612.500 165.0 9
## 40 540.200 132.0 9
## 41 547.940 148.0 23
## 42 496.100 163.0 8
## 43 479.350 219.0 15
## 44 549.300 239.0 12
## 45 863.750 223.0 23
## 46 1184.350 226.0 15
## 47 867.899 192.0 15
## 48 1007.180 226.0 6
## 49 1134.325 229.0 11
## 50 1255.374 268.0 12
## 51 999.123 381.0 42
## 52 836.037 289.0 21
## 53 644.648 259.0 13
## 54 584.343 248.0 18
## 55 679.350 283.0 19
## 56 733.838 253.0 20
## 57 1058.262 263.0 12
## 58 1129.275 287.0 25
## 59 995.150 252.2 22
## 60 1213.950 254.0 27
## 61 1660.612 295.2 23
## 62 1001.212 249.4 12
## 63 941.050 301.4 19
## 64 647.650 299.4 22
## 65 703.562 265.8 26
## 66 610.000 193.0 25
## 67 649.800 250.6 20
## 68 518.100 237.0 26
## 69 984.480 227.8 16
## 70 295.150 86.0 7model<- lm(Y~X1+X2, data = data)
summary(model)##
## Call:
## lm(formula = Y ~ X1 + X2, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -543.27 -212.40 -73.32 164.21 739.59
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 502.8442 125.1660 4.017 0.000151 ***
## X1 1.5118 0.5543 2.728 0.008138 **
## X2 -0.3689 5.5591 -0.066 0.947287
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 318.9 on 67 degrees of freedom
## Multiple R-squared: 0.1401, Adjusted R-squared: 0.1144
## F-statistic: 5.456 on 2 and 67 DF, p-value: 0.00638Model yang dihasilkan adalah \[y_i=502.84+1.5118x_1-0.3689 x_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. 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%. Selanjutnya dapat dilihat juga nilai \(R^2=0.1401\). Artinya, sebesar 14.01% keragaman nilai Y dapat dijelaskan oleh peubah X1 dan X2. 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,70,1), sisaan, col = "steelblue", pch = 20, xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,70 ,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.95227, p-value = 0.009554ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan))##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: sisaan
## D = 0.10956, p-value = 0.3448
## alternative hypothesis: two-sidedBerdasarkan uji formal Saphiro-Wilk dan 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 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.1379, p-value = 3.844e-05
## alternative hypothesis: true autocorrelation is greater than 0Berdasarkan hasil DW Test, didapatkan nilai \(DW = 1.1379\) dan p-value = \(3.844e-05\). Berdasarkan tabel Durbin-Watson diperoleh nilai \(DL = 1.5245\) dan \(DU = 1.7028\). 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, data = data)
##
## number of interaction: 12
## rho 0.482981
##
## Durbin-Watson statistic
## (original): 1.13791 , p-value: 3.844e-05
## (transformed): 2.04588 , p-value: 5.753e-01
##
## coefficients:
## (Intercept) X1 X2
## 307.834851 2.340725 -1.483302Hasil keluaran model setelah dilakukan penanganan adalah sebagai berikut. \[y_i=307.83+ 2.340725 x_1 -1.483302x_2\] Hasil juga menunjukkan bahwa nilai DW dan p-value meningkat menjadi \(2.04588\) dan \(0,5753\). Nilai DW sudah berada pada rentang DU < DW < 4-DU atau \(1.7028 < DW < 2.2972\). 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.482981\). Nilai tersebut dapat diketahui dengan syntax berikut.
Rho optimum
rho<- modelCO$rho
rho## [1] 0.4829805Selanjutnya akan dilakukan transformasi secara manual dengan syntax berikut ini.
Transformasi Manual
Y.trans<- data$Y[-1]-data$Y[-70]*rho
X1.trans<- data$X1[-1]-data$X1[-70]*rho
X2.trans<- data$X2[-1]-data$X2[-70]*rho
modelCOmanual<- lm(Y.trans~X1.trans+X2.trans)
summary(modelCOmanual)##
## Call:
## lm(formula = Y.trans ~ X1.trans + X2.trans)
##
## Residuals:
## Min 1Q Median 3Q Max
## -936.48 -146.53 -30.32 113.02 763.25
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 159.1566 81.3864 1.956 0.054752 .
## X1.trans 2.3407 0.6637 3.527 0.000771 ***
## X2.trans -1.4833 5.6869 -0.261 0.795038
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 286.7 on 66 degrees of freedom
## Multiple R-squared: 0.1998, Adjusted R-squared: 0.1755
## F-statistic: 8.239 on 2 and 66 DF, p-value: 0.0006392Hasil 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[2]
b2 <- modelCOmanual$coefficients[3]
b0## (Intercept) X2.trans
## 307.834851 -2.868948b1## X1.trans
## 2.340725b2## X2.trans
## -1.483302Hasil 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]
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]
return(lm(y~x1+x2))
}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 6274825
## 2 0.2 5886758
## 3 0.3 5618208
## 4 0.4 5465836
## 5 0.5 5427998
## 6 0.6 5506714
## 7 0.7 5709098
## 8 0.8 6046586
## 9 0.9 6531333Pertama-tama akan dicari di mana kira-kira \(ρ\) yang menghasilkan SSE minimum. Pada hasil di atas terlihat \(ρ\) minimum ketika 0.5. 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.6.
Rho optimal di sekitar 0.5
rOpt <- seq(0.3,0.6, 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
## 184 0.483 5426323
## 183 0.482 5426329
## 185 0.484 5426329
## 182 0.481 5426346
## 186 0.485 5426347
## 181 0.480 5426375#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.478, y=5431759 , labels = "rho=0.478", cex = 0.8)
Perhitungan yang dilakukan aplikasi R menunjukkan bahwa
nilai \(ρ\) optimum, yaitu saat SSE
terkecil terdapat pada nilai \(ρ=0.478\). 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.4829805, model)
summary(modelHL)##
## Call:
## lm(formula = y ~ x1 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -936.48 -146.53 -30.32 113.02 763.25
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 159.1566 81.3864 1.956 0.054752 .
## x1 2.3407 0.6637 3.527 0.000771 ***
## x2 -1.4833 5.6869 -0.261 0.795038
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 286.7 on 66 degrees of freedom
## Multiple R-squared: 0.1998, Adjusted R-squared: 0.1755
## F-statistic: 8.239 on 2 and 66 DF, p-value: 0.0006392#Transformasi Balik
cat("y = ", coefficients(modelHL)[1]/(1-0.6454578), "+", coefficients(modelHL)[2],"x1",coefficients(modelHL)[3],"x2",sep = "")## y = 448.9074+2.340725x1-1.483302x2Setelah dilakukan tranformasi balik, didapatkan model dengan metode Hildreth-Lu sebagai berikut. \[y_i=-1062.032+0.5597492x_t\]
Deteksi autokorelasi
dwtest(modelHL)##
## Durbin-Watson test
##
## data: modelHL
## DW = 2.0459, p-value = 0.5753
## alternative hypothesis: true autocorrelation is greater than 0Hasil uji Durbin-Watson juga menunjukkan bawah nilai DW sebesar \(2.0598\) berada pada selang daerah tidak ada autokorelasi, yaitu pada rentang DU < DW < 4-DU atau \(1.7028 < DW < 2.2972\). Hal tersebut juga didukung oleh p-value sebesar \(0.5834\), 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`[-1]
sseModelawal## [1] 447.7984 6812722.5017sseModelCO <- anova(modelCOmanual)$`Sum Sq`[-1]
sseModelCO## [1] 5593.285 5426323.447sseModelHL <- anova(modelHL)$`Sum Sq`[-1]
sseModelHL## [1] 5593.284 5426323.447mseModelawal <- sseModelawal/length(data$Y)
mseModelCO <- sseModelCO/length(data$Y)
mseModelHL <- sseModelHL/length(data$Y)
akurasi <- matrix(c(sseModelawal,sseModelCO,sseModelHL,
mseModelawal,mseModelCO,mseModelHL),nrow=2,ncol=3,byrow = T)## Warning in matrix(c(sseModelawal, sseModelCO, sseModelHL, mseModelawal, : data
## length differs from size of matrix: [12 != 2 x 3]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 447.7984 6812722.502 5593.285
## MSE 5426323.4465 5593.284 5426323.447Berdasarkan hasil tersebut dapat diketahui bahwa hasil penanganan autokorelasi dengan metode Cochrane-Orcutt dan Hildreth-Lu memiliki SSE yang beda dan lebih baik dibandingkan model awal ketika autokorelasi masih terjadi.
Simpulan
Autokorelasi yang terdapat pada data N02BE terjadi akibat adanya korelasi di antara unsur penyusunnya. Indikator N02BE yang erat hubungannya dengan N05B dan N05C 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. Namun, autokorelasi tersebut dapat ditangani dengan metode Cochrane-Orcutt dan Hildreth-Lu. Kedua metode menghasilkan nilai SSE yang beda, Model Hildreth-Lu baik untuk digunakan karena memiliki nilai SSE terkecil .