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)
## Warning: package 'TTR' was built under R version 4.3.1
library(forecast)
## Warning: package 'forecast' was built under R version 4.3.1
## 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(zoo)
library(orcutt) #untuk membuat model regresi Cochrane-Orcutt
## Warning: package 'orcutt' was built under R version 4.3.1
library(HoRM) #untuk membuat model regresi Hildreth-Lu
## Warning: package 'HoRM' was built under R version 4.3.1
library(rio)
## Warning: package 'rio' was built under R version 4.3.1

Input Data

Data yang digunakan dalam kesempatan kali ini adalah data Pharma Sales yang berjumlah 70 baris.

data <- import('https://raw.githubusercontent.com/divanm/mpdw/main/pertemuan%202/salesmonthlyp2.csv')
data
##      N02BA    N02BE  N05B
## 1  152.100  878.030 354.0
## 2  177.000 1001.900 347.0
## 3  147.655  779.275 232.0
## 4  130.900  698.500 209.0
## 5  132.100  628.780 270.0
## 6  122.900  548.225 323.0
## 7  129.300  491.900 348.0
## 8  123.800  583.850 420.0
## 9  122.100  887.820 399.0
## 10 191.600 1856.815 472.0
## 11 142.700  723.800 489.0
## 12 111.124 1015.660 492.0
## 13 141.000 1044.240 463.0
## 14 131.830  953.252 243.0
## 15 133.800 1084.850 208.0
## 16 122.100  940.170 192.0
## 17 136.040  765.900 194.0
## 18 145.460  746.788 217.0
## 19 125.500  708.828 203.0
## 20 133.400  790.788 265.5
## 21 110.400  852.125 243.5
## 22 146.200 1574.335 222.0
## 23 145.900 1277.725 228.0
## 24 137.000 1258.349 286.0
## 25 172.500 1476.324 248.0
## 26 134.200 1224.862 239.0
## 27 148.400 1150.700 250.0
## 28 147.700  998.337 318.0
## 29 130.550  997.150 275.0
## 30 117.750  760.050 311.0
## 31 137.900  652.362 240.0
## 32 132.700  753.050 275.5
## 33 116.700 1118.699 307.0
## 34 160.150 1617.275 312.0
## 35 133.850 1062.686 246.0
## 36 132.400 1624.335 257.0
## 37   0.000    0.000   1.0
## 38  97.000  526.350 144.0
## 39 107.350  612.500 165.0
## 40 100.500  540.200 132.0
## 41  98.950  547.940 148.0
## 42 119.600  496.100 163.0
## 43  75.200  479.350 219.0
## 44  84.400  549.300 239.0
## 45  88.150  863.750 223.0
## 46 100.400 1184.350 226.0
## 47 104.450  867.899 192.0
## 48 115.150 1007.180 226.0
## 49 101.150 1134.325 229.0
## 50 114.650 1255.374 268.0
## 51 122.300  999.123 381.0
## 52  84.600  836.037 289.0
## 53  89.400  644.648 259.0
## 54  86.800  584.343 248.0
## 55  87.200  679.350 283.0
## 56  88.250  733.838 253.0
## 57  86.500 1058.262 263.0
## 58  76.050 1129.275 287.0
## 59 102.150  995.150 252.2
## 60  84.750 1213.950 254.0
## 61  99.700 1660.612 295.2
## 62 110.200 1001.212 249.4
## 63  83.350  941.050 301.4
## 64  88.100  647.650 299.4
## 65 104.100  703.562 265.8
## 66 103.200  610.000 193.0
## 67  92.800  649.800 250.6
## 68  84.200  518.100 237.0
## 69  93.500  984.480 227.8
## 70  20.650  295.150  86.0

Eksplorasi Data

Sebelum melakukan regresi, akan diperlihatkan plot time-series dari kadar no2 dan so2 dalam air quality new delhi.

