DỰ BÁO TỈ LỆ SINH

Time Series Analysis

Chúng ta tiếp tục với một ví dụ về phân tích dữ liệu dạng Time Series để tìm hiểu những thống kê và kết quả cần đạt được trong phân tích một bộ dữ liệu Time Series, trong trường hợp này là số ca sinh, thay đổi có tính chu kì theo tháng (seasional) kết hợp với trend và random (errors).

1 Chuẩn bị dữ liệu

Đây là bộ dữ liệu về số ca sinh theo tháng ở thành phố New York từ tháng 01/1946 đến tháng 12/1959 được lấy từ file http://robjhyndman.com/tsdldata/data/nybirths.dat.

births <- scan("http://robjhyndman.com/tsdldata/data/nybirths.dat")
head(births, 10)

##  [1] 26.663 23.598 26.931 24.740 25.806 24.364 24.477 23.901 23.175 23.227

birthstimeseries <- ts(births, frequency=12, start=c(1946,1))
plot.ts(birthstimeseries)

Có 168 giá trị về số ca sinh theo tháng trong khoảng thời gian nói trên. Biểu đồ với trục x là thời gian, trục y là số ca sinh, đơn vị là ngàn (1.000). Đường biễu diễn cho thấy biến thiên có chu kì trong một năm, giảm dần đến năm 1948-1949 rồi tăng trở lại ổn định đến năm 1959.

Mục tiêu là phân tích bản chất biến thiên của dữ liệu và dự đoán phát triển số ca sinh trong thời gian tới, 1 năm hoặc 2 năm.

Trước hết, chúng ta xem biến thiên của các số sinh trung bình qua các năm và trung bình số ca sinh theo tháng trong Biểu đồ sau:

library(dplyr)
library(data.table)
# Chuyển một ts objects thành data.frame
dmn <- list(month.abb, unique(floor(time(birthstimeseries))))
birthseries <- as.data.frame(t(matrix(birthstimeseries, 12, dimnames = dmn)))
setDT(birthseries, keep.rownames = TRUE)[]

##       rn    Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep
##  1: 1946 26.663 23.598 26.931 24.740 25.806 24.364 24.477 23.901 23.175
##  2: 1947 21.439 21.089 23.709 21.669 21.752 20.761 23.479 23.824 23.105
##  3: 1948 21.937 20.035 23.590 21.672 22.222 22.123 23.950 23.504 22.238
##  4: 1949 21.548 20.000 22.424 20.615 21.761 22.874 24.104 23.748 23.262
##  5: 1950 22.604 20.894 24.677 23.673 25.320 23.583 24.671 24.454 24.122
##  6: 1951 23.287 23.049 25.076 24.037 24.430 24.667 26.451 25.618 25.014
##  7: 1952 23.798 22.270 24.775 22.646 23.988 24.737 26.276 25.816 25.210
##  8: 1953 24.364 22.644 25.565 24.062 25.431 24.635 27.009 26.606 26.268
##  9: 1954 24.657 23.304 26.982 26.199 27.210 26.122 26.706 26.878 26.152
## 10: 1955 24.990 24.239 26.721 23.475 24.767 26.219 28.361 28.599 27.914
## 11: 1956 26.217 24.218 27.914 26.975 28.527 27.139 28.982 28.169 28.056
## 12: 1957 26.589 24.848 27.543 26.896 28.878 27.390 28.065 28.141 29.048
## 13: 1958 27.132 24.924 28.963 26.589 27.931 28.009 29.229 28.759 28.405
## 14: 1959 26.076 25.286 27.660 25.951 26.398 25.565 28.865 30.000 29.261
##        Oct    Nov    Dec
##  1: 23.227 21.672 21.870
##  2: 23.110 21.759 22.073
##  3: 23.142 21.059 21.573
##  4: 22.907 21.519 22.025
##  5: 24.252 22.084 22.991
##  6: 25.110 22.964 23.981
##  7: 25.199 23.162 24.707
##  8: 26.462 25.246 25.180
##  9: 26.379 24.712 25.688
## 10: 27.784 25.693 26.881
## 11: 29.136 26.291 26.987
## 12: 28.484 26.634 27.735
## 13: 27.945 25.912 26.619
## 14: 29.012 26.992 27.897

