# Are Current Decoding Strategies Capable of Facing
the Challenges of Visual Dialogue?

Amit Kumar Chaudhary
CIMeC, University of Trento
<EMAIL_ADDRESS>
&Alex J. Lucassen
CIMeC, University of Trento
<EMAIL_ADDRESS>
Ioanna Tsani
CIMeC, University of Trento
<EMAIL_ADDRESS>
&Alberto Testoni
DISI, University of Trento
<EMAIL_ADDRESS>
###### Abstract

Decoding strategies play a crucial role in natural language generation
systems. They are usually designed and evaluated in open-ended text-only
tasks, and it is not clear how different strategies handle the numerous
challenges that goal-oriented multimodal systems face (such as grounding and
informativeness). To answer this question, we compare a wide variety of
different decoding strategies and hyper-parameter configurations in a Visual
Dialogue referential game. Although none of them successfully balance lexical
richness, accuracy in the task, and visual grounding, our in-depth analysis
allows us to highlight the strengths and weaknesses of each decoding strategy.
We believe our findings and suggestions may serve as a starting point for
designing more effective decoding algorithms that handle the challenges of
Visual Dialogue tasks.

##  1 Introduction

The last few years have witnessed remarkable progress in developing efficient
generative language models. The choice of the decoding strategy plays a
crucial role in the quality of the output (see Zarrieß et al. (2021) for an
exhaustive overview). It should be noted that decoding strategies are usually
designed for and evaluated in text-only settings. The most-used decoding
strategies can be grouped into two main classes. On the one hand, decoding
strategies that aim to generate text that maximizes likelihood (like greedy
and beam search) are shown to generate generic, repetitive, and degenerate
output. Zhang et al. (2021) refer to this phenomenon as the likelihood trap,
and provide evidence that these strategies lead to sub-optimal sequences. On
the other hand, stochastic strategies like pure sampling, top-k sampling, and
nucleus sampling Holtzman et al. (2020) increase the variability of generated
texts by taking random samples from the model. However, this comes at the cost
of generating words that are not semantically appropriate for the context in
which they appear. Recently, Meister et al. (2022) used an information-
theoretic framework to propose a new decoding algorithm (typical decoding),
which samples tokens with an information content close to their conditional
entropy. Typical decoding shows promising results in human evaluation
experiments but, given its recent release, it is not clear yet how general
this approach is.

Multimodal vision & language systems have recently received a lot of attention
from the research community, but a thorough analysis of different decoding
strategies in these systems has not been carried out. Thus, the question
arises of whether the above-mentioned decoding strategies can handle the
challenges of multimodal systems. i.e., generate text that not only takes into
account lexical variability, but also grounding in the visual modality.
Moreover, in goal-oriented tasks, the informativeness of the generated text
plays a crucial role as well. To address these research questions, in this
paper we take a referential visual dialogue task, GuessWhat?! De Vries et al.
(2017), where two players (a Questioner and an Oracle) interact so that the
Questioner identifies the secret object assigned to the Oracle among the ones
appearing in an image (see Figure 1 for an example). Apart from well-known
issues, such as repetitions in the output, this task poses specific challenges
for evaluating decoding techniques compared to previous work. On the one hand,
the generated output has to be coherent with the visual input upon which the
conversation takes place. As highlighted by Rohrbach et al. (2018); Testoni
and Bernardi (2021b), multimodal generative models often generate hallucinated
entities, i.e., tokens that refer to entities that do not appear in the image
upon which the conversation takes place. On the other hand, the questions must
be informative, i.e., they must help the Questioner to incrementally identify
the target object.

We show that the choice of the decoding strategy and its hyper-parameter
configuration heavily affects the quality of the generated output. Our results
highlight the specific strengths and weaknesses of decoding strategies that
aim at generating sequences with the highest probability vs. strategies that
randomly sample words. We find that none of the decoding strategies currently
available is able to balance task accuracy and linguistic quality of the
output. However, we also show which strategies perform better at important
challenges, such as incremental dialogue history, human evaluation,
hallucination rate, and lexical diversity. We believe our work may serve as a
starting point for designing decoding strategies that take into account all
the challenges involved in Visual Dialogue tasks.

Figure 1: Example of a GuessWhat game from De Vries et al. (2017)

##  2 Task & Dataset

GuessWhat?! De Vries et al. (2017) is a simple object identification game in
English where two participants see a real-world image from MSCOCO Lin et al.
(2014) containing multiple objects. One player (the Oracle) is secretly
assigned one object in the image (the target) and the other player (the
Questioner) has to guess it by asking a series of binary yes-no questions to
the Oracle. The task is considered to be successful if the Questioner
identifies the target. The dataset for this task was collected from human
players via Amazon Mechanical Turk. The authors collected 150K dialogues with
an average of 5.3 binary questions per dialogue. Figure 1 shows an example of
a GuessWhat game from the dataset.

##  3 Model and Decoding Strategies

We use the model and pre-trained checkpoints of the Questioner agent made
available by Testoni and Bernardi (2021c) for the GuessWhat?! task. This model
is based on the GDSE architecture Shekhar et al. (2019). It uses a ResNet-152
network He et al. (2016) to encode the images and an LSTM network to encode
the dialogue history. A multimodal shared representation is generated and then
used to train both the question generator (which generates a follow-up
question given the dialogue history) and the Guesser module (which selects the
target object among a list of candidates at the end of the dialogue) in a
joint multi-task learning fashion. Testoni and Bernardi (2021c) added an
internal Oracle module to the GDSE architecture, which guides a cognitively-
inspired beam search re-ranking strategy (Confirm-it) at inference time: this
strategy promotes the generation of questions that aim at confirming the
model’s intermediate conjectures about the target. In our work, at inference
time the Questioner agent always interacts with the baseline Oracle agent
proposed in De Vries et al. (2017).

We analyse the effect of a large number of decoding strategies as well as
hyper-parameter configuration for each strategy: as highlighted by Zhang et
al. (2021), it is crucial to evaluate different hyper-parameter configurations
when comparing multiple decoding strategies. Among the ones that maximize the
likelihood of the sequence, we consider plain beam search (with a beam size of
3) and greedy search. We also consider Confirm-it, the cognitively-inspired
beam search re-ranking strategy proposed in Testoni and Bernardi (2021c) for
promoting the generation of questions that aim at confirming the model’s
intermediate conjectures about the target. This strategy re-ranks the set of
candidate questions from beam search and selects the one that helps the most
in confirming the model’s hypothesis about the target. As for stochastic
strategies, we analyse pure sampling, top-k sampling (with different $k$
values), and nucleus sampling (with different $p$ values), a strategy proposed
in Holtzman et al. (2020) which selects the highest probability tokens whose
cumulative probability mass exceeds a given threshold $p$. We also consider
typical decoding (with different $\tau$ values), a recently proposed strategy
Meister et al. (2022) based on an information-theoretic framework. We refer to
the respective papers for additional details on decoding strategies. We let
the model generate 5 questions222Except for the accuracy per turn metric in
Section 5.3, where the dialogues consisted of 10 questions. at test time and
average the results over five random seeds.

##  4 Metrics

We are interested in evaluating different decoding strategies against a set of
metrics that reflect the complexity of the different skills required to
successfully solve multimodal referential games.

Linguistic Quality: We compute the percentage of games with at least one
repeated question, the overall number of unique words used by the model and,
in line with the observations in Testoni and Bernardi (2021a), the number of
rare words generated by the model, defined as those words that appear fewer
than 20 times in the training set.

Visual Grounding: To quantify the rate of object hallucination in the
generated dialogues, we compute the CHAIR metric Rohrbach et al. (2018);
Testoni and Bernardi (2021b). This metric, originally proposed for image
captioning, detects hallucination by checking each object mentioned in a
generated image caption against the ground-truth MSCOCO objects for that
image. The metric consists of two distinct variants: CHAIR-i, or per-instance
variant (number of hallucinated objects divided by the total number of objects
mentioned in each dialogue), and CHAIR-s, or per-sentence variant (number of
dialogues with at least one hallucination divided by the total number of
dialogues).333Testoni and Bernardi (2021b) first adapted the CHAIR metric for
Visual Dialogue. However, the authors did not investigate the effect of
different decoding strategies.

Informativeness: To study the informativeness of the generated questions, we
report the raw accuracy of the model in guessing the target object after each
dialogue turn and at the end of the dialogue. A game is considered successful
if the model identifies the target object assigned to the Oracle. Similarly,
we also report the accuracy of human annotators when guessing the target by
reading machine-generated dialogues.

##  5 Results

|  | Accuracy (%) $\uparrow$
---
CHAIR-i $\downarrow$ | CHAIR-s $\downarrow$ |  | % games with
---
repetitions $\downarrow$
| Vocabulary
---
Size $\uparrow$
| Rare
---
Words $\uparrow$
Confirm-it | 51.39 | 15.09 | 28.48 | 30.33 | 858 | 34
Beam Search (beam size = 3) | 47.05 | 18.33 | 31.08 | 38.49 | 731 | 27
Nucleus Sampling ($p$ = 0.3) | 46.92 | 17.96 | 33.60 | 32.35 | 1016 | 78
Greedy Search | 46.58 | 17.75 | 32.97 | 35.63 | 834 | 46
Typical Decoding ($\tau$ = 0.7) | 45.45 | 21.84 | 37.81 | 16.18 | 1703 | 247
Top-k Sampling ($k$ = 5) | 45.10 | 22.84 | 37.71 | 14.93 | 1462 | 171
Pure Sampling | 43.13 | 26.55 | 43.23 | 8.32 | 2609 | 793
Table 1: Comparison between decoding strategies and their best-performing (in
terms of accuracy) hyper-parameters. The decoding strategies are sorted by
accuracy.

###  5.1 Quantitative Results

Table 1 shows the performance of different decoding strategies against
accuracy and dialogue quality, as described by the metrics in Section 4.
444Here we only report the best-performing configuration for each decoding
strategy (see SM for all configurations). Confirm-it is by far the best
decoding strategy in terms of accuracy and hallucination rate. However, it
uses a restricted vocabulary compared to other strategies. A similar issue is
observed for greedy and beam search. We find nucleus sampling (with a
$p$-value of 0.3, much lower than the one used by the authors in Holtzman et
al. (2020)) to effectively increase the lexical variety compared to beam
search, without damaging accuracy and hallucination rate. Typical decoding,
top-k and pure sampling, instead, clearly decrease repetitions and increase
the vocabulary richness by generating tokens that are not related to the
source input, as indicated by the high hallucination rate. It thus looks like
there exists a trade-off between informativeness / visual grounding and
linguistic quality.

###  5.2 Effect of Hyper-Parameter Choice