#Membentuk objek time series
data1.ts<-ts(data$N02BE)
data1.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
#Membentuk objek time series
data2.ts<-ts(data$N05B)
data2.ts
## Time Series:
## Start = 1 
## End = 70 
## Frequency = 1 
##  [1] 354.0 347.0 232.0 209.0 270.0 323.0 348.0 420.0 399.0 472.0 489.0 492.0
## [13] 463.0 243.0 208.0 192.0 194.0 217.0 203.0 265.5 243.5 222.0 228.0 286.0
## [25] 248.0 239.0 250.0 318.0 275.0 311.0 240.0 275.5 307.0 312.0 246.0 257.0
## [37]   1.0 144.0 165.0 132.0 148.0 163.0 219.0 239.0 223.0 226.0 192.0 226.0
## [49] 229.0 268.0 381.0 289.0 259.0 248.0 283.0 253.0 263.0 287.0 252.2 254.0
## [61] 295.2 249.4 301.4 299.4 265.8 193.0 250.6 237.0 227.8  86.0
#Membuat plot time series
ts.plot(data1.ts, xlab="Time Period ", ylab="N02BE", main= "Time Series Plot of N02BE")
points(data1.ts)

#Membuat plot time series
ts.plot(data2.ts, xlab="Time Period ", ylab="N05B", main= "Time Series Plot of N05B")
points(data2.ts)

Regresi

#Pembuatan Model Regresi
#model regresi
model<- lm(N02BA~N02BE+N05B, data = data)
summary(model)
## 
## Call:
## lm(formula = N02BA ~ N02BE + N05B, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -51.217 -15.626   1.561  17.699  47.274 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 48.097810  10.718313   4.487 2.91e-05 ***
## N02BE        0.038773   0.009449   4.104 0.000113 ***
## N05B         0.123291   0.037639   3.276 0.001671 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24.66 on 67 degrees of freedom
## Multiple R-squared:  0.395,  Adjusted R-squared:  0.377 
## F-statistic: 21.88 on 2 and 67 DF,  p-value: 4.871e-08

Model yang dihasilkan adalah \[y_i=23.988 + 75.765 x_1-7.604 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 ketiga parameter regresi, yaitu intersep, b1, dan b2 menunjukan memiliki p-value < \(\alpha\) (5%) sehingga nyata dalam taraf 5%. Selanjutnya dapat dilihat juga nilai \(R^2=0.395\). Artinya, sebesar 39.5% keragaman nilai N02BA dapat dijelaskan oleh peubah N02BE dan N05B. 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,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.97573, p-value = 0.1914
ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan))
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.089502, p-value = 0.5974
## alternative hypothesis: two-sided

Berdasarkan 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 = 0.72364, p-value = 2.492e-10
## alternative hypothesis: true autocorrelation is greater than 0

Berdasarkan 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 = N02BA ~ N02BE + N05B, data = data)
## 
##  number of interaction: 5
##  rho 0.645458
## 
## Durbin-Watson statistic 
## (original):    0.72364 , p-value: 2.492e-10
## (transformed): 2.16141 , p-value: 7.573e-01
##  
##  coefficients: 
## (Intercept)       N02BE        N05B 
##   41.756193    0.045601    0.115288

Hasil keluaran model setelah dilakukan penanganan adalah sebagai berikut. \[y_i=41.756193+0.045601x_1+0.115288x_2\] Hasil juga menunjukkan bahwa nilai DW dan p-value meningkat menjadi \(2.16141\). Nilai DW sudah berada pada rentang DU < DW < 4-DU atau \(1.331 < DW < 2.669\). 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.645458\). Nilai tersebut dapat diketahui dengan syntax berikut.

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

Selanjutnya akan dilakukan transformasi secara manual dengan syntax berikut ini.

