Pemodelan SARIMAX (Seasonal Autoregressive Integrated and Moving Average with Exogenous Variable) dengan R

Pemodelan SARIMAX (Seasonal Autoregressive Integrated and Moving Average with Exogenous Variable) dengan R

Seasonal ARIMA dengan variabel eksogen (SARIMAX) di R

Hai teman-teman, jumpa lagi dengan blog jokoding.com, kali ini kita akan melanjutkan belajar bersama mengenai pemodelan statistik. Jenis pemodelan kita kali ini saya angkat karena ada sebuah permintaan (request) dari sahabat kita di LinkedIn beberapa waktu lalu ketika saya coba sharing terkait model Autoregressive Intergrated and Moving Average (ARIMA).

Model kita kali ini disebut sebagai SARIMAX. Apa kepanjangannya? Ya benar, Seasonal Autoregressive Integrated and Moving Average with Exogenous Regressor Model. Jadi, model ini merupakan model yang masih sekeluarga dari ARIMA. Model ini merupakan pengembangan ARIMA mengingat adanya dugaan faktor lain yang mempengaruhi sebuah variabel selain oleh masa lalu dirinya sendiri. Di dalam model SARIMAX, ini saya langsung pakai Seasonal karena data yang akan kita gunakan dalam praktik merupakan data yang memiliki pola musiman.

Setidak ada 2 istilah yang kita kenal di dalam pemodelan ini, yaitu variabel endogen dan variabel eksogen. Kendati ada istilah itu, namun esensinya sama dengan model regresi biasanya, variabel endogen mudahnya merupakan variabel dependen atau variabel yang dipengaruhi oleh variabel lain baik secara langsung maupun secara tidak langsung. Sedangkan variabel eksogen sederhananya merupakan variabel independen atau variabel yang memberikan pengaruh (memengaruhi) variabel lain secara langsung.

Kedua sebutan variabel ini, kelak akan kita bahas lebih lanjut dalam analasis Path (Path Analysis) atau bentuknya pengembangannya, yaitu Structural Equations Model (SEM).

Baik, kali ini yang akan kita praktikkan adalah pemodelan SARIMAX dengan variabel eksogen adalah jumlah kunjungan wisatawan ke Jawa Timur triwulanan sejak 2016 sampai 2021. Sedangkan variabel endogennya adalah laju pertumbuhan lapangan usaha penyediaan akomodasi dan makan minum triwulanan 2016 - 2021. Kenapa kita bahas kedua variabel tersebut? Karena selama pandemi Covid-19, sektor yang paling terdampak adalah sektor pariwisata dan untuk menguatkan pembentukan variabel eksogen dan endogennya, berikut beberapa hasil penelitian terkait dengan wisatawan dan dampaknya terhadap perekonomian suatu wilayah.

Penelitian terkait pengaruh jumlah wisatawan terhadap pertumbuhan sektor penyediaan akomodasi dan makan minum telah banyak dilakukan. Seperti hasil penelitian Yoga (2012) yang menggunakan analisis Path, ia membuktikan bahwa jumlah kunjungan wisatawan mancanegara berpengaruh positif dan signifikan terhadap Produk Domestik Regional Bruto (PDRB) Bali. Demikian halnya dengan penelitian yang dilakukan oleh Damayanti (2013), bahwa wisatawan mancanegara berpengaruh positif dan signifikan terhadap PDRB industri pariwisata. Dilanjutkan dengan penelitian Asworowati (2016) yang membuktikan bahwa pengaruh pengeluaran wisatawan mancanegara berpengaruh terhadap perekonomian daerah penelitiannya. Hasil serupa juga diperoleh dari penelitian Yoel (2008). Dengan menggunakan regresi linier dengan metode estimasi Ordinary Least Square (OLS), ia menyimpulkan bahwa jumlah wisatawan berpengaruh signifikan dan positif terhadap PDRB sektor pariwisata.

Dari beberapa hasil penelitian ini, data yang akan kita gunakan telah saya siapkan sebelumnya, jadi teman-teman yang berminat mempratikkan bisa mengunduhnya dulu pada tautan berikut. Setelah diunduh, pemodelan SARIMAX dengan menggunakan R dapat mengikuti beberapa code berikut:

Code:

#Install dan aktivasi package utama pemodelan
install.packages("forecast")
install.packages("fpp2")
library(fpp2)
library(forecast)

#import data
library(readxl)
dataku <- read_excel("C:/Users/Joko Ade/Downloads/arimaxku.xlsx")

#membentuk data Triwulanan
data_ts <- ts(dataku, start = c(2016,1), frequency = 4)

#Ringkasan data
summary(data_ts)

Hasil:

   lpeqtoq          wisatawan       
 Min.   :-18.410   Min.   :0.000000  
 1st Qu.: -1.298   1st Qu.:0.004026  
 Median :  1.670   Median :0.579726  
 Mean   :  1.099   Mean   :0.481412  
 3rd Qu.:  3.308   3rd Qu.:0.732320  
 Max.   :  9.710   Max.   :1.000000

Code:

#Plot data
autoplot(data_ts, facets = T) +
  xlab("Tahun") + ylab("") +
  ggtitle("Laju Pertumbuhan Lapus Penyediaan Akomodasi dan Makan Minum \ndan Jumlah Wisawatan Datang ke Jawa Timur Q1/2016 - Q4/2021")

Hasil:

Visualisasi 1

Code:

#Mendiagnosa pola ACF
ggAcf(data_ts) + ggtitle("Plot ACF")

#Mendiagnosa pola PACF
ggPacf(data_ts) + ggtitle("Plot PACF")

Hasil:

Plot ACF untuk menentukan order MA
 

Plot PACF untuk menentukan order AR

Code:

#Melihat pola dekomposisi masing-masing variabel
lpe_akmamin <- dataku$lpeqtoq %>% ts(., frequency = 4, start = c(2016,1))
lpe_akmamin_dekom <- stl(lpe_akmamin, s.window = "periodic")
autoplot(lpe_akmamin_dekom)

Hasil:

Dekomposisi runtun waktu untuk melihat pola data trend dan musiman

Code:

#Pemodelan SARIMAX dengan karena ada pola musiman pada dekomposisi
attach(dataku)
library(lmtest)

#Model SARIMAX1
sarimax1 <- auto.arima(lpeqtoq, xreg = wisatawan, trace = T, seasonal = TRUE, stepwise = F, approximation = F)

#Model SARIMAX2
sarimax2 <- arima(lpeqtoq, order = c(2, 0, 0), xreg = wisatawan, seasonal = list(order = c(0, 0, 1), period = 4), include.mean = T)

#Model SARIMAX3
#Model SARIMAX tanpa drift karena drift tidak signifikan
sarimax3 <- arima(lpeqtoq, order = c(2, 0, 0), xreg = wisatawan, seasonal = list(order = c(0, 0, 1), period = 4), include.mean = F)

Hasil:

ARIMA(0,0,0) with zero mean     : 152.4839
 ARIMA(0,0,0) with non-zero mean : 155.1011
 ARIMA(0,0,1) with zero mean     : 153.2841
 ARIMA(0,0,1) with non-zero mean : 156.1776
 ARIMA(0,0,2) with zero mean     : Inf
 ARIMA(0,0,2) with non-zero mean : Inf
 ARIMA(0,0,3) with zero mean     : Inf
 ARIMA(0,0,3) with non-zero mean : Inf
 ARIMA(0,0,4) with zero mean     : Inf
 ARIMA(0,0,4) with non-zero mean : Inf
 ARIMA(0,0,5) with zero mean     : Inf
 ARIMA(0,0,5) with non-zero mean : Inf
 ARIMA(1,0,0) with zero mean     : 154.7761
 ARIMA(1,0,0) with non-zero mean : 157.6731
 ARIMA(1,0,1) with zero mean     : Inf
 ARIMA(1,0,1) with non-zero mean : Inf
 ARIMA(1,0,2) with zero mean     : Inf
 ARIMA(1,0,2) with non-zero mean : Inf
 ARIMA(1,0,3) with zero mean     : Inf
 ARIMA(1,0,3) with non-zero mean : Inf
 ARIMA(1,0,4) with zero mean     : Inf
 ARIMA(1,0,4) with non-zero mean : Inf
 ARIMA(2,0,0) with zero mean     : 153.658
 ARIMA(2,0,0) with non-zero mean : 156.8756
 ARIMA(2,0,1) with zero mean     : Inf
 ARIMA(2,0,1) with non-zero mean : Inf
 ARIMA(2,0,2) with zero mean     : Inf
 ARIMA(2,0,2) with non-zero mean : Inf
 ARIMA(2,0,3) with zero mean     : Inf
 ARIMA(2,0,3) with non-zero mean : Inf
 ARIMA(3,0,0) with zero mean     : 156.2828
 ARIMA(3,0,0) with non-zero mean : 159.8826
 ARIMA(3,0,1) with zero mean     : Inf
 ARIMA(3,0,1) with non-zero mean : Inf
 ARIMA(3,0,2) with zero mean     : Inf
 ARIMA(3,0,2) with non-zero mean : Inf
 ARIMA(4,0,0) with zero mean     : 159.251
 ARIMA(4,0,0) with non-zero mean : 163.308
 ARIMA(4,0,1) with zero mean     : Inf
 ARIMA(4,0,1) with non-zero mean : Inf
 ARIMA(5,0,0) with zero mean     : 162.3275
 ARIMA(5,0,0) with non-zero mean : 166.8758

 Best model: Regression with ARIMA(0,0,0) errors