Figure 2: Different hyper-parameter values and their effect on the accuracy,
hallucinations, and repetitions in typical decoding and nucleus sampling.

We study the effect of hyper-parameter configurations in stochastic
strategies. Specifically, we try various $p$-values for nucleus sampling and
$\tau$-values for typical decoding.555Results for top-k are in SM. As shown in
Figure 2, both typical and nucleus sampling peak in accuracy with the
parameter configurations that also lead to the most repetitions and fewest
hallucinations. Conversely, both strategies show the lowest accuracy with the
highest hallucination rate. These results confirm the detrimental effect of
hallucinations on the performance of the model. It is interesting to note the
robustness of typical decoding in generating few repetitions regardless of the
$\tau$ value. In line with the findings in Zhang et al. (2021), this analysis
confirms the importance of hyper-parameter configurations and the peculiar
trade-off between informativeness, repetitions, and visual grounding: so far
it has not been possible to find a single configuration that optimizes all
three at the same time.

###  5.3 Per-turn Accuracy

Figure 3: The accuracy per dialogue turn for four different decoding
strategies for dialogues of length 10.

One crucial ability in GuessWhat?! is asking informative questions that
incrementally help in identifying the target: for this reason, we check the
accuracy of the model after each new question is asked. Figure 3 shows
accuracy per dialogue turn for a set of representative strategies: Nucleus
sampling ($p$=0.3), Typical Decoding ($\tau$=0.7), Confirm-it, and pure
sampling. To get a broader picture, we let the model generate 10 questions in
this setting. Confirm-it stands out by showing the largest incremental
increase of accuracy throughout the dialogue, indicating that it generates
more effective follow-up questions. Pure sampling, on the other hand, seems to
suffer from the very beginning of the dialogue and its accuracy stabilizes
soon. It is worth noting that the accuracy of typical decoding gets closer to
that of nucleus sampling towards the end of the dialogue, with the latter
leveling off sooner. We conjecture that Confirm-it outperforms other
techniques because it takes into account the probability of the Guesser at
inference time, so it is guided to generate questions that change these
probabilities and thus avoid generic questions.

|  | Human Accuracy (%) $\uparrow$
---
Confirm-it | 72.5
Typical Sampling ($\tau$ = 0.7) | 68.0
Nucleus Sampling ($p$ = 0.3) | 67.5
Pure Sampling | 59.5
Table 2:  Human Guess Accuracy based on dialogue generated from different
decoding strategies.

###  5.4 Human Evaluation

We asked 8 human annotators to guess the target object in a sample of
GuessWhat?! games when reading dialogues generated by our model with different
decoding strategies. Each participant annotated 100 games (25 per strategy)
and the decoding strategy was not revealed during the annotation. As shown in
Table 2, humans reach the highest accuracy when reading dialogues generated by
Confirm-it, followed by typical decoding and nucleus sampling, while pure
sampling falls behind. These results, which do not mirror the accuracy result
in Table 1, allow us to disentangle the weaknesses of the Guesser (i.e., the
classification module that predicts the target) from the actual
informativeness of the dialogues. Compared to the model, human annotators seem
to better exploit the lexical richness of typical decoding and nucleus
sampling. We refer to the SM for additional information about the annotation
procedure, in line with the best-practice guidelines in van der Lee et al.
(2021).

##  6 Related Work

In the field of multimodal NLG, Zarrieß and Schlangen (2018) propose trainable
decoding for referring expression generation. The authors propose a two-stage
optimization set-up where a small network processes the RNN’s hidden state
before passing it to the decoder, using BLEU score as a reward for the
decoder. We did not analyse this approach in our paper because we focus only
on decoding strategies that do not require any change in the architecture or
training of the model. We leave for future work an analysis of trainable
decoding approaches. Inspired by the findings in Holtzman et al. (2020),
Massarelli et al. (2020) propose a hybrid decoding strategy for open-ended
text generation which combines the non-repetitive nature of sampling
strategies with the consistency of likelihood-based approaches. The authors
show that their approach generated less repetitive and more verifiable text.
The design of hybrid decoding strategies for multimodal tasks is out of the
scope of this paper, but is an interesting subject to pursue in future work.

##  7 Discussion and Conclusion

Decoding algorithms are a key component of natural language generation
systems. They are usually designed for and evaluated in text-only tasks. We
believe multimodal (vision & language) and goal-oriented tasks pose unique and
under-studied challenges to current decoding strategies. In this paper, we ran
an in-depth analysis of several decoding strategies (and their hyper-parameter
configurations) for a model playing a referential visual dialogue game. We
found that decoding algorithms that lead to the highest accuracy in the task
and the lowest hallucination rate, at the same time generate highly repetitive
text and use a restricted vocabulary. Our analyses reveal the crucial role of
hyper-parameter configuration in stochastic strategies, an issue that poses
several questions about the trade-off between lexical variety, hallucination
rate, and task accuracy. While nucleus sampling partially balances the above-
mentioned issues, human annotators seem to better exploit the richness of the
dialogues generated by typical decoding. Finally, our results demonstrate that
a beam search re-ranking algorithm (Confirm-it) generates more effective
follow-up questions throughout the dialogue turns. We believe that taking into
account the model’s intermediate predictions about the referent, like Confirm-
it does, represents a promising direction that should be applied also to
stochastic strategies in future work, aiming at preserving their lexical
richness while reducing hallucinations.

Our results demonstrate that none of the decoding strategies currently at
disposal effectively take into account both task accuracy and dialogue quality
at the same time. We also highlight peculiar features of each strategy that
may guide future research with the goal of designing decoding strategies that
properly confront the crucial challenges of multimodal goal-oriented
dialogues.

## Acknowledgements

We kindly acknowledge the support of the Department of Information Science and
Engineering (DISI) for the arrangement of the GPUs used in our research. We
would like to thank Professor Raffaella Bernardi for the suggestions and
feedback. The first three authors are enrolled in the Erasmus Mundus European
Masters Program in Language and Communication Technologies. The first author
is supported by the Erasmus+ Programme. The third author of this scientific
paper was supported by the Onassis Foundation - Scholarship ID: F ZR
060-1/2021-2022. We thank the Master’s students of the Grounded Language
Processing course at the University of Trento for the feedback provided at an
early stage of this project.

## References

  * De Vries et al. (2017) Harm De Vries, Florian Strub, Sarath Chandar, Olivier Pietquin, Hugo Larochelle, and Aaron Courville. 2017.  GuessWhat?! Visual object discovery through multi-modal dialogue.  In _2017 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2017, Honolulu, HI, USA, July 21-26, 2017_ , pages 4466–4475. IEEE Computer Society.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. 2016.  Deep residual learning for image recognition.  In _2016 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2016, Las Vegas, NV, USA, June 27-30, 2016_ , pages 770–778. IEEE Computer Society.
  * Holtzman et al. (2020) Ari Holtzman, Jan Buys, Li Du, Maxwell Forbes, and Yejin Choi. 2020.  The curious case of neural text degeneration.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net.
  * Lin et al. (2014) Tsung-Yi Lin, Michael Maire, Serge J. Belongie, James Hays, Pietro Perona, Deva Ramanan, Piotr Dollár, and C. Lawrence Zitnick. 2014.  Microsoft COCO: common objects in context.  In _Computer Vision - ECCV 2014 - 13th European Conference, Zurich, Switzerland, September 6-12, 2014, Proceedings, Part V_ , volume 8693 of _Lecture Notes in Computer Science_ , pages 740–755. Springer.
  * Massarelli et al. (2020) Luca Massarelli, Fabio Petroni, Aleksandra Piktus, Myle Ott, Tim Rocktäschel, Vassilis Plachouras, Fabrizio Silvestri, and Sebastian Riedel. 2020.  How decoding strategies affect the verifiability of generated text.  In _Findings of the Association for Computational Linguistics: EMNLP 2020_ , pages 223–235, Online. Association for Computational Linguistics.
  * Meister et al. (2022) Clara Meister, Tiago Pimentel, Gian Wiher, and Ryan Cotterell. 2022.  Typical decoding for natural language generation.  _CoRR_ , abs/2202.00666.
  * Rohrbach et al. (2018) Anna Rohrbach, Lisa Anne Hendricks, Kaylee Burns, Trevor Darrell, and Kate Saenko. 2018.  Object hallucination in image captioning.  In _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing_ , pages 4035–4045, Brussels, Belgium. Association for Computational Linguistics.
  * Shekhar et al. (2019) Ravi Shekhar, Aashish Venkatesh, Tim Baumgärtner, Elia Bruni, Barbara Plank, Raffaella Bernardi, and Raquel Fernández. 2019.  Beyond task success: A closer look at jointly learning to see, ask, and GuessWhat.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 2578–2587.
  * Testoni and Bernardi (2021a) Alberto Testoni and Raffaella Bernardi. 2021a.  The interplay of task success and dialogue quality: An in-depth evaluation in task-oriented visual dialogues.  In _Proceedings of the 16th Conference of the European Chapter of the Association for Computational Linguistics: Main Volume, EACL 2021, Online, April 19 - 23, 2021_ , pages 2071–2082. Association for Computational Linguistics.
  * Testoni and Bernardi (2021b) Alberto Testoni and Raffaella Bernardi. 2021b.  “I’ve seen things you people wouldn’t believe”: Hallucinating entities in GuessWhat?! In _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing: Student Research Workshop_ , pages 101–111, Online. Association for Computational Linguistics.
  * Testoni and Bernardi (2021c) Alberto Testoni and Raffaella Bernardi. 2021c.  Looking for confirmations: An effective and human-like visual dialogue strategy.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_ , pages 9330–9338, Online and Punta Cana, Dominican Republic. Association for Computational Linguistics.
  * van der Lee et al. (2021) Chris van der Lee, Albert Gatt, Emiel van Miltenburg, and Emiel Krahmer. 2021.  Human evaluation of automatically generated text: Current trends and best practice guidelines.  _Computer Speech & Language_, 67:101151.
  * Zarrieß and Schlangen (2018) Sina Zarrieß and David Schlangen. 2018.  Decoding strategies for neural referring expression generation.  In _Proceedings of the 11th International Conference on Natural Language Generation_ , pages 503–512, Tilburg University, The Netherlands. Association for Computational Linguistics.
  * Zarrieß et al. (2021) Sina Zarrieß, Henrik Voigt, and Simeon Schüz. 2021.  Decoding methods in neural language generation: A survey.  _Information_ , 12(9):355.
  * Zhang et al. (2021) Hugh Zhang, Daniel Duckworth, Daphne Ippolito, and Arvind Neelakantan. 2021.  Trading off diversity and quality in natural language generation.  In _Proceedings of the Workshop on Human Evaluation of NLP Systems (HumEval)_ , pages 25–33, Online. Association for Computational Linguistics.