#Transformasi Manual
N02BA.trans <- data$N02BA[-1]-data$N02BA[-70]*rho
N02BE.trans<- data$N02BE[-1]-data$N02BE[-70]*rho
N05B.trans<- data$N05B[-1]-data$N05B[-70]*rho
modelCOmanual<- lm(N02BA.trans~N02BE.trans+N05B.trans)
summary(modelCOmanual)
## 
## Call:
## lm(formula = N02BA.trans ~ N02BE.trans + N05B.trans)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -41.131 -12.275   0.901  11.633  41.667 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 14.804333   4.239661   3.492 0.000861 ***
## N02BE.trans  0.045601   0.008079   5.644 3.77e-07 ***
## N05B.trans   0.115288   0.044074   2.616 0.011023 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.12 on 66 degrees of freedom
## Multiple R-squared:  0.5054, Adjusted R-squared:  0.4905 
## F-statistic: 33.73 on 2 and 66 DF,  p-value: 8.111e-11

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[1]
b0 <- b0bintang/(1-rho)
b1 <- modelCOmanual$coefficients[2]
b2 <- modelCOmanual$coefficients[3]
b0
## (Intercept) 
##    41.75619
b1
## N02BE.trans 
##  0.04560055
b2
## N05B.trans 
##  0.1152876

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]
  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 35519.55
## 2 0.2 31749.73
## 3 0.3 28729.02
## 4 0.4 26462.73
## 5 0.5 24957.54
## 6 0.6 24220.06
## 7 0.7 24255.36
## 8 0.8 25066.22
## 9 0.9 26653.18

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