library(reshape2)
df <- melt(birthseries)
colnames(df) <- c("Year", "Month", "Births")
t = df %>% group_by(Year) %>%
  summarise_at(vars(Births), funs(mean(., na.rm = TRUE)))

p <-  ggplot(t, aes(x = Year, y = Births, group = 1))
p + geom_line(colour = "blue", size = 1) +
  scale_y_continuous("Birhts in Thousands") +
  ggtitle("Birth Means by Year")

Biểu đồ này xác nhận lại một khuynh hướng là số ca sinh trung bình giảm đến năm 1949, sau đó tăng dần đến đỉnh cao vào năm 1958.

Còn số ca sinh trung bình tính theo các tháng trong năm thì sao.

p1 <- ggplot(df, aes(x = Month, y = Births, fill = Month ))
p1 + geom_boxplot() +
  ggtitle("Births by Months") + 
  scale_y_continuous("Birhts in Thousands") +
  guides(fill=FALSE)

Tính theo tháng thì Tháng Hai có số ca sinh trung bình thấp nhất, sau đó dao động và tháng Bảy có số ca sinh cao nhất, rồi giảm dần đến tháng 11.

Từ đó chúng ta thấy loại dữ liệu này có những đặc tính thành phần cần phải phân tích, đó là:

1. Khuynh hướng biến thiên chung (trend)
2. Chu kì biến thiên (seasional)
3. Biến thiên ngẫu nhiên (random hoặc errors)

Phân tích riêng 3 thành phần trên để có được một sự đánh giá sâu sắc hơn về biến thiên của những trường hợp sinh trong khoảng thời gian này.

library(TTR)
dc <- decompose(birthstimeseries)
plot(dc)

Hiệu chỉnh biến thiên có chu kì (theo tháng), có nghĩa sau khi loại bỏ tính chu kì của số ca sinh, chúng ta sẽ có được khuynh hướng biến thiên và thay đổi ngẫu nhiên qua các năm như sau:

dc$seasonal

##             Jan        Feb        Mar        Apr        May        Jun
## 1946 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1947 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1948 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1949 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1950 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1951 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1952 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1953 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1954 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1955 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1956 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1957 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1958 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1959 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
##             Jul        Aug        Sep        Oct        Nov        Dec
## 1946  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1947  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1948  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1949  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1950  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1951  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1952  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1953  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1954  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1955  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1956  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1957  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1958  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1959  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197

dcadjusted <- birthstimeseries - dc$seasonal
plot(dcadjusted)

Nhìn biểu đồ trên ta thấy biến thiên của số ca sinh nhiều khả năng giảm và tăng dần không hằng định, tức là non-stationary. Để xác thẩm định khả năng này, ta chuyển dữ liệu về mô hình difference = 1, tức là:

yt = xt - xt-1

birthseriesdiff1 <- diff(birthstimeseries, differences=1)
plot.ts(birthseriesdiff1)

Biểu đồ của mô hình với difference = 1, cho thấy khuynh hướng biến thiên hằng định hơn, đường biểu diễn giá trị trung bình của yt tương đối nằm ngang.

2 Dự báo bằng mô hình HoltWinters (exponential smoothing)

Dự báo biến thiên số ca sinh bằng hàm HoltWinters của package forecast. Phương pháp này tính sum of squared errors trong mẫu rồi dự đoán giai đoạn tương lai.

Biểu đồ của hàm này được vẽ lên với mức độ dao động giảm hơn (smoothing).

birthseriesforecasts <- HoltWinters(birthstimeseries, beta=FALSE, gamma=FALSE)
plot(birthseriesforecasts)

Holt-Winters exponential smoothing function dự báo trong mẫu (in-sample errors) với độ chính xác khá cao. Từ đó phân tích variance và phân phối của những errors và residuals của mô hình cho ta thấy được rõ hơn quy luật biến thiên của dữ liệu.

