Getting started with midasINLA

Model

Consider a response variable \(y_{it}\), indexed by spatial unit \(i=1,\ldots,N\) and low-frequency time point \(t=1,\ldots,T\). The predictor \(x_{i\tau}\) is observed at a higher frequency. The MIDAS framework relates the low-frequency response to multiple high-frequency observations through a weighted distributed lag:

\[ \begin{aligned} &y_{it} \sim F, \;\;\; \mathbb{E}(y_{it}) = \mu_{it} \\ &g(\mu_{it}) = \beta_0 + \beta_i\sum_{k=0}^K w_kx_{i,s(t)-k} + \epsilon_{it} \\ &w_k=h(\boldsymbol{\gamma},k) \; \text{and} \; \sum_{k=0}^K w_k=1. \end{aligned} \]

Here, \(g(\cdot)\) is the link function, \(\beta_0\) is the intercept, and \(\beta_i\) is the regression coefficient for the high-frequency predictor for the \(i^{\text{th}}\) spatial unit. The function \(h(\boldsymbol{\gamma},k)\) determines the MIDAS lag weights, where \(\boldsymbol{\gamma}\in\mathbb{R}^d\) contains the parameters governing the shape of the weighting function. The weights are constrained to sum to one, which separates the overall magnitude of the predictor effect, represented by \(\beta_i\), from the relative contributions of the individual lags.

The index \(s(t)\) denotes the cumulative number of high-frequency observations up to low-frequency time point t,

\[ s(t)=\sum_{j=1}^{t}m_j, \]

where \(m_t\) is the number of high-frequency observations associated with the \(t^{\text{th}}\) low-frequency observation. In the examples below, the same high-frequency sampling structure is assumed across spatial units.

The main flexibility of the model comes from the choice of the lag-weight function \(h(\boldsymbol{\gamma},k)\). midasINLA provides functions for constructing different MIDAS weighting schemes and incorporating them into an INLA model, while allowing the regression coefficient to be either constant or spatially varying.

Spatially varying coefficients

midasINLA allows the regression coefficient \(\beta_i\) to vary across spatial units in two ways:

Spatially structured coefficient: \[ \beta_i = \beta^* + b_i, \qquad \boldsymbol{b} \sim \operatorname{iCAR}(\tau_b,\mathbf{W}), \] where \(\beta^*\) is the overall effect and \(b_i\) represents the spatially structured deviation for unit \(i\). The deviations follow an intrinsic conditional autoregressive (iCAR) model with precision \(\tau_b\) and spatial weights matrix \(\mathbf{W}\), subject to the sum-to-zero constraint \(\sum_i b_i=0\).
Unstructured coefficient: \[ \beta_i \overset{iid}{\sim} N(0,\sigma_\beta^2), \] where the coefficients are independent across spatial units.

A constant coefficient, \(\beta_i\equiv\beta\), is also supported and corresponds to the special case in which the effect is the same across all spatial units.

Constraint functions

The lag weights are obtained by normalising a constraint function \(\psi(\boldsymbol{\gamma},k)\):

\[ w_k = h(\boldsymbol{\gamma},k) = \frac{\psi(\boldsymbol{\gamma},k)} {\sum_{j=0}^K \psi(\boldsymbol{\gamma},j)}. \]

midasINLA implements several commonly used constraint functions.

Exponential Almon polynomial (order 2):

\[ \psi(\boldsymbol{\gamma},k) = \exp\left(\sum_{j=1}^{2}\gamma_j k^j\right), \]

where \(\boldsymbol{\gamma}=(\gamma_1,\gamma_2)\).

Beta polynomial:

\[ \psi(\boldsymbol{\gamma},k) = x_k^{\gamma_1-1}(1-x_k)^{\gamma_2-1}, \]

where

\[ x_k=\xi+(1-2\xi)\frac{k}{K}, \]

with \(\xi>0\) a small fixed constant and \(\boldsymbol{\gamma}=(\gamma_1,\gamma_2)\). A one-parameter version is obtained by fixing \(\gamma_1=1\).

Hyperbolic scheme:

\[ \psi(\gamma,k) = \frac{\Gamma(k+\gamma)}{\Gamma(k+1)\Gamma(\gamma)}, \] where \(\gamma>0\).

Gaussian kernel:

\[ \psi(\boldsymbol{\gamma},k)=\exp\left\{-\frac{(k-\gamma_1)^2}{2\gamma_2}\right\}. \]

where \(\boldsymbol{\gamma}=(\gamma_1,\gamma_2)\) and \(\gamma_2>0\).

