Menamakan Peubah

y <- datapharma[,9]
x1 <- datapharma[,2]
x2 <- datapharma[,3]
x3 <- datapharma[,4]
x4 <- datapharma[,5]

data <- cbind(y, x1, x2, x3, x4)
data <- as.data.frame(data)

Ringkasan data peubah respon

summary(y)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   49.88   74.10   86.66  119.81  213.04

Histogram

# Peubah Respon dan Peubah Penjelas

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.2.3

ggplot(data, aes(x=y)) + 
  geom_histogram(fill="purple", bins=15, color='black', alpha=0.7) +
    ggtitle("Sebaran Data Peubah y") + #Title
  theme(plot.title = element_text(hjust = 0.5)) + #Title Position +
  ylab("Frekuensi")+xlab("Penjualan Obat R06")

ggplot(data, aes(x=x1)) + 
  geom_histogram(fill="purple", bins=15, color='black', alpha=0.7) +
    ggtitle("Sebaran Data Peubah x1") + #Title
  theme(plot.title = element_text(hjust = 0.5)) + #Title Position +
  ylab("Frekuensi")+xlab("Penjualan Obat M01AB")

ggplot(data, aes(x=x2)) + 
  geom_histogram(fill="purple", bins=15, color='black', alpha=0.7) +
    ggtitle("Sebaran Data Peubah x2") + #Title
  theme(plot.title = element_text(hjust = 0.5)) + #Title Position +
  ylab("Frekuensi")+xlab("Penjualan Obat M01AE")

ggplot(data, aes(x=x3)) + 
  geom_histogram(fill="purple", bins=15, color='black', alpha=0.7) +
    ggtitle("Sebaran Data Peubah x3") + #Title
  theme(plot.title = element_text(hjust = 0.5)) + #Title Position +
  ylab("Frekuensi")+xlab("Penjualan Obat N02BA")

ggplot(data, aes(x=x4)) + 
  geom_histogram(fill="purple", bins=15, color='black', alpha=0.7) +
    ggtitle("Sebaran Data Peubah x4") + #Title
  theme(plot.title = element_text(hjust = 0.5)) + #Title Position +
  ylab("Frekuensi")+xlab("Penjualan Obat N02BE")

Matriks Korelasi

library(ggplot2)
library(GGally)

## Warning: package 'GGally' was built under R version 4.2.3

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

sample_data <- data.frame(y, x1, x2, x3, x4)
ggpairs(sample_data)

Berdasarkan matriks korelasi di atas, terlihat adanya hubungan / korelasi positif antara peubah y dengan peubah x1 sebesar 0.308 dan x2 sebesar 0.029, terlihat pada titik-titik plot yang naik ke arah kanan atas. Lalu terlihat pula adanya hubungan / korelasi negatif antara peubah y dengan peubah x3 sebesar -0.034 dan x4 sebesar -0.287, terlihat pada titik-titik plot yang turun ke arah kanan bawah.

Setalah mengetahui adanya hubungan antar dua peubah, maka model regresi dapat ditentukan.

Diagnostik secara Eksploratif

#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 sisaan vs urutan (sisaan saling bebas)
plot(x = 1:dim(data)[1],
     y = model$residual,
     type = 'b', 
     ylab = "Residual",
     xlab = "Observasi")

Dua plot di samping kiri digunakan untuk melihat apakah sisaan menyebar normal. Normal Q-Q Plot dan histogram dari sisaan di atas menunjukkan bahwa sisaan cenderung menyebar normal. 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.

Diagnostik secara Uji Formal

#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.9658, p-value = 0.05275

ks.test(sisaan, "pnorm", mean=mean(sisaan), sd=sd(sisaan))

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.09172, p-value = 0.5665
## 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 tidak semua dalam rentang batas dan ada yang tidak 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

#install.packages("randtests")
library(randtests)
runs.test(model$residuals)

## 
##  Runs Test
## 
## data:  model$residuals
## statistic = -2.1672, runs = 27, n1 = 35, n2 = 35, n = 70, p-value =
## 0.03022
## alternative hypothesis: nonrandomness

#install.packages("lmtest")
library(lmtest)
dwtest(model)

## 
##  Durbin-Watson test
## 
## data:  model
## DW = 1.1122, p-value = 1.49e-05
## alternative hypothesis: true autocorrelation is greater than 0

Berdasarkan hasil DW Test, didapatkan nilai \(DW = 1.1122\) dan p-value = \(0.0000149\). Dengan nilai p-value < 0.05 dapat disimpulkan bahwa tolak H0, cukup bukti mengatakan adanya autokorelasi. Oleh karena itu, diperlukan penanganan autokorelasi. Penanganan yang akan digunakan menggunakan dua metode, yaitu Cochrane-Orcutt dan Hildret-Lu.

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
library(orcutt)
modelCO <- cochrane.orcutt(model)
modelCO

## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4, data = datapharma)
## 
##  number of interaction: 12
##  rho 0.631268
## 
## Durbin-Watson statistic 
## (original):    1.11217 , p-value: 1.49e-05
## (transformed): 1.60645 , p-value: 5.659e-02
##  
##  coefficients: 
## (Intercept)          x1          x2          x3          x4 
##   38.381612    0.418933   -0.040813    0.096879   -0.023067

