Data Kematian akibat Gagal Ginjal

Import Data

library(readxl)
dataginjal =  read_excel("DataGinjalUAS.xlsx")
dataginjal

summary(dataginjal$`Deaths - Chronic kidney disease - Sex: Both - Age: All Ages (Number)`)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   22381   27084   32357   32149   37296   42131

tsginjal = ts(dataginjal$`Deaths - Chronic kidney disease - Sex: Both - Age: All Ages (Number)`, start=c(1990,1))
tsginjal

## Time Series:
## Start = 1990 
## End = 2019 
## Frequency = 1 
##  [1] 22381 22809 23412 24052 24677 25314 26018 26855 27770 28605 29306 29721
## [13] 30449 31235 31985 32728 33682 34407 35225 35997 36670 37234 37317 37624
## [25] 37822 38126 39406 40286 41234 42131

library(fpp)
library(DT)
datatable(dataginjal)

library(ggplot2)
plotly::ggplotly(ggplot(data = dataginjal, aes(x = Year, y = `Deaths - Chronic kidney disease - Sex: Both - Age: All Ages (Number)`))+geom_line()+ggtitle("Plot Chronic Kidney Disease 1990 - 2019"))

Exponential Smoothing Methods

Untuk lebih meyakinkan metode apa yang akan kita ambil, kita dapat membentuk visualisasi dengan menggunakan decompose untuk memastikan data ini sebenarnya masuk ke pola data yang mana.

library(tidyverse)
glimpse(dataginjal)

## Rows: 30
## Columns: 5
## $ Entity                                                                 <chr> …
## $ Code                                                                   <chr> …
## $ Year                                                                   <dbl> …
## $ Time                                                                   <dbl> …
## $ `Deaths - Chronic kidney disease - Sex: Both - Age: All Ages (Number)` <dbl> …

period = ts(dataginjal$`Deaths - Chronic kidney disease - Sex: Both - Age: All Ages (Number)`, start=c(1990,1),frequency = 10)
period

## Time Series:
## Start = c(1990, 1) 
## End = c(1992, 10) 
## Frequency = 10 
##  [1] 22381 22809 23412 24052 24677 25314 26018 26855 27770 28605 29306 29721
## [13] 30449 31235 31985 32728 33682 34407 35225 35997 36670 37234 37317 37624
## [25] 37822 38126 39406 40286 41234 42131

periode = ts(dataginjal$`Deaths - Chronic kidney disease - Sex: Both - Age: All Ages (Number)`, start = range(dataginjal$Time)[1],frequency = 10)
periode

## Time Series:
## Start = c(1, 1) 
## End = c(3, 10) 
## Frequency = 10 
##  [1] 22381 22809 23412 24052 24677 25314 26018 26855 27770 28605 29306 29721
## [13] 30449 31235 31985 32728 33682 34407 35225 35997 36670 37234 37317 37624
## [25] 37822 38126 39406 40286 41234 42131

decompose(periode) %>%
  autoplot()

Dari visualisasi yang terbentuk, dapat dilihat bahwa pada data memiliki trend yang naik.

fcsingle = ses(tsginjal,h=3)
summary(fcsingle)

## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
##  ses(y = tsginjal, h = 3) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
## 
##   Initial states:
##     l = 22382.5457 
## 
##   sigma:  737.5851
## 
##      AIC     AICc      BIC 
## 502.1690 503.0921 506.3726 
## 
## Error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE     ACF1
## Training set 658.3447 712.5749 658.4477 2.082823 2.083283 0.9668346 0.439506
## 
## Forecasts:
##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2020       42130.91 41185.66 43076.16 40685.27 43576.55
## 2021       42130.91 40794.19 43467.63 40086.57 44175.25
## 2022       42130.91 40493.79 43768.03 39627.15 44634.67

plot(fcsingle, col="red", ylab = "Chronic kidney disease" , xlab = "Year")

data = window(tsginjal, start=1990)
fc = holt(data, h=3)
autoplot(fc)

summary(fc)