Code:

#Ringkasan Model
summary(sarimax1)
summary(sarimax2)
summary(sarimax3)

Hasil:

Series: lpeqtoq
Regression with ARIMA(0,0,0) errors

Coefficients:
        xreg
      2.4191
s.e.  1.8351

sigma^2 estimated as 29.01:  log likelihood=-73.96
AIC=151.91   AICc=152.48   BIC=154.27

Training set error measures:
                      ME     RMSE      MAE      MPE     MAPE      MASE       ACF1
Training set -0.06541324 5.272847 3.399757 79.78266 115.5155 0.6443709 -0.1138176

Call:
arima(x = lpeqtoq, order = c(2, 0, 0), seasonal = list(order = c(0, 0, 1), period = 4),
    xreg = wisatawan, include.mean = T)

Coefficients:
         ar1      ar2     sma1  intercept  wisatawan
      0.0603  -0.9668  -1.0000    -0.4202     2.5906
s.e.  0.0927   0.0776   0.2336     0.7965     1.4145

sigma^2 estimated as 14.55:  log likelihood = -69.38,  aic = 150.77

Training set error measures:
                   ME     RMSE      MAE     MPE    MAPE      MASE       ACF1
Training set 0.167009 3.815091 2.996547 39.8852 130.295 0.5679487 -0.2962499

Call:
arima(x = lpeqtoq, order = c(2, 0, 0), seasonal = list(order = c(0, 0, 1), period = 4),
    xreg = wisatawan, include.mean = F)

Coefficients:
         ar1      ar2     sma1  wisatawan
      0.0648  -0.9609  -0.9998     1.8707
s.e.  0.0950   0.0834   0.2384     0.3762

sigma^2 estimated as 14.71:  log likelihood = -69.52,  aic = 149.04

Training set error measures:
                    ME     RMSE      MAE      MPE     MAPE      MASE       ACF1
Training set 0.1524059 3.835937 2.949479 36.03903 128.7071 0.5590278 -0.2812725

Code:

#Uji Signifikansi
coeftest(sarimax2)
rbind(coeftest(sarimax1), coeftest(sarimax2), coeftest(sarimax3))

#Uji asumsi nonautokorelasi Ljung-Box
rbind(checkresiduals(sarimax1),checkresiduals(sarimax2),checkresiduals(sarimax3))

#Uji Normalitas Residual model SARIMAX
jarque.bera.test(sarimax2$residuals)
rbind(jarque.bera.test(sarimax1$residuals),jarque.bera.test(sarimax2$residuals),jarque.bera.test(sarimax3$residuals))

#AIC dan BIC
cbind(rbind(Model1 = AIC(sarimax1),Model2 = AIC(sarimax2), Model3 = AIC(sarimax3)),
      rbind(Model1 = BIC(sarimax1),Model2 = BIC(sarimax2), Model3 = BIC(sarimax3)))

Hasil:

             Estimate Std. Error     z value     Pr(>|z|)
xreg       2.41909032 1.83512782   1.3182135 1.874322e-01
ar1        0.06029365 0.09272973   0.6502084 5.155576e-01
ar2       -0.96678932 0.07759928 -12.4587406 1.253342e-35
sma1      -0.99999027 0.23360054  -4.2807703 1.862475e-05
intercept -0.42015063 0.79645614  -0.5275251 5.978290e-01
wisatawan  2.59062370 1.41445695   1.8315324 6.702112e-02
ar1        0.06477122 0.09498716   0.6818944 4.953057e-01
ar2       -0.96092361 0.08337186 -11.5257549 9.785080e-31
sma1      -0.99982246 0.23835525  -4.1946736 2.732649e-05
wisatawan  1.87072751 0.37623269   4.9722620 6.617618e-07


    Ljung-Box test

data:  Residuals from Regression with ARIMA(0,0,0) errors
Q* = 6.0169, df = 4, p-value = 0.1979

Model df: 1.   Total lags used: 5


    Ljung-Box test

data:  Residuals from ARIMA(2,0,0)(0,0,1)[4] with non-zero mean
Q* = 7.5277, df = 3, p-value = 0.05685

Model df: 5.   Total lags used: 8


    Ljung-Box test

data:  Residuals from ARIMA(2,0,0)(0,0,1)[4] with zero mean
Q* = 6.3942, df = 3, p-value = 0.09393

Model df: 4.   Total lags used: 7

     statistic parameter p.value    method          
[1,] 6.016899  4         0.1978898  "Ljung-Box test"
[2,] 7.52767   3         0.05685173 "Ljung-Box test"
[3,] 6.394197  3         0.09392982 "Ljung-Box test"
     data.name                                                 
[1,] "Residuals from Regression with ARIMA(0,0,0) errors"      
[2,] "Residuals from ARIMA(2,0,0)(0,0,1)[4] with non-zero mean"
[3,] "Residuals from ARIMA(2,0,0)(0,0,1)[4] with zero mean"

Uji normalitas

    statistic parameter p.value      method             data.name           
[1,] 22.8638   2         1.084396e-05 "Jarque Bera Test" "sarimax1$residuals"
[2,] 6.134182  2         0.0465564    "Jarque Bera Test" "sarimax2$residuals"
[3,] 5.410167  2         0.06686473   "Jarque Bera Test" "sarimax3$residuals"

AIC dan BIC

            AIC   BIC

           [,1]     [,2]
Model1 151.9124 154.2685
Model2 150.7675 157.8359
Model3 149.0419 154.9322

Code:

#Model Terpilih SARIMAX3
#Model SARIMAX tanpa drift karena drift tidak signifikan
sarimaxku3 <- arima(lpeqtoq, order = c(2, 0, 0), xreg = wisatawan, seasonal = list(order = c(0, 0, 1), period = 4), include.mean = F)

#Peramalan
datax <- dataku$wisatawan
datay <- dataku$lpeqtoq
xreg1 <- rep(mean(datax),4)
fc <- forecast(sarimaxku3, xreg = xreg1, h = 4)
poin_fc <- fc$mean
pred <- sarimaxku3$fitted
pred <- ts(pred, start = c(2016,1), frequency = 4)

pesimis <- as.data.frame(fc$lower)[2]
optimis <- as.data.frame(fc$upper)[2]
poin_fc <- ts(fc$mean, start = c(2022,1), frequency = 4)
line_ts <- ts(dataku$lpeqtoq, start = c(2016,1), frequency = 4)
a <- plot(line_ts, col = "blue", type = "l", xlim = c(2016.1, 2023.2), pch = 19,
          ylab = "Pertumbuhan (persen)", xlab = "Periode", ylim = c(-20,20))
b <- lines(pred, col = "red", type = "b")
c <- lines(poin_fc, col = "steelblue", type = "b")
d <- lines(ts(pesimis$`95%`, start = 2022.1, frequency = 4), col = "grey", type = "b")
e <- lines(ts(optimis$`95%`, start = 2022.1, frequency = 4), col = "grey", type ="b")
text(2022.1,15, c("Skenario\n Optimis"), cex = 0.8)
text(2022.1,-15, c("Skenario\n Pesimis"), cex = 0.8)
abline(v =2020.2, col = 5, lty=2)
text(2020.1, 15, c("Kebijakan PSBB"), cex = 0.8)
legend(2016.1,-5, legend = c("Data Aktual", "Prediksi dari SARIMAX(2,0,0)[4]", "Forecast", "Upper Forecast 95%",
                              "Lower Forecast 95%"), col = c("blue", "red", "steelblue", "grey", "grey"),
       lty = 1, cex = 0.8, box.col = "white", horiz = F)

Hasil:

Plot hasil peramalan dan memprediksi skenario optimis dan pesimis