Các thông số thêm vào như alpha, beta và gamma chính là intercept, độ dốc của biểu đồ và của chu kì (seasional). Nếu set TRUE thì biểu đồ thường overfitting, vì vậy set FALSE sẽ cho biểu đồ ít phức tạp hơn, giảm variance.

library(forecast)
birthseriesforecasts <-HoltWinters(birthstimeseries)
predict(birthseriesforecasts,n.ahead=24)

##           Jan      Feb      Mar      Apr      May      Jun      Jul
## 1960 27.30020 25.92822 29.04779 27.56001 28.88793 28.50749 30.56943
## 1961 27.80419 26.43221 29.55178 28.06400 29.39192 29.01148 31.07342
##           Aug      Sep      Oct      Nov      Dec
## 1960 30.55133 29.94243 29.63266 27.53065 28.36129
## 1961 31.05533 30.44642 30.13665 28.03464 28.86528

plot(birthstimeseries,xlim=c(1946,1962))
lines(predict(birthseriesforecasts,n.ahead=24),col=2)

Với h = 24, chúng ta yêu cầu hàm fredict() dự báo cho chúng ta trong 24 thời điểm sắp tới (24 tháng). Giá trị dự báo của hàm này biến thiên tương đồng với các giá trị quan sát được trong vài năm trước đó, tuy nhiên không có các khoảng tin cậy 80% hay 95%.

Chúng ta sử dụng residuals của hàm này để khảo sát tính tương quan với nhau (autocorrelation) của dữ liệu ở thời điểm hiện tại so với các thời điểm lùi về trước (lag).

Hai hàm ACF và PACF sẽ được sử dụng để đánh giá tương quan (autocorrelation) của các giá trị xt giữa các thời điểm.

acf(birthstimeseries,lag.max=20)

pacf(birthstimeseries, lag.max=20)

Hai biểu đồ trên cho thấy mô hình AR(1) phù hợp với dữ liệu này.

Chúng ta khảo sát thêm ACF và PACF của residuals của mô hình này.

birthseriesforecasts2 <- forecast:::forecast.HoltWinters(birthseriesforecasts)
reds <- birthseriesforecasts2$residuals[-c(1:12)]
acf(reds,lag.max=20)

pacf(reds, lag.max=20)

Từ hai biểu đồ ACF, PACF này cho thấy hầu như không có tương quan có ý nghĩa sau lag 2 cho đến lag 13, vì thế mô hình Moving Average (MA) cũng phù hợp. Từ đó một mô hình phối hợp AR và MA có thể phù hợp với bộ dữ liệu này.

Chúng ta tiếp tục phân tích mean, variance và distribution của các errors để biết về qui luật biến thiên của dữ liệu.

# Function for forecasterror
plotForecastErrors <- function(forecasterrors)
{
  # make a histogram of the forecast errors:
  mybinsize <- IQR(forecasterrors)/4
  mysd <- sd(forecasterrors)
  mymin <- min(forecasterrors) - mysd*5
  mymax <- max(forecasterrors) + mysd*3
  # generate normally distributed data with mean 0 and standard deviation mysd
  mynorm <- rnorm(10000, mean=0, sd=mysd)
  mymin2 <- min(mynorm)
  mymax2 <- max(mynorm)
  if (mymin2 < mymin) { mymin <- mymin2 }
  if (mymax2 > mymax) { mymax <- mymax2 }
  # make a red histogram of the forecast errors, with the normally distributed
  #data overlaid:
  mybins <- seq(mymin, mymax, mybinsize)
  hist(forecasterrors, col="red", freq=FALSE, breaks=mybins)
  # freq=FALSE ensures the area under the histogram = 1
  # generate normally distributed data with mean 0 and standard deviation mysd
  myhist <- hist(mynorm, plot=FALSE, breaks=mybins)
  # plot the normal curve as a blue line on top of the histogram of forecast
  #errors:
  points(myhist$mids, myhist$density, type="l", col="blue", lwd=2)
}
plotForecastErrors(reds)

Biểu đồ plotForecastErrors cho thấy các errors có phân phối chuẩn với mean = 0.

3 Dựng mô hình ARIMA

