\documentclass[12pt]{article}

\usepackage[dvipsnames]{xcolor}

\usepackage{amsmath, amssymb}

\usepackage{hyperref}

\usepackage{setspace}\doublespacing

\usepackage{graphicx}

\usepackage{float}

\usepackage{comment}

\title{Response to reviewers on `New tools for network time
series \ldots' ID JRSSA-Dec-2023-0316}

\author{Nason, Salnikov and Cortina Borja}

\date{30th April 20224}


\begin{document}


\begin{singlespace}

\maketitle

\noindent
{\color{red} Guy comment in RED}\\
{\color{ForestGreen} Daniel comment in GREEN}\\
{\color{blue} Mario comment in BLUE}\\
{\color{brown} Marker for COMPLETE OR INCOMPLETE}

{\Large \sl Please use [] when making comments to each other. Don't use them when making comments to the reviewers/AE, etc}

\end{singlespace}

We would like to thank the editorial team and the two reviewers for
their hard work, comprehensive and positive reviews and suggestions
for improvement.

\section{Response to Editor}

{\color{red} [We might put something in here to respond to Adam]}

We noticed that we had used $d$ for the number of nodes/time series and $d(i, j)$ as the
number of hops (distance) between two nodes. Since this is confusing, we replaced the
latter $d(i, j)$ by $\delta(i, j)$.

\section{Response to AE}

{\em Associate Editor comments:

Following the advice of the reviewers, my recommendation is for major revision. It is a very well-written and important paper introducing very relevant new methodologies for network time series analysis, although, it could be more reader-friendly, despite their attempt to move the proofs to supplemental material. However, I share the same concern as Reviewer 2 that the used example reflects the data structure of spatio-temporal data more than network time series data, and unfortunately, the authors did not provide the R-package GNAR and scripts containing code to perform the diagnostic methods described in the article as a GNU zipped tar file (which supposedly contain the datasets used in the analysis), or the COVID-19 Series data set analysed in Section 5 as a .txt file. The authors refer to these in page 27. 
It is important to ensure that, if this paper ends being discussed, that the methodologies being introduced are the star of the show, with the alternative being a potentially heated discussion around the example used to demonstrate the methods. However, it would be recommendable to allow the authors to reply to the reviewers’ feedback before making a final decision on it being accepted for Discussion.}

{\color{red} We would like to thank the Associate Editor very much for their consideration of our
work and, in particular, suggesting that our work was `well-written' and `important'. We have tried hard
in the revision to be more `reader-friendly' and the review team will notice extensive revisions of varying
lengths because of this.

On whether the COVID NHS Trust data is really `classical' spatio-temporal or network spatio-temporal,
this is discussed by Reviewer~2's point~5 and
we believe we provide several compelling reasons in response to that
point why we believe the network
version is more appropriate (although, granted, such data can be shoehorned into all kinds of models).

We apologise for not supplying the GNAR package and associated scripts directly last 
time. However, for the last version our
code for producing the plots and the related functions was already in the GNAR package that was
deposited on the CRAN archive, but not the scripts to produce the figures.
We apologise for not making this clearer. To remedy the lack of
scripts we have, this time, produced an extensive RMarkdown document (linked below and in the paper) that
contains code to produce plots and an extensive set of extra and alternative plots, tables and analyses. The time series data and network are (and were) available in the 
freely available GNAR package (at CRAN) as the objects {\tt logMVbedMVC.vts} and
{\tt NHSTrustMVCAug120.net}.

Reviewer~1 suggests that the paper has a `split' nature, but Reviewer~2 seems to think
we have the balance right (see reply to Reviewer~1 below). However, we have worked hard
to rewrite sections, improve the flow and cross-connections (and, frankly, a few
mistakes, which could not have helped). We are happy to look at the balance again
should the review team and/or Discussion Meeting Committee wish.}



 
\section{Response to Reviewer~1}

{\em General

On balance I would be quite strongly supportive of some version of this paper being
read at a Discussion Meeting.

Why "on balance," you might ask?

On the positive side, the methodological area is important and the subject of quite a bit of current work, and the application area (covid) ought to be continue to be highlighted because the epidemic is not over, and even if it doesn't recrudesce the wealth of data now available on the epidemic has much to tell us about how to, and how not to, run public health monitoring and surveillance systems.

On the negative side, the paper has something of a split personality. Apart from
introducing the corbit plot (on which, see below for more comment) the first 4 sections
are highly technical and, as a non-specialist, I struggled at times to distinguish between
what was new and what was a review of existing results. Section 5 has obvious scope for discussion, especially if the authors could make the data-set directly available
with the paper in advance of the meeting (again, see below for more comment). A harsh judgement of the application as currently presented is that it is primarily a vehicle
for showing that the authors' methodology out-performs several rival methods in a very
specific predictive setting.}

{\color{brown} INCOMPLETE. We need to
perhaps think if there is anything more to do
here.}

{\color{red} We would like to thank the Reviewer for their work and we are heartened that they `{\em would be strongly
supportive of some version of this paper being read at a Discussion Meeting}' However, we also understand the Reviewer's
comments and we have worked hard to address them. We can see where the Reviewer is
coming from in referring to what they say as the paper's `split' nature and respect
that view. We would
politely disagree, as we see the `theoretical' part
(Section~4) as a key aspect of the paper, supporting
 other parts of the paper, not least the proper examination and interpretation
 of the COVID data.
We have tried to improve the flow and cross-connection in the paper, so perhaps this
will help? Having said this, the main theoretical meat is in Section~4, whereas
we believe Sections~1 to~3 would be accessible and familiar (at some level) to those
who have seen regular time series and/or linear modelling --- although, of course
the methods and plots are completely new, different to their time series counterparts,
but inspired by them.

In particular, we have tried to signpost what is new and what was already established
by better referencing and signposting. The basic GNAR model is not new, but our
new formulation (using Hadamard products, and matrix form), the NACF, PNACF, Corbit
and R-Corbit plots are all new and very useful (in our opinion).

We are very glad that you see that {\em Section~5 has obvious scope for discussion} and
the data set is available in the GNAR package from CRAN, as mentioned above and
(now mentioned) in the paper. Our work is bolstered by the very good predictive
performance, but that should be a compelling feature to encourage people to
use these methods. Parameter inference is fine, and we do it and discuss them, but
you never really know whether (a) the parameters are the right ones and (b) whether
they are well-estimated, whereas predictive performance is absolute and can
be tested on the actual data. Having said this, we provide more discussion on the 
parameters we estimate and their interpretation in Section~5.5.
}


{\color{ForestGreen} [Indeed GNAR outperforms rival methods in this specific predictive setting. However, one of the main goals of the analysis is studying the correlation structure between NHS Trusts during the COVID-19 pandemic, this can be done from a spatio-temporal point of view by modelling the data as a spatio-temporal conditional autoregressive model. However, this method implicitly assumes that only one-stage neighbours affect a node and does not provide a tool for assessing autocorrelation in the whole network. Another possibility is choosing a known covariance function, e.g., a Matern kernel, and modelling the data as a random field, this approach does not provide a concise summary of autocorrelation in the network, is an approximate solution to a stochastic differential equation, where further approximations are required for computing estimates, and assumes that spatial effects spread continuously in space, which might not be sensible given that the data are related to hospitals and not viral spread in space-time; see response to Reviewer 2 point 5. Also, both of these methods are computationally cumbersome, constrained to one-lag temporal effects (i.e., the autoregressive order is assumed to be equal to one), and are difficult to interpret since most effects are modelled as random. 

VAR models are well-known and provide a quick assessment tool for modelling multivariate time series. However, model order selection is challenging, the number of unknown parameters increases rapidly (exponentially) as the number of variables (time series) increases, may not capture spatial or some other effects and are difficult to estimate. This makes estimating and comparing many different VAR models challenging. 

GNAR models enable us to study the autocorrelation in the network rather than between pair-wise interactions, can be estimated quicker and compared faster than the models above. Also, they provide a tool for assessing whether or not the spread of the virus depends mostly on a Trust's one-stage neighbours and its condition the previous day rather than implicitly assuming it. Put another way, the PNACF allows us to select lagged and/or spatial effects based on the observed data rather than by assumption. At the present moment GNAR models are linear, however, a research avenue we are excited about is extending the model to nonlinear settings, incorporating penalty terms, mixed type data and extensions to multivariate nodal processes. These extensions are likely to be conceptually and computationally harder for models based on strong distributional assumptions that already require cumbersome computation in practice, e.g., CAR. We hope that the analysis in Section 5 motivates the use and development of parsimonious models by highlighting some attractive GNAR properties. 

Finally, we have compared GNAR to VAR methods in other settings and observed similar results. Our main intention is showcasing the increased interpretability and flexibility a parsimonious network based model provides. Interestingly, GNAR outperforms VAR models in this important setting.]}

{\em Details}

{\em  1.  What does the circular presentation of the corbit plot add? The circularity seems artificial, with no obvious advantage (and some loss of resolution) by comparison with superimposed conventional line-graphs of correlation vs lag, one per stage. Similarly, I thought the ``Wagner'' reference was somewhat contrived to the point of being irritating.}

{\color{brown} COMPLETE}


{\color{red} We agree that the Wagner name was
    contrived. We also discovered it was probably
    not a good name for other reasons. Hence,
    we have changed it to R-Corbit, which stands
    for `ring Corbit'.

    We believe that the circular presentation adds
    a great deal. The inner ring corresponds to 
    nearest neighbours and although they are at
    different time lags, they are all at the same
    distance from the centre. You don't get the
    same kind of effect from a conventional
    line-graph plot=`rectangular plot'. E.g.\ across the lags, where
    is `spatial zero'? 

    Further, we believe that for larger networks,
    where there are many \mbox{$r$-stages}, the circular plot
    might be easier to read than a rectangular plot,
    whose lines might overlap.

    Similarly, a rectangular line-graph plot of
    the equivalent of a R-Corbit (Wagner) plot could appear to
    be very complex visualiation. The equivalent would either need
    to plot an extra line for each covariate (or
    time period) or plot extra symbols at each point
    of a Corbit-type plot. The R-Corbit plot spreads out
    these points and they are much easier to see.
    
     We also think that the visual appeal of the
     Corbit plot
    should not be underestimated. There is an 
    aesthetic component that might attract analysts
    and researchers. Further, the Corbit
    and R-Corbit plots are distinctive and one is
    constantly grounded in the network time series
    world, which separates the modelling and analysis
    from both time and spatial data analysis.
    
    However, we like your suggestion and would like
    to give users a choice. So, we have implemented
    the `rectangular' line-graph plot and made it an
    option for the user to select, if they wish.
    You can see some examples of this kind of plot in
    the our new associated RMarkdown document in Section 2.3 there:
    
    
    \href{https://dansal182.github.io/new_tools_for_network_time_series_with_application_to_COVID_19_hospitalisations/network_autocorrelation_notebook.html}{R Markdown link.}

    In addition, we have added a small blank gap between
    the largest lag spoke and that at lag zero, as
    these two lags are not adjacent in time.
    }



{\em 2. The application has great potential for discussion but the link to coronavirus.data.gov.uk is not very helpful. After some digging, I found what I take to be the data-source under
\begin{verbatim}
https://api.coronavirus.data.gov.uk/v2/data?
areaType=overview&metric=hospitalCases&format=csv
\end{verbatim}
but I could not see how to download the data without doing this one NHS Trust at a time. Since this is publicly available data, the authors could do potential discussants (and, indirectly, themselves) a favour by giving a single link to a 452 by 140 data-array. I would make this a requirement for acceptance.}

{\color{brown} COMPLETE}

{\color{red} The full data set (multivariate time series and network)
    is contained with our GNAR package. The package was mentioned in the last submission in the Supplementary
    Material, but we did not include the Zip file of the package. 
    However, the data and associated functions that contain the functionality explained in our article
    is freely available from the CRAN R software archive as the GNAR package and this is probably
    the easiest way to get hold of the data.
    The multivariate time series is stored conveniently as a matrix.
    The NHS Trust network is stored in the GNAR network format.
    The GNAR package has functions that convert
    to and from our GNAR network format and the popular {\tt igraph} network format.
    
    The incantation we
    used to obtain the data from the Coronavirus website in June 2021 was:
    \begin{verbatim}
    https://api.coronavirus.data.gov.uk/v2/data?
    areaType=nhsTrust&metric=covidOccupiedMVBeds&format=json    
    \end{verbatim}
    However, it looks like the website {\tt api.coronavirus.data.gov.uk} has moved to
    {\tt api.ukhsa-dashboard.data.gov.uk}. From a search of their metrics
    it looks like {\tt covidOccupiedMVBeds} is no longer stored. Possibly it is
    no longer being collected and perhaps this is
    because few patients are currently mechanically ventilated
    because of COVID?
    
    We appreciate that replicating the results is an important part of our work. Thus, we have created a Markdown document that demonstrates how to load the NHS COVID-19 into an $\texttt{R}$ session and reproduce our results using the $\texttt{GNAR}$ package. The Markdown document is available at:

       \href{https://dansal182.github.io/new_tools_for_network_time_series_with_application_to_COVID_19_hospitalisations/network_autocorrelation_notebook.html}{R Markdown link.}

}

{\em 3. The authors describe the data as giving numbers "transferred to mechanical ventilation beds for day t" which sounds like incidence, whereas the data I found on the web seemed to be referring to numbers on ventilators on day t, i.e. prevalence. Maybe I've got the wrong data-set, which rather reinforces my previous point.}

{\color{brown} COMPLETE}

{\color{red} Thank you this. It is definitely prevalence. However, you are right, the wording was either
ambiguous or wrong in one or two places. We have been through and checked all references to the data are
about occupancy, or prevalence.}

{\em 4. Several statements seem to me a bit opaque and/or unjustified. For example:}

{\em a) ``each trust is connected to trusts within a 120km radius of itself.'' Why this choice
of metric and value?}

{\color{brown} COMPLETE}

{\color{red} We have explained in more detail on how we obtained our `120km' 
network (at the beginning
of Section 5), why this value was chosen and its relationship with other possible choices. In essence, it is a network that results
in a very good fit,
but, of course, there are many
other possibilities. The choice of such networks itself can be viewed as a (separate)
statistical modelling problem, consideration of which would require further investigation
and papers. We have discussed a tessellation
alternative in response to Reviewer~2's
comment~7, below. However, we have
found our $k$}