Demikian sedikit ulasan mengenai pemodelan SARIMAX dengan R. Semoga sedikit banyak membantu pemahaman teman-teman, jangan lupa share dan terus ikuti unggahan selanjutnya. Selamat mempraktikkan!

Pemodelan Regresi Ridge (Ridge Regression Model) dengan R

Pemodelan Regresi Ridge (Ridge Regression Model) dengan R

Regresi Ridge dengan R

Halo teman-teman, kali ini kita akan melanjutkan belajar bersamanya mengenai pemodelan karena tadi pagi kita telah belajar bersama bagaimana menerapkan fungsi do.call() di R. Model yang kita bahas kali ini adalah model yang agaknya masih keluarga dekat dengan model regresi linier sederhana atau berganda, namanya adalah regresi ridge (ridge regression model).

Acapkali data yang kita gunakan untuk pemodelan regresi linier tidak berjalan mulus dan lancar-lancar saja. Entah dari hubungannya sebenarnya tidak linier, atau yang lebih sering mengalami gangguan asumsi klasik tertentu. Regresi ridge ini merupakan bentuk lain dari regresi yang mampu mengakomodir adanya bias akibat adanya multikolinearitas di antara variabel independen di dalam model. Adapun sifat dari penduga parameter model ini adalah bias namun konsisten karena memiliki kemampuan untuk menurunkan Mean Square Error (MSE). Kendati demikian, pada praktiknya, model regresi ridge ini masih debatable atau masih menjadi bahan perdebatan, terutama dari aspek efektivitas modelnya.

Selain itu, karakteristik dari model regresi ridge ini adalah uji signifikansi koefisien modelnya tidak dapat dilakukan secara langsung karena di dalam koefisiennya mengandung bias yang dalam formulasinya dinotasikan sebagai c. Demikian halnya dengan uji analisis varians dari model, kita tidak bisa melakukannya secara langsung. Perdebatan terhadap model ini juga terlihat akibat bias yang dimasukkan ke dalam model digunakan untuk meningkatkan prediksi.

Lantas bagaimana pengaplikasian regresi ridge dengan R? Kali ini kita akan mempraktikkan dengan menggunakan data dummy yang kita bangkitkan secara manual. Adapun langkah-langkah pemodelannya dapat mengikuti beberapa code berikut:

Code:

y <- c(1, 2, 3, 5, 1, 3, 7, 4, 9, 8)
x1 <- c(3, 2, 5, 7, 3, 5, 9, 6, 12, 10)
x2 <- c(1, 2, 3, 5, 1, 3, 6, 4, 8, 8)

df = data.frame(x1, x2, y)

#regresi OLS
ols <- lm(y~x1+x2, data = df)

#ringkasan hasil
summary(ols)

Hasil:

Call:
lm(formula = y ~ x1 + x2, data = df)

Residuals:
     Min       1Q   Median       3Q      Max
-0.31623 -0.16148 -0.09776  0.11088  0.49885

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -0.4905     0.2416  -2.030  0.08191 .
x1            0.2577     0.1290   1.998  0.08587 .
x2            0.7787     0.1633   4.768  0.00204 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3087 on 7 degrees of freedom
Multiple R-squared:  0.991,    Adjusted R-squared:  0.9884
F-statistic: 385.2 on 2 and 7 DF,  p-value: 6.925e-08

Code:

#Uji nonmultikol
library(olsrr)
ols_vif_tol(ols)


#uji homos
library(lmtest)
bptest(ols)

#uji nonautokorel
dwtest(ols)

#uji normalitas
ks.test(ols$residuals, ecdf(ols$residuals))

Hasil:

  Variables Tolerance      VIF
1        x1 0.0586826 17.04083
2        x2 0.0586826 17.04083

    studentized Breusch-Pagan test
data:  ols
BP = 2.5441, df = 2, p-value = 0.2803

    Durbin-Watson test
data:  ols
DW = 2.6447, p-value = 0.8619
alternative hypothesis: true autocorrelation is greater than 0

    One-sample Kolmogorov-Smirnov test

data:  ols$residuals
D = 0.1, p-value = 0.9996
alternative hypothesis: two-sided

Terlihat bahwa seluruh uji asumsi klasik memenuhi kecuali uji non-multikolinearitas, karena nilai dari VIF > 10, maka alternatifnya menggunakan model regresi ridge

Code:

#PEMODELAN RIDGE REGRESSION
library(glmnet)

#Definisikan variabel dependen
y <- df$y

#Definisikan variabel independen
x <- data.matrix(df[1:2])

#Model ridge
ridge <- glmnet(x, y, alpha = 0)

#ringkasan model ridge
summary(ridge)

Hasil:

          Length Class     Mode   
a0        100    -none-    numeric
beta      200    dgCMatrix S4     
df        100    -none-    numeric
dim         2    -none-    numeric
lambda    100    -none-    numeric
dev.ratio 100    -none-    numeric
nulldev     1    -none-    numeric
npasses     1    -none-    numeric
jerr        1    -none-    numeric
offset      1    -none-    logical
call        4    -none-    call   
nobs        1    -none-    numeric

Code:

#mencari lamda optimal dengan K-fold
cv_ridge <- cv.glmnet(x, y, nfolds = 3)

lamop <- cv_ridge$lambda.min

#Plot lambda
plot(cv_ridge)

Hasil:

Plot lambda untuk mendapatkan lambda optimal

Code:

#Model Final
ridgemod <- glmnet(x, y, alpha = 0, lambda = lamop)
coef(ridgemod)

Hasil:

3 x 1 sparse Matrix of class "dgCMatrix"
                    s0
(Intercept) -0.3570000
x1           0.3762883
x2           0.5668323

#Model regresi ridgenya: -0,3570000 + 0,3762883x1 + 0,5668323x2

Code:

#Menggunakan model fit untuk prediksi
y_predicted <- predict(ridgemod, s = lamop, newx = x)

#menghitung nilai Sum Square of Total dan Sum Square of Error
sst <- sum((y - mean(y))^2)
sse <- sum((y_predicted - y)^2)

#Menghitung nilai R Square model
rsq <- 1 - sse/sst
rsq
#atau bisa juga dilihat dari dev.ratio
ridgemod$dev.ratio

Hasil:

[1] 0.9863052

[1] 0.9863052

Code:

#Nilai MSE OLS
olspred <- predict(ols, newdata = df)
mean((olspred-df$y)^2)

#Nilai MSE Ridge
ridgepred <- predict(ridgemod, newx = x)
mean((ridgepred-df$y)^2)

Hasil:

OLS: [1] 0.06671445

Rigde: [1] 0.1014786

Terlihat bahwa MSE dari model ridge juga tidak lebih rendah daripada regresi OLS. Sedangkan kalau dilihat dari R square-nya relatif sama nilainya. Baik, demikian sekilas belajar bersama kita mengenai regresi ridge dengan R. Jangan lupa untuk terus menyimak unggahan berikutnya dan jangan lupa share. Selamat mempraktikkan!

Penerapan Fungsi do cal (do.call function) dengan R

Penerapan Fungsi do cal (do.call function) dengan R

Penerapan fungsi do call di R

Halo teman-teman, kali ini kita masih break ya dengan membahas beberapa fungsi ringan dan renyah dulu saja. Kali ini kita akan belajar bersama bagaimana fungsi dan manfaat do.call() di R.

Biasanya, dalam bahasa pemrograman, kita biasa melakukan instruksi untuk dieksekusi dengan menyiapkan fungsinya terlebih dahulu. Tapi di R, ada fungsi do.call() yang prinsip kerjanya terbalik, jadi kita menyiapkan bahan atau datanya terlebih dahulu, setelah siap kita langsung menggunakan fungsi ini untuk melakukan penugasan tertentu pada data yang telah kita siapkan, misalkan untuk mentransformasi data, menggabungkan data, atau perintah lainnya.

Untuk kali ini, kita mengenerate data secara manual saja untuk kemudahan praktikumnya. Kita dapat mengikuti  beberapa code berikut:

Code:

#Penerapan do.call
# do.call hanya bisa diterapkan untuk data list

#Membuat Data Frame
df1 <- data.frame(tim=c('A', 'B', 'C'),
                  point=c(22, 27, 38))

df2 <- data.frame(tim=c('D', 'E', 'F'),
                  point=c(22, 14, 20))

df3 <- data.frame(tim=c('G', 'H', 'I'),
                  point=c(11, 15, 18))

Hasil:

 tim point
1   A    22
2   B    27
3   C    38

  tim point
1   D    22
2   E    14
3   F    20

  tim point
1   G    11
2   H    15
3   I    18

Code:

#Transformasi dari data frame ke bentuk list
df_list <- list(df1, df2, df3)

