Prak02
Angga Fathan Rofiqy
2023-09-12
I. Setting
Abaikan bagian ini. Lansung saja loncat ke
1. Tentang Data . Function biar gak perlu
ganti backslash (\) jadi slash
(/).
path <- function() gsub ( "\\\\", "/", readClipboard () )
#Copy path, Panggil function di console
#Copy r path, paste ke var yang diinginkanFunction biar gak perlu repot buat
install("") atau load() package.
# -=( Install & Load Package Function )=-
install_load <- function (package1, ...) {
# convert arguments to vector
packages <- c(package1, ...)
# start loop to determine if each package is installed
for(package in packages){
# if package is installed locally, load
if(package %in% rownames(installed.packages()))
do.call('library', list(package))
# if package is not installed locally, download, then load
else {
install.packages(package)
do.call("library", list(package))
}
}
}Variabel untuk warna pada function box
plot.
# -=( Default Color Palette )=-
colord <- data.frame(
f = c('#75bfc8', 'Fill color'), c = c('#444442', 'Stroke color'),
oc = c('#0b5f6a', 'Outlier color'), med = c('#c3ac8b', 'Median value color'),
iqr = c('#c3ac8b', 'IQR3 Value'), mean = c('#608981', 'Mean Marker'),
v = c('#c3eca2', 'Violin'),
stringsAsFactors = FALSE
)Variabel untuk ggplot theme.
install_load('ggplot2')## Warning: package 'ggplot2' was built under R version 4.2.3theme1 <- list(
guides(fill="none") , #No Legends
theme(
text = element_text(size = 66),
axis.title = element_text(),
axis.text.x = element_text(hjust = 0.5, face = "bold", margin = margin(b=50)),
plot.title = element_text(hjust = 0.5, face = "bold", margin = margin(b=50)),
plot.subtitle = element_text(hjust = 0.5),
panel.background = element_rect(fill = 'transparent'),
plot.background = element_rect(fill='transparent', color=NA),
panel.grid.major = element_line(colour = "grey50"),
axis.line = element_line(linewidth = 5, colour = "grey80"))
)Box Plot Function.
# -=( Box Plot Function )=-
boxp <- function(dt, cp=colord, labs=lab1){ #data, color palette, label
install_load('ggplot2')
b <- ggplot(data=dt, aes(y="", x=y)) +
#Violin Plot
geom_violin(scale="count", fill=cp$v[1], color=NA, alpha=0.5,
trim = FALSE) +
#Box Plot
geom_boxplot(fill=cp$f[1], color=cp$c[1],
outlier.size=15, outlier.color=cp$oc[1],
alpha=0.7,notch=T, width = 0.25) +
#Median Value
stat_summary(geom="text", fun=median,
aes(label=sprintf("%1.1f", after_stat(x))),
position=position_nudge(y=0.425), size=16, color=cp$med[1]) +
#IQR1 Value
stat_summary(geom="text", fun="quantile", fun.args=list(probs=0.25),
aes(label=sprintf("%1.1f", after_stat(x))),
position=position_nudge(y=-0.425), size=16, color=cp$iqr[1]) +
#IQR3 Value
stat_summary(geom="text", fun="quantile", fun.args=list(probs=0.75),
aes(label=sprintf("%1.1f", after_stat(x))),
position=position_nudge(y=-0.425), size=16, color=cp$iqr[1]) +
#Mean Marker
stat_summary(fun='mean', geom="point",colour=cp$mean[1],
shape=18, size=16, alpha=0.75) +
#Title
ggtitle(labs$title) +
theme(plot.title = element_text(hjust = 0.5)) + #Title Position
#Labels
ylab(labs$ylab) + xlab(labs$xlab) +
#Mean Value
geom_text(data = data.frame(mean = round(mean(dt$y), 2))
, aes(label = paste0("Mean : ",mean), x = mean),
angle=45, size=16, vjust=-0.5, hjust=0.5, color=cp$mean[1]) +
#Theme
theme1
return(b)
}Default label untuk function box
plot.
# -=( Default Label/Title )=-
lab1 <- data.frame(title='Distribution of Data X', ylab='', xlab='Data X')A. Tentang Data
Dataset ini dibangun dari dataset awal yang terdiri dari \(600.000\) data transaksi yang dikumpulkan dalam waktu 6 tahun (periode 2014-2019), yang mencatat tanggal dan waktu penjualan, nama merek obat farmasi, dan jumlah yang terjual. Data ini diekspor dari sistem Point-of-Sale di apotek individu. Kelompok obat yang dipilih dari dataset (57 obat) diklasifikasikan ke dalam kategori berikut dalam Sistem Klasifikasi Anatomi Terapeutik Kimia (ATC):
- M01AB - Produk antiinflamasi dan antirheumatik, non-steroid, Turunan asam asetat dan substansi terkait
- M01AE - Produk antiinflamasi dan antirheumatik, non-steroid, Turunan asam propionat
- N02BA - Analgesik dan antipiretik lainnya, Asam salisilat dan turunannya N02BE/B -
- Analgesik dan antipiretik lainnya, Pirazolona dan Anilida
- N05B - Obat psikoleptik, Obat ansiolitik
- N05C - Obat psikoleptik, Obat hipnotik dan sedatif
- R03 - Obat untuk penyakit saluran napas obstruktif
- R06 - Antihistamin untuk penggunaan sistemik
Data penjualan diambil ulang dalam periode per jam, per hari, per minggu, dan per bulan. Data sudah diproses sebelumnya, termasuk deteksi dan penanganan outlier serta pengisian data yang hilang.
Data yang akan saya gunakan adalah data dalam periode per bulan. Obat
M01AB dan M01AE merupakan obat yang sama,
namun dengan bahan dasar yang berbeda. Saya penasaran “Apakah
penjualan obat M01AB dipengaruhi oleh penjualan obat
M01AB dan obat-obat lainnya?“ Sehingga saya memilih
Peubah M01AB sebagai peubah respon Y. Dan
peubah lainnya sebagai peubah penjelas X.
Import Data
install_load('rio')
raw.data <- import("https://raw.githubusercontent.com/Zen-Rofiqy/STA1341-MPDW/main/Pertemuan%202/salesmonthly.csv")Assign Variable
install_load('dplyr')## Warning: package 'dplyr' was built under R version 4.2.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, uniondata <- raw.data[,2:9] %>%
rename_all(~ c('y', 'x1', 'x2', 'x3', 'x4', 'x5', 'x6', 'x7'))B. Explorasi
Sebaran Peubah Respon
boxp( dt = data,
labs = data.frame(title='Sebaran Peubah y', ylab='', xlab='Obat M01AB')
)
Matriks Korelasi
install_load('psych')## Warning: package 'psych' was built under R version 4.2.3##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alphapairs.panels(data,
method = "pearson", # correlation method
hist.col = "#00AFBB", #Coloring histogram
density = TRUE, # show density plots
ellipses = TRUE, # show correlation ellipses
smooth = TRUE, #show loess smooths
pch = 20, #Scatter = cirlce / dot
rug = TRUE, #Rug under histogram
stars = TRUE #Significance of corr
)
C. Regresi
Model Regresi
model <- lm(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7, data)
summary(model)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -55.649 -11.046 0.109 13.309 62.473
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 51.26828 15.04508 3.408 0.00116 **
## x1 0.45722 0.16379 2.792 0.00696 **
## x2 0.16184 0.14071 1.150 0.25451
## x3 0.01414 0.01764 0.802 0.42574
## x4 -0.02988 0.04765 -0.627 0.53296
## x5 -0.07883 0.43916 -0.179 0.85813
## x6 0.01912 0.05446 0.351 0.72666
## x7 0.23362 0.07338 3.184 0.00228 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.05 on 62 degrees of freedom
## Multiple R-squared: 0.4758, Adjusted R-squared: 0.4166
## F-statistic: 8.04 on 7 and 62 DF, p-value: 6.328e-07Model yang dihasilkan adalah \(\hat{y}= 51.27 + 0.46x_1 + 0.16x_2 + 0.01x_3 - 0.03x_4 - 0.08x_5 + 0.02x_6 + 0.023x_7\) . Berdasarkan ringkasan model dapat diketahui bahwa hasil uji F memiliki \(p.value ( 6.328 \times10^{-7} ) < \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.4758\). Artinya, sebesar \(47.58\%\) keragaman penjualan obat
M01AB dapat dijelaskan oleh penjualan 7 obat lainnya.
Selanjutnya kita perlu melakukan uji terhadap sisaannya seperti berikut
ini.
Plot Sisaan & Q-Q Plot
sisaan<- residuals(model)
fitValue<- predict(model)
#Diagnostik dengan eksploratif
par(mfrow = c(2,2)) #Buat 4 plot
#Q-Q Plot
qqnorm(sisaan)
qqline(sisaan, col = "steelblue", lwd = 2)
#Sisaan vs Fitted Values
plot(fitValue, sisaan, col = "steelblue", pch = 20,
xlab = "Sisaan", ylab = "Fitted Values", main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)
#Histogram
hist(sisaan, col = "steelblue")
#Sisaan vs Order
plot(1:length(sisaan), sisaan, col = "steelblue", pch = 20,
xlab = "Sisaan",
ylab = "Order", main = "Sisaan vs Order")
lines(1:length(sisaan), 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 dan histogram dari sisaan menunjukkan hal yang serupa. Selanjutnya, dua plot di samping kanan digunakan untuk melihat autokorelasi. Terlihat ada pola pada sisaan yang ditunjukkan oleh Plot Sisaan vs Fitted Value dan terutama pada Plot Sisaan vs Order yang terlihat dengan garis warna merahnya yang membentuk pola. Ini merupakan tanda adanya autokorelasi.
Dengan ggplot2
install_load('ggplot2','cowplot')## Warning: package 'cowplot' was built under R version 4.2.3# Membuat QQ Plot
qq_plot <- ggplot(data, aes(sample = y)) +
geom_qq(fill=NA, alpha=0.5) +
geom_qq_line(color = "steelblue") +
labs(x = "Theoretical Quantiles", y = "Sample Quantiles") +
ggtitle("QQ Plot") + theme_minimal()
# Membuat plot sisaan vs Fitted Values dengan ggplot2
plot1 <- ggplot(data, aes(x = fitted(model), y = residuals(model))) +
geom_point(col = "steelblue", pch = 20) +
geom_abline(intercept = 0, slope = 0, color = "black",
size = 0.75, alpha=0.5) +
labs(x = "Fitted Values", y = "Residuals") +
ggtitle("Residuals vs Fitted Values") + theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.# Membuat histogram sisaan dengan ggplot2
plot2 <- ggplot(data, aes(x = residuals(model))) +
geom_histogram(binwidth = 10, fill = "steelblue", color = "black") +
labs(x = "Residuals", y = "Frequency") +
ggtitle("Histogram of Residuals") + theme_minimal()
# Membuat plot sisaan vs Order dengan ggplot2
plot3 <- ggplot(data, aes(x = 1:length(residuals(model)), y = residuals(model))) +
geom_point(col = "steelblue", pch = 20) +
geom_abline(intercept = 0, slope = 0, color = "black",
size = 1, alpha=0.5) +
geom_line(aes(x = 1:length(residuals(model)),
y = residuals(model)), color = "red", alpha=0.5) +
labs(x = "Order", y = "Residuals") +
ggtitle("Residuals vs Order") + theme_minimal()
# Menggabungkan semua plot dengan cowplot
plot_grid(qq_plot, plot1, plot2, plot3, nrow = 2, ncol = 2)
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
\(H_0\) : Sisaan mengikuti sebaran normal
\(H_1\) : Sisaan tidak mengikuti sebaran normal
shapiro.test(sisaan)##
## Shapiro-Wilk normality test
##
## data: sisaan
## W = 0.99006, p-value = 0.8574ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan))##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: sisaan
## D = 0.083832, p-value = 0.6774
## 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 hanya ada satu yang signifikan (terindikasi adanya autokorelasi) pada plot ACF dan ada satu pada Plot PACF Namun, untuk lebih memastikan akan dilakukan uji formal dengan uji Durbin Watson.
Deteksi autokorelasi dengan uji-Durbin Watson
\(H_0\) : Tidak ada autokorelasi
\(H_1\) : Ada autokorelasi
install_load('randtests')
runs.test(model$residuals)##
## Runs Test
##
## data: model$residuals
## statistic = -1.6856, runs = 29, n1 = 35, n2 = 35, n = 70, p-value =
## 0.09188
## alternative hypothesis: nonrandomnessinstall_load('lmtest')## Warning: package 'lmtest' was built under R version 4.2.3## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericdwtest(model)##
## Durbin-Watson test
##
## data: model
## DW = 1.2495, p-value = 0.0001048
## alternative hypothesis: true autocorrelation is greater than 0Berdasarkan hasil DW Test, didapatkan nilai \(DW = 1.2495\) dan \(p.value = 0.0001048\). Dengan nilai \(p.value < 0.05\) dapat disimpulkan bahwa tolak \(H_0\), cukup bukti mengatakan adanya autokorelasi. Oleh karena itu, diperlukan penanganan autokorelasi. Penanganan yang akan digunakan menggunakan dua metode, yaitu Cochrane-Orcutt dan Hildret-Lu.
D. 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.
install_load('orcutt')## Warning: package 'orcutt' was built under R version 4.2.3modelCO <- cochrane.orcutt(model)
modelCO## Cochrane-orcutt estimation for first order autocorrelation
##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7, data = data)
##
## number of interaction: 28
## rho 0.623643
##
## Durbin-Watson statistic
## (original): 1.24950 , p-value: 1.048e-04
## (transformed): 1.91190 , p-value: 3.785e-01
##
## coefficients:
## (Intercept) x1 x2 x3 x4 x5
## 26.693862 0.318248 0.378026 -0.000166 0.097630 -0.029786
## x6 x7
## 0.020880 0.175696Hasil keluaran model setelah dilakukan penanganan adalah sebagai berikut. \(y_i= 26.693862 + 0.318248x_1 + 0.378026 - 0.000166x_3 + 0.097630x_4 - 0.029786x_5 + 0.020880x_6 + 0.175696x_7\) Hasil juga menunjukkan bahwa nilai DW dan p-value meningkat dari \(1.048 \times 10^{-4}\) menjadi \(0.3785\), dengan nilai \(p.value > 0.05\), artinya belum cukup bukti menyatakan bahwa sisaan terdapat autokorelasi pada taraf nyata \(5\%\).
Rho optimum
rho <- modelCO$rho
rho## [1] 0.6236433Untuk nilai \(\hat{\rho}\) optimum yang digunakan adalah \(0.6236433\).
Transformasi Manual
Selanjutnya akan dilakukan transformasi secara manual dengan :
data tanpa amatan pertama -
data tanpa amatan terakhir \(\times\) rho
trans <- data[-1,] - data[-nrow(data),] * rho
modelCOmanual <- lm(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7, trans)
summary(modelCOmanual)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7, data = trans)
##
## Residuals:
## Min 1Q Median 3Q Max
## -46.159 -13.079 -1.617 12.868 47.550
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.0464144 5.8607533 1.714 0.0916 .
## x1 0.3182475 0.1376392 2.312 0.0242 *
## x2 0.3780261 0.1476258 2.561 0.0129 *
## x3 -0.0001661 0.0148848 -0.011 0.9911
## x4 0.0976296 0.0611365 1.597 0.1155
## x5 -0.0297857 0.4130362 -0.072 0.9427
## x6 0.0208799 0.0531617 0.393 0.6959
## x7 0.1756964 0.0785588 2.236 0.0290 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.52 on 61 degrees of freedom
## Multiple R-squared: 0.5823, Adjusted R-squared: 0.5344
## F-statistic: 12.15 on 7 and 61 DF, p-value: 1.33e-09Hasil model transformasi bukan merupakan model sesungguhnya. Koefisien regresi masih perlu dicari kembali mengikuti \({\beta_0}^*=\beta_0+\hat{\rho}\beta_0\) dan untuk koefisen selanjutnya mengikuti \(\beta_1^*=\beta_1\), \(\beta_2^*=\beta_2\), \(\beta_3^*=\beta_3\), \(\beta_4^*=\beta_4\), \(\beta_5^*=\beta_5\), \(\beta_6^*=\beta_6\), \(\beta_7^*=\beta_7\).
Mencari Penduga Koefisien Regresi setelah Transformasi ke Persamaan Awal
#Mencari Penduga Koefisien Regresi setelah Transformasi ke Persamaan Awal
b0bintang <- coef(modelCOmanual) %>%
as.data.frame() %>%
t() %>%
as_tibble() %>%
setNames(c("b0", "b1", "b2", "b3", "b4", "b5", "b6", "b7"))
b0bintang$b0 <- b0bintang$b0/(1 - rho)
install_load('knitr','kableExtra')##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rowskable(b0bintang) %>% kable_styling()| b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 |
|---|---|---|---|---|---|---|---|
| 26.69386 | 0.3182475 | 0.3780261 | -0.0001661 | 0.0976296 | -0.0297857 | 0.0208799 | 0.1756964 |
Hasil perhitungan koefisien regresi tersebut akan menghasilkan hasil yang sama dengan model yang dihasilkan menggunakan packages.
E. 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.func<- function(r, model){
x1 <- model.matrix(model)[,2]
x2 <- model.matrix(model)[,3]
x3 <- model.matrix(model)[,4]
x4 <- model.matrix(model)[,5]
x5 <- model.matrix(model)[,6]
x6 <- model.matrix(model)[,7]
x7 <- model.matrix(model)[,8]
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]
x4 <- x4[t]-r*x4[t-1]
x5 <- x5[t]-r*x5[t-1]
x6 <- x6[t]-r*x6[t-1]
x7 <- x7[t]-r*x7[t-1]
return(lm(y~x1+x2+x3+x4+x5+x6+x7))
}Pencarin 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)) ) )
kable(round(tab, 4)) %>% kable_styling()| rho | SSE |
|---|---|
| 0.1 | 33474.50 |
| 0.2 | 31622.80 |
| 0.3 | 30219.01 |
| 0.4 | 29213.27 |
| 0.5 | 28563.52 |
| 0.6 | 28265.17 |
| 0.7 | 28402.90 |
| 0.8 | 29189.53 |
| 0.9 | 30918.32 |
min(tab$SSE)## [1] 28265.17Pertama-tama akan dicari di mana kira-kira \(\rho\) yang menghasilkan SSE minimum. Pada hasil di atas terlihat \(\rho\) minimum ketika \(0.6\). Namun, hasil tersebut masih kurang teliti sehingga akan dicari kembali \(\rho\) yang lebih optimum dengan ketelitian yang lebih. Jika sebelumnya jarak antar \(\rho\) yang dicari adalah \(0.1\), kali ini jarak antar \(\rho\) 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))}))
kable(head(tabOpt[order(tabOpt$SSE),])) %>% kable_styling()| rho | SSE | |
|---|---|---|
| 101 | 0.500 | 28563.52 |
| 100 | 0.499 | 28568.33 |
| 99 | 0.498 | 28573.17 |
| 98 | 0.497 | 28578.04 |
| 97 | 0.496 | 28582.95 |
| 96 | 0.495 | 28587.89 |
min(tabOpt$SSE)## [1] 28563.52Grafik 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.500, y=28563.52, labels = "rho=0.500", cex = 0.8)
Perhitungan yang dilakukan aplikasi R menunjukkan bahwa
nilai \(\rho\) optimum, yaitu saat SSE
terkecil terdapat pada nilai \(\rho=0.500\) Hal tersebut juga ditunjukkan
pada plot. Selanjutnya, model dapat didapatkan dengan mengevaluasi nilai
\(\rho\) 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.500, model)
summary(modelHL)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.322 -11.988 -0.674 11.744 43.963
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.276334 7.546434 2.157 0.0350 *
## x1 0.351409 0.143260 2.453 0.0170 *
## x2 0.337828 0.147281 2.294 0.0253 *
## x3 0.001195 0.015469 0.077 0.9387
## x4 0.057134 0.058611 0.975 0.3335
## x5 0.012208 0.435671 0.028 0.9777
## x6 0.032262 0.053251 0.606 0.5469
## x7 0.189946 0.077786 2.442 0.0175 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.64 on 61 degrees of freedom
## Multiple R-squared: 0.5548, Adjusted R-squared: 0.5037
## F-statistic: 10.86 on 7 and 61 DF, p-value: 8.287e-09Transformasi Balik
cat("y = ", coef(modelHL)[1]/(1-0.500), "+", coef(modelHL)[2],"x", sep = "")## y = 32.55267+0.351409xSetelah 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 = 1.7683, p-value = 0.1576
## alternative hypothesis: true autocorrelation is greater than 0Hasil uji Durbin-Watson juga menunjukkan bawah nilai DW sebesar \(1.7683\) berada pada selang daerah tidak ada autokorelasi. Hal tersebut juga didukung oleh p-value sebesar \(0.1576\) , di mana \(p.value > \alpha (5\%).\) Artinya tak tolak \(H_0\) atau belum cukup bukti menyatakan bahwa ada autokorelasi dalam data nilai penjualan obat 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:7)]
sseModelCO <- anova(modelCOmanual)$`Sum Sq`[-c(1:7)]
sseModelHL <- anova(modelHL)$`Sum Sq`[-c(1:7)]
mseModelawal <- sseModelawal/nrow(data)
mseModelCO <- sseModelCO/nrow(data)
mseModelHL <- sseModelHL/nrow(data)
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 35854.4575 28252.6084 28563.5224
## MSE 512.2065 403.6087 408.0503Berdasarkan hasil tersebut dapat diketahui bahwa hasil penanganan autokorelasi dengan metode Cochrane-Orcutt dan Hildreth-Lu memiliki SSE yang sama, sebesar \(28252.6084\) dan lebih baik dibandingkan model awal ketika autokorelasi masih terjadi, yaitu sebesar \(28563.5224\).
II. Kesimpulan
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 sama, artinya keduanya baik untuk digunakan.