#Rho optimal di sekitar 0.4
rOpt <- seq(0.5,0.7, by= 0.001)
tabOpt <- data.frame("rho" = rOpt, "SSE" = sapply(rOpt, function(i){deviance(hildreth.lu.func(i, model))}))
tabOpt[order(tabOpt$SSE),]
##       rho      SSE
## 146 0.645 24140.09
## 147 0.646 24140.09
## 145 0.644 24140.16
## 148 0.647 24140.17
## 144 0.643 24140.31
## 149 0.648 24140.33
## 143 0.642 24140.54
## 150 0.649 24140.57
## 142 0.641 24140.85
## 151 0.650 24140.88
## 141 0.640 24141.23
## 152 0.651 24141.27
## 140 0.639 24141.69
## 153 0.652 24141.74
## 139 0.638 24142.23
## 154 0.653 24142.28
## 138 0.637 24142.85
## 155 0.654 24142.91
## 137 0.636 24143.54
## 156 0.655 24143.61
## 136 0.635 24144.31
## 157 0.656 24144.38
## 135 0.634 24145.16
## 158 0.657 24145.24
## 134 0.633 24146.09
## 159 0.658 24146.17
## 133 0.632 24147.09
## 160 0.659 24147.18
## 132 0.631 24148.17
## 161 0.660 24148.27
## 131 0.630 24149.33
## 162 0.661 24149.44
## 130 0.629 24150.57
## 163 0.662 24150.68
## 129 0.628 24151.88
## 164 0.663 24152.00
## 128 0.627 24153.27
## 165 0.664 24153.40
## 127 0.626 24154.74
## 166 0.665 24154.87
## 126 0.625 24156.28
## 167 0.666 24156.43
## 125 0.624 24157.91
## 168 0.667 24158.06
## 124 0.623 24159.61
## 169 0.668 24159.76
## 123 0.622 24161.38
## 170 0.669 24161.55
## 122 0.621 24163.24
## 171 0.670 24163.41
## 121 0.620 24165.17
## 172 0.671 24165.35
## 120 0.619 24167.18
## 173 0.672 24167.37
## 119 0.618 24169.27
## 174 0.673 24169.47
## 118 0.617 24171.43
## 175 0.674 24171.64
## 117 0.616 24173.67
## 176 0.675 24173.89
## 116 0.615 24175.99
## 177 0.676 24176.22
## 115 0.614 24178.39
## 178 0.677 24178.62
## 114 0.613 24180.86
## 179 0.678 24181.11
## 113 0.612 24183.41
## 180 0.679 24183.67
## 112 0.611 24186.04
## 181 0.680 24186.31
## 111 0.610 24188.75
## 182 0.681 24189.02
## 110 0.609 24191.53
## 183 0.682 24191.82
## 109 0.608 24194.39
## 184 0.683 24194.69
## 108 0.607 24197.33
## 185 0.684 24197.63
## 107 0.606 24200.34
## 186 0.685 24200.66
## 106 0.605 24203.44
## 187 0.686 24203.76
## 105 0.604 24206.61
## 188 0.687 24206.95
## 104 0.603 24209.85
## 189 0.688 24210.20
## 103 0.602 24213.18
## 190 0.689 24213.54
## 102 0.601 24216.58
## 191 0.690 24216.95
## 101 0.600 24220.06
## 192 0.691 24220.45
## 100 0.599 24223.61
## 193 0.692 24224.01
## 99  0.598 24227.25
## 194 0.693 24227.66
## 98  0.597 24230.96
## 195 0.694 24231.38
## 97  0.596 24234.75
## 196 0.695 24235.19
## 96  0.595 24238.61
## 197 0.696 24239.07
## 95  0.594 24242.55
## 198 0.697 24243.02
## 94  0.593 24246.57
## 199 0.698 24247.06
## 93  0.592 24250.67
## 200 0.699 24251.17
## 92  0.591 24254.85
## 201 0.700 24255.36
## 91  0.590 24259.10
## 90  0.589 24263.43
## 89  0.588 24267.83
## 88  0.587 24272.32
## 87  0.586 24276.88
## 86  0.585 24281.52
## 85  0.584 24286.23
## 84  0.583 24291.02
## 83  0.582 24295.89
## 82  0.581 24300.84
## 81  0.580 24305.86
## 80  0.579 24310.97
## 79  0.578 24316.14
## 78  0.577 24321.40
## 77  0.576 24326.73
## 76  0.575 24332.14
## 75  0.574 24337.63
## 74  0.573 24343.19
## 73  0.572 24348.84
## 72  0.571 24354.56
## 71  0.570 24360.35
## 70  0.569 24366.22
## 69  0.568 24372.17
## 68  0.567 24378.20
## 67  0.566 24384.31
## 66  0.565 24390.49
## 65  0.564 24396.75
## 64  0.563 24403.08
## 63  0.562 24409.50
## 62  0.561 24415.99
## 61  0.560 24422.56
## 60  0.559 24429.20
## 59  0.558 24435.92
## 58  0.557 24442.72
## 57  0.556 24449.60
## 56  0.555 24456.55
## 55  0.554 24463.58
## 54  0.553 24470.69
## 53  0.552 24477.87
## 52  0.551 24485.13
## 51  0.550 24492.47
## 50  0.549 24499.89
## 49  0.548 24507.38
## 48  0.547 24514.95
## 47  0.546 24522.60
## 46  0.545 24530.32
## 45  0.544 24538.12
## 44  0.543 24546.00
## 43  0.542 24553.95
## 42  0.541 24561.99
## 41  0.540 24570.09
## 40  0.539 24578.28
## 39  0.538 24586.54
## 38  0.537 24594.88
## 37  0.536 24603.30
## 36  0.535 24611.79
## 35  0.534 24620.36
## 34  0.533 24629.01
## 33  0.532 24637.74
## 32  0.531 24646.54
## 31  0.530 24655.42
## 30  0.529 24664.37
## 29  0.528 24673.41
## 28  0.527 24682.52
## 27  0.526 24691.70
## 26  0.525 24700.97
## 25  0.524 24710.31
## 24  0.523 24719.72
## 23  0.522 24729.22
## 22  0.521 24738.79
## 21  0.520 24748.44
## 20  0.519 24758.16
## 19  0.518 24767.97
## 18  0.517 24777.84
## 17  0.516 24787.80
## 16  0.515 24797.83
## 15  0.514 24807.94
## 14  0.513 24818.13
## 13  0.512 24828.39
## 12  0.511 24838.73
## 11  0.510 24849.15
## 10  0.509 24859.64
## 9   0.508 24870.21
## 8   0.507 24880.86
## 7   0.506 24891.59
## 6   0.505 24902.39
## 5   0.504 24913.26
## 4   0.503 24924.22
## 3   0.502 24935.25
## 2   0.501 24946.36
## 1   0.500 24957.54
#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.6454578, y=26672.12, labels = "rho=0.6454578", cex = 0.8)