INLA implementation

The MIDAS lag structure is incorporated into the latent Gaussian model through INLA’s rgeneric interface. The MIDAS constraint functions define the lag weights as a function of a low-dimensional parameter vector, while the resulting weighted high-frequency covariates are represented as part of the latent model.

The functions in midasINLA construct the required rgeneric model components and interface them with INLA::inla(). This allows the MIDAS lag-weight parameters and regression coefficients to be estimated within the INLA framework, while retaining the spatial structure specified for the regression coefficients.

Users do not need to construct the rgeneric model directly; this is handled internally by the package functions demonstrated below.

Simulated spatial Poisson example

We show an example with two high-frequency covariates with different lag-weight constraints. The first covariate has a spatially varying coefficient, whereas the second has a constant coefficient.

The example illustrates the main functions in midasINLA:

prepare_Minla_spatial() prepares high-frequency covariates for inclusion in the model;
fit_Minla_spatial() fits the resulting model using INLA;
compute_beta_spatial() obtains posterior summaries of the regression coefficients;
compute_weights() obtains posterior summaries of the MIDAS lag weights;
predict_midas() generates posterior predictions.

We consider an outcome \(y_{it}\) observed at 16 spatial locations and 192 time points. Two high-frequency covariates, \(x_{1it}\) and \(x_{2it}\), are available for each location. There are 30 high-frequency observations corresponding to each response time point.

The first covariate uses a hyperbolic lag-weight constraint with \(\gamma = 0.9\) and 29 lags. Its regression coefficient varies spatially according to an intrinsic conditional autoregressive (iCAR) model.

The second covariate uses a Gaussian lag-weight constraint with \(\gamma_1 = 10\) and \(\sqrt{\gamma_2} = 12\) and 45 lags. Its regression coefficient is constant across locations.

The simulated data are included with the package and can be loaded using:

data("data_spatialpoisson_example")

The dataset is provided as a list containing the response, two high-frequency covariates, the spatial polygons, and the true parameter values used to generate the data. The available components can be inspected with:

names(data_spatialpoisson_example)
#>  [1] "data_x1"  "data_x2"  "data_y"   "weights1" "weights2" "eta"     
#>  [7] "beta0"    "beta1"    "beta2"    "icar"     "tau"      "grid_sf"

The response data are stored in data_y and contain the outcome together with the spatial and temporal indices:

head(data_spatialpoisson_example[["data_y"]])
#>    y loc Time
#> 1 NA   1    1
#> 2 75   1    2
#> 3 21   1    3
#> 4 51   1    4
#> 5 63   1    5
#> 6  9   1    6

The high-frequency covariates are stored in data_x1 and data_x2. Each contains the covariate values together with their corresponding spatial indices:

head(data_spatialpoisson_example$data_x1)
#>           x1 loc Time
#> 1 5.42730258   1    1
#> 2 0.02233117   1    2
#> 3 0.56899585   1    3
#> 4 0.36948942   1    4
#> 5 2.36128703   1    5
#> 6 0.19844343   1    6
head(data_spatialpoisson_example$data_x2)
#>           x2 loc Time
#> 1 -0.8913253   1    1
#> 2  0.8279532   1    2
#> 3  1.7493964   1    3
#> 4  0.1564567   1    4
#> 5 -1.3907526   1    5
#> 6 -0.1627427   1    6

The neighbourhood structure used to generate the spatially varying coefficient is stored as an inla.graph object. It can be loaded from the package using:

g <- INLA::inla.read.graph(
  filename = system.file("map.adj", package = "midasINLA")
)

The data were generated using a Poisson model of the form

\[ y_{it} \sim \operatorname{Poisson}(\mu_{it}), \]

with

\[ \log(\mu_{it}) = \beta_0 + (\beta_1^* + b_i) \sum_{k=0}^{29} w_{1k}x_{1i,s(t)-k} + \beta_2 \sum_{k=0}^{45} w_{2k}x_{2i,s(t)-k}, \]

where \(\boldsymbol{b}\) follows an iCAR model. The true values of the parameters are \(\beta_0 = 1\), \(\beta_1^* = 1.1\), and \(\beta_2 = -2\).

Preparing the MIDAS covariates

Before fitting the model, each high-frequency covariate is prepared using prepare_Minla_spatial().

For the first covariate, we use the hyperbolic constraint and allow the coefficient to vary spatially according to an iCAR model:

