 
\section{Modelling COVID--19 ventilation bed prevalence}
\label{sec: application}
Define the multivariate time series
$\Tilde{X}_{i, t}$ 
to be the recorded
number of COVID--19 patients  occupying  mechanical ventilation beds
for day $t$, for NHS Trust $i$, across $d=140$ trusts in England. 
The $T=452$ \mbox{$d$-dimensional}
observations were obtained from the UK Government Coronavirus Dashboard
({\tt coronavirus.data.gov.uk}) between April 2020 and July 2021.
Patients who  occupy  mechanical ventilation beds
are a subset of those COVID-positive patients who are, or become, seriously ill and require
artificial ventilation. Some of these patients are directly transferred
in on ventilation from other hospitals and care settings.
\par
The $\Tilde{X}_{i, t}$ observations are count data and it is not appropriate to use regular GNAR
models directly and so, to stabilise variance and bring the data closer to normality, we study
the well--known and understood transformed version:
$X_{i, t} := \log\, ( 1 + \Tilde{X}_{i, t})$. We stack the individual trust $X_{i, t}$
time series and denote $\boldsymbol{X}_t = (X_{1, t}, \cdots, X_{140, t} ),$
for $t=1, \ldots, T=452$.

\par

In reality, there is no tangible single network connecting NHS Trusts.
There are many ways a reasonable network could be constructed for modelling purposes
and we outline one approach here. 
Our network was built using  geographical `as the crow flies'
distances between NHS trusts.  We define a $D$--connected network to be one where Trusts are
directly connected, if they are less than $D$~km apart. To find a good value of $D$ we
computed several $D$-connected networks with different values of $D >0$. For each one
we computed the residual sum of squares using a
full fitted GNAR model, GNAR$(12, [4^{(12)}])$. The network that achieved the smallest residual sum of
squares was with $D=120$~km (approx 75~miles). A full GNAR model was used to capture any structure within a large
possible number of lags and neighbour stages, any of which might be statistically significant in later model
fitting with our chosen network (i.e.,\ as we do not want to inadvertently rule out
 detailed  structure at this point). There are many reasonable networks that could be used
instead, including those that delineate communications/transmission
links more comprehensively. Our choice is intuitively attractive
and can be seen to be a generalisation of partial correlation networks,
see \cite{partialcorrnetworks}, for example.


The $X_{i, t}$ series for Trust $i$ is associated
with vertex $i$ of our Trust network, which is depicted in Figure \ref{fig: rstage neighbors}a.
%Figure~\ref{fig: correlation structure}a shows strong correlations between series at different Trusts.
Our underlying assumption is that the behaviour of COVID-19 infection does not change drastically across
location, and that the needs of similar hospitals can be described by the parsimonious GNAR model, which enables us to exploit the underlying structure.
Hence, we propose modelling $\boldsymbol{X}_t$ as a stationary {\em global-$\alpha$} GNAR process~\eqref{eq: vector wise}. 

GNAR parsinomity allows us to use significantly
more observations per parameter when estimating the autoregressive coefficients rather than separately for each node,
as would be required
for a (overparameterized) VAR model or  even competitor CAR-like models . 
For instance, fitting a VAR$(p)$ to $\boldsymbol{X}_t$ requires estimating (or, at the very least, having to consider) $p\, \times\,140^2=19600\,\times\,p$ unknown parameters, whereas, fitting a $\GNAR(p, [s_1, \dots, s_p])$ requires estimating $p + \sum_{k= 1}^{p} s_k$ parameters, which is upper--bounded by $p\,(1 + 6)$ given that the maximum $\rstage$ per lag is six and we are modelling the data as a global--$\alpha$ GNAR process.

%Before fitting the global-$\alpha$ GNAR model and selecting the order, we pause and mention that the data under study are counts which are likely to violate our assumptions of stationarity, however, we model the log-transformed count data as a stationary process in the interest of illustrating GNAR robustness to model misspecification.

\subsection{GNAR Model Selection}
Using the ideas from Section~\ref{sec:corbit}, we focus on the correlation structure in $\{ X_{i, t} \}$ and employ Corbit plots as graphical aids for assisting model selection. The observed NACF and PNACF values for our COVID-19 network time series are shown in Figures~\ref{fig: corbits}(a) and (b). 
 \begin{figure}
        \centering
        \begin{subfigure}{0.7\textwidth}
            \includegraphics[width=\textwidth]{images/covid_nacf.pdf}
            \caption{COVID-19 series NACF Corbit plot.}
            \label{fig: corbit nacf}
        \end{subfigure}
        \hfill
        \begin{subfigure}{0.7\textwidth}
            \includegraphics[width=\textwidth]{images/covid_pnacf.pdf}
            \caption{COVID-19 series PNACF Corbit plot.}
            \label{fig: corbit pnacf}
        \end{subfigure}
        \caption{Corbit plots for (a) observed NACF, and (b) for observed PNACF with respect to the COVID-19 series; maximum lag is equal to 20 and maximum $\rstage$ depth is equal to six. Plot (b) suggests a strong autoregressive behavior, which cuts-off after the first lag.}
        \label{fig: corbits}
    \end{figure}
Figure~\ref{fig: corbit pnacf} shows that the partial network autocorrelation cut--offs after the first time lag
and sharply decays after the first stage at the first time lag.
It remains roughly constant at the first lag for stages larger than or equal to two at lag~$1$, and cuts-off sharply at all stages for all lags $h \geq 2$.
Hence, our Corbit plot suggests looking at $1$--lag models for which there does not appear to be a large contribution to autocorrelation from $\rstage$ neighbours beyond the first stage. Looking at Corbit plots is considerably faster than a full evaluation of information criteria values such as
AIC or BIC, as mentioned earlier.
\par
We proceed by fitting a $\GNAR(1, [1])$ model and then compare it to $\GNAR(1, [6])$ and   $\GNAR(6, [6^{(6)}])$  models,
the latter being the one with the largest possible number of active $\rstage$ regressions with the same number of lags. Table~\ref{tab: gnar comparison} compares results of different GNAR models by looking at
one- and two-step predictive performance. For the two-step
predictive performance each model is fitted to the first $447$ observations of the COVID--19 network time series, and two--step prediction is performed by using the one-step prediction as a 
pseudo--observation. Table~\ref{tab: gnar comparison} shows that the order-6 model is the least--best choice based on one-- and two-step predictive performance. The predictive errors for the remaining first order models are very similar with the $\GNAR(1, [6])$ best for one--step predictive error, very closely followed by the $\GNAR(1, [1])$ model, which is best for
the two-step predictive performance, but there is little to choose between here. The AIC and BIC values are almost identical and, interestingly, select the model with the largest MSPE and number of parameters. However, the difference is less than 10$\%$ between the largest and smallest values for both criteria.
%GNAR models comparison
\begin{table}

\centering
\begin{tabular}{ccccc}
 Model & AIC & BIC & One--step MSPE & Two--step MSPE \\ 
\hline
GNAR(1, [1]) & -452.39 & -452.37 & 5.09 & 8.81 \\
 GNAR(1, [6]) & -452.53 & -452.47 & 5.06 & 8.8 \\
 GNAR(6, [6, 6]) & -467.35 & -466.96 & 5.24 & 9.88
\end{tabular}
\caption{AIC, BIC, One-- and two--step mean--squared prediction error (MSPE)
performance for different GNAR model order choices.}
\label{tab: gnar comparison}

\end{table}

% Table~\ref{tab: gnar comparison} shows that the order-6 model is the least-best choice based on one- and two-step predictive
% performance. The predictive errors for the remaining first order models are very similar with the
% $\GNAR(1, [6])$ best for one-step predictive error, very closely followed by the $\GNAR(1, [1])$ model, which is best for
% the two-step predictive performance, but there is little to choose between here.

\subsection{Comparison to Other Models}
\label{sec: comparison}
The underlying network structure aids not only in proposing a parsimonious GNAR model, it also permits estimation of the autoregressive coefficients more accurately and, we conjecture, reduces generalization error. In the interests of interpretability and illustration we select the  $\GNAR(1, [1])$ model
and compare it to other 

popular models.
Table~\ref{tab: prediction comparison} compares
the predictive performance of the $\GNAR(1, [1])$ fit  to VAR$(1)$, sparse VAR$(1)$, restricted VAR$(1)$, decoupled AR$(1)$ models, and to a conditional autoregressive first order temporal autoregressive (CARar) model, which is fitted using the \texttt{CARBayesST} package; see \cite{carbayes} based on \cite{RUSHWORTH201429}. The comparison performs one--step prediction over the ten most recent times of the COVID--19 (network) time series. See the
RMarkdown document (Section 3)  for further details.

%% Forecasting comparison table
\begin{table}[H]
    \centering
    \begin{tabular}{crr}
        Model & MSPE(sd) & Mean No. Parameters \\
        \hline
         GNAR(1, [1])+ & 7.2(2.82) & 2 \\
         GNAR(1, [1]) &  7.4(2.98) & 2 \\
         CARar(1) & 7.5(3.10) & 4(+140) \\
         Naive & 7.6(3.12) & 0 \\
         Res. VAR(1) & 10.6(2.48) & 3773 \\
         Sparse VAR(1) & 11.4(2.80) & 3089 \\
         VAR(1) & 12.3(3.18) & 19600 \\
         AR(1) & 87.4(5.15) & 140
    \end{tabular}
    \caption{Mean number of parameters and one-step predictive performance for the COVID-19 (network) time series for five different
        time series models over ten predictions. The mean-squared
        prediction error (MSPE) standard deviation is shown within parenthesis. Naive is the observed value at the previous time step. $\GNAR(1, [1])+$ denotes a GNAR$(1, [1])$ model estimated on centred data (i.e., subtracting column means before model fitting). The $(+140)$ for the CARar model corresponds to the
        random effects in that model.}
    \label{tab: prediction comparison}
\end{table}
The models are fitted using the \texttt{GNAR} package \citep{gnar_package} for the $\GNAR(1, [1])$, \texttt{VARS} package \citep{varspack} for fitting the VAR$(1)$ and the restricted VAR$(1)$, \texttt{forecast} package \citep{fcst_pack} for the $140$
individual AR$(1)$ models, and \texttt{sparsevar} \citep{sparsevar} for the Sparse VAR$(1)$. We restricted each model to one lag
for a fair comparison based on the autocorrelation analysis the Corbit plot in Figure \ref{fig: corbit pnacf} provides.
\par
Table~\ref{tab: prediction comparison}  shows that
GNAR's MSPE is about $35\%$ smaller than the nearest time series models Sparse VAR$(1)$ and almost $30\%$ smaller than restricted VAR$(1)$. The centred GNAR is $4\%$ better than CARar$(1)$. In this particular situation CARar and GNAR are similar because CARar is structurally equivalent to a restricted GNAR with correlated error
{\em limited to first stage neighbours}, but effectively having to also estimate $d = 140$ random effects `pseudo--parameters'. However, it should be noted that CARar is significantly more computationally expensive than GNAR and much less flexible. 
We remark that CARar software does not provide functionality for forecasting, and so we  created an {\em ad hoc} procedure that recasts the full Bayesian forecasting problem into a generalised least squares problem, where the correlated errors have the same covariance structure as the random effects. The forecasting performance of CARar using a full Bayesian method might be better or worse, but certainly enormously more computationally expensive, as it requires multiple code invocation  equivalent to the parameter estimation, which is already expensive. 
See Appendix~\ref{supp_sec: GNAR CARar comparison} in the supplementary material for details.