Để có được mô hình dự báo số ca sinh tương lai khác hơn hàm Holt-Winters exponential smoothing function, chúng ta xây dựng mô hình ARIMA.

Trước hết, chúng ta sử dụng hàm auto.arima() để tìm ra các thông số p, d, q và seasional của mô hình này.

auto.arima(birthstimeseries)

## Series: birthstimeseries 
## ARIMA(2,1,2)(1,1,1)[12]                    
## 
## Coefficients:
##          ar1      ar2      ma1     ma2     sar1     sma1
##       0.6539  -0.4540  -0.7255  0.2532  -0.2427  -0.8451
## s.e.  0.3004   0.2429   0.3228  0.2879   0.0985   0.0995
## 
## sigma^2 estimated as 0.4076:  log likelihood=-157.45
## AIC=328.91   AICc=329.67   BIC=350.21

Các thông số của hàm ARIMA là ARIMA(2,1,2)(1,1,1)[12], với p = 2, d = 1, q = 2. Seasional là c(1,1,1) và chu kì là 12. So sánh với kết quả từ difference = 1 và ACF cho biết p = 2, khá thống nhất.

Code cho mô hình này sẽ là:

fit <- arima(birthstimeseries, c(2, 1, 2),seasonal = list(order = c(1, 1, 1), period = 12))
birthseriesforecasts <- forecast(fit, h=24)
birthseriesforecasts

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 1960       27.69056 26.88679 28.49433 26.46130 28.91982
## Feb 1960       26.07680 24.98034 27.17326 24.39991 27.75369
## Mar 1960       29.26544 28.04020 30.49069 27.39160 31.13929
## Apr 1960       27.59444 26.29165 28.89723 25.60199 29.58689
## May 1960       28.93193 27.54860 30.31527 26.81630 31.04757
## Jun 1960       28.55379 27.07347 30.03411 26.28983 30.81774
## Jul 1960       29.84713 28.26538 31.42888 27.42806 32.26620
## Aug 1960       29.45347 27.77916 31.12778 26.89284 32.01410
## Sep 1960       29.16388 27.40777 30.91999 26.47814 31.84962
## Oct 1960       29.21343 27.38167 31.04519 26.41200 32.01486
## Nov 1960       27.26221 25.35695 29.16746 24.34837 30.17604
## Dec 1960       28.06863 26.09098 30.04628 25.04408 31.09318
## Jan 1961       27.66908 25.63754 29.70062 24.56210 30.77606
## Feb 1961       26.21255 24.12791 28.29719 23.02437 29.40073
## Mar 1961       29.22612 27.08633 31.36591 25.95360 32.49865
## Apr 1961       27.58011 25.38474 29.77547 24.22259 30.93762
## May 1961       28.71354 26.46431 30.96277 25.27363 32.15344
## Jun 1961       28.21736 25.91651 30.51820 24.69852 31.73620
## Jul 1961       29.98728 27.63644 32.33812 26.39198 33.58258
## Aug 1961       29.96127 27.56138 32.36117 26.29095 33.63160
## Sep 1961       29.56515 27.11690 32.01339 25.82088 33.30941
## Oct 1961       29.54543 27.04965 32.04121 25.72846 33.36240
## Nov 1961       27.57845 25.03603 30.12088 23.69015 31.46676
## Dec 1961       28.40796 25.81976 30.99616 24.44965 32.36627

plot(birthseriesforecasts)

h = 24 có nghĩa là ta muốn mô hình ARIMA dự báo cho chúng ta số ca sinh trong 24 tháng tới (24 time points).

Biểu đồ được vẽ trên đây kèm theo giai đoạn 24 tháng từ tháng 01/1960 đến tháng 12/1961, đi kèm khoảng tin cậy 80% và 95%.

Để vẽ chi tiết giai đoạn 1960-1961 trên, chúng ta xử lí số liệu như sau:

library(data.table)
births6061 <- as.data.frame(birthseriesforecasts)
setDT(births6061, keep.rownames = TRUE)[]