##  Appendix A Supplementary Material

###  A.1 Effect of Hyper-parameters

Figures 4, 5, and 6 illustrate how hyper-parameter choice affects the
accuracy, the hallucinations, and the repetitions. Top-k sampling (Figure 4)
shows decreased accuracy and repetitions, and increased hallucinations, as the
$k$-value gets higher. The same general pattern can be observed with the
gradual increase of the $p$-value in nucleus sampling (Figure 6). On the other
hand, typical decoding accuracy peaks at $\tau$ = 0.7 (Figure 5). This is also
the point at which the repetitions are at their highest and the hallucinations
are at their lowest. Both very high and very low $\tau$-values cause lower
accuracy, fewer repetitions, and an increase of hallucinations.

Figure 4: Hyper-parameter choices’ effect on the accuracy, hallucinations, and repetitions in top-k sampling. Figure 5: Hyper-parameter choices’ effect on the accuracy, hallucinations, and repetitions in typical decoding. Figure 6: Hyper-parameter choices’ effect on the accuracy, hallucinations, and repetitions in nucleus sampling. |  | % Accuracy $\uparrow$
---
CHAIR-i $\downarrow$ | CHAIR-s $\downarrow$ |  | % games with
---
repetitions $\downarrow$
| Vocabulary
---
Size $\uparrow$
| Rare
---
Words $\uparrow$
| mean | std | mean | std | mean | std | mean | std | mean | std | mean | std
CI | 51.39 | 0.00 | 15.09 | 0.00 | 28.48 | 0.00 | 30.33 | 0.00 | 858 | 0.0 | 34 | 0.0
BS (beam = 3) | 47.05 | 0.00 | 18.33 | 0.00 | 31.08 | 0.00 | 38.49 | 0.00 | 731 | 0.0 | 27 | 0.0
NS ($p$ = 0.3) | 46.92 | 0.17 | 17.96 | 0.09 | 33.60 | 0.14 | 32.35 | 0.21 | 1016 | 8.3 | 78 | 5.5
NS ($p$ = 0.2) | 46.90 | 0.11 | 17.65 | 0.06 | 34.79 | 0.12 | 46.41 | 0.18 | 919 | 7.1 | 59 | 4.5
NS ($p$ = 0.4) | 46.73 | 0.27 | 18.38 | 0.17 | 34.27 | 0.28 | 29.16 | 0.25 | 1097 | 11.9 | 88 | 6.4
NS ($p$ = 0.1) | 46.63 | 0.02 | 17.76 | 0.01 | 32.98 | 0.01 | 35.66 | 0.06 | 839 | 1.6 | 46 | 1.9
NS ($p$ = 0.5) | 46.62 | 0.17 | 19.10 | 0.09 | 35.03 | 0.10 | 25.73 | 0.26 | 1192 | 18.0 | 103 | 3.4
GS | 46.58 | 0.00 | 17.75 | 0.00 | 32.97 | 0.00 | 35.63 | 0.00 | 834 | 0.0 | 46 | 0.0
NS ($p$ = 0.6) | 46.13 | 0.38 | 20.04 | 0.15 | 35.69 | 0.35 | 22.35 | 0.26 | 1303 | 12.7 | 126 | 9.4
NS ($p$ = 0.7) | 45.85 | 0.19 | 21.14 | 0.11 | 36.82 | 0.31 | 19.11 | 0.26 | 1451 | 9.0 | 162 | 9.5
TD ($\tau$ = 0.7) | 45.45 | 0.32 | 21.84 | 0.15 | 37.81 | 0.23 | 16.18 | 0.29 | 1703 | 13.0 | 247 | 12.6
NS (p = 0.8) | 45.38 | 0.14 | 22.20 | 0.16 | 37.97 | 0.30 | 15.80 | 0.19 | 1643 | 23.0 | 219 | 12.9
TD ($\tau$ = 0.6) | 45.29 | 0.28 | 22.08 | 0.16 | 38.17 | 0.30 | 15.75 | 0.17 | 1723 | 21.9 | 248 | 18.3
TD ($\tau$ = 0.8) | 45.16 | 0.18 | 22.21 | 0.20 | 37.93 | 0.29 | 15.32 | 0.28 | 1712 | 10.8 | 244 | 13.6
TD ($\tau$ = 0.5) | 45.12 | 0.15 | 22.60 | 0.17 | 38.69 | 0.36 | 14.89 | 0.22 | 1745 | 7.3 | 262 | 8.7
Top-k ($k$ = 5) | 45.10 | 0.27 | 22.84 | 0.21 | 37.71 | 0.26 | 14.93 | 0.10 | 1462 | 12.6 | 171 | 5.2
TD ($\tau$ = 0.4) | 44.83 | 0.17 | 23.11 | 0.13 | 39.11 | 0.44 | 13.94 | 0.24 | 1755 | 19.0 | 265 | 12.2
TD ($\tau$ = 0.9) | 44.61 | 0.16 | 23.74 | 0.18 | 39.59 | 0.19 | 12.41 | 0.25 | 1919 | 13.5 | 334 | 9.1
Top-k ($k$ = 10) | 44.61 | 0.24 | 24.03 | 0.29 | 39.62 | 0.26 | 11.96 | 0.16 | 1692 | 13.8 | 235 | 10.5
NS ($p$ = 0.91) | 44.60 | 0.15 | 23.92 | 0.13 | 39.79 | 0.23 | 12.21 | 0.13 | 1948 | 22.3 | 342 | 14.0
NS ($p$ = 0.9) | 44.56 | 0.23 | 23.82 | 0.07 | 39.62 | 0.17 | 12.44 | 0.10 | 1912 | 20.2 | 332 | 13.6
TD ($\tau$ = 0.91) | 44.45 | 0.27 | 23.92 | 0.14 | 39.82 | 0.31 | 12.18 | 0.19 | 1945 | 11.7 | 345 | 16.1
TD ($\tau$ = 0.3) | 44.38 | 0.31 | 24.07 | 0.21 | 39.88 | 0.22 | 13.14 | 0.23 | 1791 | 14.8 | 278 | 13.7
NS ($p$ = 0.95) | 44.38 | 0.12 | 24.93 | 0.15 | 41.08 | 0.24 | 10.56 | 0.05 | 2129 | 11.3 | 438 | 11.4
TD ($\tau$ = 0.95) | 44.09 | 0.21 | 24.83 | 0.24 | 40.82 | 0.28 | 10.51 | 0.20 | 2117 | 18.9 | 435 | 17.1
Top-k ($k$ = 20) | 43.89 | 0.10 | 25.12 | 0.30 | 41.00 | 0.39 | 10.39 | 0.17 | 1879 | 23.1 | 305 | 17.5
TD ($\tau$ = 0.2) | 43.36 | 0.20 | 25.19 | 0.16 | 40.82 | 0.34 | 12.11 | 0.09 | 1815 | 21.5 | 287 | 10.7
PS | 43.13 | 0.28 | 26.55 | 0.25 | 43.23 | 0.36 | 8.32 | 0.17 | 2609 | 9.3 | 793 | 11.4
TD ($\tau$ = 0.1) | 42.85 | 0.15 | 26.25 | 0.14 | 41.54 | 0.09 | 11.12 | 0.11 | 1825 | 18.4 | 286 | 13.5
Table 3: Comparison between decoding strategies and their hyper-parameters (CI
= Confirm-it, BS = Beam Search, NS = Nucleus Sampling, GS = Greedy Search, TD
= Typical Decoding, Top-k = Top-k Sampling, PS = Pure Sampling).

###  A.2 Experiments

Table 3 presents our results in detail for all the parameter configurations we
considered. We have computed accuracy percentage, CHAIR-i, CHAIR-s, percentage
of games with repeated questions, vocabulary size and number of rare words for
each decoding method and its respective hyper-parameter configurations. These
results are sorted by decreasing accuracy. The 3 best results of each metric
are in bold.

###  A.3 Human Annotation Details

Figure 7: Example of the games displayed to the participants for the
annotation task. Participants had to select one target object among the list
of candidate objects on the right. The machine-generated dialogue is in the
red box.

The annotation was done by 8 human annotators on a sample of GuessWhat?!
games. They were recruited within our organization on a voluntary basis and
they did not receive any payment for the annotation. Written informed consent
was obtained from all the participants. Participants were 4 males and 4
females with high educational level and from different ethnic groups. Before
the beginning of the annotation task, each annotator was briefed on the
GuessWhat?! gameplay and purpose, and was asked to annotate some sample games
in order to get familiar with the annotation process. We used the makesense.ai
online software for image recognition. Each image had a minimum of 3 and a
maximum of 6 candidate objects. The annotators could see both the bounding box
and the category for each candidate object in the image. They could also see
the full dialogue between the Questioner and the Oracle. The annotators then
had to pick the object they believed was the right one, based on the
information given by the dialogue. Figure 7 provides an example of the games
we asked the participants to annotate. Overall, we extracted 200 images from
the GuessWhat?! test set and generated one dialogue per decoding strategy per
image. In total, we thus generated 800 dialogues. Each human participant
annotated 25 images per decoding strategy. To prevent biases, the participants
were not given any information about the decoding strategy that each image
represented.

###  A.4 Comparison of Decoding Strategies

Table 4 shows a comparison between generated dialogues of various decoding
strategies (see Figure 8 for the image the dialogues refer to).

Figure 8: The image for game 147344. Turn | Speaker | Confirm-it | Nucleus Sampling | Typical Decoding | Pure Sampling
---|---|---|---|---|---
0 | Questioner | is it a bird? | is it a bird? | is it a bird? | is it a bird?
Oracle | yes | yes | yes | yes
1 | Questioner | is it on the left side? | is it on the left side? | is it on the air? | is it close?
Oracle | yes | yes | yes | yes
2 | Questioner |  | is it on the left side
---
of the group?
| is it one of the
---
first birds?
is it 1 of the 3 on the right? | is it squatting?
Oracle | no | yes | no | yes
3 | Questioner | is it in the middle? |  | is it the first one
---
from the left?
| is it one of the three birds
---
that are next to each other?
| is it one of the three
---
all the way to the right?
Oracle | yes | no | yes | yes
4 | Questioner |  | is it one of the three
---
on the right?
| is it the second one
---
from the left?
is it further left? | is it 1st?
Oracle | no | no | yes | no
Table 4: Comparison between dialogues generated by various decoding strategies
for the image reported in Figure 8.


# VARX Granger Analysis: Modeling, Inference, and Applications

Lucas C Parra
Biomedical Engineering
City College of New York
<EMAIL_ADDRESS>
&Aimar Silvan
Biomedical Engineering
City College of New York
<EMAIL_ADDRESS>
&Maximilian Nentwich
The Feinstein Institutes for Medical Research
Northwell Health
<EMAIL_ADDRESS>
&Jens Madsen
Biomedical Engineering
City College of New York
<EMAIL_ADDRESS>
&Behtash Babadi
Electrical and Computer Engineering
University of Maryland
<EMAIL_ADDRESS>
###### Abstract