## 
## Forecast method: Holt's method
## 
## Model Information:
## Holt's method 
## 
## Call:
##  holt(y = data, h = 3) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.5923 
## 
##   Initial states:
##     l = 21447.3492 
##     b = 606.8282 
## 
##   sigma:  271.0375
## 
##      AIC     AICc      BIC 
## 443.8783 446.3783 450.8843 
## 
## Error measures:
##                    ME     RMSE      MAE        MPE      MAPE     MASE
## Training set 16.83478 252.3222 149.9195 0.05030664 0.4751712 0.220135
##                    ACF1
## Training set 0.02502777
## 
## Forecasts:
##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2020       43036.95 42689.60 43384.30 42505.72 43568.17
## 2021       43942.89 43289.82 44595.96 42944.11 44941.67
## 2022       44848.84 43847.73 45849.94 43317.78 46379.89

tise = ts(dataginjal$`Deaths - Chronic kidney disease - Sex: Both - Age: All Ages (Number)`, start = range(dataginjal$Time)[1],frequency = 10)

hw<-HoltWinters(tise)

predic<-forecast(hw,h=3)
summary(predic)

## 
## Forecast method: HoltWinters
## 
## Model Information:
## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = tise)
## 
## Smoothing parameters:
##  alpha: 1
##  beta : 0
##  gamma: 1
## 
## Coefficients:
##           [,1]
## a   41841.4400
## b     732.5194
## s1    254.4600
## s2    -84.4400
## s3   -117.2400
## s4    -81.5900
## s5    -73.9400
## s6   -176.5900
## s7   -164.4400
## s8    -24.8900
## s9    179.1100
## s10   289.5600
## 
## Error measures:
##                    ME   RMSE     MAE       MPE     MAPE       MASE     ACF1
## Training set 122.4574 842.83 353.148 0.4611765 1.080599 0.05194995 0.047564
## 
## Forecasts:
##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 4.00       42828.42 41731.99 43924.85 41151.57 44505.27
## 4.10       43222.04 41671.45 44772.63 40850.62 45593.46
## 4.20       43921.76 42022.68 45820.83 41017.38 46826.14

SSE<-hw$SSE
SSE

## [1] 14207247

plot(hw)

Metode Box Jenkins

Membuat Plot Time Series

plot.ts(tsginjal, ylab = "Chronic kidney disease", xlab = "Year", lwd =2 , col = 'red')

Berdasarkan time series plot diatas diketahui bahwa jumlah kematian akibat Chronic kidney disease d mengalami peningkatan dari tahun ke tahun atau memiliki pola trend naik.

Melakukan Uji Stasioneritas

Uji stasioner dapat dilihat dengan melalui plot Autocorrelation Function dan uji Augmented Dickey-Fuller. Data yang belum stasioner merupakan data yang terlihat meluruh secara lambat menuju 0

Membuat Plot ACF dan PACF

par(mfrow = c(1,2))
acf(tsginjal, main = "Plot ACF")
pacf(tsginjal, main = "Plot PACF")

Jika kita melihat plot ACF maka dapat disimpulkan bahwa grafik terlihat meluruh secara lambah menuju 0 artinya bahwa data tersebut belum stasioner.

Melakukan Uji ADF

adf.test(tsginjal)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  tsginjal
## Dickey-Fuller = -2.7924, Lag order = 3, p-value = 0.2683
## alternative hypothesis: stationary

Pada hasil uji ADF, dengan menggunakan taraf signifikansi 5% data yang tersedia tidak mendukung hipotesis nol yang berarti data tersebut tidak stasioner. Maka dari itu diperlukan proses differencing.

dataginjal_diff1 = diff(tsginjal)
ts.plot(dataginjal_diff1, main = "Differensi Orde 1")

adf.test(dataginjal_diff1)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dataginjal_diff1
## Dickey-Fuller = -3.8578, Lag order = 3, p-value = 0.03039
## alternative hypothesis: stationary

dataginjal_diff2 = diff(tsginjal, differences = 2)
ts.plot(dataginjal_diff2, main = "Differensi Orde 2")

adf.test(dataginjal_diff2)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dataginjal_diff2
## Dickey-Fuller = -3.839, Lag order = 3, p-value = 0.03215
## alternative hypothesis: stationary

Estimasi Model

Selanjutnya dalam mengestimasi model dapat dilihat dari plot ACF dan PACF sebagai berikut.

par(mfrow = c(1,2))
Acf(dataginjal_diff1)
pacf(dataginjal_diff1)

Berdasarkan pada plot ACF dan PACF, dapat dilihat bahwa terdapat lag yang keluar dari limit maka akan didaptkan AR (1) dan MA(1). Sehingga didapatkan model ARIMA (1,1,1) yang berarti bahwa model mengandung auto regressive AR(1) dan Moving Average MA(1). Maka unTk mendapatkan model yang terbaik, akan dilakukan overfitting model disekitar model utama. yaitu