##           rn Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
##  1: Jan 1960       27.69056 26.88679 28.49433 26.46130 28.91982
##  2: Feb 1960       26.07680 24.98034 27.17326 24.39991 27.75369
##  3: Mar 1960       29.26544 28.04020 30.49069 27.39160 31.13929
##  4: Apr 1960       27.59444 26.29165 28.89723 25.60199 29.58689
##  5: May 1960       28.93193 27.54860 30.31527 26.81630 31.04757
##  6: Jun 1960       28.55379 27.07347 30.03411 26.28983 30.81774
##  7: Jul 1960       29.84713 28.26538 31.42888 27.42806 32.26620
##  8: Aug 1960       29.45347 27.77916 31.12778 26.89284 32.01410
##  9: Sep 1960       29.16388 27.40777 30.91999 26.47814 31.84962
## 10: Oct 1960       29.21343 27.38167 31.04519 26.41200 32.01486
## 11: Nov 1960       27.26221 25.35695 29.16746 24.34837 30.17604
## 12: Dec 1960       28.06863 26.09098 30.04628 25.04408 31.09318
## 13: Jan 1961       27.66908 25.63754 29.70062 24.56210 30.77606
## 14: Feb 1961       26.21255 24.12791 28.29719 23.02437 29.40073
## 15: Mar 1961       29.22612 27.08633 31.36591 25.95360 32.49865
## 16: Apr 1961       27.58011 25.38474 29.77547 24.22259 30.93762
## 17: May 1961       28.71354 26.46431 30.96277 25.27363 32.15344
## 18: Jun 1961       28.21736 25.91651 30.51820 24.69852 31.73620
## 19: Jul 1961       29.98728 27.63644 32.33812 26.39198 33.58258
## 20: Aug 1961       29.96127 27.56138 32.36117 26.29095 33.63160
## 21: Sep 1961       29.56515 27.11690 32.01339 25.82088 33.30941
## 22: Oct 1961       29.54543 27.04965 32.04121 25.72846 33.36240
## 23: Nov 1961       27.57845 25.03603 30.12088 23.69015 31.46676
## 24: Dec 1961       28.40796 25.81976 30.99616 24.44965 32.36627
##           rn Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95

colnames(births6061) <- c("month.year", "births", "Lo.80", "Hi.80", "Lo.95", "Hi.95")
births6061$id <- NA
births6061$id <- c(1:24)

Và vẽ lại biểu đồ dự báo số ca sinh trong giai đoạn này kèm theo khoảng tin cậy 80% và 95%.

library(ggplot2)
g <- ggplot(births6061, aes(x = id, group = 1))
g + geom_line(aes(y=Lo.80), colour="orange", linetype = "dotted", size=0.1) +
  geom_line(aes(y=Hi.80), colour="orange", linetype = "dotted", size=0.1) +
  geom_line(aes(y=Lo.95), colour="green", linetype = "dotted", size=0.1) +
  geom_line(aes(y=Hi.95), colour="green", linetype = "dotted", size=0.1) +
  ggtitle("Forecasts for Births in 1960 - 1961")+
  scale_y_continuous("Births in Thousands") +
  scale_x_continuous("24 Months from Jan. 1960 to Dec. 1961") +
  theme(axis.line = element_line(colour="black"))+
  geom_ribbon(aes(x=id, ymax=Hi.95, ymin=Lo.95), fill="green", alpha=.1)+
  geom_ribbon(aes(ymax=Hi.80, ymin=Lo.80), fill="pink", alpha=.7) +
  geom_point(aes(y=births)) +
  geom_line(aes(y=births), colour = "blue", size=1) 

Thời điểm tương lai càng xa thì khoảng tin cậy 80% và 95% của giá trị dự báo càng lớn, giảm mức độ chính xác.

Đây là bài viết về thực hành phân tích và dự báo một dạng dữ liệu Time Series thường gặp trong nghiên cứu y học với dạng dữ liệu biến thiên có chu kì (seasional), khuynh hướng (trend) và thay đổi ngẫu nhiên (random).

Cám ơn các anh chị, các bạn đã dành thời gian đọc bài viết. Mong nhận được những góp ý, bàn luận để hiểu biết hơn về chủ đề này.