Midas_x1 <- prepare_Minla_spatial(
  x = data_spatialpoisson_example$data_x1$x1,
  loc_x = data_spatialpoisson_example$data_x1$loc,
  constraint = "hyperbolic",
  K = 0:29,
  m = 30,
  svc = TRUE,
  svc_prior = "icar",
  g = g
)

Midas_x2 <- prepare_Minla_spatial(
  x = data_spatialpoisson_example$data_x2$x2,
  loc_x = data_spatialpoisson_example$data_x2$loc,
  constraint = "gaussian",
  K = 0:45,
  m = 30,
  svc = FALSE
)

The resulting objects contain the information required by fit_Minla_spatial() to construct the MIDAS components of the model.

Fitting the model

The response data are stored in data_y. To illustrate prediction, we reserve the final 10 response time points at each location as a test set.

response_data <- data_spatialpoisson_example[["data_y"]]

response_data$y_all <- response_data$y

response_data[which(response_data[["Time"]] %in% 183:192),"y"] <- NA

The model can then be fitted using fit_Minla_spatial():

fit_res <- fit_Minla_spatial(
  formula = y ~ 1,
  data = response_data,
  loc_var = "loc",
  time_var = "Time",
  family = "poisson",
  hf_input = list(Midas_x1, Midas_x2),
  inla_options = list(verbose = FALSE,
                      control.predictor = list(
                        compute = TRUE,link = 1)))

A summary of the fitted model can be obtained using the standard summary() method:

summary(fit_res[["res"]])
#> Time used:
#>     Pre = 0.947, Running = 20.1, Post = 0.242, Total = 21.3 
#> Fixed effects:
#>              mean   sd 0.025quant 0.5quant 0.975quant  mode kld
#> (Intercept) 0.997 0.02      0.957    0.997      1.037 0.997   0
#> 
#> Random effects:
#>   Name     Model
#>     hf_idx_1 RGeneric2
#>    hf_idx_2 RGeneric2
#> 
#> Model hyperparameters:
#>                      mean    sd 0.025quant 0.5quant 0.975quant  mode
#> Theta1 for hf_idx_1  2.17 0.076       2.02     2.16       2.32  2.15
#> Theta2 for hf_idx_1  2.28 0.389       1.53     2.28       3.06  2.26
#> Theta1 for hf_idx_2 -1.31 0.035      -1.38    -1.31      -1.24 -1.31
#> Theta2 for hf_idx_2  2.50 0.023       2.46     2.50       2.55  2.50
#> Theta3 for hf_idx_2 -1.99 0.023      -2.04    -1.99      -1.95 -1.99
#> 
#> Marginal log-Likelihood:  -8839.62 
#>  is computed 
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

The fitted model contains two MIDAS components, corresponding to the two high-frequency covariates. The model also estimates the spatially varying coefficient associated with the first covariate. The parameters governing the MIDAS lag-weight functions are reported with names beginning with hf_idx_.

Posterior summaries of regression coefficients

Posterior summaries of the MIDAS regression coefficients can be obtained using compute_beta_spatial().

beta_results <- compute_beta_spatial(
  model = fit_res,
  n_loc = 16
)

The returned object contains results for each high-frequency covariate, including marginal distributions and posterior summaries.

For the first covariate, the posterior summaries of the spatially varying component \(b_i\) can be accessed using:

beta_results$hf_index_1$summary.icar.beta
#>             Mean          SD         2.5%          50%       97.5%
#> b1   0.144486652 0.004972584  0.134365723  0.144546535  0.15422010
#> b2   0.137016169 0.005420709  0.125940897  0.137002802  0.14782478
#> b3   0.006070001 0.006197002 -0.005814864  0.006151536  0.01798475
#> b4  -0.223484377 0.007895843 -0.238952783 -0.223766772 -0.20783057
#> b5   0.137748078 0.005092015  0.127834969  0.137790783  0.14741063
#> b6   0.033968095 0.006056428  0.022138423  0.033911399  0.04574504
#> b7  -0.105208311 0.006596310 -0.118723432 -0.105307289 -0.09293775
#> b8   0.040061614 0.006011510  0.028007704  0.040201841  0.05252542
#> b9   0.125352055 0.005744286  0.113941288  0.125605270  0.13634697
#> b10  0.019680601 0.005916555  0.008263325  0.019732012  0.03122072
#> b11 -0.059567582 0.007141423 -0.073331629 -0.059312934 -0.04635210
#> b12 -0.170726979 0.007140385 -0.185612277 -0.170432548 -0.15741375
#> b13 -0.140813435 0.007064259 -0.155242567 -0.141043980 -0.12677287
#> b14 -0.132180330 0.007114006 -0.146294796 -0.132004402 -0.11810291
#> b15 -0.059499088 0.006540695 -0.072403182 -0.059393586 -0.04611928
#> b16  0.245896571 0.004794204  0.236833718  0.245748877  0.25544649