{\em b) ``Figure 6a shows strong correlations between series at different nodes'' ...but the maximum correlation seems to be about 0.2?}

{\color{brown} COMPLETE}

{\color{red} The reviewer is most perceptive and has observed the phenomenon
where correlations in our kinds of network time series tend to be smaller than one
might experience in, e.g., regular time series analysis. The highest is a bit
more than 0.2. However,
such smaller correlations are still influential and {\em detectable}, partly
because they are the result of smoothing over a
large number of pairs (much larger than is usually formed in regular univariate time series).
The detectability of small, but influential, correlations improves for
larger sample sizes and, indeed, some of them would not be detected with reasonably large
probability, if the sample
size were not large enough. (People perhaps
underestimate this issue in regular time
series too). Subtle effects can only be detected with larger sample sizes.
For example, in white noise testing, to achieve a statistical
power of $80\%$, a sample size of only 32 is required for an AR$(1)$ process
for AR parameters $\alpha=\pm0.9$, but a size of 256 is 
required when $\alpha=0.3$ [see, e.g.\ Savchev, D. and Nason, G.P.
(2014) White noise testing using wavelets, {\em Stat}, {\bf 3}, 351--362 and references
therein]).

It's also the case that our correlations tend to be smaller as they are uniformly
reduced by the $\lambda$ introduced in Definition~3 to ensure that the
overall correlations lie in $[-1, 1]$, which is not a feature of regular autocorrelation.}



