Initial model: Total hospital admissions
We start the analysis with a model for the total number of hospital admissions over time in Brazil, i.e., a temporal model that does not consider the hospital admissions on each state.
From the figure above we can see that there is a consistent trend of (log) linear decay of the rate of hospital admissions over time, until April of 2020, when there is an abrupt reduction of hospital admission due to the pandemic of COVID-19 (Ribeiro et al., 2022), which is then followed by what seems to be a return to the previous level, although more observations would be necessary to be sure. We can also note that the data has a clear seasonal pattern, with a period of \(12\) months.
We begin with a very simply model. Let \(Y_t\) be the total number of hospital admissions on Brazil at time \(t\). Assume that:
\[ \begin{aligned} Y_t|\eta_t &\sim Poisson(\eta_t)\\ \ln{\eta_t}&=\lambda_t=\theta_{1,t}\\ \theta_{1,t}&=\theta_{1,t-1}+\theta_{2,t-1}+\omega_{1,t}\\ \theta_{2,t}&=\theta_{2,t-1}+\omega_{2,t}.\\ \end{aligned} \]
First we define the model structure:
Then we define the outcome:
Then we fit the model:
Finally, we can see how our model performed with the and methods:
Fitted DGLM with 1 outcomes.
distributions:
Series.1: Poisson
---
No static coeficients.
---
See the coef.fitted_dlm for the coeficients with temporal dynamic.
One-step-ahead prediction
Log-likelihood : -24560.92
Interval Score : 52492.80142
Mean Abs. Scaled Error: 1.41602
---
Clearly the model described above is too simple to describe the data. In particular, it does not take into account any form of seasonal pattern. Let us proceed then by assuming the following model:
\[ \begin{aligned} Y_t|\eta_t &\sim Poisson(\eta_t)\\ \ln{\eta_t}&=\lambda_t=\theta_{1,t}+\theta_{3,t}\\ \theta_{1,t}&=\theta_{1,t-1}+\theta_{2,t-1}+\omega_{1,t}\\ \theta_{2,t}&=\theta_{2,t-1}+\omega_{2,t},\\ \begin{bmatrix}\theta_{3,t}\\\theta_{4,t}\end{bmatrix}&=R\begin{bmatrix}\theta_{3,t}\\\theta_{4,t}\end{bmatrix}+\begin{bmatrix}\omega_{3,t}\\\omega_{4,t}\end{bmatrix}\\ R&=\begin{bmatrix} \cos(\frac{2\pi}{12}) &\sin(\frac{2\pi}{12})\\ -\sin(\frac{2\pi}{12}) & \cos(\frac{2\pi}{12})\end{bmatrix} \end{aligned} \]
Where \(R\) is a rotation matrix with angle \(\frac{2\pi}{12}\), such that \(R^{12}\) is equal to the identity matrix.
To define the structure of that model we can use the
harmonic_block function alongside the
polynomial_block function:
structure <- polynomial_block(
rate = 1, order = 2, D = c(0.95, 0.975),
name = "Trend"
) +
harmonic_block(
rate = 1, period = 12, D = 0.98,
name = "Season"
)Then we fit the model (using the previously defined outcome):
Fitted DGLM with 1 outcomes.
distributions:
Series.1: Poisson
---
No static coeficients.
---
See the coef.fitted_dlm for the coeficients with temporal dynamic.
One-step-ahead prediction
Log-likelihood : -17339.1
Interval Score : 41481.46809
Mean Abs. Scaled Error: 0.70882
---
Notice that this change significantly improves all metrics provided by the model summary, which indicates that we are going in the right direction. We encourage the reader to test different orders for the harmonic block.
The previous model could capture the mean behavior of the series reasonably well. However, two deficiencies of that model standout: first, the overconfidence in the predictions, evidenced by the particularly thin credibility interval; nnd second, the difficulty the model had to adapt to the pandemic period.
The first problem comes from the fact that we are using a Poisson model, which implies that \(Var[Y_t|\eta_t]=\mathbb{E}[Y_t|\eta_t]\), which means that \(Var[Y_t]=\mathbb{E}[Var[Y_t|\eta_t]]+Var[\mathbb{E}[Y_t|\eta_t]]=\mathbb{E}[\eta_t]+Var[\eta_t]\). For latter observations we expect \(Var[\eta_t]\) to be relatively small; as such, the variance of \(Y_t\) should be very close to its mean after a reasonable amount of observations. In this scenario, the coefficient of variation, defined as \(\frac{\sqrt{Var[Y_t]}}{\mathbb{E}[Y_t]}\) goes to \(0\) as \(\mathbb{E}[Y_t]\) grows, in particular, for data in the scale we are working with in this particular problem, we would expect a very low coefficient of variation if the Poisson model were adequate, but that is not what we observe. This phenomena is called and is a well known problem in literature . To solve it, we can include a block representing a white noise that is added to the linear predictor at each time, but does not affect previous or future observation, so as to capture the overdipersion. In this case, we will assume the following model:
\[ \begin{aligned} Y_t|\eta_t &\sim Poisson(\eta_t)\\ \ln{\eta_t}&=\lambda_t=\theta_{1,t}+\theta_{3,t}+\epsilon_t\\ \theta_{1,t}&=\theta_{1,t-1}+\theta_{2,t-1}+\omega_{1,t}\\ \theta_{2,t}&=\theta_{2,t-1}+\omega_{2,t},\\ \begin{bmatrix}\theta_{3,t}\\\theta_{4,t}\end{bmatrix}&=R\begin{bmatrix}\theta_{3,t}\\\theta_{4,t}\end{bmatrix}+\begin{bmatrix}\omega_{3,t}\\\omega_{4,t}\end{bmatrix}\\ \epsilon_t & \sim \mathcal{N}(0,\sigma_t^2) \end{aligned} \]
This structure can be defined using the function, alongside the previously used functions:
structure <- polynomial_block(
rate = 1, order = 2, D = c(0.95, 0.975),
name = "Trend"
) +
harmonic_block(
rate = 1, period = 12, D = 0.98,
name = "Season"
) +
noise_block(rate = 1, name = "Noise")For the second problem, that of slow adaptation after the start of the pandemic. The ideal approach would be to make an intervention, increasing the uncertainty about the latent states at the beginning of the pandemic period and allowing our model to quickly adapt to the new scenario (see West and Harrison, 1997, Chapter 11). We recommend this approach when we already expect a change of behavior in a certain time, even before looking at the data (which is exactly the case). Still, for didactic purposes, we will first present how the automated monitoring can also be used to solve this same problem. In general, we recommend the automated monitoring approach when we do not known if or when a change of behavior happened before looking at the data, i.e., we do not known of any particular event that we expect to impact our outcome.
Following what was presented in the Subsection Intervention and monitoring, we can use the following code to fit our model:
structure <- polynomial_block(
rate = 1, order = 2, D = c(0.95, 0.975),
name = "Trend", monitoring = c(TRUE, TRUE)
) +
harmonic_block(
rate = 1, period = 12, D = 0.98,
name = "Season"
) +
noise_block(rate = 1, name = "Noise")Notice that we set the monitoring of the
polynomial_block to c(TRUE,TRUE). By default,
the polynomial_block function only activates the monitoring
of its first component (the level), but, by the visual analysis made at
the beginning, it is clear that the pandemic affected both the level and
the slope of the average number of hospital admissions, as such, we
would like to monitor both parameters.
# To activate the automated monitoring it is enough to set the p.monit argument to a valid value
fitted.model <- fit_model(structure, outcome, p.monit = 0.05)
summary(fitted.model)Fitted DGLM with 1 outcomes.
distributions:
Series.1: Poisson
---
No static coeficients.
---
See the coef.fitted_dlm for the coeficients with temporal dynamic.
One-step-ahead prediction
Log-likelihood : -1283.634
Interval Score : 14965.54610
Mean Abs. Scaled Error: 0.71239
---
The summary presented above shows a massive improvement in the comparison metrics with the new changes introduced. Moreover, we can see that the automated monitoring detected the exact moment where the series \(Y_t\) changed behavior, which allowed the model to immediately adapt to the pandemic period.
One aspect of the model that may bother the reader is the exceedingly
high uncertainty at the first observations. This behavior is duo to our
approach to the estimation of the variance of the white noise introduced
by the noise_block function (see dos Santos et al., 2024 and the
associated documentation for details), which can be a bit too
sensitive to bad prior specification at the initial steps. As such, we
highly recommend the user to perform a sensitivity analysis to choose
the initial variance of the white noise:
structure <- polynomial_block(
rate = 1, order = 2, D = c(0.95, 0.975),
name = "Trend"
) +
harmonic_block(
rate = 1, period = 12, D = 0.98,
name = "Season"
) +
noise_block(rate = 1, R1 = "H", name = "Noise") # Setting the initial variance as a unknown parameter
structure <- structure |>
intervention(time = 124, var.index = c(1:2), D = 0.005)
search.model <- fit_model(
structure, outcome,
lag = -1, # Using the model likelihood (f(y|M)) as the comparisson metric.
H = seq.int(0, 0.1, l = 101)
)
fitted.model <- search.model$modelNotice that, this time around, we chose to make an intervention at the beginning of the pandemic, instead of an automated approach. As mentioned before, this approach is preferable in this scenario, since we were aware that the pandemic would affect our outcome before even looking at the data.
Fitted DGLM with 1 outcomes.
distributions:
Series.1: Poisson
---
No static coeficients.
---
See the coef.fitted_dlm for the coeficients with temporal dynamic.
One-step-ahead prediction
Log-likelihood : -1220.102
Interval Score : 7470.8794
Mean Abs. Scaled Error: 0.5993
---Again, the new changes improve the comparison metrics even further, leading to the conclusion that our last model is the best among those presented until now. We highly encourage the reader to run this example and experiment with some of the options the kDGLM package offers, but that were not explored, such as changing the discount factors used in each block, the order of the blocks, adding/removing structural components, etc..
As a last side note, the user may not like the approach of choosing a
specific value for the initial variance of the white noise introduced by
the noise_block. Indeed, one may wish to define a prior
distribution for this parameter and estimate it along with the others.
While we will not detail this approach for the sake of brevity (since it
is not directly supported), we would like to point out that we do offer
tools to facilitate this procedure:
search.result <- search.model$search.data[order(search.model$search.data$H), ]
H.vals <- search.result$H
log.prior <- dgamma(H.vals, 1, 1, log = TRUE)
log.like <- search.result$log.like
l.fx <- log.prior + log.like
pre.fx <- exp(l.fx - max(l.fx))
fx <- pre.fx / sum(pre.fx * (H.vals[2] - H.vals[1]))
plot(H.vals, fx,
type = "l", xlab = "H", ylab = "Density",
main = "Posterior density for the unknown hyperparameter H"
)
Since there are