str(df_list)

Hasil:

List of 3
 $ :'data.frame':    3 obs. of  2 variables:
  ..$ tim  : chr [1:3] "A" "B" "C"
  ..$ point: num [1:3] 22 27 38
 $ :'data.frame':    3 obs. of  2 variables:
  ..$ tim  : chr [1:3] "D" "E" "F"
  ..$ point: num [1:3] 22 14 20
 $ :'data.frame':    3 obs. of  2 variables:
  ..$ tim  : chr [1:3] "G" "H" "I"
  ..$ point: num [1:3] 11 15 18

Code:

#Menggabungkan baris df1, df2, dan df3
do.call(rbind, df_list)

Hasil:

   tim point
1   A    22
2   B    27
3   C    38
4   D    22
5   E    14
6   F    20
7   G    11
8   H    15
9   I    18

Baik, demikian sedikit ulasan mengenai fungsi dan manfaat do.call() di R. Terus simak unggahan berikutnya. Selamat mempraktikkan!

Visualisasi Matriks Korelasi dan Uji Signifikansi Korelasi antar Variabel dengan R

Visualisasi Matriks Korelasi dan Uji Signifikansi Korelasi antar Variabel dengan R

Visualisasi matriks korelasi dan uji signifikansi korelasi variabel

Hai teman-teman, kembali lagi kita akan belajar bersama-sama. Kalau kemarin kita sedikit serius membahas bagaimana pemodelan ARIMA dengan menggunakan R secara cukup urut dan rinci, kali ini kita break dulu untuk membahas visualisasi data.

Visualisasi yang akan kita ulas kali ini adalah hasil korelasi antar variabel yang kita gunakan dalam penelitian. Secara sederhana, korelasi sendiri merupakan besarnya hubungan keeratan antara suatu variabel dengan variabel lainnya yang ditunjukkan oleh besar dan arah. Mirip dengan pengertian besaran vektor yang juga ditunjukkan oleh besar dan arahnya, korelasi biasanya digunakan sebagai awalan deskripsi variabel penelitian untuk melihat seberapa besar keterkaitan antar variabel. Aspek yang perlu ditekankan sebagai pembeda antara korelasi dan regresi yaitu tujuannya. Kalau regresi menunjukkan hubungan sebab-akibat variabel yang digunakan, sedangkan kalau korelasi tidak menunjukkan sebab-akibat.

Seperti yang kita ketahui, korelasi yang populer hingga saat ini adalah korelasi Spearman atau disebut pula korelasi Rank Spearman dan korelasi Pearson atau disebut juga korelasi moment product. Perbedaan dari dua jenis korelasi ini umumnya pada jenis penelitiannya, kalau korelasi Pearson biasa digunakan untuk data-data parametrik atau data kontinu serta diksrit. Sedangkan korelasi Spearman biasanya digunakan untuk data-data nonparametrik baik kontinu maupun diskrit.

Adapun persamaan dari kedua jenis korelasi tersebut adalah pada besarnnya. Baik korelasi Spearman maupun korelasi Pearson, keduanya memiliki nilai korelasi pada selang interval 0 sampai 1. Semakin menuju 1, korelasi antara dua variabel dikatakan semakin kuat dan dalam praktiknya, para ahli statistika kemudian melakukan pengelompokan kekuatan korelasi menjadi beberapa tingkatan.

Untuk korelasi yang bernilai 0 dikatakan bahwa kedua variabel tidak berkorelasi. Bila nilai korelasi 0,01 - 0,20 dikatakan berkorelasi sangat rendah, korelasi 0,21 - 0,40 dikatakan berkorelasi rendah, 0,41 - 0,60 dikatakan agak rendah, 0,61 - 0,80 dikatakan korelasinya cukup, 0,81 - 0,99 dikatakan korelasinya tinggi, dan 1 berkorelasi sangat tinggi (sempurna). Selain pengelompokkan seperti ini, pengelompokan lainnya juga banyak dilakukan namun esensinya sama, semakin menuju 1 maka korelasi dua variabel dikatakan semakin kuat.

Perihal arah korelasi, sebenarnya ia hanya menunjukkan arah kecenderungan dalam interpretasinya. Bukan menunjukkan sebab dan akibat. Dalam hukum permintaan (product demand) misalkan, ketika harga barang semakin tinggi (mahal), permintaan konsumen akan cenderung menurun. Pada kasus ini dikatakan korelasinya memiliki arah negatif. Atau pada kasus lain misalkan, semakin tinggi lama pendidikan seseorang, maka ada kecenderungan konsumsi mereka semakin tinggi (banyak). Pada kasus ini dikatakan korelasinya memiliki arah yang positif.

Lantas? Bagaimana penerapan visualisasi sekaligus uji signifikansi korelasi banyak variabel dengan R? Untuk datanya kali ini kita akan men-generate data sendiri kemudian kita lakukan visualisasi matriks korelasi berikut dengan uji signifikansi korelasinya dengan menggunakan beberapa code berikut:

Code:

#Membuat Matriks Korelasi
#Membuat Data Frame
df <- data.frame(lamapendidikan=c(4, 5, 5, 6, 7, 8, 8, 10),
                 pengeluaran=c(12, 14, 14, 13, 15, 16, 15, 12),
                 pendapatan=c(22, 25, 26, 27, 29, 32, 39, 40))

#Menampilkan Data Frame
df

Hasil:

  lamapendidikan pengeluaran pendapatan
1              4          12         22
2              5          14         25
3              5          14         26
4              6          13         27
5              7          15         29
6              8          16         32
7              8          15         39
8             10          12         40

Code:

#Membuat Matriks Korelasi
cor(df)

Hasil:

              lamapendidikan pengeluaran pendapatan
lamapendidikan      1.0000000   0.1780213  0.9404385
pengeluaran         0.1780213   1.0000000  0.1646659
pendapatan          0.9404385   0.1646659  1.0000000

Terlihat bahwa korelasi antara pendapatan dan lama pendidikan arahnya positif dan cukup tinggi 0,9404385 sedangkan korelasi antara pengeluaran dan lama pendidikan sangat rendah dan positif 0,1780213, demikian halnya dengan korelasi pendapatan dan pengeluaran.

Code:

#Plot Matriks Korelasi
library(corrplot)
corrplot(cor(df))

Hasil:

visualisasi 1

Code:

#Mengcustome bentuk lain matriks korelasi
par(mfrow=c(2, 2))
corrplot(cor(df), method = "circle")
corrplot(cor(df), method = "pie")
corrplot(cor(df), method = "color")
corrplot(cor(df), method = "number")

Hasil:

Visualisasi 2

Code:

#Plot matriks korelasi dengan package ggcorrrplot
par(mfrow=c(1, 1))
library(ggcorrplot)
ggcorrplot(cor(df))

Hasil:

Visualisasi 3

Code:

#Menguji signifikansi korelasi
attach(df)
cor.test(lamapendidikan, pengeluaran)
cor.test(lamapendidikan, pendapatan)
cor.test(pengeluaran, pendapatan)

Hasil:

    Pearson's product-moment correlation

data:  lamapendidikan and pengeluaran
t = 0.44314, df = 6, p-value = 0.6732
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.6021950  0.7843056
sample estimates:
      cor
0.1780213

    Pearson's product-moment correlation

data:  lamapendidikan and pendapatan
t = 6.776, df = 6, p-value = 0.0005049
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.6989819 0.9894208
sample estimates:
      cor
0.9404385

    Pearson's product-moment correlation

data:  pengeluaran and pendapatan
t = 0.40893, df = 6, p-value = 0.6968
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.6108923  0.7789526
sample estimates:
      cor
0.1646659

Terlihat bahwa korelasi variabel yang signifikan hanya lama pendidikan dan pendapatan, selain tinggi dan positif korelasinya juga secara statistik signifikan.

Baik, demikian sedikit pembahasan kita terkait visualisasi korelasi antar variabel dan uji signifikansinya menggunakan R. Jangan lupa share dan bertanya bila ada di kolom komentar. Sampai jumpa pada unggahan selanjutnya. Selamat mempraktikkan!


Pemodelan Autoregressive Integrated Moving Average (ARIMA Model) dengan R

Pemodelan Autoregressive Integrated Moving Average (ARIMA Model) dengan R

ARIMA dengan R

Jumpa lagi teman-teman, sebelumnya saya mohon maaf karena kemarin tidak sempat membuat unggahan terbaru di blog ini. Baik, sebelumnya kita telah mengulas tentang pemodelan Geographically Weigthed Regression (GWR) dengan R. Kali ini, kita akan melanjutkan belajar bersama mengenai pemodelan yang tak asing lagi dan populer hingga kini, yaitu pemodelan Autoregressive Integrated Moving Average (ARIMA).