Hasil keluaran model setelah dilakukan penanganan adalah sebagai berikut. \[y=38.3816+0.418933x_1-0.040813x_2+0.096879x_3-0.023067x_4\] Hasil juga menunjukkan bahwa nilai DW dan p-value meningkat menjadi \(1.60645\) dan \(0.05659\). Nilai p-value > 0.05, artinya belum cukup bukti menyatakan bahwa sisaan terdapat autokorelasi pada taraf nyata 5%. Autokorelasi telah berhasil tertangani.

Untuk nilai \(ρ ̂\) optimum yang digunakan adalah \(0.6312679\). Nilai tersebut dapat diketahui dengan syntax berikut.

#Rho optimum
rho <- modelCO$rho
rho

## [1] 0.6312679

Selanjutnya akan dilakukan transformasi secara manual dengan syntax berikut ini.

#Transformasi Manual
y.trans <- y[-1]-y[-12]*rho
x1.trans<- x1[-1]-x1[-12]*rho
x2.trans<- x2[-1]-x2[-12]*rho
x3.trans<- x3[-1]-x3[-12]*rho
x4.trans<- x4[-1]-x4[-12]*rho
modelCOmanual<- lm(y.trans~x1.trans+x2.trans+x3.trans+x4.trans)
summary(modelCOmanual)

## 
## Call:
## lm(formula = y.trans ~ x1.trans + x2.trans + x3.trans + x4.trans)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -39.739 -11.327  -3.699   8.793  41.516 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  3.74164    9.81426   0.381  0.70428   
## x1.trans     0.47776    0.19661   2.430  0.01792 * 
## x2.trans     0.56453    0.26606   2.122  0.03773 * 
## x3.trans    -0.08618    0.21724  -0.397  0.69291   
## x4.trans    -0.05759    0.02029  -2.839  0.00606 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.82 on 64 degrees of freedom
## Multiple R-squared:  0.1917, Adjusted R-squared:  0.1412 
## F-statistic: 3.795 on 4 and 64 DF,  p-value: 0.007864

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[-2]
b0 <- b0bintang/(1-rho)
b1 <- modelCOmanual$coefficients[-1]
b0

## (Intercept)    x2.trans    x3.trans    x4.trans 
##  10.1473257   1.5310121  -0.2337204  -0.1561872

b1

##    x1.trans    x2.trans    x3.trans    x4.trans 
##  0.47775775  0.56453329 -0.08618020 -0.05759123

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]
  x3 <- model.matrix(model)[,4]
  x4 <- model.matrix(model)[,5]
  
  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]
  
  
  return(lm(y~x1+x2+x3+x4))
}

#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 90341.01
## 2 0.2 83717.37
## 3 0.3 78330.97
## 4 0.4 74150.66
## 5 0.5 71257.69
## 6 0.6 69853.34
## 7 0.7 70221.24
## 8 0.8 72654.65
## 9 0.9 77384.21

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.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))}))
head(tabOpt[order(tabOpt$SSE),])

##       rho      SSE
## 101 0.500 71257.69
## 100 0.499 71279.67
## 99  0.498 71301.80
## 98  0.497 71324.08
## 97  0.496 71346.51
## 96  0.495 71369.09

#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.500, y=71257.69    , labels = "rho=0.500", cex = 0.8)

Perhitungan yang dilakukan aplikasi R menunjukkan bahwa nilai \(ρ\) optimum, yaitu saat SSE terkecil terdapat pada nilai \(ρ=0.500\). 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.500, model)
summary(modelHL)

## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -53.517 -22.325  -3.658  15.815  99.520 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 17.67857   11.24554   1.572   0.1209  
## x1           0.47628    0.18522   2.571   0.0125 *
## x2           0.01142    0.22957   0.050   0.9605  
## x3           0.08331    0.22653   0.368   0.7143  
## x4          -0.03437    0.01989  -1.728   0.0888 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.37 on 64 degrees of freedom
## Multiple R-squared:  0.1451, Adjusted R-squared:  0.0917 
## F-statistic: 2.716 on 4 and 64 DF,  p-value: 0.03736

#Transformasi Balik
cat("y = ", coef(modelHL)[1]/(1-0.500), "+", coef(modelHL)[2],"x", sep = "")

## y = 35.35713+0.4762813x

Setelah dilakukan tranformasi balik, didapatkan model dengan metode Hildreth-Lu sebagai berikut. \[y_i=17.67857+0.47628x_1+0.01142x_2+0.08331x_3--0.03437x_4\]

#Deteksi autokorelasi
library(lmtest)
dwtest(modelHL)

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

Hasil uji Durbin-Watson juga menunjukkan bawah nilai DW sebesar \(1.4925\) berada pada selang daerah tidak ada autokorelasi, yaitu pada rentang DU < DW < 4-DU. Hal tersebut juga didukung oleh p-value sebesar \(0.01578\), di mana p-value > \(\alpha\)=5%. Artinya tak tolak \(H_0\) atau belum cukup bukti menyatakan bahwa ada autokorelasi dalam data 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,2,3,4)]
sseModelCO <- anova(modelCOmanual)$`Sum Sq`[-c(1,2,3,4)]
sseModelHL <- anova(modelHL)$`Sum Sq`[-c(1,2,3,4)]
mseModelawal <- sseModelawal/length(y)
mseModelCO <- sseModelCO/length(y)
mseModelHL <- sseModelHL/length(y)
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  98760.831            20334.5928         71257.686
## MSE   1410.869              290.4942          1017.967

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