Vector Autoregressive models with exogenous input (VARX) provide a powerful
framework for modeling complex dynamical systems like brains, markets, or
societies. Their simplicity allows us to uncover linear effects between
endogenous and exogenous variables. The Granger formalism is naturally suited
for VARX models, but the connection between the two is not widely understood.
We aim to bridge this gap by providing both the basic equations and easy-to-
use code. We first explain how the Granger formalism can be combined with a
VARX model using deviance as a test statistic. We also present a bias
correction for the deviance in the case of L2 regularization, a technique used
to reduce model complexity. To address the challenge of modeling long
responses, we propose the use of basis functions, which further reduce
parameter complexity. We demonstrate that p-values are correctly estimated,
even for short signals where regularization influences the results.
Additionally, we analyze the model’s performance under various scenarios where
model assumptions are violated, such as missing variables or indirect
observations of the underlying dynamics. Finally, we showcase the practical
utility of our approach by applying it to real-world data from neuroscience,
physiology, and sociology. To facilitate its adoption, we make Matlab, Python,
and R code available here: https://github.com/lcparra/varx

##  1 Introduction

The Granger formalism relies on temporal precedence to establish the direction
of effects: It tests whether the linear prediction of a variable $y(t)$ from
its past can be improved by adding the history of another variable $x(t)$ as
predictor [1]. If it does, we say that $x$ has an "effect" on $y$.111In this
work we avoid calling the effect ”causal”, as there are too many caveats to
this interpretation; see Discussion. This basic idea had originally been
proposed by Wiener [2] leading to the term "Wiener-Granger Causality" [3]. For
multiple variables, this linear prediction naturally fits in the formalism of
vectorial autoregressive (VAR) models [4]. The VAR-Granger formalism is
appealing to neuroscientists, economists, and sociologists because it allows
one to quantify the strength and direction of effects in interactive dynamical
systems, such as brains [3, 5], markets [6], and societies [7]. When such
systems are subject to an external stimulus, the appropriate model is a VAR
model with eXogenous input – the VARX model.

Exogenous variables were incorporated into the Granger formalism as
conditional dependence by Geweke [8]. In practice, this "conditional
causality" has been used to control for spurious correlations due to common
causes. Code implementations of this idea [9, 10] only use exogenous variables
to remove confounds. In contrast, we propose to model the total system
response as a combination of both exogenous effects and endogenous dynamics.
In this view, the exogenous effects are not a nuisance, but an important
component of the model to be estimated, encompassing multiple time delays.
Although economists have employed VARX models to capture exogenous effects,
the use of the Granger formalism to establish the effects of individual
variables is not as widely used. Indeed, the correspondence of the
"conditional Granger-causality" [3, 11] with VARX models is not well known.
While the VARX model is common in statistics toolboxes, we are not aware of
any implementation of the VARX model with the Granger-Geweke test to assess
the effects.

Here we apply the Granger-Geweke formalism in the language of VARX models,
capturing both internal dynamics and exogenous inputs. We summarize the basic
equations needed to estimate model parameters, effect size, and statistical
significance. To cope with high-dimensional datasets and longer prediction
filters, we use L2 regularization and basis functions, effectively reducing
the number of parameters. Notably, we derive the de-biased estimates for L2
regularization for the first time. We also illustrate instances where
interpreting the Granger formalism as ’causal’ may be misleading, such as in
cases of missing variables or colliders. We offer examples showcasing the
application of this formalism to neural signals, highlighting the differences
between the VARX models and "temporal response functions" [12] often used in
neuroscience. We also show examples with physiological and sociological data.
Finally, we close with a discussion on the specifics of implementation in
code, and some caveats on the interpretation of model results.

##  2 Methods

###  2.1 VARX model

Consider the vectorial “input” signal ${\bf x}(t)$ and the vectorial “output”
signal ${\bf y}(t)$ of dimensions $d_{x}$ and $d_{y}$ respectively, with both
assumed to be observable (lower-case bold characters represent vectors). In
the present context, the input may be multiple features of a continuous
natural stimulus, say luminance and sound volume of a movie. The output could
be neural activity recorded at multiple locations in the brain. The simplest
model we can envision is one where the current signal ${\bf y}(t)$ can be
predicted linearly from the input ${\bf x}(t)$ and also linearly from the
preceding output ${\bf y}(t-1)$.

${\bf y}(t)={\bf A}*{\bf y}(t-1)+{\bf B}*{\bf x}(t)+{\bf e}(t),$ (1)

where ${\bf A}$ and ${\bf B}$ are filter matrices of dimensions
$[d_{y},d_{y}]$ and $[d_{y},d_{x}]$, with filters of length $n_{a}$ and
$n_{b}$, respectively. The additive term ${\bf e}(t)$ represents an unobserved
"innovation" that introduces an error in the prediction. In linear systems it
is called innovation because it injects novelty into the recurrent dynamic.
Here we have used a compact formulation of a multi-input multi-output
convolution, which for the auto-regressive and moving average filters reads:

$\displaystyle{\bf A}*{\bf y}(t-1)$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{n_{a}}{\bf A}(l)\cdot{\bf y}(t-l),$ (2)
$\displaystyle{\bf B}*{\bf x}(t)$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{n_{b}-1}{\bf B}(l)\cdot{\bf x}(t-l).$ (3)

Equation (1) is called a VARX model or “equation error” model [13]. It is
different from an “output error” model where the recursion has no error, but
the recursive signal ${\bf z}(t)$ is hidden and only observed with additive
noise (see Fig. 1):

$\displaystyle\begin{split}{\bf z}(t)&={\bf A}*{\bf z}(t-1)+{\bf B}*{\bf
x}(t),\\\ {\bf y}(t)&={\bf z}(t)+{\bf e}(t).\end{split}$ (4)

Note that in both cases the total response of the dynamical system (the
impulse response) can be written in the $Z$ or Fourier domains simply as

$\mathcal{H}=\mathcal{Y}\mathcal{X}^{-1}=({\bf
I}-\mathcal{A})^{-1}\mathcal{B}.$ (5)

where cursive upper-case corresponds to transform-domain variables. In this
view, what we are proposing is to model the total system response as a
combination of a Moving Average (MA) filters ${\bf B}$ and Auto-Regressive
(AR) filters $({\bf I}-{\bf A})^{-1}$. In the time domain, the total system
response can be simply computed by passing impulses in each input variable
through th systems, while setting the error/innovation to zero, ${\bf
e}(t)=0$.

Figure 1: Model schematic: In the output error model the internal state is not
observable. The gray box represents the overall system response $\mathcal{H}$.

Given observed ${\bf x}(t),{\bf y}(t)$ one can estimate the parameters ${\bf
A}$ and ${\bf B}$ by minimizing the mean of square error:

${\boldsymbol{\sigma}}^{2}=\frac{1}{T}\sum_{t=1}^{T}{\bf e}^{2}(t).$ (6)