Kita akan membahas ARIMA secara langsung tanpa membahas AR dan MA secara tersendiri mengingat pada dasarnya ARIMA adalah model perpaduan antara model AR dengan order p, MA dengan order q dan aspek differencing dengan order d. Artinya, ketika kita mendengar istilah AR(1), maka sebenarnya itu adalah ARIMA(1, 0, 0), ketika kita mendengar ARI(1,1), maka aslinya itu ARIMA(1, 1, 0), atau bila mendengar MA(3), itu sebenarnya ARIMA(0, 0, 3) atau IMA(2,1) sebenarnya adalah ARIMA(0, 2, 1).

Data runtun waktu atau time series merupakan salah satu jenis data yang hingga kini banyak digunakan dalam penelitian. Popularitas jenis data ini disebabkan karena manfaat yang dapat kita peroleh selain mampu mengamati pola data serta historis data, kita dapat pula melakukan estimasi atau peramalan menggunakan data tersebut yang dipengaruhi oleh masa lalu data itu sendiri. Beberapa contoh kasus nyata bahwa suatu data dapat dipengaruhi oleh dirinya sendiri pada periode sebelumnya atau masa lalunya adalah inflasi, indeks harga saham, sentimen, atau yang akan kita gunakan datanya dalam praktik kali ini yaitu pertumbuhan ekonomi.

Bagaimana dengan batasan titik waktu yang harus digunakan dalam pemodelan ARIMA? Sebenarnya untuk jumlah titik waktu minimal yang dapat digunakan pemodelan ARIMA ini hingga kini masih dalam perdebatan. Menurut Pyndick dalam bukunya Makroekonomi minimal 30 titik waktu, namun kita ambil pilihan amannya, semakin banyak titik waktu yang digunakan akan semakin baik model ARIMA yang kita peroleh, berapa jumlah titik waktunya? Ukuran relevannya adalah dengan data runtun waktu yang kita gunakan, kita telah mampu menangkap pola dari data tersebut secara visual, misalkan pola musiman, pola harian, pola tahunan, pola trend, atau pola campuran. Ketika dari sekian periode waktu data itu, kita telah mampu mendiagnosa pola datanya, sebenarnya kita telah cukup untuk digunakan dalam pemodelan ARIMA. Namun jika dengan data yang tersedia, terdapat waktu-waktu di mana data kita menunjukkan pola yang tidak pasti atau dalam ketidakpastian, maka ada baiknya memang jumlah titik waktu data atau series datanya ditambah.

Kita kembali ke ARIMA, jadi di dalam pemodelan ARIMA ini, kita akan dihadapkan dengan istilah-istilah unik, misalnya stasioneritas data. Sebuah data dikatakan stasioner bila rata-rata dan varians data tersebut tidak beruba-ubah (berfluktuasi) seiring dengan perubahan waktu atau periodenya. Kondisi ini kemudian memunculkan istilah stasioner dalam rata-rata (mean) dan varians (ragam). Dalam penerapannya di R, uji stasioneritas biasanya menggunakan uji Augmented Dicky Fuller (ADF) yang merupakan uji pengembangan dari Dicky Fuller (DF). Fungsi uji ADF dalam R terdapat di dalam package tseries.

Bila sebuah data tidak stasioner pada level atau data aslinya (d = 0), maka kita akan menggunakan bentuk differencing data aslinya untuk kemudian diuji kembali stasioneritasnya. Bila didapatkan nilai p-value < alpha penelitian, maka dikatakan bahwa data telah stasioner. Bila data stasioner pada level maka nilai  = 0, jika stasioner pada differencing 1, maka = 1, dan seterusnya. Khusus untuk mengantisipasi data yang kita gunakan tidak stasioner dalam varians, kita dapat menggunakan data hasil transformasi Box-Cox sebagaimana yang telah kita pelajari pada unggahan transformasi variabel di blog ini. Kendati demikian, transformasi Box-Cox ini hanya mampu mengakomodir data yang bernilai positif saja, sedangkan untuk data pertumbuhan yang mengandung nilai negatif tidak bisa dan bahkan saya menemukan eror dalam eksekusi di R. Maka solusi alternatifnya bisa kita gunakan transfromasi invers Box-Cox.

Lantas bagaimana cara mendiagnosa order AR dan MA? Untuk mendiagnosa order AR kita mengacu pada sebaran data kelemabaman dari Partial Autocorrelation Function (PACF), sedangkan order MA dapat kita diagnosa berdasarkan sebaran data kelemabaman Autocorrelation Function (ACF). Berapa nilai order yang menjadi kandidat model AR dan MA? kita amati setiap garis yang melewati batas garis Bartlett sebaran data kelambaman ACF dan PACF tadi.

Setelah ordernya kita tentukan, tahapan selanjutnya adalah mengestimasi model ARIMA yang memungkinkan. Dalam R kita dapat menggunakan fungsi auto.arima() atau pemodelan secara manual dengan fungsi Arima() di dalam package forecast. Setelah seluruh kandidat model telah dimodelkan, berikutnya kita cek performa model dan akurasinya untuk mendapatkan model terbaik (best model). Adapun kriteria atau ukuran pemilihan model terbaik bisa berdasarkan nilai RMSE, MAE, MAPE, R Square, Akaike Information Criteria (AIC), atau Bayessian Information Criteria (BIC), atau ukuran lainnya. Namun, dalam praktik kali ini kita pilih model terbaiknya berdasarkan nilai RMSE terkecil.

Ketika model terbaik telah terpilih, maka selanjutnya kita dapat melakukan peramalan untuk beberapa periode selanjutnya dengan titik waktu sesuai kebutuhan analisis dengan menggunakan fungsi forecast() dalam package forecast.

Oke, demikian sekilas pengantar ARIMA, selanjutnya kita akan coba mempraktikkan bagaimana pemodelan ARIMA dengan R. Terlebih dulu kita akan siapkan data yang akan kita gunakan, yaitu data laju pertumbuhan ekonomi Jawa Timur Triwulanan, mulai Triwulan I-2011 hingga Triwulan IV-2021 pada tautan berikut. Setelah datanya telah kita unduh, pemodelan ARIMA dengan R dapat mengikuti beberapa code berikut:

Code:

#Import data laju pertumbuhan ekonomi Jatim Triwulan 1/2011- Triwulan IV/2021
library(readxl)
lpejatim <- read_excel("C:/Users/Joko Ade/Downloads/lpejatim2021.xlsx")

#Mengattach data
attach(lpejatim)

#Melihat sekilas data
str(lpejatim)
unlist(lpejatim)
as.numeric(lpejatim$Y)
head(lpejatim)

Hasil:

tibble [44 x 1] (S3: tbl_df/tbl/data.frame)
 $ Y: chr [1:44] "6.89" "6.40" "6.06" "6.43" ...

    Y1      Y2      Y3      Y4      Y5      Y6      Y7      Y8      Y9     Y10     Y11     Y12     Y13
 "6.89"  "6.40"  "6.06"  "6.43"  "6.49"  "6.62"  "6.86"  "6.59"  "6.17"  "5.84"  "5.72"  "6.59"  "5.85"
    Y14     Y15     Y16     Y17     Y18     Y19     Y20     Y21     Y22     Y23     Y24     Y25     Y26
 "5.93"  "6.14"  "5.52"  "5.20"  "5.41"  "5.43"  "5.71"  "5.58"  "5.81"  "5.64"  "5.27"  "5.37"  "5.05"
    Y27     Y28     Y29     Y30     Y31     Y32     Y33     Y34     Y35     Y36     Y37     Y38     Y39
 "5.64"  "5.76"  "5.41"  "5.56"  "5.39"  "5.53"  "5.60"  "5.80"  "5.33"  "5.40"  "2.90" "-5.86" "-3.47"
    Y40     Y41     Y42     Y43     Y44
"-2.65" "-0.44"  "7.07"  "3.27"  "4.59"

[1]  6.89  6.40  6.06  6.43  6.49  6.62  6.86  6.59  6.17  5.84  5.72  6.59  5.85  5.93  6.14  5.52  5.20
[18]  5.41  5.43  5.71  5.58  5.81  5.64  5.27  5.37  5.05  5.64  5.76  5.41  5.56  5.39  5.53  5.60  5.80
[35]  5.33  5.40  2.90 -5.86 -3.47 -2.65 -0.44  7.07  3.27  4.59

Code:

#konversi ke data runtun waktu
#frekuensi = 4 karena triwulanan
lpe_y <- ts(Y, start = c(2011,1), frequency = 4)