{\em c) ``We interpret the gap period (July 2020 to September 2020) as an absorption state,  which is characterised by constantly zero-valued observations'' So what are you modelling?}

{\color{brown} COMPLETE}

{\color{red} Thanks for this. This sentence
was poorly worded.
We have rewritten it in the paper and
reproduced it
immediately below for convenience. The words `absorbtion state' was unfortunate in this description as we don't mean an absorbtion state in a stochastic process sense.
We're merely saying that the pandemic was two waves with a quiet state in between
and then studies modelling on the second wave in a bit more detail. We also
remark that the gap will affect the `full-data' modelling, but effectively inflate
the variance (``lag/stage'' zero correlations). We now write:
\begin{quote}
``The gap period of July 2020 to September 2020
is characterised by constantly zero-valued observations, 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
small than would be the case should the two main waves be analysed separately.''
\end{quote}}

{\em d) ``The location of these influential Trusts might have implications for future epidemic mitigations that could be taken to minimize viral spread.'' what might these implications be?}

{\color{brown} COMPLETE}

{\color{red} Thank you for this comment. We have added further explanation of
what we mean by this in the last eight lines of the penultimate paragraph of Section~5.4.
We describe a tentative nuanced local strategy compared to the
blanket movement restrictions enacted during the 2001 foot-and-mouth epidemic.
We are not saying this is perfect, or will definitely work, but the concept is
suggestive of what might be done and/or further researched.}