For zero-mean signals this is also the variance, hence the conventional symbol
${\boldsymbol{\sigma}}^{2}$, which is vectorial here as it is computed and
minimized for each dimension in ${\bf y}$ individually. For the equation error
or VARX model, the linear predictors ${\bf y}(t-1)$ and ${\bf x}(t)$ are
observable, and parameter estimation results in a simple linear least-squares
problem with a well-established closed form solution. Equation (1 can be
rewritten as

${\bf Y}={\bf X}\cdot\bf H+{\bf E},$ (7)

where ${\bf X}$ is a block-Toeplitz matrix of the predictors with dimensions
$[T,N]$. Here $N$ is the total number of free parameters for each predicted
dimension in $y(t)$, i.e. $N=d_{y}n_{a}+d_{x}n_{b}$ , ${\bf Y}$ is the output
signal ${\bf y}(t)$ arranged as a matrix of dimensions $[T,d_{y}n_{a}]$, and
${\bf H}=[{\bf A},{\bf B}]^{\top}$ is a matrix of dimensions $[d_{y}n_{a},N]$
combining the AR and MA filters. The least-squares estimate is then simply:

$\hat{\bf H}={\bf R}_{xx}^{-1}{\bf R}_{xy}.$ (8)

Matrices ${\bf R}_{xx}$ and ${\bf R}_{xy}$ are block-toeplitz capturing cross-
correlations:

$\displaystyle{\bf R}_{xx}={\bf X}^{\top}\cdot{\bf X},$ (9) $\displaystyle{\bf
R}_{xy}={\bf X}^{\top}\cdot{\bf Y}.$ (10)

The estimated output is:

$\hat{\bf Y}={\bf X}*\hat{\bf H}.$ (11)

The residual errors of this model prediction for each output channel are the
diagonal elements of the correlation matrix for the errors:

$\displaystyle{\bf R}_{ee}$ $\displaystyle=$ $\displaystyle({\bf Y}-\hat{\bf
Y})^{\top}\cdot({\bf Y}-\hat{\bf Y})=\hat{\bf H}^{\top}\cdot{\bf
R}_{xx}\cdot\hat{\bf H}-2\hat{\bf H}^{\top}\cdot{\bf R}_{xy}+{\bf R}_{yy},$
(12) $\displaystyle\boldsymbol{\sigma}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{T}\mathbf{diag}({\bf R}_{ee}).$ (13)

For the output error model, the square error is not a quadratic function of
the parameters ${\bf A}$ and ${\bf B}$, and a few different iterative
optimization approaches have been proposed, such as an expectation
maximization (EM) algorithm [14], gradient back-propagation through time [15,
16], or “pseudo regression” [13]. The pros and cons of the equation-error
versus the output-error models are elaborated in the Discussion section.

Note that Equations (2)-(11) are identical in the “temporal response function”
(TRF) formalism [12], whereby matrix ${\bf X}$ only contains the input ${\bf
x}$ and with $N=d_{x}n_{b}$, and the impulse response ${\bf H}$ is then called
a multivariate TRF.

###  2.2 Granger formalism

To establish if any of the channels in filters ${\bf A}$ and ${\bf B}$
significantly improve predictions – i.e. have an "effect" – one can use a
likelihood-ratio test [17]. In this formalism, one uses deviance as the test
statistics to quantify the contribution of a given predictor in ${\bf X}$ for
each output in ${\bf Y}$. The approach consists of estimating filter
parameters ${\bf H}$ for ${\bf X}$ intact, which is referred to as the “full”
model, and then again with one of the predictors removed, which is referred to
as the “reduced” model. We compute the resulting square error
$\hat{\boldsymbol{\sigma}}_{f}^{2}$ and $\hat{\boldsymbol{\sigma}}_{r}^{2}$,
for the full and reduced models, and obtain the _deviance_ between the two
models as test statistics (there is one deviance value for each dimension in
${\bf y}$):

$\mathscrbf{D}:=T\log\left(\hat{\boldsymbol{\sigma}}^{2}_{r}/\hat{\boldsymbol{\sigma}}^{2}_{f}\right),$
(14)

where the division of the two variance vectors and the log operator are
interpreted element-wise. For normal, independent, and identically distributed
error, the vector $\mathscrbf{D}$ contains the log-likelihood ratios (times a
factor of 2), with each element following a chi-square distribution [18].
Notice that the test statistics vector $\mathscrbf{D}$ is formed by computing
the log-likelihood ratio for each output dimension, and for each predictor
dimension that is removed in the reduced model. Thus, one can estimate the
statistical significance of each channel in ${\bf A}$ and ${\bf B}$ by
computing the full model once, and then removing each predictor variable
individually from the full model. The statistical significance for a non-zero
contribution from a particular predictor to a particular output is then given
by an element of the “p-value” vector computed with the corresponding deviance
vector:

$\mathbf{p}=\mathbf{1}-F(\mathscrbf{D},n),$ (15)

Here $F$ is the cumulative distribution function for the chi-square
distribution and $n$ is the number of parameters that were removed in the
reduced model, i.e $n_{a}$ or $n_{b}$ depending if an element of ${\bf
y}(t-1)$ or ${\bf x}(t)$ was removed. The operation of $F$ on a vector is
interpreted element-wise.

###  2.3 Debiased Deviance for L2 Regularization

To avoid overtraining on small sample sizes, i.e. where $T$ is not much larger
than $N$, we decided to use an L2 penalty, with Tikhonov regularization. The
advantage over other forms of regularization, such as L1 [19, 20] or a state
space model [10, 10] is computational efficiency thanks to the closed-form
solution:

$\displaystyle\hat{\bf H}=({\bf R}_{xx}+\gamma{\bf\Gamma})^{-1}{\bf R}_{xy},$
(16)

where we selected ${\bf\Gamma}={\rm diag}({\bf R}_{xx})$ so that all variables
are regularized equally regardless of their scale. The choice of $\gamma$ is
discussed in the results section. This regularization introduces
(purposefully) a bias in the estimate, and the deviance estimate has to be
corrected to account for this bias [21]. The term that corrects the log-
likelihood in case of L2 regularization is (derived in the Appendix):

$\mathbf{b}=\frac{1}{2}\mathbf{diag}({\bf R}_{xe}^{\top}\cdot{\bf
R}^{-1}_{xx}\cdot{\bf R}_{xe})/\mathbf{diag}({\bf R}_{ee}),$ (17)

where the division between the two diagonal vectors is element-wise and ${\bf
R}_{xe}={\bf R}_{xy}-{\bf R}_{xx}\cdot\hat{\bf H}$. This bias term has to be
computed for the full and reduced models, giving $\mathbf{b}_{f}$ and
$\mathbf{b}_{r}$ respectively. The corresponding de-biased deviance is then:

$\mathscrbf{D}^{de-
biased}=T\log\left(\hat{\boldsymbol{\sigma}}^{2}_{r}/\hat{\boldsymbol{\sigma}}^{2}_{f}\right)-\mathbf{b}_{r}+\mathbf{b}_{f},$
(18)

and can be used to compute the p-values as before. We have found empirically
that we obtain a better (conservative) estimate of p-values if we use
$T^{\prime}=T-N$ instead of $T$ in this calculation of the de-biased deviance.
$T^{\prime}$ represents the effective degrees of freedom of the full model and
converges to $T$ in the asymptotic limit for which the de-biased deviance
formula was derived.

###  2.4 Effect size

Note that deviance increases linearly with $T$, that is, the statistical
evidence increases with the length of the signals. The actual benefit of
adding a link to the model is given by the log-ratio of the residual
variances, and that remains constant with increasing signal duration. This
log-ratio can be seen as a measure of effect size. Regularization can make
this log-ratio negative and this is what the bias terms correct for. In this
work, we take $\mathscrbf{D}^{de-biased}/T$ as a metric of effect size (e.g.
in Fig. 9). A traditional definition of effect size in the context of reduced
and full linear models is the coefficient of determination, or generalized
R-square [22]:

$R^{2}:=\mathbf{1}-\exp\left(-\mathscrbf{D}^{de-biased}/T\right),$ (19)

where the exponential of a vector is interpreted element-wise.

###  2.5 Basis functions for the moving average filters

The filter length (number of parameters) used in AR filters is typically kept
relatively short, to avoid over-fitting, reduce the odds of instability in the
recursion, and because even a single delay can already represent an infinite
impulse response. This is not the case for MA filters, where longer responses
have to be modeled explicitly, resulting in a relatively large number of
parameters, with a risk of over-fitting. We have found empirically that the
corrections we introduced in the Deviance estimate for short signals (Eq. 18)
do not work well when the filter lengths for ${\bf A}$ and ${\bf B}$ are very
different. A solution to both these problems (imbalance in number of
parameters and filter length) is to use basis functions for the ${\bf B}$
filters, following an approach used previously for TRFs [23, 24, 25]. In this
formalism, we have:

${\bf B}=\underline{\bf B}\circ{\bf W},$ (20)

where the inner product $\circ$ is along the lag-axis of the filter matrix
${\bf B}$, and the goal now is to find the optimal $\underline{\bf B}$. The
matrix ${\bf W}$ has dimensions $[n_{b},\underline{n}]$ so that the number of
parameters per filter is reduced from $n_{b}$ to $\underline{n}$. The linear
least-squares problem remains unchanged with the closed-form solutions using
now $\underline{\bf X}={\bf W}{\bf X}$. In the equations above this can be
implemented by replacing ${\bf R}_{xx}$ and ${\bf R}_{xy}$ with:

$\displaystyle{\bf R}_{\underline{xx}}$ $\displaystyle=$ $\displaystyle{\bf
W}\circ{\bf R}_{xx}\circ{\bf W}^{\top},$ (21) $\displaystyle{\bf
R}_{\underline{x}y}$ $\displaystyle=$ $\displaystyle{\bf W}\circ{\bf R}_{xy}.$
(22)

Note that the new ${\bf R}_{\underline{xx}}$ and ${\bf R}_{\underline{x}y}$
are no longer Toeplitz matrices. The Granger formalism applies without change.

Here we implemented Gaussian windows as shown in Figure 5A. With this, we are
not only reducing the number of parameters, i.e. regularizing the solutions,
but also selecting among a set of smooth filters ${\bf B}$.

##  3 Results

Details on all results provided next can be found in the accompanying code
repository https://github.com/lcparra/varx.

###  3.1 Test of model estimation on known model

To validate the estimation algorithm and code, we simulated a simple VARX
model with two outputs and one input ($d_{y}=2,d_{x}=1$). The algorithm
correctly recovers the AR and MA parameters (Fig.2). VARX model estimation is
available as part of the econometric toolbox in Matlab, but is limited to
instantaneous input $n_{b}=1$, i.e. no filtering of the input. When limiting
the simulation to this case, the algorithms obtain similar results. Small
variations are expected based on how the initial conditions (boundaries) are
handled and numerical differences. In our implementation we omit from the
estimate all samples that do not have a valid history. The code handles
missing values (NaN) in the same way.

Figure 2: Comparison of estimated parameters to true parameters in a simple
toy example. Here $x_{1}\rightarrow y_{2}$ and $y_{1}\rightarrow y_{2}$ were
zero in the simulated data. In this simulation signals are generated with
i.i.d innovation. We used $T=1000,\lambda=0$ and $n_{a}=3,n_{b}=1$ in both
simulation and estimation.

###  3.2 Validation of p-values

To determine if the estimation computes correct $p$-values we simulated and
VARX model with some channels set to zero, and others with non-zero values. We
did this with a small and large simulated dataset, generated with normal
i.i.d. innovations. We repeat the simulation 1000 times and determine how many
times the zero channels report a $p<0.05$. We find a false discovery rate of
approximately 0.05 for that specific null-channel, suggesting that $p$-values
are correctly estimated (Fig. 3).

Figure 3: Validation of $p$-values in both a small and a larger model
($d_{y}=6$ and $d_{y}=60$). Significance is set of p<0.05, so we expect a
false discovery rate of 0.05. Simulation here used with
$d_{x}=1,n_{a}=2,n_{b}=2,T=1000,\lambda=0.0$. Filter coefficients for ${\bf
B}$ were selected at random from unit variance normal, and ${\bf A}$ values
were set to be $\pm 0.05$ with sign selected at random (this insured stable
recursion in practice). Only two links are set to zero
$A(:,2,2)=0,B(:,5,1)=0$, for these two, the $p$-value is correctly estimated
at approximately $0.05$ (green).

###  3.3 Validation of p-values with L2 regularization

As $T$ increases the need for regularization decreases, so we scale the
regularization factor as $\gamma=\lambda/\sqrt{df}$, where $\lambda$ is a
choice of regulation in the order of 1. This particular scaling was
established empirically to accomplish two things. One, if $\gamma$ is kept
large for large $T$, then the asymptotic approximation used to compute the de-
biased deviance is no longer correct and $p$-values are miss-estimated.
Second, the specific scaling with the square root is an empirical finding that
places the optimal $\gamma$ at similar values of $\lambda$ in simulation with
different signal lengths $T$. In Fig. 4 we simulated a simple model for
varying $\lambda$ and $T$. In this example we find, as expected, that training
error increases with regularization (Fig. 4a), but test error decreases with
increasing regularization (Fig. 4b). This effect becomes clear with larger
noise (note that the explained variance here is relatively small) and with
increasing correlation in input (here random normal inputs were drawn
independently). The second observation is that with increasing sample size
$T$, test error drops, and regularization becomes less important. Finally, we
see that the bias correction provides conservative control of the false
discovery rate (Fig. 4c). The expected value in this example is $0.05$ and at
small $T$ we are well below that. Indeed, for small $T$ we seem to be over-
correcting, which limits the power to detect a true effect Fig. 4d).

Figure 4: Effect of L2 regularization on small toy examples $d_{y}=2,d_{x}=2$.
We ran the model and fit the data computing analysis of p-values in 1000
simulations with normal distributed error (zero mean, std of 2) and inputs x
(zero mean, std of 1). (a) the training error is the relative error, i.e.
$\sigma_{e}^{2}/std(y)^{2}$ computed on the same data that the model was fit
on (length $T$ samples). Error bars indicate SEM over the 1000 simulations.
(b) test error is the same relative error but computed on newly simulated data
with the estimated model. (c) False discovery is the fraction of times where
the null path is p<0.05. Solid lines are for false discovery in ${\bf A}$, and
dotted lines are for false discovery in ${\bf B}$ in this simulation. (d)
Statistical power is the fraction of times that the $p>0.05$ for the path that
has non-zero coefficients. The specific models tested here had effects of
$y_{1}\rightarrow y_{2}$ and $x_{2}\rightarrow y_{2}$ set to zero and
otherwise $n_{a}=2,n_{b}=10$ coefficients with non-zero values.