#melihat data runtun waktu
lpe_y
str(lpe_y)

Hasil:

     Qtr1  Qtr2  Qtr3  Qtr4
2011 6.89  6.40  6.06  6.43
2012 6.49  6.62  6.86  6.59
2013 6.17  5.84  5.72  6.59
2014 5.85  5.93  6.14  5.52
2015 5.20  5.41  5.43  5.71
2016 5.58  5.81  5.64  5.27
2017 5.37  5.05  5.64  5.76
2018 5.41  5.56  5.39  5.53
2019 5.60  5.80  5.33  5.40
2020 2.90  -5.86 -3.47 -2.65
2021 -0.44 7.07  3.27  4.59 

 Time-Series [1:44] from 2011 to 2022: 6.89 6.40 6.06 6.43 ...

Code:

#Visualisasi grafik garis
plot(lpe_y, col = "blue", lwd = 1, main = "Laju Pertumbuhan Ekonomi Jawa Timur \nTw I/2011 - Tw IV/2021 year on year", xlab = "Periode",
     ylab = "persen")
#Membubuhi grid pada grafik
grid()

#memberi garis vertikal warna putus-putus pada grafik
abline(v =2020.2, col = "red", lty=2)

#menambahkan teks keterangan pada garis vertikal merah
text(2019.2, 0, c("Efek Kebijakan PSBB"))

Hasil:

Visualisasi 1

Code:

#Uji Stasioneritas pada mean
library(tseries)
adf.test(lpe_y)


#uji stasioneritas pada varians (ragam)
require(car)
lamdafm1 <- boxCox(lpe_y~1, family = "yjPower")
lamdamod <- lamdafm1$x[which.max(lamdafm1$y)]
library(forecast)
lpe_box <- forecast(as.numeric(lpejatim$Y), lambda = lamdamod)
lpe_bc <- lpe_box$fitted

Hasil:

    Augmented Dickey-Fuller Test

data:  lpe_y
Dickey-Fuller = -4.3265, Lag order = 3, p-value = 0.01
alternative hypothesis: stationary

Nilai Lambda Box-Cox

Yang digunakan selanjutnya adalah lpe_bc

Code:

#Identifikasi Model Arima dan Plot ACF dan PACF
par(mfrow=c(2, 1))
acf(lpe_bc)
pacf(lpe_bc)

Hasil:

Diagnosa order ARIMA dari plot ACF dan PACF
Terlihat bahwa model ARIMA yang mungkin adalah ARIM(1, 1, 3), ARIMA(1, 1, 2), dan ARIMA(1, 1, 1) baik dengan drift (konstanta) maupun tanpa konstanta atau setidaknya ada 6 kandidat model ARIMA

Code:

#Memodelkan arima
arima1 <- Arima(lpe_bc, order = c(1, 1, 3), include.drift = F)
arima2 <- Arima(lpe_bc, order = c(1, 1, 2), include.drift = F)
arima3 <- Arima(lpe_bc, order = c(1, 1, 1), include.drift = F)
arima4 <- Arima(lpe_bc, order = c(1, 1, 3), include.drift = T)
arima5 <- Arima(lpe_bc, order = c(1, 1, 2), include.drift = T)
arima6 <- Arima(lpe_bc, order = c(1, 1, 1), include.drift = T)

#Performa model
library(performance)
cbind(Model = c("arima1_113WD", "arima2_112WD", "arima3_111WD", "arima4_113D", "arima5_112D", "arima6_111D"),
      rbind(performance(arima1), performance(arima2), performance(arima3), performance(arima4),
            performance(arima5), performance(arima6)))

#Akurasi Model
cbind(Model = c("arima1_113WD", "arima2_112WD", "arima3_111WD", "arima4_113D", "arima5_112D", "arima6_111D"),
      rbind(accuracy(arima1), accuracy(arima2), accuracy(arima3), accuracy(arima4),
            accuracy(arima5), accuracy(arima6)))

Hasil:

         Model      AIC      BIC     RMSE    Sigma
1 arima1_113WD 172.4274 181.2334 1.551617 1.648081
2 arima2_112WD 171.0174 178.0622 1.561728 1.637954
3 arima3_111WD 169.0502 174.3338 1.549629 1.605322
4  arima4_113D 171.1483 181.7155 1.456645 1.567431
5  arima5_112D 170.0821 178.8881 1.468487 1.559783
6  arima6_111D 170.9937 178.0385 1.548613 1.624199

             Model          ME                    RMSE               MAE                
Training set "arima1_113WD" "-0.259609486763692"  "1.55161669987552" "0.671955748856711"
Training set "arima2_112WD" "-0.241673526778238"  "1.56172802583685" "0.65957405897493"
Training set "arima3_111WD" "-0.0552082070266032" "1.54962865789642" "0.692723965049203"
Training set "arima4_113D"  "0.0803514319937035"  "1.45664538456452" "0.727631232089374"
Training set "arima5_112D"  "0.0992586023896908"  "1.46848703826404" "0.756376385151208"
Training set "arima6_111D"  "0.00023447190214626" "1.54861284512487" "0.688015537688339"
             MPE                  MAPE               MASE               ACF1                 
Training set "2.78490489815687"   "17.1241960247521" "1.06356838866507" "-0.0895988988457491"
Training set "3.61167160071564"   "16.8526596491993" "1.04397070834323" "-0.0275484078350913"
Training set "-0.591718127659727" "17.4012301341839" "1.09644022325965" "-0.0799684281753089"
Training set "9.4241238117586"    "17.3589874339763" "1.15169128082078" "-0.0614280088800428"
Training set "11.0168221505381"   "17.6050937493751" "1.19718897345296" "-0.0430505221739149"
Training set "0.113333421180444"  "17.3254527504112" "1.0889877466496"  "-0.0800315290870742"

Code:

#Uji Signifikasi Model Arima
library(lmtest)
rbind(cbind(Model = c("arima1_113WD"),rbind(coeftest(arima1))),
cbind(Model = c("arima2_112WD"),rbind(coeftest(arima2))),
cbind(Model = c("arima3_111WD"),rbind(coeftest(arima3))),
cbind(Model = c("arima4_113D"),rbind(coeftest(arima4))),
cbind(Model = c("arima5_112D"),rbind(coeftest(arima5))),
cbind(Model = c("arima6_111D"),rbind(coeftest(arima6))))

Hasil:

      Model          Estimate             Std. Error           z value              Pr(>|z|)              
ar1   "arima1_113WD" "0.636286741337564"  "0.449053441049423"  "1.41695104228705"   "0.156497241726144"   
ma1   "arima1_113WD" "-0.419573044318619" "0.460601485844076"  "-0.910924209351452" "0.362335308848335"   
ma2   "arima1_113WD" "-0.800891249391554" "0.193247621699613"  "-4.14437829737679"  "3.40736855884907e-05"
ma3   "arima1_113WD" "0.352319589531121"  "0.352257530521445"  "1.00017617511139"   "0.317225256934432"   
ar1   "arima2_112WD" "0.0378948644684753" "0.295683694715215"  "0.12816014256374"   "0.898022241521018"   
ma1   "arima2_112WD" "0.143113146939234"  "0.281474408184658"  "0.508441061701587"  "0.61114406248608"    
ma2   "arima2_112WD" "-0.706199595456142" "0.169172931135828"  "-4.17442430485015"  "2.98740524659347e-05"
ar1   "arima3_111WD" "-0.634075276768642" "0.127015437988826"  "-4.99211187874991"  "5.97226418857177e-07"
ma1   "arima3_111WD" "0.999999877748336"  "0.067114488567726"  "14.8999105720551"   "3.30020911121018e-50"
ar1   "arima4_113D"  "0.580617375518388"  "0.326779107733077"  "1.77678854546801"   "0.0756030372083661"  
ma1   "arima4_113D"  "-0.493567262310261" "0.349083581049896"  "-1.4138942336555"   "0.157392959282349"   
ma2   "arima4_113D"  "-0.862153783283784" "0.180703835364512"  "-4.77108735154772"  "1.83234049334131e-06"
ma3   "arima4_113D"  "0.355725447777676"  "0.333809752697306"  "1.06565324980257"   "0.286580418226939"   
drift "arima4_113D"  "-0.121719569017145" "0.0444373506687291" "-2.73912749489811"  "0.00616024792593575"
ar1   "arima5_112D"  "0.171766595217114"  "0.186618023289883"  "0.920418039956951"  "0.357354345010246"   
ma1   "arima5_112D"  "-0.113718270169818" "0.167655681566651"  "-0.678284619448518" "0.497591248730932"   
ma2   "arima5_112D"  "-0.886272428131124" "0.16250746519766"   "-5.45373363034822"  "4.93230691444996e-08"
drift "arima5_112D"  "-0.124376284025448" "0.0384252558337676" "-3.23683684927214"  "0.00120862489519324"
ar1   "arima6_111D"  "-0.634080545480926" "0.126927940941883"  "-4.99559467187176"  "5.86547403550869e-07"
ma1   "arima6_111D"  "0.999999571312214"  "0.0673260320919646" "14.8530893658824"   "6.64351382618105e-50"
drift "arima6_111D"  "-0.069295087857588" "0.291670145614196"  "-0.237580324553502" "0.812206612354163"