The spatially varying coefficient at location \(i\) is defined as

\[ \beta_{1,i} = \beta_1^* + b_i. \] Posterior summaries of the resulting total coefficient can be accessed using:

beta_results$hf_index_1$summary.total.beta
#>             Mean         SD      2.5%       50%     97.5%
#> beta1  1.2457216 0.01048782 1.2237387 1.2456787 1.2655496
#> beta2  1.2375718 0.01067897 1.2157503 1.2372824 1.2584179
#> beta3  1.1076913 0.01101020 1.0862391 1.1076086 1.1289906
#> beta4  0.8769760 0.01248564 0.8503544 0.8775068 0.9001637
#> beta5  1.2385937 0.01046528 1.2180243 1.2383180 1.2601783
#> beta6  1.1353647 0.01135307 1.1128136 1.1350674 1.1568555
#> beta7  0.9958926 0.01155036 0.9722900 0.9961221 1.0180031
#> beta8  1.1406232 0.01155987 1.1170668 1.1404468 1.1636177
#> beta9  1.2260177 0.01085150 1.2050063 1.2260251 1.2472006
#> beta10 1.1211220 0.01141409 1.0982131 1.1215215 1.1432208
#> beta11 1.0412863 0.01190198 1.0182609 1.0411389 1.0650395
#> beta12 0.9302476 0.01206977 0.9060055 0.9300331 0.9532307
#> beta13 0.9596178 0.01191715 0.9357719 0.9594845 0.9815521
#> beta14 0.9687006 0.01224003 0.9446940 0.9685399 0.9931769
#> beta15 1.0410059 0.01176284 1.0171099 1.0406877 1.0654608
#> beta16 1.3467105 0.01041105 1.3274772 1.3469468 1.3665332

For the second covariate, which has a constant regression coefficient, posterior summaries of \(\beta_2\) are available using:

beta_results$hf_index_2$summary.beta
#>           Mean         SD      2.5%      50%     97.5%
#> beta -1.991056 0.02222116 -2.033759 -1.99112 -1.948287

Estimating the MIDAS lag weights

The posterior distributions of the lag weights can be obtained using compute_weights():

res_weights <- compute_weights(fit_res)

The result is a list containing one data frame for each high-frequency covariate. Each data frame contains the lag, posterior mean, and lower and upper posterior quantiles.

For example, the estimated weights for the first covariate are:

head(res_weights$hf_1)
#>   lag       mean       q2.5      q97.5
#> 1   0 0.04548848 0.04367942 0.04749251
#> 2   1 0.04081357 0.03977345 0.04194673
#> 3   2 0.03871741 0.03799511 0.03949764
#> 4   3 0.03739207 0.03686256 0.03796024
#> 5   4 0.03643224 0.03603847 0.03685206
#> 6   5 0.03568417 0.03539393 0.03599141

The weights for the second covariate can be inspected similarly:

head(res_weights$hf_2)
#>   lag       mean       q2.5      q97.5
#> 1   0 0.03021630 0.02896381 0.03163707
#> 2   1 0.03212713 0.03091376 0.03361481
#> 3   2 0.03392939 0.03262788 0.03538489
#> 4   3 0.03559200 0.03429011 0.03701866
#> 5   4 0.03708518 0.03567537 0.03855939
#> 6   5 0.03838124 0.03695189 0.03988035

The posterior summaries can be used to visualise the estimated lag-weight functions. Because the data are simulated, the true lag weights are also available for comparison.

ggplot(res_weights$hf_1, aes(x = lag, y = mean)) +
  geom_errorbar(
    aes(
      ymin = q2.5,
      ymax = q97.5
    ),
    width = 0.2,
    colour = "grey30"
  ) +
  geom_point(
    aes(
      colour = "Posterior mean"
    ),
    size = 2
  ) +
  geom_point(
    aes(
      y = data_spatialpoisson_example$weights1,
      colour = "True value"
    ),
    size = 2
  ) +
  scale_colour_manual(
    name = NULL,
    values = c(
      "Posterior mean" = "red",
      "True value" = "blue"
    )
  ) +
  labs(
    x = "Lag",
    y = "Lag weight"
  ) +
  theme_bw() +
  theme(
    legend.position = "bottom"
  )