###  3.4 Example of using basis functions

For large $T$ there is minimal regularization and regularization can not be
used to “smooth” the estimated filters. For that, we developed the basis
functions decomposition using smooth basis functions. A second motivation was
also to make the number of parameters of ${\bf A}$ and ${\bf B}$ filters
comparable, as we noticed empirically that an imbalance results in quite
skewed $p$-values. We simulated examples with very different filter lengths
$n_{a}=2,n_{b}=60,n=6$ (Fig. 5). We find that introducing basis functions
improves the estimates of false discovery rate while providing smooth filter
estimates. We noticed empirically that these basis functions, for slow
signals, made the signals significantly larger in power. As a result,
conventional ridge regression regularized different signals differently,
resulting in a miss-estimation of values. This motivated the adoption of
Tikhonov regularization so that all variables are regularized by the same
amount.

Figure 5: Estimate MA filter coefficients ${\bf B}$ with and without basis
functions (basis vs free, respectively). (A) The Gaussian basis functions used
here, both for the generation of the data and for estimation. (B)
Corresponding filter estimates and true filters used during data generation.
Here $d_{y}=2$ and $d_{x}=1$, with $y_{1}\rightarrow y_{2}$ and $x\rightarrow
y_{2}$ channels set to zero. Data was simulated with $T=300$ and $n_{a}=3$,
$n_{b}=60$, $n_{b}ar=4$. The total number of free parameters was
$4*n_{a}+2*n_{b}=132$ (free) and $4*na+2*\underline{n}=20$ (basis). No
regularization was used. (C) False discovery rate for AR and MA channels with
and without basis functions. False discovery (bars) indicates how many times
in 5000 random simulations these missing links reported a $p<0.05$, therefore
we expect discovery rate to be 0.05 (dotted line).

###  3.5 Example: Intracranial recordings in human

In our view, one of the main advantages of a VARX model is the factorization
of the overall system response into an AR and MA portion, which separates the
effect of the endogenous variables into an initial response that then
reverberates in the recurrence of the dynamical system. To demonstrate this we
analyze an example of an intracranial EEG recording in a 49-year-old male
patient, included in [26]. The patient watched a total of 43.6 minutes of
movie clips while recording iEEG and eyetracking. We focus on a subset of 50
electrodes covering visual brain areas, including the occipital cortex,
fusiform face area and parahippocampal cortex. The neural data was
preprocessed as described in [26] before extracting high-frequency broadband
power (70-150Hz and downsampled to 60Hz). The exogenous "input" is a pulse
train indicating the start of fixations, i.e. the moments when the subject
receives new visual input. What is most evident is that the ${\bf B}$ response
is shorter than the full impulse response ${\bf H}$. When estimating it
directly or as ${\bf H}=(1-{\bf A})^{-1}{\bf B}$ we see that the two are very
similar. This suggests that the VARX model decomposes the total response into
a fast response followed by a prolonged response due to the recurrence in the
brain network. The VARX model appears to be therefore a valid factorization of
the full system response estimated traditionally.

Figure 6: Example of intracranial recording in humans. Here a VARX model is
fit to Broad-band high-frequency activity during free viewing of 7 videos
recorded in ($d_{y}=50$) electrodes. A total of 43.6 minutes of data were used
at a sampling rate of 60Hz ($T=156,955$) from a single patient. (A) Effect
size $R$ for the recurrent connectivity ${\bf A}$ between recording electrodes
– in the language of neuroscience, this could be called "functional
connectivity". (B) Effect size $R$ of fixation onset as exogenous variable on
different electrodes. (C) ${\bf A}$ filter coefficients ($n_{a}=4$). Diagonal
terms dominate. (D) ${\bf B}$ filter coefficients $n_{b}=30,\underline{n}=20$.
$\lambda=0$

(E) System response estimated directly – in the language of neuroscience, this
is the multivariate "temporal response function" (mTRF). (F) System response
resulting from the VARX model estimate (Eq. 5). Data from [26].

###  3.6 Example: Physiological signals in human

Another dynamical system is the body, with multiple interconnected signals
such as respiration, heart, pupil, and brain. Here we are analyzing data from
[27] which showed bidirectional effects for several of these signals during
various controlled conditions. We estimated a VARX model for this data
suggesting possible directional effects between various physiological
variables (Fig. 7). Pupil size and heart rate were measured in the experiment
as they are conventional metrics of physiological arousal. As variables are
added to the VARX model the connectivity structure is typically preserved. In
this specific example using the controlled breathing task, we initially
observe a bidirectional link between pupil size and heart rate, however this
disappears once respiration is taken into account. Instead, it is explained by
an effect of respiration on pupil size, together with the well-established
bidirectional link between respiration and heart rate that is recovered in
this data. The effect of saccades on pupil identified here is also well
established. In general, adding variables can remove links – if the addition
is a common cause, or add links – if the addition is a "collider". This is
well established for i.i.d samples [28], and is no different for temporally
correlated time sequence data.

Figure 7: Example of physiological signals in human. This data was collected
while people carried out a controlled breathing intervention. In this case
there was no exogenous stimulus, so we only fit a VAR model. Links in ${\bf
A}$ are show if $p<0.001$. Here we had 26 minutes of data from compiled across
multiple subjects samples at 25Hz ($T=26*60*25$). Data from [27].

###  3.7 Example: Sociology

As a final example, we present an analysis from the area of sociology. Here we
analyze the history of union participation, and how it relates to strikes
(Fig. 8A). The hypothesis we want to test is if strikes increase union
membership in subsequent years. We assume that the unemployment rate can not
be influenced by any of the union variables. It therefore is modeled as an
exogenous input variable, whereas the number of workers in unions and on
strike, as well as the number of strikes, can in principle all affect each
other. VARX Granger analysis suggests (Fig. 8B) suggests that unemployment
affects unionization, which in turn affects number of strikes, which obviously
affects the number of workers on strike. These results depended on the choice
of parameters $n_{a},n_{b},\lambda$. Only the effect of NumberOfStrikes
$\rightarrow$ WorkersOnStrike was robust to parameter choice. What did not
robustly emerge from this data is evidence for the initial hypothesis that
strikes lead to an increase in union participation.

Figure 8: Example on union participation and strikes: (A) Historical data from
the US. We treated the unemployment rate as an exogenous input in the VARX
model, and the others as endogenous variables. (B) Significant effects in
${\bf A}$ and ${\bf B}$ are indicated in blue and red, respectively ($p<0.05$,
$n_{a}=n_{b}=3,T=91,\lambda=0$). Note missing data around 1980, which was
omitted during the estimation, including a 3-year history.

###  3.8 The case of a missing and superfluous variable

Finally, we wanted to evaluate the case where the model does not match the
data generation process. We simulated three possible data generation processes
with one input and two outputs ($d_{x}=1,d_{y}=2$). We simulated three cases,
where the exogenous input $x$ (conditionaing variable) will be either a common
cause, a collider, or an independent variable (see Figure 9). In all cases,
the simulation implements a one-directional effect $y_{1}\rightarrow y_{2}$.
We then measure how frequently we find $p<0.05$ for this path, i.e. the power
of the test to identify a correct path, and how frequently we find $p<0.05$
for $y_{1}\leftarrow y_{2}$, i.e. the rate of false discovery. Note that
conditioning on a collider is known to introduce spurious correlations [28].
We test this for data generated with both the equation (VARX) and output error
models. The results in Figure 9 indicate the false discovery rate is correctly
estimated at 0.05 in most scenarios, i.e. we are not finding causal effects
above chance where there were none. This result holds regardless of whether
$x$ was included as input (i.e. as a control variable with instant or delayed
effect) or whether it did or did not have a true effect on $y_{1}$ and $y_{2}$
(common cause vs independent). Only when incorrectly modeling a collider as
input, did we obtain spurious effects. Statistical power was improved when
including the input to the model. In summary, there is no risk of false
discovery when including input variables, even if they don’t have a true
effect, except if they are actually affected by the internal variables $y$.

Despite using the VARX model, when the data was generated with an output error
model, the false discovery remains limited at the target of 0.05. However, a
common unobserved cause can generate spurious effects $y1\leftarrow y_{2}$. It
has been suggested that running a Granger causality model on time-reversed
data provides a control for this situation [29, 30]. We have found that
running the model on time-reversed data results in spurious effects in all
conditions tested here, so it is not clear to us how this can provide a
remedy, and the caveat of an unobserved common cause remains when we are not
directly observing endogenous variables, but only a noisy version of the
internal dynamic.

Figure 9: The VARX (equation error) model (1) requires a larger $T$ to obtain
similar power as the output error model (4). Simulation here used
$n_{a}=n_{b}=3,T=5000$ and normal i.i.d. error.

##  4 Discussion

### Equation error versus output error models

The estimation of model parameters for a VARX model has a closed-form
solution, thus being much faster than finding parameters for an output error
model, which requires iterative algorithms. The gain in computational
efficiency results from the assumption that $y(t)$ is observable. This may not
be a good assumption in the case of brain signals measured across the skull,
such as EEG/MEG where only a linear mixture, possibly with added noise, is
observed. In that case, iterative algorithms are needed, but the Granger
formalism can still be used with some effort [20]. In the VARX model, however,
we do not need to assume that all internal activity is directly observable.
Any unobserved activity is captured as innovation $e(t)$. We only need to be
aware that any recurrent connectivity may be due to those unobserved common
"causes". In particular symmetric effect sizes $R$ will be suggestive of such
a missing variable. The role of the error is quite different in the two
models. In the VARX (equation error) the error is an internal source of
innovation driving the recurrent dynamic similar to the drive that comes from
the input. The internal states are fully observable. In the output error
model, the input entirely drives the system, and the error only affects the
observations and is not injected into the dynamic.

#### Sensitivity to parameters

A caveat to all results above is that individual links can be sensitive to the
model assumptions, namely, which variables are selected as endogenous (and can
be affected by all others), and which variables are selected as exogenous (and
cannot be affected). An example of that was the choice of the unemployment
rate as exogenous to the dynamic of unions. The results can also depend on
which endogenous variables are included, as we saw in the example of
physiological signals. Results can also depend on the number of parameters
$n_{a}$ and $n_{b}$ and regularization factor $\lambda$ (we saw this in the
Union example). Although we did demonstrate this here, these parameters could
be established with cross-validation.