# Model 1 MA (1)
model1=Arima(tsginjal,order = c(0,1,1)) 
model1

## Series: tsginjal 
## ARIMA(0,1,1) 
## 
## Coefficients:
##          ma1
##       0.8490
## s.e.  0.0865
## 
## sigma^2 = 204438:  log likelihood = -218.58
## AIC=441.17   AICc=441.63   BIC=443.9

# Model 2 MA (2)
model2=Arima(tsginjal,order = c(0,1,2)) 
model2

## Series: tsginjal 
## ARIMA(0,1,2) 
## 
## Coefficients:
##          ma1     ma2
##       0.9620  0.9999
## s.e.  0.4707  0.9684
## 
## sigma^2 = 139154:  log likelihood = -214.71
## AIC=435.42   AICc=436.38   BIC=439.52

# Model 3 AR (1)
model3=Arima(tsginjal,order = c(1,1,0)) 
model3

## Series: tsginjal 
## ARIMA(1,1,0) 
## 
## Coefficients:
##          ar1
##       0.9368
## s.e.  0.0536
## 
## sigma^2 = 65280:  log likelihood = -202.44
## AIC=408.88   AICc=409.34   BIC=411.62

# Model 4 ARMA (1,1)
model4=Arima(tsginjal,order = c(1,1,1)) 
model4

## Series: tsginjal 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.9683  -0.2618
## s.e.  0.0383   0.2063
## 
## sigma^2 = 63706:  log likelihood = -201.64
## AIC=409.28   AICc=410.24   BIC=413.38

# Model 5 ARMA (2,2)
model5=Arima(tsginjal,order = c(2,1,2)) 
model5

## Series: tsginjal 
## ARIMA(2,1,2) 
## 
## Coefficients:
##          ar1     ar2     ma1      ma2
##       0.2850  0.6675  0.4553  -0.2821
## s.e.  0.3591  0.3377  0.3996   0.2140
## 
## sigma^2 = 67541:  log likelihood = -201.4
## AIC=412.81   AICc=415.42   BIC=419.65

fitmodel= auto.arima(tsginjal)
summary(fitmodel)

## Series: tsginjal 
## ARIMA(1,1,0) with drift 
## 
## Coefficients:
##          ar1     drift
##       0.4697  679.9663
## s.e.  0.1644   74.2785
## 
## sigma^2 = 51294:  log likelihood = -197.49
## AIC=400.98   AICc=401.94   BIC=405.08
## 
## Training set error measures:
##                   ME     RMSE      MAE        MPE      MAPE      MASE
## Training set 5.65271 214.8585 140.7229 0.01195042 0.4183572 0.2066311
##                      ACF1
## Training set -0.002900604

checkresiduals(model3)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,0)
## Q* = 6.5286, df = 5, p-value = 0.2581
## 
## Model df: 1.   Total lags used: 6

tsdiag(model3)

Berdasarkan plot gambar diatas dapat dilihat pada plot ACF residual data merupakan model white noise ditandai dengan tidak adanya lag yang keluar dari garis batas interval. Sedangkan dilihat dari plot p-value for Ljung-Box statistic, nilai p-value ada yang lebih dari α=0,05 yang artinya tidak terdapat autokorelasi pada data sehingga asumsi no-autokorelasi terpenuhi. Kemudian dilihat dari normaliatas data sebagai berikut.

Selanjutnya adalah Peramalan

pred= predict(model3,n.ahead = 3)
pred

## $pred
## Time Series:
## Start = 2020 
## End = 2022 
## Frequency = 1 
## [1] 42971.35 43758.63 44496.18
## 
## $se
## Time Series:
## Start = 2020 
## End = 2022 
## Frequency = 1 
## [1] 255.4996 556.9285 909.5533

summary(model3)

## Series: tsginjal 
## ARIMA(1,1,0) 
## 
## Coefficients:
##          ar1
##       0.9368
## s.e.  0.0536
## 
## sigma^2 = 65280:  log likelihood = -202.44
## AIC=408.88   AICc=409.34   BIC=411.62
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE      MAPE      MASE      ACF1
## Training set 61.05644 246.8361 160.0461 0.2046482 0.4844479 0.2350044 -0.253209