Code:

#Model terbaik berdasarkan RMSE terkecil
arima <- arima5

#Melakukan Peramalan h = berapa titik yang akan diramal
fc <- forecast(arima, h = 4)

#Ploting secara keseluruhan termasuk hasil forecast
par(mfrow=c(1,1))
plot(forecast(arima, h = 4))

Hasil:

Plot hasil peramalan dengan model ARIMA terbaik

Code:

#Mendiagnosa hasil pemodelan
tsdiag(arima)

#Mengenerate Resid2 setiap model
resid2 <- arima$residuals

#Melakukan uji nilai tengah residuals dengan uji t
t.test(resid2, mu = 0, alternative = "two.sided")

Kalau p-value uji t > alpha penelitian maka aman

#Melihat Akurasi Peramalan
accuracy(arima)

Hasil:

Hasil diagnosa kenormalan residual ARIMA terpilih

    One Sample t-test

data:  resid2
t = 0.44425, df = 43, p-value = 0.6591
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.3513310  0.5498482
sample estimates:
mean of x
0.0992586 

Terlihat residual terlihat normal

Code:

#Melihat R Square Model Fit
library(caret)
pred <- arima$fitted
R2(pred, as.numeric(lpejatim$Y))

Hasil:

[1] 0.2345608

nilai R square  ARIMA terpiliha kita sebesar 23,45 persen

Code:

#Menampilkan Prediksi dan nilai asli
library(ggfortify)
lpepred <- as.numeric(pred)
lpeaktual <- as.numeric(lpejatim$Y)
predtrueval <- data.frame(lpepred, lpeaktual)
predtrueval.ts <- ts(predtrueval)
theme_set(theme_bw())
autoplot(predtrueval.ts,xlab = "Periode", ylab = "Laju Pertumbuhan Jawa Timur yoy") +
  ggtitle("Plot Prediksi dan Data Aktual Laju Pertumbuhan Ekonomi Jawa Timur") +
  theme(plot.title = element_text(hjust = .5))

Hasil:

Plot hasil peramalan dan prediksi dengan model ARIMA terpilih

Demikian sekilas pemodelan ARIMA dengan R. Jangan lupa share dan bertanya bila ada di kolom komentar. Selamat mempraktikkan!

Pemodelan Geographically Weigthed Regression (GWR Model) menggunakan R

Pemodelan Geographically Weigthed Regression (GWR Model) menggunakan R

Geographically Weigthed Regression (GWR) dengan R

Halo teman-teman, kemarin kita telah sejenak break dari pemodelan dengan membahas fungsi sprintf() yang ada di R. Kali ini kita akan melanjutkan kembali pembahasan model-model statistik dengan R. Yang akan kita bahas kali ini adalah bagaimana menerapkan model Geographically Weigthed Regression atau biasa disingkat model GWR.

GWR merupakan sebuah model alternatif bila di dalam pemodelan regresi linier (dalam parameter) kita terganggu oleh asumsi heteroskedastisitas spasial atau asumsi homoskedastisitas residual model tidak terpenuhi. Mengapa disebut heteroskedastisitas spasial? Ya karena di dalam data kita terdiri atas amatan-amatan lokasi yang masing-masing memiliki variabel pembobot spasial dan di dalam praktiknya masing-masing lokasi ini ditunjukkan oleh posisi astronomis longitute dan latitude.

Kalau di model regresi Robust yang sebelumnya kita bahas, ia merupakan model alternatif bila residual modelnya terganggu asumsi normalitas dan atau heteroskedastisitas, tapi setiap amatannya bukan terdiri atas lokasi-lokasi amatan. Bila ternyata di dalam data kita memiliki lokasi-lokasi amatan yang ditunjukkan oleh longitude dan latitude, maka model GWR bisa juga menjadi alternatif model untuk analisis.

Baik, kita kembali ke model GWR. Bila di dalam model regresi dengan metode estimasi parameter Ordinary Least Square (OLS) hanya melihat pengaruh variabel independen terhadap variabel dependennya baik secara parsial maupun simultan, maka di model GWR kita dapat melihat seberapa besar pengaruh variabel independen terhadap variabel dependen di sebuah wilayah atau lokasi amatan dan seberapa besar pengaruh variabel independen seluruh lokasi terhadap variabel dependennya. Kalau diistilahkan, di dalam model GWR ini kita akan mengenal model lokal dan model global. Artinya, pada masing-masing lokasi atau wilayah akan memiliki model regresi GWR yang berbeda dan secara global juga memiliki sebuah model GWR.

Secara mendasar, karakteristik dari model GWR ini ditunjukkan oleh matriks pembobot dengan menggunakan beberapa teknik dengan menggunakan fungsi Kernel, yaitu fixed kernel Gaussian, adaptive kernel Gaussian, fixed kernel Bisquare, adaptive kernel Bisquare atau lainnya. Menurut Fotheringham et.al (2002), fungsi matriks pembobot ini memiliki bandwith (h) yang berbeda untuk setiap lokasi. Bandwith ini sebenarnya merupakan parameter penghalus non-negatif berbentuk lingkaran yang memiliki radius sebesar b dari titik pusat lokasi sebagai dasar penentuan bobot model regresi di setiap lokasi amatan. Dengan demikian, untuk pengamatan yang dekat dengan lokasi atau wilayah X misalkan, maka ia akan lebih berpengaruh terhadap pembentukan parameter lokasi atau wilayah X tersebut.

Dalam kasus model GWR, nilai bandwith yang kecil menyebabkan varians membesar sehingga diperlukan metode-metode bagaimana memilih bandwith sedemikian rupa sehingga mampu memperkecil varians. Sebab varians yang besar (tidak homogen) menyebabkan penduga parameter yang dihasilkan tidak akurat. Adapun metode pemilihan bandwith yang optimum ini ada beberapa, yaitu dengan berdasarkan nilai Cross Validation (CV), Akaike Information Criteria (AIC), General Cross Validation (GCV), atau Bayessian Information Criteria (BIC).

Secara matematis, bentuk model GWR secara umum dapat dituliskan sebagai berikut:

Bentuk umum model Geographically Weigthed Regression (GWR)

Baik, itu sekilas teori mengenai GWR teman-teman. Selanjutnya kita akan mempraktikkan bagaimana memodelkan GWR dengan R. Data yang akan kita gunakan dalam praktik kali ini adalah data mengenai jumlah kematian bayi, indeks ketersediaan pangan, dan jumlah persalinan yang dibantu tenaga kesehatan. Data tersebut bersumber dari Profil Kesehatan Jawa Timur dan kemeterian Pertanian tahun 2020 yang dapat teman-teman unduh pada tautan berikut. Setelah datanya telah siap, pemodelan GWR dapat mengikuti beberapa code berikut:

Code:

#Mengimport Data
library(readxl)
dataku <- read_excel("C:/Users/Joko Ade/Downloads/datagwr.xlsx")
#Sesuaikan dengan letak penyimpanan file

#melihat sekilas struktur data
str(dataku)

Hasil:

tibble [38 x 6] (S3: tbl_df/tbl/data.frame)
 $ Wilayah  : chr [1:38] "Pacitan" "Ponorogo" "Trenggalek" "Tulungagung" ...
 $ long     : num [1:38] 111 111 112 112 112 ...
 $ lat      : num [1:38] -8.13 -7.87 -8.18 -8.91 -8.1 ...
 $ JKB      : num [1:38] 62 139 38 158 115 163 86 174 338 117 ...
 $ IP       : num [1:38] 74.1 76.2 71.6 76.5 76.1 ...
 $ Salinakes: num [1:38] 6.55 10.62 8.99 14.3 15.52 ...

Code:

#Melakukan Seleksi data dan melihat ringkasan data
dataku <- dataku[,c(-1)]
summary(dataku)

#Mengattach data
attach(dataku)