Perhitungan yang dilakukan aplikasi R menunjukkan bahwa nilai \(ρ\) optimum, yaitu saat SSE terkecil terdapat pada nilai \(ρ=0.341\). 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.6454578, model)
modelHL
## 
## Call:
## lm(formula = y ~ x1 + x2)
## 
## Coefficients:
## (Intercept)           x1           x2  
##     14.8043       0.0456       0.1153
summary(modelHL)
## 
## Call:
## lm(formula = y ~ x1 + x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -41.131 -12.275   0.901  11.633  41.667 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 14.804332   4.239661   3.492 0.000861 ***
## x1           0.045601   0.008079   5.644 3.77e-07 ***
## x2           0.115288   0.044074   2.616 0.011023 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.12 on 66 degrees of freedom
## Multiple R-squared:  0.5054, Adjusted R-squared:  0.4905 
## F-statistic: 33.73 on 2 and 66 DF,  p-value: 8.111e-11
#Transformasi Balik
cat("y = ", coefficients(modelHL)[1]/(1-0.6454578), "+", coefficients(modelHL)[2],"x1","+",coefficients(modelHL)[3],"+","x2",sep = "")
## y = 41.75619+0.04560055x1+0.1152876+x2

Setelah dilakukan tranformasi balik, didapatkan model dengan metode Hildreth-Lu sebagai berikut. \[y_i=14.804332+0.045601x_1+0.115288x_2\]

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

Hasil uji Durbin-Watson juga menunjukkan bawah nilai DW sebesar \(2.1614\) berada pada selang daerah bukan autokorelasi, yaitu pada rentang DU < DW < 4-DU atau \(1.7028 < DW < 2.2972\). Hal tersebut juga didukung oleh p-value sebesar \(0.7573\), di mana p-value <>\(\alpha\)=5%. Artinya tak tolak \(H_0\) atau tidak cukup bukti menyatakan bahwa ada autokorelasi dalam data nilai N02BA 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]  6526.339 40752.559
sseModelCO <- anova(modelCOmanual)$`Sum Sq`[-1]
sseModelCO
## [1]  2502.62 24140.08
sseModelHL <- anova(modelHL)$`Sum Sq`[-1]
sseModelHL
## [1]  2502.62 24140.08
mseModelawal <- sseModelawal/length(data$N02BA)
mseModelCO <- sseModelCO/length(data$N02BA)
mseModelHL <- sseModelHL/length(data$N02BA)
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   6526.339              40752.56           2502.62
## MSE  24140.079               2502.62          24140.08

Berdasarkan hasil tersebut dapat diketahui bahwa hasil penanganan autokorelasi dengan metode Cochrane-Orcutt dan Hildreth-Lu memiliki SSE yang berbeda, sebesar \(40752.56\) dan \(2502.62\).

Simpulan

Autokorelasi yang terdapat pada data N02BA terjadi akibat adanya korelasi di antara unsur penyusunnya. Indikator N02BA yang erat hubungannya dengan N02BE dan N05B 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. Namun, autokorelasi tersebut dapat ditangani dengan metode Cochrane-Orcutt dan Hildreth-Lu. Kedua metode menghasilkan nilai SSE yang berbeda dengan metode Hildreth-Lu menghasilkan SSE lebih kecil sehingga lebih baik digunakan.