The corresponding lag-weight function for the second covariate can be visualised in the same way:

ggplot(res_weights$hf_2, aes(x = lag, y = mean)) +
  geom_errorbar(
    aes(
      ymin = q2.5,
      ymax = q97.5
    ),
    width = 0.2,
    colour = "grey30"
  ) +
  geom_point(
    aes(
      colour = "Posterior mean"
    ),
    size = 2
  ) +
  geom_point(
    aes(
      y = data_spatialpoisson_example$weights2,
      colour = "True value"
    ),
    size = 2
  ) +
  scale_colour_manual(
    name = NULL,
    values = c(
      "Posterior mean" = "red",
      "True value" = "blue"
    )
  ) +
  labs(
    x = "Lag",
    y = "Lag weight"
  ) +
  theme_bw() +
  theme(
    legend.position = "bottom"
  )

Prediction for held-out observations

The final 10 time points were withheld from the model fit and are used here to illustrate posterior prediction.

Posterior predictions can be generated using predict_midas():

pred_res <- predict_midas(
  model = fit_res,
  family = "poisson",
  Ntrials = NULL,
  nsamples = 1000
)

The returned object contains posterior summaries of the predicted outcome, including posterior means and 95% credible intervals, as well as posterior samples of the latent predictor. For example:

str(pred_res, max.level = 2)
#> List of 2
#>  $ computed_y:List of 4
#>   ..$ mean : num [1:3056] 91.8 20.4 54.6 77.1 14.7 ...
#>   ..$ sd   : num [1:3056] 9.91 4.43 7.41 8.88 3.85 ...
#>   ..$ q2.5 : num [1:3056] 73 12 41 60 8 ...
#>   ..$ q97.5: num [1:3056] 112 29 70 95 23 ...
#>  $ samples   :List of 1
#>   ..$ latent_predictor: num [1:3056, 1:1000] 4.51 3.01 3.99 4.33 2.69 ...

The posterior summaries of the predicted outcome can be accessed using:

head(pred_res$computed_y$mean)
#> [1] 91.794 20.450 54.627 77.073 14.741 11.797
head(pred_res$computed_y$q2.5)
#> [1] 73.000 12.000 40.975 60.000  8.000  6.000
head(pred_res$computed_y$q97.5)
#> [1] 112.000  29.000  70.000  95.025  23.000  19.000

The predicted and observed outcomes can also be compared graphically. The following example shows the results for the first four spatial locations. The dashed vertical line indicates the boundary between the training and held-out prediction periods.

plot_data <- data.frame(
  observed = fit_res$data_final$y_all,
  predicted = pred_res$computed_y$mean,
  lower = pred_res$computed_y$q2.5,
  upper = pred_res$computed_y$q97.5,
  loc = fit_res$data_final$loc,
  Time = fit_res$data_final$Time
)

plot_long <- plot_data |>
  dplyr::filter(loc %in% 1:4) |>
  tidyr::pivot_longer(
    cols = c(observed, predicted),
    names_to = "series",
    values_to = "value"
  )

# Determine the training/held-out boundary for each location
non_na <- !is.na(fit_res$data_final$y)

segment <- cumsum(
  non_na != dplyr::lag(non_na, default = TRUE)
)
segment[!non_na] <- NA

rel_idx <- ave(
  seq_along(non_na),
  segment,
  FUN = seq_along
)

first_na <- which(
  diff(c(FALSE, is.na(fit_res$data_final$y))) == 1
)

vlines <- data.frame(
  cut = rel_idx[first_na - 1] + 1,
  loc = seq_len(16)
) |>
  dplyr::filter(loc %in% 1:4)

ggplot(plot_long, aes(x = Time, y = value, colour = series)) +
  geom_ribbon(
    data = plot_data |>
      dplyr::filter(loc %in% 1:4),
    aes(
      x = Time,
      ymin = lower,
      ymax = upper
    ),
    inherit.aes = FALSE,
    fill = "red",
    alpha = 0.2
  ) +
  geom_line() +
  geom_vline(
    data = vlines,
    aes(xintercept = cut),
    colour = "black",
    linetype = "dashed"
  ) +
  facet_wrap(
    ~loc,
    ncol = 2,
    labeller = labeller(
      loc = function(x) paste("Loc =", x)
    )
  ) +
  scale_colour_manual(
    values = c(
      "observed" = "blue",
      "predicted" = "red"
    )
  ) +
  labs(
    x = "Time",
    y = "Outcome",
    colour = NULL
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom"
  )