Hasil:

      long            lat              JKB              IP          Salinakes     
 Min.   :111.1   Min.   :-8.912   Min.   :  6.0   Min.   :58.40   Min.   : 2.103  
 1st Qu.:112.0   1st Qu.:-7.971   1st Qu.: 44.5   1st Qu.:67.11   1st Qu.: 8.729  
 Median :112.5   Median :-7.712   Median : 96.0   Median :75.11   Median :14.629  
 Mean   :112.6   Mean   :-7.691   Mean   :101.8   Mean   :73.58   Mean   :15.062  
 3rd Qu.:113.1   3rd Qu.:-7.463   3rd Qu.:148.8   3rd Qu.:79.41   3rd Qu.:17.820  
 Max.   :114.4   Max.   :-6.895   Max.   :338.0   Max.   :90.11   Max.   :43.534

Code:

#Pemodelan regresi OLS
ols <- lm(JKB~IP+Salinakes)

#Ringkasan model OLS
summary(ols)

#Melihat AIC
AIC(ols)

Hasil:

Call:
lm(formula = JKB ~ IP + Salinakes)

Residuals:
     Min       1Q   Median       3Q      Max
-123.989  -32.292    1.592   29.144  105.588

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 210.3187    76.3647   2.754  0.00927 **
IP           -2.5151     1.0633  -2.365  0.02368 *  
Salinakes     5.0789     0.8212   6.184 4.42e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 49.12 on 35 degrees of freedom
Multiple R-squared:  0.5261,    Adjusted R-squared:  0.499
F-statistic: 19.42 on 2 and 35 DF,  p-value: 2.114e-06 

AIC: [1] 408.6804

Code:

#Uji Asumsi Normalitas Residual
ks.test(ols$residuals, ecdf(ols$residuals))
shapiro.test(residuals(ols))

#Uji Non Autokorelasi Residual
dwtest(ols)

#Uji Non Multikolinieritas Variabel Bebas Model
ols_vif_tol(ols)

#Uji Homoskedastisitas Residual Model
bptest(ols)

Hasil:

    One-sample Kolmogorov-Smirnov test

data:  ols$residuals
D = 0.026316, p-value = 1
alternative hypothesis: two-sided

    Shapiro-Wilk normality test

data:  residuals(ols)
W = 0.99049, p-value = 0.9836

    Durbin-Watson test

data:  ols
DW = 1.6189, p-value = 0.09199
alternative hypothesis: true autocorrelation is greater than 0

  Variables Tolerance      VIF
1        IP 0.9317424 1.073258
2 Salinakes 0.9317424 1.073258

    studentized Breusch-Pagan test

data:  ols
BP = 11.072, df = 2, p-value = 0.003943

Interpretasi: terlihat bahwa model regresi memenuhi semua uji asumsi klasik kecuali uji homoskedastisitas karena nilai p-value Breush-Pagan 0,0003943 atau < 0,05 (alpha penelitian), maka model GWR dapat menjadi alternatif bila setiap amatan kita adalah lokasi dengan memiliki data longitude dan latitude

Code:

#Karena Terganggu Heteroskedastisitas Opsi Lain Bisa Menggunakan GWR
#Install dan aktivasi beberapa package untuk pemodelan GWR
install.packages("spgwr")
install.packages("sp")
install.packages("spData")
library(spgwr)
library(sp)
library(spData)

#Pemodelan GWR
#Mencari Bandwith, misalkan Bisquare Adaptive
band_bidap <- gwr.sel(formula = JKB~IP+Salinakes, data = dataku,
                      coords = cbind(dataku$long, dataku$lat), adapt = TRUE, gweight = gwr.bisquare)

Hasil:

Adaptive q: 0.381966 CV score: 71986.42
Adaptive q: 0.618034 CV score: 85144.32
Adaptive q: 0.236068 CV score: 79166.29
Adaptive q: 0.3985704 CV score: 72813.99
Adaptive q: 0.3493846 CV score: 71763.77
Adaptive q: 0.3061015 CV score: 78009.31
Adaptive q: 0.3617677 CV score: 71361.09
Adaptive q: 0.3639216 CV score: 71345.02
Adaptive q: 0.3650074 CV score: 71341.32
Adaptive q: 0.3658301 CV score: 71340.32
Adaptive q: 0.3659504 CV score: 71340.29
Adaptive q: 0.3659911 CV score: 71340.29
Adaptive q: 0.3660318 CV score: 71340.3
Adaptive q: 0.3659911 CV score: 71340.29 ----> nilai bandwith yang digunakan 0,3659911

Code:

#pemodelan GWR
modelfix <- gwr(formula = JKB~IP+Salinakes, data = dataku,
                coords = cbind(dataku$long, dataku$lat), bandwidth = band_bidap,
                hatmatrix = TRUE)

#Pengecekan koefisien model
modelfix$lm$coefficients

#Mengecek signifikansi variabel dalam model
LMZ.F3GWR.test(modelfix)

Hasil:

(Intercept)          IP   Salinakes
 210.318714   -2.515098    5.078869

Leung et al. (2000) F(3) test

            F statistic Numerator d.f. Denominator d.f.     Pr(>)    
(Intercept)      5.9564        13.6357           24.734 6.602e-05 ***
IP               5.7204        12.0966           24.734 0.0001238 ***
Salinakes        3.4095        13.6284           24.734 0.0039280 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Interpretasi: terlihat seluruh variabel independen berpengaruh signifikan statistik terhadap variabel dependen

Code:

#Melihat AIC
modelfix$results[8]

Hasil:

$AICh
[1] 355.3489

Code:

#cek Anova modelnya perbandingan OLS dan GWR
#bila signifikan berarti model GWR lebih baik dari OLS
anova(modelfix)

#mencari nilai F tabel
qf(p = 0.05, df1 = anova(modelfix)[2,1], df2 = anova(modelfix)[3,1],
   lower.tail = FALSE)
#Jika nilai |F hitung| > F tabel tolak H0 artinya GWR lebih baik dari OLS
# Atau Mencari p-value nilai F Anova, jika < 0.05 tolak H0 artinya GWR lebih baik dari OLS
pf(q = anova(modelfix)[3,4], df1 = anova(modelfix)[2,1],
   df2 = anova(modelfix)[3,1], lower.tail = FALSE)

Hasil:

Analysis of Variance Table
                    Df Sum Sq Mean Sq F value
OLS Residuals    3.000  84451                
GWR Improvement 16.169  67176  4154.7        
GWR Residuals   18.831  17275   917.4  4.5289

Nilai F-Tabel: [1] 2.217443

Nilai p-value: [1] 0.00118195

Code:

#Selain itu uji model Global dapat dilakukan dengan uji F LMZ
LMZ.F1GWR.test(modelfix)

Hasil:

    Leung et al. (2000) F(1) test

data:  modelfix
F = 0.38019, df1 = 24.734, df2 = 35.000, p-value = 0.007164
alternative hypothesis: less
sample estimates:
SS OLS residuals SS GWR residuals
        84451.27         17274.95

Code:

#Melihat Intersep model lokal
modelfix$SDF$`(Intercept)`

#melihat Koefisien IP lokal
modelfix$SDF$IP

Hasil:

[1]  -165.493877  -109.437424  -158.089894 -1046.972846   -67.279399    20.705466   212.146159   483.523836
 [9]   555.629387   836.886416   552.258880   659.191891   487.971285   158.777521   137.298122   120.523552
[17]   111.351302     6.488024   -57.322304   -62.434970   -23.977988    68.010960   395.432418   244.031644
[25]   201.522895   150.840227   124.501058   259.844691   353.003696   -26.771901   -67.279399    86.833596
[33]   343.234161   156.035701   146.394134   -66.493632   162.328836    58.899471


[1]   2.00291333   1.42004407   2.23550425  16.07065781   1.33894888  -0.04054907  -2.36473580  -6.71981226
 [9]  -7.97684097 -11.37571776  -7.54229252  -8.63917298  -7.19002370  -1.92665900  -1.74410450  -1.46782969
[17]  -1.33347504   0.06484369   0.80541727   0.76526999   0.29533918  -0.72940634  -4.83500655  -3.08973403
[25]  -2.62090249  -2.14313147  -2.06405256  -4.44891112  -5.27246568   0.56831402   1.33894888  -0.77049139
[33]  -4.85706069  -1.97571435  -1.83141695   0.88599258  -2.16050614  -0.45529833


Oke, demikian sedikit pembahasan mengenai penerapan model Geographically Weigthed Regression (GWR) dengan R. Jangan lupa share dan bisa komentar bila ada pertanyaan di kolom komentar. Dan pastinya, simak dan ikuti terus unggahan-unggahan berikutnya dalam blog ini. Selamat mempraktikkan!