#### Causality

In Granger’s original formalism [1], the error of the full and reduced model
refers to one-dimensional signals where $y(t-1)$ is used in both cases and
$x(t)$ is either used or omitted. If the error is significantly reduced by
including $x(t)$ in the model, Granger argues that $x$ "causes" $y$. This
interpretation is inaccurate for several reasons [31]. As we saw, when common
causes are not observed (either as external input or internal variables) they
can generate spurious links [32]. Similarly, including colliders can cause
spurious links. All this is well explained by Pearl’s approach to causal
inference [28]. Therefore one should not think of the Granger formalism as
serious evidence for a causal graph without a well-justified prior graphical
model [33]. In particular for large dimensional datasets such as brain data,
where we only observe a tiny fraction of all the variables, the risk of
unobserved common causes is much too large to take the resulting graph
seriously as a causal graph. Nevertheless, asymmetries in the ${\bf A}$ matrix
can be seen as evidence of temporal precedence suggestive of an asymmetric
"information flow".

#### Non-stationarity

Deviance makes a statistical judgment for the entire channel, not individual
delays, as is common we simply treat each delay as a new predictor with its
statistical test (this is the approach of the Matlab VARX). There are multiple
methods under the umbrella of "Granger causality" that attempt to decide on
how many tabs or which delays to use. By collapsing statistical evidence into
a single test statistic, Deviance, this approach has greater statistical
power. This is reflected in the linearity of D with the number of samples T.
The flip side is that this statistic is very sensitive to violations of its
assumptions. For instance, it assumes that all $T$ samples of the innovation
process are independent and identically distributed. The AR portion of the
model assures that the linear-fit residual errors are uncorrelated in time,
however, if there is any non-stationarity, this will no longer be the case.
Therefore, non-stationarity will cause spurious correlations [34]. In
particular, any transient will cause larger deflection and correlation across
samples. In particular, transient that affect several signals, say a common
edge at the start or end of the signal may appear to behave like a common
drive with high amplitude that results in a spurious link. Therefore, in the
present approach, one has to treat edges and transients with utmost care to
avoid spurious links.

Some have argued that issues with non-linearity and non-stationarity can be
addressed [35]. Barnett et al. proposed a State Space model that can cope to
some degree with missing variables, does not need to compute a reduced model,
and can deal with non-linearity and non-stationarity [35]. However, Stokes and
Purdon showed that even the state-space Granger is not immune to confounding
effects, non-stationarities, etc. The topic remains a matter of debate [36].

An alternative is to avoid using analytic expressions for the $p$-value, Eq.
(15), and instead use standard non-parametric statistics. For time series, the
simplest is to randomly time-delay channels relative to one-another,
potentially with a circular wraparound. All else in the model identification,
i.e. estimates ${\bf A}$, ${\bf B}$ and effect size $R^{2}$ remain valid
estimates of linear predictions even in the presence of non-stationarity and
non-linearity.

#### Toolboxes

There are several software tools to estimate VARX models. We are not aware of
one that provides the Granger-Geweke test for significance. The MVGC toolbox
[10] does not report results on exogenous variables,
https://github.com/SacklerCentre/MVGC1. Implementations of VARX models in
Matlab Econometrics Toolbox and in SAS for instance, make significance
statements for individual delays but do not allow for the exogenous variable,
i.e. $nb=1$. We welcome feedback if there other tools that that should be
mentioned here. Our code emphasized computational efficiency to handle large
datasets, with comparatively long ${\bf B}$ filters.

#### Stability