%Further simulations
 \begin{table}[H]
        \begin{tabular}{crrrrrr}
            St. Ahead  & GNAR(1, [1]) & GNAR(1, [1])* & Res. VAR(1) & VAR(1) & Sp.VAR(1) \\ 
            \hline
            one & 2.25 & 2.66 &  6.15 & 10.30 &  6.07 \\
            two & 5.05 & 4.15 &  9.96 & 13.53 & 10.36 \\
            three & 7.98 & 6.26 & 13.80 & 14.84 & 16.52 
        \end{tabular}
        \caption{Three-step ahead squared prediction error model comparison. GNAR$(1, [1])$ denotes global-$\alpha$ and GNAR$(1, [1])*$ denotes local-$\alpha$. The realisation end point taken is $T = 300$, i.e., at the start of the second wave's peak. AR(1) were over 800 in all cases.}
        \label{tab: 3 step ahead at second wave peak}
\end{table}

Table~\ref{tab: 3 step ahead at second wave peak} shows one--, two-- and three--step ahead comparisons for some of the methods described earlier but for an earlier end point at the start of the second wave's peak. It can be seen that GNAR reduces the mean squared prediction error by at least $40\%$ in all cases. We did not include CARar as there is no forecasting functionality in the \texttt{CARBayesST} package, and we have not validated our
{\em ad hoc} method further into the future.

\par
Hence, for the right kind of data, i.e.,\ those for which the GNAR model is (at least approximately) valid, then GNAR forecasting can be extremely competitive. Moreover, model parsinomity eases interpretation of the results: our analysis provides evidence that mechanical ventilation  bed occupancies  `spreads' mostly through direct trust neighbours,
and also on the previous number of beds occupied at the trust. The parameter estimates for the $\GNAR(1, [1])$ were
$\hat{\alpha}_1 \approx 0.95$ and $\hat{\beta}_{1, 1} \approx 0.043$, and they were both statistically significant at the $0.1\%$ level.

For completeness, we also tried fitting local--$\alpha$ models and
predictive power improvement is negligible. The results are shown
in the Appendices as Tables~\ref{tab: pred 443}  to~\ref{tab: pred 452},
where even with an extra 139 parameters the improvement in predictive
error is never more than 3\%, often nearer 1 or 2\%, and in two cases out
of the ten tables the simpler GNAR does better. This modelling   exercise 
provides further validation for our initial global--$\alpha$ choice.

\subsection{Model Interpretation and Analysis}
\label{sec: model interpretation application}
Another attractive property of GNAR models is that they can handle missing data and zero-values by weight readjustment; see \cite{gnar_paper}. We extend this idea to study the correlation structure of the two main COVID--19 outbreaks in England, the first outbreak was from April 2020 to July 2020 and the second one was from September 2020 to July 2021. We visualise the differences in the correlation structure during these different time periods
by looking at R--Corbit plots, such as the one for NACF in Figure \ref{fig: R-Corbits}, the equivalent
plot for PNACF is shown in Figure~\ref{fig: R-Corbit pnacf} in the Appendix.
 %\begin{figure}[H]
 %      \centering
 %       \begin{subfigure}{0.65\textwidth}
 %           \includegraphics[width=\textwidth]{images/R-Corbit_nacf_covid.pdf}
 %           \caption{COVID-19 series NACF R-Corbit plot.}
 %           \label{fig: R-Corbit nacf}
 %       \end{subfigure}
 %       \hfill
%      \begin{subfigure}{0.65\textwidth}
%          \includegraphics[width=\textwidth]{images/R-Corbit_pnacf_covid.pdf}
%            \caption{COVID-19 series PNACF R-Corbit plot.}
%            \label{fig: R-Corbit pnacf}
%        \end{subfigure}
%        \caption{Figure \ref{fig: R-Corbit nacf} shows the R-Corbit plot.}
%        \label{fig: R-Corbits}
%    \end{figure}
\begin{figure}[H]
\centering
   \hspace{-3.5cm} \includegraphics[scale=0.43]{images/R-Corbit_nacf.pdf}
        \caption{COVID--19 series NACF R-Corbit plot. The plot compares the NACF values between the two COVID-19 outbreaks and the gap between them.
        The maximum lag is ten and the maximum $\rstage$ depth is equal to six.  See text for description.}
    \label{fig: R-Corbits}
\end{figure}
R--Corbit plots allow us to compare the NACF and PNACF values for different time-slices and/or covariates, we read the plot in the same manner as the Corbit plot and look at the legend on the right--hand side for distinguishing between covariates and/or time--slices. The point at the centre is the mean value of the NACF or PNACF values arising from the time--slices and/or covariate data splits. Essentially, if $c \in \{1, \dots, C\}$, where $C \in \mathbb{N}$ is the number of  levels in the  covariate or time-slices, then the value at the centre is
$\nacf(h, r) = C^{-1} \sum_{c = 1}^{C} \nacf_c (h, r),$ where $\nacf_c(h, r)$ is the NACF value corresponding to covariate--level/time--slice $c$.  
\par
The R--Corbit plot in Figure \ref{fig: R-Corbits} suggests that the correlation structure is different between the three time-slices, largest for the second outbreak and smallest for the gap.
 The gap period of July 2020 to September 2020
is characterised by constantly zero-valued observations for many Trusts (but not all of them), which, in the full--data
models above, contribute to the overall variance, but not to the non--zero
lag/stage correlations and related parameters and perhaps making them appear
smaller than would be the case should the two main waves be analysed separately.


 We  focus on model comparisons for the second outbreak in the table below. Prediction error is calculated using the last five realisations as test observations for the model fitted using the rest of the data.
\begin{table}[H]
    \centering
    \begin{tabular}{ccccc}
    & \multicolumn{2}{c}{Mean-squared Prediction Error} &&\\
     Model & One-Step & Two-Step & No. Parameters & Sig. Parameters \\ 
    \hline
    GNAR(1, [1]) & 5.43 & 12.74  & 2 & 2 \\ 
    GNAR(1, [6]) & 5.46 & 12.73 & 7 & 7 \\
    GNAR(6, [$6^{(6)}$]) & 5.87 & 12.99 & 42 & 20 
    \end{tabular}
    \caption{Comparison of different GNAR model orders for the series corresponding to the number of mechanical ventilation beds needed during the second COVID-19 wave for hospital trusts in the network shown by Figure \ref{fig: 1 stage neighbors}.}
    \label{tab: second outbreak comparison}
\end{table}
Table \ref{tab: second outbreak comparison} reflects the information that the Corbit plots in Figure \ref{fig: corbits} provide by noting that models with more than one lag should be discarded based on prediction error and the number of statistically significant parameters at the $1\%$ level. Moreover, it shows the slight performance improvements increasing $r^*$ from one to six has. Therefore, we, as previously, select GNAR(1, [1]), and exploit the ideas in Section \ref{sec: graph model connection} by introducing the following measures of relevance. 

\subsection{Node Relevance}
\label{sec: node influence}
Once we fix an estimate of $r^*$, we can find the sparse weight matrix $\bmat{W} \odot \bmat{S}_{(r^*)}$, where $\bmat{S}_{(r^*)} = \sum_{r = 1}^{r^*} \bmat{S}_r$, from which it is possible to compare the relevance each node has. Motivated by (\ref{eq: nacf inequality}) we define the global relevance index as 
\begin{equation} \label{eq: global index}
    \globindex(X_{i, t}) := \bigg (\sum_{j = 1}^{d} [\bmat{W} \odot \bmat{S}_{(r^*)}]_{j i} \bigg )
   \bigg \{ \underset{l \in \mathcal{K}}{\max} \bigg ( \sum_{j = 1}^{d} [\bmat{W} \odot \bmat{S}_{(r^*)}]_{j l} \bigg ) \bigg \}^{-1},
\end{equation}
which computes the ratio between the largest column sum for active nodes and a particular node. We interpret this as the relevance each node has globally on the correlation structure. 
\par
%\newpage
Next, we define a measure of the strength of pairwise interactions as the weighted contribution from all estimated $\hat{\beta}_{kr}$ across active stages (i.e., $1 \leq r \leq r^*)$, formally,
\begin{equation} \label{eq: local influence}
    \local_{i} (j) := \bigg ( w_{ij} \sum_{k = 1}^{p} |\hat{\beta}_{kr}| \bigg ) \bigg \{ \sum_{l \in \mathcal{N} (i)} \sum_{r = 1}^{r^*} \sum_{k = 1}^{p}  w_{il} |\hat{\beta}_{kr}| \bigg) \bigg \}^{-1},
\end{equation}
%Alternative, 
%\begin{equation*}
%    \local (i, j) := \bigg ( w_{ij} \sum_{k = 1}^{p} |\hat{\beta}_{kr}| \bigg ) \bigg \{ \underset{{l \in \mathcal{N} (i) }}{\max} \bigg ( w_{il} \sum_{r = 1}^{r^*} \sum_{k = 1}^{p}  |\hat{\beta}_{kr}| \bigg) \bigg \}^{-1},
%\end{equation*}
where $\mathcal{N} (i) = \cup_{r = 1}^{r^*} \mathcal{N}_r (i)$, and $i$ and $j$ can only be \mbox{$r$-stage} neighbours for one $r$. Note that this results in $w_{ij}$ if we restrict the above to a specific $\rstage$. This index is larger as $w_{ij}$ gets closer to one. Hence, the index computes the percentage of the neighbourhood coefficient across all lags corresponding to the specific pair (i.e., how strong, relative to other active nodes, $j$ is for forecasting $i$). Thence, $\local_{i} (j)$ is not a symmetric function of $i$ and $j$.
\par
Also, by Theorem \ref{th: pop_sparsity gnar}, we can plot the correlation structure of our selected model based on the distances in the underlying network, which we plot by distinguishing between active nodes: $\delta(i, j) \leq r^*$, colliding nodes: $r^* < \delta(i,j) \leq 2 r^*$ and conditionally uncorrelated nodes: $\delta(i,j ) \geq 2r^* + 1$. Formally, we can also compute the relative strength of conditionally correlated nodes as
\begin{equation} \label{eq: conditional corr index}
    \rscc (i, j) := \mathbb{I} \{ \delta(i, j) \leq r^* \} \{ \delta(i,j) \}^{-1}+ \mathbb{I} \{ r^* < \delta(i, j) \leq 2 r^* \} \{2 \delta(i,j) \}^{-1},
\end{equation}
where the distances are computed on the underlying network, and we divide by two the distance for nodes that are conditionally correlated despite not being active in their respective neighbourhood regressions.
%Mention the change to the RMarkdown document
\par
\cbstart
These measures of global relevance, $\globindex(X_{i, t})$ and local relevance, $\local_{i} (j)$, are shown in detail in Appendix \ref{app: node relevance}. Below, we plot a coloured version of the NHS network, where colours are assigned using the values 
 of $\globindex(X_{i, t})$ as
given by~\eqref{eq: global index} for $r^{\ast} = 1$. 
We remark that these diagnostic statistics depend on the time series values through the choice of $r^*$ and maximum lag, as well as on network connectivity and edge weights. 
\cbend

%\begin{subfigure}{0.48\textwidth}
%    \centering
%    \includegraphics[width=\textwidth]{images/cross_corr_map.jpeg}
%    \caption{}
%    \label{fig: cross-corr map}
%\end{subfigure}
%\hfill
%\begin{subfigure}{0.48\textwidth}
%    \centering
%    \includegraphics[width=\textwidth]{images/local_relevance_r=1.jpeg}
%%    \caption{}
%   \label{fig: local relevance}
%\end{subfigure}
%\hfill
\begin{figure}[H]
    \centering
    \includegraphics[width=0.95\linewidth]{images/node_relevance_network_plot.pdf}
    \caption{NHS Trusts network coloured by global node relevance values, i.e., $\globindex(X_{i, t})$ given by~\eqref{eq: global index} for $r^* = 1$. Nodes that have a larger effect on network autocorrelation are brighter ({\color{yellow} yellow}) and less relevant ones darker ({\color{violet} violet}). See the RMarkdown document mentioned in the Conclusion (link provided in the Appendix) for $r^* = 3$ and $r^* = 6$, and Figure~\ref{fig: node relevance network A4} in
    Appendix \ref{app: node relevance} for a full page size version of this plot.
%{\cbsta
}
    \label{fig: global relevance}
\end{figure}

% \begin{figure}
% \centering
% \begin{subfigure}{0.65\textwidth}
%     \includegraphics[width=\textwidth]{images/node_relevance_large_legend.pdf}
%     \caption{}
%     \label{fig: global relevance}
% \end{subfigure}
% \hfill
% \begin{subfigure}{0.65\textwidth}
%     \includegraphics[width=\textwidth]{images/local_node_relevance_large_legend.pdf}
%     \caption{}
%     \label{fig: active nodes}
% \end{subfigure}
% \caption{ Figure~\ref{fig: global relevance} shows the global node relevance, $\globindex(X_{i, t})$ in~\eqref{eq: global index}, which distinguishes nodes that have a large effect on network autocorrelation 
% %{ 
% (a network coloured version of node relevance can be found in the R Markdown document mentioned in the conclusion).
% %} 
% Figure~\ref{fig: active nodes} shows the function $\local_{i} (j)$ from \eqref{eq: local influence}, which shows active neighbourhood regressions by highlighting the strength of pairwise neighbourhood regression coefficients. In other words, $\local_{i} (j)$ quantifies the portion of the neighbourhood effect (connection weight, $w_{ij}$) weighted according to its network parameters $\beta$. E.g., the largest (yellow) coefficient on the right-hand side of
% the plot corresponds to transmission from the University Hospitals Plymouth to the Royal
% Cornwall Hospitals NHS Trusts.
% }
% \label{fig: correlation structure}
% \end{figure}
% Figure~\ref{fig: cross-corr map}  and \ref{fig: local relevance} compare the observed correlation with the
% conditional correlation.
% Both show clustering across the network in a self-similar manner, which is reflected by the strong autoregressive behaviour and dependence on one-stage neighbours.
Figure~\ref{fig: global relevance} shows the
$\globindex(X_{i, t})$ values for all NHS Trusts.
Individual Trust
identities are hard to discern from this
plot, but it is instructive
to examine them in a bit more detail.
The top ten
Trusts for $\globindex(X_{i, t})$ in decreasing order
are 1.\ Buckinghamshire,
2.\ East Cheshire,
3.\ Stockport,
4.\ Royal Berkshire,
5.\ Harrogate District, 
6.\ Tameside Glossop, 
7.\ Great Western (Swindon),
8.\ University Hospitals of North Midlands, 
9.\ Oxford Health and 
10.\ Oxford University Hospitals
with $\globindex$ relevance ranging 
from $1.000$ to $0.952$.
Figure~\ref{fig: conduitTowns} shows the
most relevant $\globindex(X_{i, t})$ Trusts for both the
first ten and first sixty Trusts in decreasing order.
It is interesting that the most relevant Trusts
relating to the network time series
dynamics 
are located in two clear clusters: one positioned to the
north-west of London (between London, the West Midlands [Birmingham],
Bristol and the Southampton/Portsmouth urban centres)
and the other between and around the urban centres of the
West Midlands
[Birmingham], Manchester, Sheffield, Nottingham, West
Yorkshire and Liverpool. It is noticeable
that in the top ten $\globindex(X_{i, t})$ relevance
list there are {\em no} Trusts in urban centres and {\em only} eighteen in the top sixty list. From the point of view
of network time series dynamics, the most relevant trusts
are not urban centres {\em per se}, but intermediately-located Trusts. Having said this, it is intriguing
that there do not seem to be relevant Trusts
immediately to
the south and east of Birmingham, nor between
Yorkshire and Tyneside, nor in the densely populated
area south and east of London.
The location of these relevant Trusts might have
implications
for future epidemic mitigation
that could be taken to minimise
viral spread.

For example, rather than national or regional lockdown the aim would be
to impede movement between areas as a firebreak. Such a firebreak is
a larger, regional, version of, e.g., firebreak culls that existed as a control measure during the 
2001 foot--and--mouth epidemic \citep{fmdnature} and part of
the current foot--and--mouth control strategy for Great Britain,
see \cite{defra}.

Knowledge of such
relevant locations might enable a more targeted and efficient set of
restrictions to be implemented, rather than a blanket approach, which
imposes costs on unafflicted areas. Obviously, such mitigations require
further investigation and research.


The least relevant $\globindex(X_{i, t})$ Trusts are those in the extremities of the network, typically coastal towns, and are shown as black squares
in Figure~\ref{fig: conduitTowns} (the exception
being the Wye Valley NHS Trust, which is an extremity
of the English NHS Trust network on the border with Wales).
That the least relevant Trusts are located at network extremities
is perhaps not surprising as there are no susceptibles located in the sea and
the epidemic can only be reflected back to the region(s) that initially
`infected' the least
relevant regions. 
%local node influence brief description

\cbstart    
Figure~\ref{fig: REF local relevance} in Appendix \ref{app: node relevance} shows neighbourhood regressions sparseness by highlighting the strength of pairwise neighbourhood regression coefficients measured by $\local_{i} (j)$ in \eqref{eq: local influence} for NHS Trust REF (located in Cornwall). Looking at similar plots of local relevance for different nodes suggests that effects of one--stage neighbours are evenly spread, which highlights parsimonious properties of GNAR models. Further, it appears that Trusts connected to many one-stage neighbours behave in a group--like manner, i.e., there is no single node that appears to be more relevant. Whereas, if a Trust is in a more isolated community, e.g., nodes corresponding to Trusts in Cornwall, then the Trust is affected more by each of its one-stage neighbours. The largest ({\color{yellow}yellow}) coefficient on the right-hand side of the plot corresponds to transmission from the University Hospitals Plymouth to the Royal Cornwall Hospitals NHS Trusts, which can be thought of as the main conduit town into Cornwall. Thus, the plot suggests, unsurprisingly, that communities at the network's edge have different characteristics than 
 the ones closer to the main population centres. %Figure~\ref{fig: REF local relevance} in Appendix \ref{app: node relevance} further points out that GNAR can manage this differing community dynamics by adjusting the apportionment of the neighbourhood regression coefficient $\beta$ accordingly.
\cbend

\begin{figure}
\hspace{-2cm}
\begin{subfigure}{0.76\textwidth}
    \centering
    \includegraphics[width=\textwidth]{images/conduitTowns10}
    \caption{}
    \label{fig: conduitTowns10}
\end{subfigure}
\hspace{-5cm}
\begin{subfigure}{0.76\textwidth}
    \centering
    \includegraphics[width=\textwidth]{images/conduitTowns60.pdf}
    \caption{}
    \label{fig: conduitTowns60}
\end{subfigure}
\caption{Blue dots show 
top eleven urban centres in the UK.
Blue circles indicate urban centre population size with
London as approximately 9.78 million people.
Black squares are ten least `relevant' $\globindex(X_{i, t})$ Trusts.
Red dots: highest relevant NHS Trusts according to
$\globindex(X_{i, t})$. (a) top ten $\globindex(X_{i, t})$ Trust
locations (b) top sixty $\globindex(X_{i, t})$ Trust locations. \label{fig: conduitTowns}}
\end{figure}

\subsection{Discussion} \label{sub: application discussion}
An attractive property of GNAR parsinomity is that it enables us to interpret a large number of interactions with a relatively small number of parameters, in our case we can interpret the interactions between 140 NHS trusts with just two parameters
and {\em still} have superior forecasting performance, as shown
by Table~\ref{tab: gnar comparison}. The Corbit plots in Figure \ref{fig: corbits} show the strong autoregressive structure of the series $\boldsymbol{X}_t$, highlighting that correlation cuts-off after lag one. 
\par

Our chosen model is a global--$\alpha \enspace \mathrm{GNAR}(1, [1])$, which mandates that the expected value for each node is a linear combination of its value at time
lag one (with coefficient $\hat{\alpha}_1 \approx = 0.95$ from Section~\ref{sec: comparison}) and its $1$--stage neighbours at time lag one
(with coefficient $\hat{\beta}_{1, 1} \approx 0.043$). The transfer from the same trust at the previous
time step at $0.95$ is considerably larger than the $0.043$ neighbor contribution, but there are typically many neighbors, making the neighbour contribution non--neglible.

Further, the expected number of beds at trust $i$ at time $t$ increases,
if either the number of beds at trust $i$ increases at time $t-1$ and/or the number of transfers at its $1$-stage neighbours increases at $t-1$. Hence,
the outbreak acts in a `continuous' matter, i.e.,
bed demand spreads between $1$-stage neighbours. Crucially,
we concluded that demand does not spread via 
from higher-stage neighbours and we inferred this as a result of {\em statistical model
selection}, not by an {\em a priori} assumption, as is implicitly, and perhaps unthinkingly
in some simpler network time series
and conditional autoregressive models.
\par
Relating to our discussion on mitigations earlier,
if local communities can block the spread to $1$--stage neighbours of their $1$--stage neighbours (i.e., their $2$--stage neighbours), then, on average, this would also create a firewall between them and other communities in the network.
This obvious--sounding statement is encapsulated by Theorem 2, which reflects this by noting that two nodes are uncorrelated at all lags given their $2$--stage neighbours, i.e., the outbreak cannot jump between nodes separated by at least two edges. It must travel through each nodes $1$--stage neighbours. 
\par
%  Table \ref{tab: prediction comparison} shows the excellent predictive performance of our chosen GNAR(1, [1]) model, which suggests that the outbreak is mostly dependent on the current node's past and gets worse by absorbing cases from one-stage neighbours, and, by Theorem \ref{th: pop_sparsity gnar}, nodes are conditionally uncorrelated from all other nodes given their two-stage neighbours --- i.e., communities can protect themselves by isolation.
% \par
Further, and interestingly, the R--Corbit plot in Figure \ref{fig: R-Corbits} shows that there is an increase in correlation during the second outbreak, which we attribute to less stringent measures and an increase in the number of interactions between nodes (i.e., more people travelling across different areas of England); see \cite{owidcoronavirus}. These results reflect the node relevance in Figure \ref{fig: global relevance}, which shows that the nodes with the most relevance are close to London, and that the less relevant nodes are at the extremities,  e.g., NHS Trust  REF, i.e., Royal Cornwall Hospitals NHS Trusts. Nevertheless, all nodes are similarly important given that we constrain the model to 1-stage neighbourhood regression, which shows the global properties of how demand for mechanical ventilation beds due to COVID--19 cases spread throughout England during the pandemic.