{\em e) ``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 2, nodes are conditionally uncorrelated from all other nodes given their two-stage neighbours — i.e. communities can protect themselves by isolation.'' I'm having trouble parsing this text}

{\color{brown} COMPLETE}

{\color{red} Thank you. Yes, this sentence was compressed and
garbled. We have replaced it with two, hopefully clearer, paragraphs that expand on what
we meant (new second and third paragraphs of Section~5.5). We interpret
the GNAR model, explain how the statistical modelling process enables us to conclude
scientifically about the spread of bed occupancy (rather than just assuming) and how this
information can be used to suggest a local mitigation strategy.}


{\em f) ``The least influential globindex Trusts are those at the network extremities'' Is this coincidental, or at least in part a consequence of the geometry - all other things being equal, Trusts near the coast will automatically be less connected than Trusts away from the coast?}

{\color{brown} COMPLETE}

{\color{red} As you have alluded to, our comment about the least influential is
 not surprising given their location at network extremes. We have added
 a sentence of two to reinforce that statement. The more interesting point is,
 perhaps, that we initially thought
 that the most influential Trusts would be those in densely-populated urban
 areas, but that turned out not to be case. We were surprised when it turned out
 to be `conduit towns' and the text does reflect these two aspects. A full readable
 list of all the Trust's and their influence is provided in our new RMarkdown document,
 which is part of our supplementary materials.

 Also, the influence metric is not
 solely dependent on the design/network
 topology/position, but also on the
 time series values, through the value
 of $r^\ast$. This is important as
 the influence list could have been different,
 if $r^\ast$ had of been different, which
 it could have been for different time
 series values. We say more on this
 in the response to Reviewer~2 point~10 below.
 }



\section{Response to Reviewer~2}

{\em This is a contribution to the network autoregression literature, which has attracted quite some attention in recent years. It addresses an interesting problem, provides very practical modelling and diagnostic methodology and gives a topical and substantial application. The paper is well written, with focus on the important parts of the methods and with technical details sensibly placed in supplementary materials. Most of my comments are minor. They are all intended to be helpful.}

{\color{brown} POSSIBLY INCOMPLETE}

{\color{red} Thank you very much for you comprehensive and supportive response. We have
 addressed your constructive points fully as described below.

[Note that R2 clearly says that the {\em technical details sensibly placed in supplementary
materials}. This is some evidence that perhaps the paper is not so `split personality'. So, perhaps
we can address R1's general comments by better signposting (e.g. if you don't want to read the
technical bits, skip to section XXX below) and perhaps reducing the technicalities in the paper?)]

}