A word about ${\bf A}$ is in order. The AR filter $1/(1-{\bf A})$ can be
unstable. We have not implemented any mechanism for this vectorial AR filter
to remain stable. Lack of stability only manifests in the systems when
computing the overall systems response $\mathcal{\bf H}$, which is not
necessary during estimation for the ${\bf B},{\bf A}$ nor the calculation of
statistical significance of each path (contrary to the output error model,
where the recurrence has to be run back in time to estimated gradients,
risking issues of stability (is this true?). We rarely encountered unstable AR
estimates, and where we did, L2 regularization addressed the issues. But
again, there is nothing in our formalism to ensure the stability of
$\mathcal{\bf H}$.

### 1/f spectrum

Another word about ${\bf A}$. The diagonal elements of ${\bf A}$ in practice
will always be high-pass filters, as we saw in the example of intra-cranial
recordings. We advise not to take individual delays in the diagonal terms
literally. The reason for this is that the innovation process is assumed to be
white (constant spectrum), whereas all natural signals tend to have a 1/f
spectrum. As a result, $1/(1-{\bf A})$ has to have a 1/f spectrum, and ${\bf
A}$ has to scale with f, i.e. be high-pass. In practice, we find that this is
entirely accomplished by the diagonal elements of ${\bf A}$. But the caveat in
principle applies also to the off-diagonal elements. Future work could
consider a VARMAX model where the innovation is first filtered and then
injected into the recurrent dynamic [13]. However, estimation of VARMAX model
parameters is a non-convex optimization problem with similar complications to
the output-error model.

##  5 Conclusion

The "multivariate temporal response function" (TRF) approach for modeling the
response of dynamical systems to exogenous variables essentially boils down to
a vectorial MA model. And, unlike the VAR model, TRF is not commonly examined
in the Granger formalism. Although Geweke introduced exogenous variables into
the Granger formalism [8], the clear connection to the VARX model is not yet
widely recognized. We hope to have bridged this gap. Different models vary in
their assumptions about how to break down the system’s overall response. When
estimating with the VARX model, the data itself dictates how to factor the
overall response into MA versus AR dynamics.

##  6 Acknowledgement

We would like to thank Vera E. Parra<EMAIL_ADDRESS>for posing the
question on union participation and compiling the corresponding datasets.

##  7 Appendix: Derivation of the De-biased Deviance

Consider the linear model in Eq. (7) for $n_{y}=1$:

$\mathbf{y}=\mathbf{X}\cdot\mathbf{h}+\mathbf{e},$ (23)

for simplicity of exposition. The case of $n_{y}>1$ is a natural
generalization, as we will explain below. A ridge regressor for $\mathbf{h}$
is obtained by:

$\hat{\mathbf{h}}:=\underset{\mathbf{h}}{\operatorname{argmin}}\quad\|\mathbf{y}-\mathbf{X}\cdot\mathbf{h}\|^{2}+\gamma\|\mathbf{h}\|^{2},$
(24)

where $\gamma$ is the ridge parameter. To perform de-biasing, stationarity
conditions for $\hat{\mathbf{h}}$ give:

$-\mathbf{X}^{\top}(\mathbf{y}-\mathbf{X}\hat{\mathbf{h}})+\lambda\hat{\mathbf{h}}=\mathbf{0}.$
(25)

Following the arguments of Eqs. (3)–(5) in [37], the following estimator is
the de-biased estimator:

$\hat{\mathbf{h}}^{de-
biased}=\hat{\mathbf{h}}+(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\left(\mathbf{y}-\mathbf{X}\hat{\mathbf{h}}\right)=\hat{\mathbf{h}}+\mathbf{R}_{xx}^{-1}\mathbf{r}_{xe}.$
(26)

Now, consider the GC inference setting where we have full and reduce models
given by:

$\textrm{Full
model:}\quad\mathbf{y}=\mathbf{X}^{f}\mathbf{h}^{f}+\mathbf{e}^{f},\quad\textrm{Reduced
model:}\quad\mathbf{y}=\mathbf{X}^{r}\mathbf{h}^{r}+\mathbf{e}^{r}.$ (27)

The log-likelihoods of the full and reduced models are given by:

$\ell^{f}:=-\frac{T}{2}\ln(2\pi)-\frac{T}{2}\ln\sigma^{2}_{f}-\frac{1}{2\sigma_{f}^{2}}\left\|\mathbf{y}-\mathbf{X}^{f}\mathbf{h}^{f}\right\|^{2},$
(28)
$\ell^{r}:=-\frac{T}{2}\ln(2\pi)-\frac{T}{2}\ln\sigma^{2}_{r}-\frac{1}{2\sigma_{r}^{2}}\left\|\mathbf{y}-\mathbf{X}^{r}\mathbf{h}^{r}\right\|^{2},$
(29)

and the deviance difference is given by
$\mathscr{D}:=2\big{(}\hat{\ell}^{f}-\hat{\ell}^{r}\big{)}$, where
$\hat{\ell}^{f}$ and $\hat{\ell}^{r}$ are, respectively, the full and reduced
log-likelihoods evaluated at the ridge estimates of full and reduced
parameters:

$\displaystyle\hat{\mathbf{h}}^{f}$
$\displaystyle:=\underset{\mathbf{h}}{\operatorname{argmin}}\quad\|\mathbf{y}-\mathbf{X}^{f}\mathbf{h}\|^{2}+\lambda\|\mathbf{h}\|^{2},$
(30) $\displaystyle\hat{\mathbf{h}}^{r}$
$\displaystyle:=\underset{\mathbf{h}}{\operatorname{argmin}}\quad\|\mathbf{y}-\mathbf{X}^{r}\mathbf{h}\|^{2}+\lambda\|\mathbf{h}\|^{2}.$
(31)

If the log-likelihoods were evaluated at the ML estimators, then $\mathscr{D}$
would be asymptotically chi-square/non-central chi-square, allowing to perform
precise statistical tests. But, the ridge estimator is different than the ML
estimator. So, we need to de-bias the deviance difference.

Following Eq. (20) and Eqs. (B.2)–(B.6) in [20], the de-biased deviance
difference is given by:

$\mathscr{D}^{de-
biased}=2\big{(}\hat{\ell}^{f}-\hat{\ell}^{r}\big{)}-B_{r}+B_{f},$ (32)

where the bias function $B(\cdot)$ is defined as:

$\displaystyle b_{f}$
$\displaystyle:=\frac{1}{2\hat{\sigma}^{2}_{f}}\big{(}\mathbf{y}-\mathbf{X}^{f}\hat{\mathbf{h}}^{f}\big{)}^{\top}\mathbf{X}^{f}\left(\mathbf{X}^{f\top}\mathbf{X}^{f}\right)^{-1}\mathbf{X}^{f\top}\big{(}\mathbf{y}-\mathbf{X}^{f}\hat{\mathbf{h}}^{f}\big{)}=\frac{1}{2r^{f}_{ee}}\mathbf{r}_{xe}^{f\top}\mathbf{R}_{xx}^{f\,\,-1}\mathbf{r}^{f}_{xe},$
(33) $\displaystyle b_{r}$
$\displaystyle:=\frac{1}{2\hat{\sigma}^{2}_{r}}\big{(}\mathbf{y}-\mathbf{X}^{r}\hat{\mathbf{h}}^{r}\big{)}^{\top}\mathbf{X}^{r}\left(\mathbf{X}^{r\top}\mathbf{X}^{r}\right)^{-1}\mathbf{X}^{r\top}\big{(}\mathbf{y}-\mathbf{X}^{r}\hat{\mathbf{h}}^{r}\big{)}=\frac{1}{2r^{r}_{ee}}\mathbf{r}_{xe}^{r\top}\mathbf{R}_{xx}^{r\,\,-1}\mathbf{r}^{r}_{xe}.$
(34)

Simplifying the deviance difference $\mathscr{D}$ for the linear model, we can
express $\mathscr{D}^{de-biased}$ as:

$\mathscr{D}^{de-
biased}=T\ln\frac{\hat{\sigma}^{2}_{r}}{\hat{\sigma}^{2}_{f}}-b_{r}+b_{f},$
(35)

where

$\hat{\sigma}_{f}^{2}:=\frac{1}{T}\left\|\mathbf{y}-\mathbf{X}^{f}\hat{\mathbf{h}}^{f}\right\|^{2}=r_{ee}^{f},\qquad\hat{\sigma}_{r}^{2}:=\frac{1}{T}\left\|\mathbf{y}-\mathbf{X}^{r}\hat{\mathbf{h}}^{r}\right\|^{2}=r_{ee}^{r},$
(36)

are the full and reduced prediction variances. For $n_{y}>1$, given that each
channel of $\mathbf{Y}$ is treated separately, the full and reduced bias
vectors are given by:

$\mathbf{b}_{f}=\frac{1}{2}\mathbf{diag}({\bf R}_{xe}^{f\top}\cdot{\bf
R}^{f\,\,-1}_{xx}\cdot{\bf R}_{xe}^{f})/\mathbf{diag}({\bf
R}_{ee}^{f}),\qquad\mathbf{b}_{r}=\frac{1}{2}\mathbf{diag}({\bf
R}_{xe}^{r\top}\cdot{\bf R}^{r\,\,-1}_{xx}\cdot{\bf
R}_{xe}^{r})/\mathbf{diag}({\bf R}_{ee}^{r}),$ (37)

and the de-biased deviance vector is given by:

$\mathscrbf{D}^{de-
biased}=T\log\left(\hat{\boldsymbol{\sigma}}^{2}_{r}/\hat{\boldsymbol{\sigma}}^{2}_{f}\right)-\mathbf{b}_{r}+\mathbf{b}_{f}.$
(38)

Under the null hypothesis, the de-biased deviance difference can be shown to
be asymptotically, as $T\rightarrow\infty$, chi-square distributed with $n$
degrees of freedom, where $n$ is the number of parameters removed in the
reduced model and $T$ is the number of observations (i.e., column dimension of
$\mathbf{Y}$) [18]. This asymptotic property requires that ridge estimator is
consistent. To guarantee this, $\gamma$ must grow slower than
$\mathcal{O}(T)$. In practice, fixing $\gamma$ using cross-validation
typically guarantees this.

## References

  * [1] Clive WJ Granger.  Investigating causal relations by econometric models and cross-spectral methods.  Econometrica: journal of the Econometric Society, pages 424–438, 1969.
  * [2] Norbert Wiener.  The Theory of Prediction.  In Edwin F. Beckenbach, editor, Modern Mathematics of the Engineer: First series, pages 165–190. Dover Publications, Inc., 1956.
  * [3] Steven L Bressler and Anil K Seth.  Wiener–granger causality: a well established methodology.  Neuroimage, 58(2):323–329, 2011.
  * [4] Ali Shojaie and Emily B Fox.  Granger causality: A review and recent advances.  Annual Review of Statistics and Its Application, 9:289–319, 2022.
  * [5] Anil K Seth, Adam B Barrett, and Lionel Barnett.  Granger causality analysis in neuroscience and neuroimaging.  Journal of Neuroscience, 35(8):3293–3297, 2015.
  * [6] Duo Qin.  Rise of var modelling approach.  Journal of Economic Surveys, 25(1):156–174, 2011.
  * [7] John R Freeman.  Granger causality and the times series analysis of political relationships.  American Journal of Political Science, pages 327–358, 1983.
  * [8] John F Geweke.  Measures of conditional linear dependence and feedback between time series.  Journal of the American Statistical Association, 79(388):907–915, 1984.
  * [9] Shuixia Guo, Anil K Seth, Keith M Kendrick, Cong Zhou, and Jianfeng Feng.  Partial granger causality—eliminating exogenous inputs and latent variables.  Journal of neuroscience methods, 172(1):79–93, 2008.
  * [10] Lionel Barnett and Anil K Seth.  The MVGC multivariate granger causality toolbox: a new approach to granger-causal inference.  J. Neurosci. Methods, 223:50–68, February 2014.
  * [11] Mingzhou Ding, Yonghong Chen, and Steven L Bressler.  Granger causality: basic theory and application to neuroscience.  Handbook of time series analysis: recent theoretical developments and applications, pages 437–460, 2006.
  * [12] Michael J Crosse, Giovanni M Di Liberto, Adam Bednar, and Edmund C Lalor.  The multivariate temporal response function (mTRF) toolbox: A MATLAB toolbox for relating neural signals to continuous stimuli.  Front. Hum. Neurosci., 10:604, November 2016.
  * [13] L Ljung.  System identification-theory for the user 2nd edition ptr prentice-hall.  Upper Saddle River, NJ, 1999.
  * [14] R H Shumway and D S Stoffer.  An approach to time series smoothing and forecasting using the em algorithm.  J. Time Ser. Anal., 3(4):253–264, July 1982.
  * [15] J J Shynk.  Adaptive IIR filtering.  IEEE ASSP Mag., 6(2):4–21, April 1989.
  * [16] Carles Igual, Jorge Igual, Janne M Hahne, and Lucas C Parra.  Adaptive auto-regressive proportional myoelectric control.  IEEE Transactions on Neural Systems and Rehabilitation Engineering, 27(2):314–322, 2019.
  * [17] John Geweke.  Measurement of linear dependence and feedback between multiple time series.  J. Am. Stat. Assoc., 77(378):304–313, June 1982.
  * [18] R R Davidson and W E Lever.  The limiting distribution of the likelihood ratio statistic under a class of local alternatives.  Sankhyā Indian J. Stat. Ser, 32:209–224, 1961.
  * [19] Proloy Das and Behtash Babadi.  Non-asymptotic guarantees for reliable identification of granger causality via the LASSO.  IEEE Trans. Inf. Theory, 69(11):7439–7460, November 2023.
  * [20] Behrad Soleimani, Proloy Das, I M Dushyanthi Karunathilake, Stefanie E Kuchinsky, Jonathan Z Simon, and Behtash Babadi.  NLGC: Network localized granger causality with application to MEG directional functional connectivity analysis.  Neuroimage, 260(119496):119496, October 2022.
  * [21] Alireza Sheikhattar, Sina Miran, Ji Liu, Jonathan B Fritz, Shihab A Shamma, Patrick O Kanold, and Behtash Babadi.  Extracting neuronal functional network dynamics via adaptive granger causality analysis.  Proc. Natl. Acad. Sci. U. S. A., 115(17):E3869–E3878, April 2018.
  * [22] Lonnie Magee.  R2Measures based on wald and likelihood ratio joint significance tests.  Am. Stat., 44(3):250–253, August 1990.
  * [23] Sahar Akram, Jonathan Z Simon, and Behtash Babadi.  Dynamic estimation of the auditory temporal response function from MEG in competing-speaker environments.  IEEE Trans. Biomed. Eng., 64(8):1896–1905, August 2017.
  * [24] Sina Miran, Sahar Akram, Alireza Sheikhattar, Jonathan Z Simon, Tao Zhang, and Behtash Babadi.  Real-time tracking of selective auditory attention from M/EEG: A bayesian filtering approach.  Front. Neurosci., 12, May 2018.
  * [25] I M Dushyanthi Karunathilake, Jason L Dunlap, Janani Perera, Alessandro Presacco, Lien Decruy, Samira Anderson, Stefanie E Kuchinsky, and Jonathan Z Simon.  Effects of aging on cortical representations of continuous speech.  J. Neurophysiol., 129(6):1359–1377, June 2023.
  * [26] Maximilian Nentwich, Marcin Leszczynski, Brian E Russ, Lukas Hirsch, Noah Markowitz, Kaustubh Sapru, Charles E Schroeder, Ashesh D Mehta, Stephan Bickel, and Lucas C Parra.  Semantic novelty modulates neural responses to visual change across the human brain.  Nat. Commun., 14(1):2910, May 2023.
  * [27] Jens Madsen and Lucas C Parra.  Bidirectional brain-body interactions during natural story listening.  Cell Reports, 43(4), 2024.
  * [28] Judea Pearl.  Linear models: A useful “microscope” for causal analysis.  J. Causal Inference, 1(1):155–170, May 2013.
  * [29] Martin Vinck, Lisanne Huurdeman, Conrado A Bosman, Pascal Fries, Francesco P Battaglia, Cyriel M A Pennartz, and Paul H Tiesinga.  How to detect the granger-causal flow direction in the presence of additive noise?  Neuroimage, 108:301–318, March 2015.
  * [30] Irene Winkler, Danny Panknin, Daniel Bartz, Klaus-Robert Muller, and Stefan Haufe.  Validity of time reversal for testing granger causality.  IEEE Trans. Signal Process., 64(11):2746–2760, June 2016.
  * [31] Mariusz Maziarz.  A review of the granger-causality fallacy.  The journal of philosophical economics: Reflections on economic and social issues, 8(2):86–105, 2015.
  * [32] Helmut Lütkepohl.  Non-causality due to omitted variables.  Journal of econometrics, 19(2-3):367–378, 1982.
  * [33] J. Pearl.  Causality.  Causality: Models, Reasoning, and Inference. Cambridge University Press, 2009.
  * [34] Clive WJ Granger and Paul Newbold.  Spurious regressions in econometrics.  Journal of econometrics, 2(2):111–120, 1974.
  * [35] Lionel Barnett and Anil K Seth.  Granger causality for state-space models.  Phys. Rev. E Stat. Nonlin. Soft Matter Phys., 91(4):040101, April 2015.
  * [36] Lionel Barnett, Adam B Barrett, and Anil K Seth.  Misunderstandings regarding the application of granger causality in neuroscience.  Proceedings of the National Academy of Sciences, 115(29):E6676–E6677, 2018.
  * [37] Sara van de Geer, Peter Bühlmann, Ya’acov Ritov, and Ruben Dezeure.  On asymptotically optimal confidence regions and tests for high-dimensional models.  The Annals of Statistics, 42(3):1166–1202, June 2014.  Publisher: Institute of Mathematical Statistics.