{\color{ForestGreen} [I believe that the technical split is Section 4, however, I agree with Guy that it's an important part of the paper and shows how GNAR extends partial correlation graph models (i.e., from edges to paths). Also, Section 4.2 shows some form of transforming a VAR into a GNAR by exploiting the network and suggests some interesting properties. Perhaps, the best course of action is to mention that Section 4 should be skipped by readers interested in the methodology and application, and not so much in the models construction and properties?]}

{\em 1. I really like the ideas and the basic simplicity of the model and the estimation procedure. But it would be good to bring out more clearly what exactly is the new methodology in the paper. So far as I can tell the network and partial network autocorrelation functions and the various plots are new, but not the underlying model or estimation procedure. For example Equations (1) and (2) together in this paper are the same as Equation (1) of Knight, Leeming, Nason and Nunes (2020). Obviously the paper is intended for a general readership and needs to introduce GNAR for the non-specialists (which it does, successfully). Still, it would be helpful to flag more often what is new and what is extant, without changing the paper to any material extent.}

{\color{brown} COMPLETE}

{\color{red} We are very happy to do this and we describe below how we have
    improved the signposting in the article to make it more clear what is new.
    R2 is correct in that the network and partial network autocorrelation are
    new. The Corbit and R-Corbit plots are also brand new. The estimation
    procedure is not new and is, as previously, least squares. However, we
    have introduced a new way of writing the model down, using Hadamard notation,
    which, at the very least, provides a new way to understand the model and, perhaps
    makes it clearer to some. Since putting out original submission on arXiv, we have
    noticed another paper `borrowing' our Hadamard notation already (also on arXiv).
    The material in Section 4 is all new: a new interpretation of GNAR as
    a generalised graphical model (which are well-known) and, perhaps more
    significantly, a new explanation on how GNAR processes, with known network
    structure, can be viewed through the massively popular lens of variable
    selection and shrinkage. The COVID-19 Section~5 is also completely new.

    Specifically, to improve the signposting as requested,
    we have modified the abstract, four places in the Introduction,
    Section~2.2, Section~3 and~4 and these are indicated with (small)
    changebars. In the `2.3 GNAR Model Estimation' section, we have added two
    sentences that makes clear that the least squares estimation was used previously.
    If the reviewer thinks that any of this is overdone, please do let us know, and
    we can rein back if necessary.} 



{\em 2. More intuition around nacf in Definition 3 might be helpful.}

{\color{brown} COMPLETE}

{\color{red} Good idea, thanks. We have rewritten and re-ordered content in Section~3.1,
    which contains the NACF
    definition in Definition 3. We swapped the mention of the network autocorrelation with
    the regular autocorrelation, so that the reader is grounded in the familiar autocorrelation
    first. We then explain how the NACF is similar but has two main differences, which we explain in
    general terms and then we explain the mathematical purpose of those extra terms
    directly after the definition.}

{\em 3. The Corbit plot is very nice but the actual correlations are quite low so there isn’t a huge contrast between the highest and lowest. Maybe a different example with a more dramatic range of correlations would bring out the merits of the plot more clearly? This figure is just to illustrate the merits of the plot after all, so use data that shows the plot to best effect.}

{\color{brown} COMPLETE}

{\color{red}
    Again, thank for the compliment and your
     perceptive comment. The largest correlation in, e.g., Figure~2
    is $0.14$. Yes, in general the correlations are
    smaller. This is for two reasons. 1.\ In GNAR
    the correlations are formed from weighted values
    for each neighbour-pair. In related time series
    models, e.g.\ VAR, there isn't (generally) a weight.
    The weight, $w_{i, j}$
    is an important feature in GNAR
    as two time series are associated by a
    (sometimes virtual) edge: the
    correlation should be smaller if the time
    series are further apart (in the network structure)
    and, conversely, larger if they are closer together.
    We use inverse distance weights by
    default, so the effect of the
    network parameter (one of the $\beta$s)
    per node is
    fractionally distributed across the neighbours,
    because the inverse distance weights always
    sum to one. This depresses the values.
    2.\ To ensure the correlation is between
    $-1$ and $1$ we introduce the correlation bound
    $\lambda$, although this typically does not 
    reduce the correlations by much.

    We have inserted a remark into
    Section~3.2,
    which introduces the Corbit
    plot, to briefly explain
    the why correlations tend to
    be smaller.
    
    Having said all this, and following your suggestion,
    our new RMarkdown document includes examples with larger NACF values.}



{\em 4. Knight et al (2020) suggest using BIC to inform model selection (with an AIC option in their package). Perhaps we should see BIC (or AIC) alongside Figures 2 and 3 for a range of models for the fiveNet network. Also for the COVID example later.}

{\color{brown} COMPLETE}

{\color{red}  A clear rationale
    for our work is to move away from AIC/BIC as it
    is highly computationally intensive when being
    used to compare large numbers of models, as
    of often occurs during high-dimensional time series.
    Our proposal is to use NACF/PNACF instead, which
    is quicker and also very familiar to those
    who have undertaken even a basic time series course.
    We are also guided strongly by predictive performance.
    
    Having said this, AIC or BIC can, of course, be
    used for model inference in our context and its use,
    over or instead of NACF/PNACF or predictive performance
    can depend on ones' intended aims and or personal
    preference. We have included AIC/BIC values in
    Table~1 of the paper, which compares different
    GNAR models.

    Another potential disadvantage of AIC/BIC in practice
    is that it can be sometimes difficult to compare
    models of different types (e.g.\ ARIMA models with
    different numbers of differences) as there is not
    uniformity over how the constant in the log-likelihood
    component is handled, or even present at all in
    different software implementation.
    }



{\em 5. I am not fully convinced that the COVID data represent a network problem: these are surely spatio-temporal data. If the authors want to apply a network autoregressive model then it would be better to compare results with that of a spatio-temporal analysis. We could easily take a conditional autoregressive model (given information to just before t) combined with a parsimonious spatially-referenced error structure, such as Mat\'{e}rn with known or unknown shape parameter. This would be a fairer comparison than with an unstructured VAR. And having correlation depend on distance avoids the fair arbitrary selection of a connection if trusts are less than 120km apart.}

{\color{brown} COMPLETE}

{\color{red} We have confidence in modelling the COVID data using a network based model, and believe it is more
appropriate than a spatial data-based
model (which could admittedly be used also).
Why? The hospitals can be seen as nodes and the patients do not necessarily all come from a spatial area
surrounding a hospital. Many patients do arise from the area surrounding the hospital, but many come as
transfers from other hospitals. Some patients were already in the hospital and catch COVID. Some of the
staff can catch COVID in the hospital (and outside too, but the chances of catching it for them in the hospital
are quite large due to poor PPE provision, for example). Also, it is surely more accurate to talk about
the number of cases at the hospital on a given day and the hospital is a relatively small geographical area.
There are no humans on mechanical ventilation beds in any geographical area outside of the hospital grounds.
This would be different if the variable was `COVID caught' as that could be anywhere in the geographic area or,
with other diseases, such as foot and mouth, you could consider number of animals infected in a wider
spatial area, particularly as the wind spreads the vector in that instance.


Many spatial models rely on observations that occur over a fixed spatial area. The COVID data won't be like that necessarily.
For example, patients often have a choice over which hospital to attend. So, if your live equidistant between two or
three hospitals, or you might think a certain hospital is better even though it is further away. Hence, there
is not necessarily a coherent mapping between Trust and a spatial area. However, the definition is clear
with a particular Trust.

Having said all this, we decided to follow your suggestion and apply a Conditional Autogressive Model to our
data. We chose the {\tt CARBayesST} package of Lee, Rushworth, Napier and Pettersson (2023), which
is a software implementation of Lee D., Rushworth A. and Napier G. (2018) Spatio-Temporal Areal Unit Modeling in R with Conditional Autoregressive Priors Using the CARBayesST Package.
{\em Journal of Statistical Software}, {\bf 84}, Number 9. This does not use a Matern kernel, but the
network structure is encoded into the model.

We have discovered that, for the simplest GNAR$(1, [1])$ model the two models have similarities and their
forecasting performance is also similar,
see Table~\ref{tab: prediction comparison with CAR}.
(We have updated Table~2 in our paper to include the Naive and
CARar results).
The comparable GNAR$(1, [1])+$ is a little
better (MSPE=7.2) than CARar$(1)$,
but not by much.
One reason, therefore surely, is to prefer GNAR due to Occam's razor as it is a bit better, but also
a simpler explanation.
However, there are other
reasons to prefer GNAR, as described below.
  \begin{table}[H]
  {\color{red}
    \centering
    \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.}
    \begin{tabular}{crr}
        Model & MSPE(sd) & Mean No. Parameters \\
        \hline
         GNAR(1, [1]) &  7.4(2.98) & 2 \\
         Sparse VAR(1) & 11.4(2.80) & 3089 \\
         Res. VAR(1) & 10.6(2.48) & 3773 \\
         VAR(1) & 12.3(3.18) & 19600 \\
         AR(1) & 87.4(5.15) & 140 \\
         CARar(1) & 7.5(3.10) & 5 \\
         GNAR(1, [1])+ & 7.2(2.82) & 2 \\
         Naive & 7.6(3.12) & 0
    \end{tabular}
    
    \label{tab: prediction comparison with CAR}}
\end{table}

Statistical inference in the two software packages is carried out
differently. GNAR uses least-squares and CARBayesST is fully Bayesian. The latter means that
CARBayesST does adopt and use a coherent
Bayesian approach. However, it is \underline{extremely
computationally expensive}. Table~\ref{tab: execution time comparison} shows that for similar
problems CARar$(1)$ is over one thousand
times slower than GNAR.
\begin{table}[H]
{\color{red} 
\centering
\caption{Execution time comparison between global-$\alpha$ GNAR models and CARar models. The total execution time is the time computing one-step ahead forecasts ten times (i.e., the forecasts used for producing Table \ref{tab: prediction comparison with CAR}) takes in seconds. The average time is the total time divided by ten.}
\begin{tabular}{crr}
     Model & total & average \\
     \hline
     GNAR(1, [1]) & 6.2 & 0.62 \\
     CARar(1) & 7821.9 & 782.19
\end{tabular}

\label{tab: execution time comparison}}
\end{table}    
One of the reasons for
the relatively slow computation is that
CARBayesST uses MCMC, which is widely
acknowledged to be slow, especially as the
dimensionality of the problem increases.
The other, perhaps insurmountable, problem
is that the CARar model requires the inversion
(or, at least, equation solving) of a square
matrix of dimension of the number of nodes at
each step: $d \times d$.
So, for the COVID data, this is
a matrix of dimension $140\times 140$.
In the GNAR problem a square
matrix of dimension `number of parameters'
or $q \times q$, using the notation of
the paper, is inverted in the $(X^T X)^{-1}$ part of the parameter estimation. An advertised benefit of GNAR
models is their parsimonity, so for
the GNAR$(1, [1])$ model, e.g., we have
$q=2$ and the dimension of the matrix
to invert in the least squares problem
$X^T X$ is only $2 \times 2$, which is
rapidly computable.

Another reason to prefer GNAR for
network time series problems is that easily
permits the modelling of associations at different time lags,
different neighbourhood stages and a combination of both.
This aligns with the usual time series
paradigm of considering and selecting
lags and more general neighbourhood,
compared to the CARar we have seen 
that only considers two time lags and
nearest neighbours (just adjacency matrix). So, GNAR has the
extra flexibility to choose an appropriate
lagged model, $p$, and number $r$ of
\mbox{$r$-stage} neighbours, but CARar,
at the present time, does not. One
could envisage extending CARar, perhaps,
but then, given the onerous computational
effort for the current model, would it
be practically feasible to compute for
the equivalent or more lags and neighbours and,
more to the point, perform variable selection
for $p$ and $r$? 

GNAR is also capable of handling networks that are not fully connected, whereas the model behind
CARar does not seem able to be. This is because the $Q$ precision matrix is given by
\begin{equation}
   Q(S_1, \rho_S) = \rho_S \left\{ \operatorname{diag} \left( \sum_{j=1}^d S_{1, i, j}\right) - S_1 \right\} + (1 - \rho_S) I_d,
\end{equation}
where $S_1$ is the usual network adjacency matrix in the notation of our paper and $\operatorname{diag}(v_1, \ldots, v_d)$ is the
diagonal matrix with $v_1, \ldots, v_d$ on the main diagonal. Here, $\rho_S$ is a hyperparameter, which acts as 
a kind of ridge parameter.

GNAR and Bayesian CAR-type models handle missing data differently, with different interpretations. GNAR
is especially tolerant of large amounts of missing data due to its parsimonious nature --- large, sometimes
huge, numbers
of (pairs of) observations contribute to the estimation of parameters, so even sizeable quantities of missing
data do not affect the parameter estimates unduly. However, we know that for some specialised variants
within the {\tt CARBayesST} package will not operate with missing data. For the CARar used here it can be
accommodated in the usual nice coherent Bayesian manner of `integrating out', but we believe it increases
the already significant computation cost further. It would also be possible to embed GNAR in a Bayesian
context (indeed Matt Nunes at Bath wrote a version using dynamic linear models via
the {\tt dlm} R package written by Petris and Gilks (2022)) and similar `integrating out'
could be used to deal with the missing values as a coherent Bayesian alternative.
   
Another key GNAR advantage is that the weights matrix does not have to be positive definite, it's just a form of ``projecting'' data in a network informed manner, not a covariance matrix.].

We believe there are several other advantages of GNAR over CARar-type models, which perhaps we should not
go into too much detail here, but merely mention a few points:
\begin{itemize}
    \item CARar difficult to extend, e.g.\ adding a functional mean or penalty term, which would be easy for GNAR.
    %\item Does not adapt easily to tensor-linked data settings
    \item {\bf Forecasting difficulty:} No clear $n$-step ahead forecast computation, not even for $n = 1$.
    \item Difficult to compare/mix multiple CAR models (e.g.\ with GNAR we can imagine multiple networks for the
        same vector time series and then do bagging).
%    \item Does not extend easily to mixed type data
    \item Random effects in CARar are perhaps difficult to interpret
    \item Constrained to one-stage neighbours (in current formulations)
    \item Software limited to lag-two models (might not capture higher lag seasonality)
    %\item Assumes static network and network weights, i.e., difficult to extend to dynamic network settings
    %\item Assumes that spatial autocorrelations are stationary
    %\item {\bf MCMC estimates:} could be difficult to reproduce
    %\item Not as easy to interpret as GNAR
\end{itemize}

All in all, we think that the insertion of the CARar model
and results is interesting, although GNAR is demonstrably more flexible
and faster. We think your suggestion of looking at CARar was
fascinating and  perhaps \underline{enhances the potential for interesting
discussion} of this paper.
}

{\em 6. On this point, the 120km selection is based on best MSE predictive distance for a one-stage neighbourhood regression. Please explain more. Prediction at which t? Is there a potential problem in using the data twice, thus biasing comparisons with VAR or other methods?}

{\color{brown} COMPLETE}

{\color{red} We use MSE for the best fitting model over all $t$ to choose the 120km and the explanation
    has been expanded, and hopefully clarified, in the third paragraph of Section~5. Probably, if
    we then went on to do `traditional' model fitting with either looking at \mbox{$p$-values}
    in an ANOVA table or AIC/BIC values, then there is a potential issue over 
    `double-dipping'. However, with the 120km selection, we then go on to choose 
    models either through NACF/PNACF or predictive performance assessment. It is not clear to us
    how or whether the `double-dipping' is a problem, especially in the latter context?
    However, you raise a timely
    and
    interesting point and so we have inserted a short paragraph about this in the (renamed) Section~6
    on `Conclusion and Further Work'.}



{\em 7. A tessellation of the country based on trust locations, with neighbours defined by shared edges, seems a natural way to create a network. Using a fixed distance to- tally ignores population density. If you used a tessellation then presumably first-order
1
neighbours in urban areas would be geographically close, but first-order neighbours in rural areas would be more spread out. This seems sensible.}

{\color{brown} COMPLETE}

{\color{red} The NHS Trusts are already roughly spread according to population density.
We have looked, but we do not think that this is due to any  rule or guideline, but it is
partly historical and responding to needs of population centres as they grew.

Thank you for your suggestion about tessellation. We did consider the Dirichlet tessellation
and using points as immediate neighbours if they shared a Dirichlet cell boundary.
Unfortunately, such a tessellation is not great for this purpose because of points near
to the `boundary' of England, which results in (the well-known phenomenon of)
long thin triangles. So, for example, we did try it, but it ends up with long connections
between, e.g. James Paget and East Kent (195 road miles) or Northern Devon and Wirral
Community (275 road miles). These connections, and others like them, do not seem realistic or
useful for the COVID problem.

We have added a short paragraph in the further work section~6 that mentions alternative possibilities for
networks and disadvantages of the Dirichlet tessellation.
}


{\em 8. It is impressive to see such a parsimonious model do so well in prediction. Unless I’ve missed something, I think the points used for prediction are towards the end of the time series. Are values low at that time? Maybe also look at predictions at earlier times during the epidemic, when case numbers are high? Using only data before that time for estimation of course.}


{\color{brown} COMPLETE}

{\color{red}  The point you raise
is an excellent one --- it might be 
easy to forecast something `low' when the
past has been `low'. The RMarkdown document,
mentioned above, carries out a parallel
analysis, but chooses the end point
for training at an
earlier time at the peak of the second
wave (at $t=300$).
The forecasting
results from this second example are
presented in the table at the end of
Section~3.2 (the
original ones appear in the middle of that
section). The RMarkdown document also
contains full code that demonstrates
what we did, and also shows the individual
mean-squared prediction errors for each
of the ten separate rolling forecasts.
The GNAR results for are similarly superior
for high values as for low values and 
provide some evidence of reproducibility
across the series.

In addition, a multi-step ahead forecasts from $t=300$ is
now included as Table~3 in Section~5.2.
}

{\em 9. Should we have more on the data? Perhaps mean and standard deviation over time? Not too important.}

{\color{brown} COMPLETE}

{\color{red} Certainly, yes. We have added several plots relating to these
    items at the beginning of Section~3 in the RMarkdown document.}

{\em 10. I suggest that ``influence'' might be a misleading term when used in the globindex discussion. In standard regression, ``influence'' usually means the effect a particular case has on the parameter estimate, depending on both design x and response y. The globindex measure (17) depends on design only, and effectively simply measures how many other trusts are within 75 miles. London and Birmingham are about 120 miles apart and so trusts roughly midway will presumably be within 75 miles of the majority of both conurbations. Similarly Leeds and Birmingham are roughly 120 miles apart so trusts between them will be neighbours of many other trusts, including Liverpool and Manchester. Hence no surprise in Figure 7. The role of these (or other) trusts in viral spread isn’t known.}



{\color{brown} COMPLETE}

{\color{red} The globindex measure [was formula~(17), now is formula~(18)] does not just depend on the design. The formula is
    subtle, but it depends on $r^{\ast}$ \underline{and} chosen lag $p$, both of
    which are estimates that depend on the multivariate time series
    values {\em and} the network design.
    The pictures in the paper, particularly Figures~6c and~7, are based on $r^{\ast} = 1$. However, on reading our previous text, we did
    not explicitly state the value of $r^{\ast}$ in this section, but assumed its value as selected from the end of the previous
    section. We have added a small piece of text to state the selected value of $r^{\ast}$ at about lines 4--5 on page~24. To confirm
    that the influence pictures look different, for different values of $r^{\ast}$, the RMarkdown document in sections 3.4.1 and 3.4.2
    there show quite different influential Trust locations compared to Figure~7 in the paper, so choice of $r^{\ast}$ does make a difference
    and it does not just depend on network design.}


{\em 11. Finally some minor points. You can have also a few typos and trivia that I spotted while reading.}

{\color{brown} INCOMPLETE (Daniel/Guy to sit down and improve size of smaller labels)}

{\em i Just under Theorem 1: ``contrains''.}

{\color{brown} COMPLETE}


{\color{red} Done.}

{\em ii Perhaps think about changing the colour scale in Figure 2 and later. Yellow corresponds to the highest values but it looks the weakest in the plot, so the high values don’t stand out. A very bold colour for large values would be better.}

{\color{brown} COMPLETE}


{\color{red} Changing the colour map is very easy by using the \verb+viridis_color_option+
    in the \verb+corbit_plot+ function. The RMarkdown document shows several options
    with different colour maps, e.g.\ in that document's Section~2.1 and~2.2.}

{\em iii The captions to Figures 2 and 3 say ``200 realizations''. Does this mean that the points are average values, so masking sample to sample variation?}

{\color{brown} COMPLETE}


{\color{red} Apologies, this was poor phrasing on our part. We mean $T=200$ samples
and have changed the captions accordingly.}

{\em iv The figure labels and legends should be larger. Remove the vertical axis label in Figure 6c.}

{\color{brown} COMPLETE}


{\color{red}
We have taken a look at all figures and changed/improved them as follows.
We have increased the overall size of Figure~1. In principle, we could have inflated the
labels themselves, but this would make the 1-stage network picture even more confused with its already
overlapping labels. It is very difficult to plot labels on nodes when many of the nodes are located so
close to each other, for example, in the north west or London. For most
of the other figures
we have increased the global size of the plots, but if you feel the labels inside any plots need to be bigger
still, please let us know. With Figure~6, we have removed 6a and 6b.
For 6c (now called 6a)
we increased the colour bar legend text size as suggested.
An alternative plot to 6c are the relevances plotted according
to node location in the R Markdown document and a reference to it
has been inserted into the caption.}


{\em v Section 4.1 doesn’t contribute much and could be removed.}

{\color{brown} COMPLETE}


{\color{red}
We would prefer to keep this section for the following reasons.
We are often asked and challenged on the differences between
our work on GNAR and graphical models as there are
superficial similarities. However, there are genuine
differences not least that GNAR models can be constructed with
a genuine ``hard'' network (such as a railway), but also can be
constructed via partial correlation graph and, of course, GNAR models
are extended in terms of not only considering first stage neighbours.
As well as this, in terms of the current submission, this section
provides background and buildup to the node relevance
statistic defined in Section~5.4 (we have added a forward reference
to that section in Section~4.1).

If, of course, you still would like us to remove this section, then
we could move it to the Supplementary Material and then insert
appropriate cross-references from the main text to it.}

{\em vi When mentioning Figure 6 in para 3, say something about it coming later.}

{\color{brown} COMPLETE}


{\color{red} We have added the word ``below'' to indicate this.}

{\em vii Figures 4 and 5 are nice. The reader, or at least this one, could have more help with each of the four panels in Figure 6. What messages are we supposed to get?}

{\color{brown} COMPLETE}


{\color{red} Agreed, we have moved Figures 6a and 6b to the RMarkdown document. We have given further guidance and explanation of Figures
6c and 6d [now they are called 6a and b]. In particular, your
comment influenced us to change our notation from
$\text{local}(i, j)$ to $\text{local}_i (j)$ as it is not
a symmetric function. E.g.\ University Hospitals Plymouth heavily influences
the Royal Cornwall Hospitals Trust, but not vice versa. 
}

{\color{brown} COMPLETE}


{\em viii For calibration purposes, you might add to the table in 5.2 the MSPE using simple last observation carried forward? And perhaps an extrapolation of the last two observations.}

{\color{red} We have updated Table~2 in Section~5.2 to include the
simple last observation (which we refer to as Naive) and also the CARar
results.}


{\em ix Top of page 20: ``Table demonstrates the superlative performance...'' Don’t use superlatives like ``superlative'' about your own work. Let the reader decide.}

{\color{brown} COMPLETE}


{\color{red} We have removed the word `superlative' and rewritten the sentence. We also removed the
comparative word `superior' in the abstract,
and replaced by `excellent' (as it demonstrably
is).}

{\em x The same para mentions ``simulations''. Are these simulations?}

{\color{brown} COMPLETE}


{\color{red} We replaced the word ``simulation'', which was wrong.}



{\em xi Last para of Section 5.2: “Tables 6 thru 15”. Standard spelling please!}

{\color{brown} COMPLETE}

{\color{red} We have changed `thru' to `to'.}


{\em xii Is (18) presented correctly? Should there be a sum over r in the first term on the right?}

{\color{brown} COMPLETE}


{\color{red} The formula was correct (it is now (19)), but we did
    not make clear the range of $r$, which we have added as
    well as a clarifying remark afterwards that says that
    $i$ and $j$ can only be \mbox{$r$-stage} neighbours for one $r$.}

{\em xiii In looking into the literature a came across this paper: A general modeling frame- work for network autoregressive processes, by Yin, Safikhani and Michailidis (Technometrics, 2023). It looks worth a mention.}


{\color{brown} COMPLETE}


{\color{red} Thank you very much for this useful contribution and we have added the
paper to our literature review. It introduces some neat ideas
such as per-node network parameter estimation, testing for their
equality and related asymptotic theory. However, they stick with the
nearest-neighbour only regime and hence don't deal with a key topic 
of this paper, which is model selection.}





\begin{comment}
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.4]{images/rectangular_corbit_plot.pdf}
    \caption{Rectangular Corbit, first version. NACF values.}
    \label{fig:rect corbit one}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[scale=0.4]{images/nacf_sim1_rectangular_version_two.pdf}
    \caption{Rectangular Corbit, second version. NACF values.}
    \label{fig:rect corbit two}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[scale=0.4]{images/pnacf_sim1_rectangular_version_two.pdf}
    \caption{Rectangular Corbit, second version. PNACF values.}
    \label{fig: rect corbit 3}
\end{figure}
\end{comment}

\begin{figure}[H]
    \centering
    \includegraphics[scale=0.50]{images/node_relevance_network_plot.pdf}
    \caption{Node relevance hierarchy imposed on the network. }
    \label{fig: node relevance network}
\end{figure}


\end{document}