(i) $S_{C}=\\{u_{2}=u_{1}\\}$, (ii) $S_{U}=\\{g(u_{1})=u_{2}\\}$, and (iii)
$S_{f}=\\{f(a,c)=a,f(b,c)=u_{1},f(c,u_{2})=b\\}$.

Different total orderings on constants are used to illustrate how different
canonical forms can be generated. Consider a total ordering $f\gg g\gg a\gg
b\gg u_{2}\gg u_{1}$. $R_{C}=\\{1.~{}u_{2}\rightarrow u_{1}\\}$ normalizes the
uninterpreted equation and it is oriented as:
$R_{U}=\\{2.~{}g(u_{1})\rightarrow u_{1}\\}$.

To generate a reduced canonical rewrite system for AC $f$-terms, an admissible
ordering on $f$-terms must be chosen. The degree-lexicographic ordering on
monomials will be used for simplicity: $f(M_{1})\gg_{f}f(M_{2})$ iff (i)
$|M_{1}|>|M_{2}|$ or (ii) $|M_{1}|=|M_{2}|\land M_{1}-M_{2}$ includes a
constant bigger than every constant in $M_{2}-M_{1}$. $f$-equations are
normalized using rules 1 and 2, and are oriented as:
$\\{3.~{}f(a,c)\rightarrow a,~{}4.~{}f(c,u_{1})\rightarrow
b,~{}5.~{}f(b,c)\rightarrow u_{1}\\}$.

Applying the SingleACCompletion algorithm, the superposition between rules
$3,5$ is $f(a,b,c)$ with the critical pair: $(f(a,b),f(a,u_{1}))$, leading to
a rewrite rule $6.~{}f(a,b)\rightarrow f(a,u_{1})$; the superposition between
the rules $3,4$ is $f(a,c,u_{1})$ with the critical pair:
$(f(a,u_{1}),f(a,b))$ which is trivial by rule 6. The superposition between
the rules $4,5$ is $f(b,c,u_{1})$ with the critical pair:
$(f(b,b),f(u_{1},u_{1}))$ giving: $7.~{}f(b,b)\rightarrow f(u_{1},u_{1})$. The
rewrite system $R_{f}=\\{3,4,5,6,7\\}$ is a reduced canonical rewrite system
for $S_{f}$. $R_{S}=\\{1,2\\}\cup R_{f}$ is a reduced canonical rewrite system
associated with the AC congruence closure of the input and serves as its
decision procedure.

In the original signature, the above rewrite system $R_{S}$ is:
$\\{g(f(b,c)))\rightarrow f(b,c),~{}f(a,c)\rightarrow a,~{}f(b,c,c)\rightarrow
b,~{}f(a,b)\rightarrow f(a,b,c),~{}f(b,b)\rightarrow f(b,b,c,c)\\}$ with 5
being trivial. The reader would observe this rewrite system is locally
confluent but not terminating.

Consider a different ordering: $f\gg g\gg a\gg b\gg u_{1}\gg u_{2}$.

This gives rise to a related rewrite system in which $u_{1}$ is rewritten to
$u_{2}$: $\\{u_{1}\rightarrow u_{2},g(u_{2})\rightarrow
u_{2},f(a,c)\rightarrow a,f(c,u_{2})\rightarrow b,f(b,c)\rightarrow
u_{2},f(a,b)\rightarrow f(a,u_{2}),f(b,b)\rightarrow f(u_{2},u_{2})\\}.$ In
the original signature, the rewrite system is:
$R^{\prime}_{S}=\\{f(b,c)\rightarrow g(f(b,c)),g(g(f(b,c)))\rightarrow
g(f(b,c)),f(a,c)\rightarrow a,f(c,g(f(b,c)))\rightarrow b,f(a,b)\rightarrow
f(a,g(f(b,c))),f(b,b)\rightarrow f(g(f(b,c)),g(f(b,c)))\\}$. The reader should
compare $R^{\prime}_{S}$ which seems more complex using complicated terms,
with the above system $R_{S}$ with different orientation between
$u_{1},u_{2}.$

Example 10: This example illustrates interaction between congruence closures
over uninterpreted symbols and AC symbols. Let
$S=\\{g(b)=a,g(d)=c,a*c=c,b*c=b,a*b=d\\}$ where $g$ is uninterpreted and $*$
is AC. Let $*\gg g\gg a\gg b\gg c\gg d$. Applying the steps of the algorithm,
$R_{U}$, the congruence closure over uninterpreted symbols, is
$\\{g(b)\rightarrow a,g(d)\rightarrow c\\}$, Completion on the $*$-equations
using degree lexicographic ordering, oriented as $\\{a*c\rightarrow
c,b*c\rightarrow b,a*b\rightarrow d\\}$, generates an implied constant
equality $b=d$ from the critical pair of $a*c\rightarrow c,b*c\rightarrow b$.
Using the rewrite rule $b\rightarrow d$, the AC rewrite system reduces to:
$R_{*}=\\{a*c\rightarrow c,c*d\rightarrow d,a*d\rightarrow d\\}$, which is
canonical.

The implied constant rule $b\rightarrow d$ is added to $R_{C}:\\{b\rightarrow
d\\}$ and is also propagated to $R_{U}$, which makes the left sides of
$g(b)\rightarrow a$ and $g(d)\rightarrow c$ equal, generating another implied
constant rule $a\rightarrow c$ which is added to $R_{C}$: $\\{a\rightarrow
c,~{}b\rightarrow d\\}$. $R_{U}$ becomes $\\{g(d)\rightarrow c\\}.$ The
implied constant rule is propagated to the AC rewrite system on $*$. $R_{C}$
normalizes $R_{*}$ to $\\{c*c\rightarrow c,~{}c*d\rightarrow d\\}.$

The output of the algorithm is a canonical rewrite system: $\\{g(d)\rightarrow
c,~{}b\rightarrow d,~{}a\rightarrow c,~{}c*c\rightarrow c,~{}c*d\rightarrow
d\\}$. In general, the propagation of equalities can result in the left sides
of the rules in AC subsystems change, generating new superpositions.

Example 11: Consider an example with multiple AC symbols and a single
uninterpreted symbol from [NR91]: $S=\\{a+b=a*b,a*c=g(d),d=d^{\prime}\\}$ with
AC symbols $+,*$ and an uninterpreted symbol $g$. This example also illustrate
flexibility in choosing orderings in the proposed framework.

Purification leads to introduction of new constants :
$\\{1.~{}a+b=u_{0},~{}2.~{}a*b=u_{1},~{}3.~{}a*c=u_{2},4.~{}g(d)=u_{2},5.~{}d=d^{\prime},6.~{}u_{0}=u_{1}\\}$.
Depending upon a total ordering on constants $a,b,c,d,d^{\prime}$, there are
thus many possibilities depending upon the desired canonical forms.

Recall that $a+b$ cannot be compared with $a*b$, however $u_{0}$ and $u_{1}$
can be compared. If canonical forms are desired so that the canonical form of
$a+b$ is $a*b$, the ordering should include $u_{0}~{}>~{}u_{1}$. Consider a
total ordering $a\gg b\gg c\gg d\gg d^{\prime}\gg u_{2}\gg u_{0}\gg u_{1}$.
Degree-lexicographic ordering is used on $+$ and $*$ monomials. This gives
rise to: $R_{C}=\\{6.~{}u_{0}\rightarrow u_{1},5.~{}d\rightarrow
d^{\prime}\\},~{}R_{U}=\\{4^{\prime}.~{}g(d^{\prime})\rightarrow u_{2}\\}$
along with $R_{+}=\\{1.~{}a+b\rightarrow u_{1}\\}$ and
$R_{*}=\\{2.~{}a*b\rightarrow u_{1},3.~{}a*c\rightarrow u_{2}\\}$. The
canonical rewrite system for $*$ is: $\\{2.~{}a*b\rightarrow
u_{1},3.~{}a*c\rightarrow u_{2},7.~{}b*u_{2}\rightarrow c*u_{1}\\}.$

The rewrite system $R_{S}=\\{1.~{}a+b\rightarrow u_{1},~{}2.~{}a*b\rightarrow
u_{1},~{}3.~{}a*c\rightarrow u_{2},~{}4^{\prime}.~{}g(d^{\prime})\rightarrow
u_{2},~{}5.~{}d\rightarrow d^{\prime},~{}6.~{}u_{0}\rightarrow
u_{1},~{}7.~{}b*u_{2}\rightarrow c*u_{1}\\}$ is canonical. Both $a+b$ and
$a*b$ have the same normal form $u_{1}$ standing for $a*b$ in the original
signature.

With $u_{1}\gg u_{0}$, another canonical rewrite system
$R_{S}=\\{1.~{}a+b\rightarrow u_{0},~{}2.~{}a*b\rightarrow
u_{0},~{}6^{\prime}.~{}u_{1}\rightarrow u_{0},~{}3.~{}a*c\rightarrow
u_{2},~{}4.~{}g(d^{\prime})\rightarrow u_{2},~{}5.~{}d\rightarrow
d^{\prime},~{}7.~{}b*u_{2}\rightarrow c*u_{0}\\}$ in which $a+b$ and $a*b$
have the normal form $u_{0}$ standing for $a+b$, different from the one above.

##  11\. A Gröbner Basis Algorithm as an AC Congruence Closure

Buchberger’s Gröbner basis algorithm when extended to polynomial ideals over
integers [KK88] can be interestingly viewed as a special congruence closure
algorithm, extending the congruence closure algorithm over multiple AC symbols
$+$ and $*$ which in addition, satisfy the properties of a commutative ring
with unit, especially the distributivity property. $+$ satisfies the
properties of an Abelian group with $0$ as its identity; $*$ is AC with
identity 1 and additionally with the distributivity property: $\forall
x,y,z,~{}x*(y+z)=(x*y)+(x*z)$; $x*0=0$ as well. To incorporate the
distributivity property in the combination, additional superpositions and
associated critical pairs must be considered. Relationship between Gröbner
basis algorithm and the Knuth-Bendix completion procedure has been
investigated in [KK83, Win84, KKW89, Mar96], but the proposed insight is
novel. Below, we illustrate this new insight using an example from [KK88].

Abusing the notation, $x+-(y)$ will be written as $x-y$. An additive monomial
$c~{}t$, where $c\neq 0$, is an abbreviation of repeating the monomial $t$,
$c$ times; if $c$ is positive, say $3$, then $3~{}t$ is an abbreviation for
$t+t+t$; if $c$ is negative, say $-2$, then $-2~{}t$ is an abbreviation for
$-t-t$. A $*$-monomial with the coefficient 1 is a pure term expressed in $*$,
but a monomial $3~{}y*y*y$ is a mixed term $y*y*y+y*y*y+y*y*y$ with $+$ as the
outermost symbol which has $*$ monomials as subterms.

It can be easily proved that a ground term in $+,*$ and constants can be
brought into a polynomial form using the rules of a commutative ring with
unity and by collecting identical $*$ monomials; for example, $y*(y+y)+y*y$ is
equivalent to its canonical form $3~{}y*y*y$. Similarly, a ground equation
built using these function symbols can be transformed to a polynomial
equation.

Associated with a finite set $S$ of ground equations in $+,*$ and constants,
is a basis $B$ of polynomial equations over the constants. It can be proved
that the congruence closure $CR1CC(S)$ of $S$ using the properties of a
commutative ring with 1 corresponds to the ideal $I(B)$ generated by the
polynomials in the basis $B$; the membership test in $CR1CC(S)$ is precisely
the membership test in the associated ideal $I(B)$.

Because of the distributivity property, a mixed term of the form
$m*(t_{1}+t_{2})$ can possibly be rewritten by the rules in a rewrite system
of $*$-monomials as well as by the rules in a rewrite system of $+$ monomials;
further it can be expanded using the distributivity property as well to be
$m*t_{1}+m*t_{2}$ which can be further rewritten. This is captured below using
new superpositions and the associated critical pairs due to the distributivity
property. This is first illustrated using an example.

Example 12: Consider an example [KK88]. The input basis is:
$B=\\{7~{}x*x*y=3~{}x,4~{}x*y*y=x*y,3~{}y*y*y=0\\}$ of a polynomial ideal over
the integers [KK88]. Ground terms in $B$ can be written in many different
ways, resulting in different sets of ground equations on terms built using
$+,0,-,*,1$; the above basis $B$ is one such instance.

None of these equations is pure in $+$ or $*$; in other words, they are mixed
terms both in $*$ and $+$. Purification of the above equations using new
constants leads to:
$\\{1.~{}7~{}u_{1}=3~{}x,2.~{}4~{}u_{2}=u_{3},3.~{}3~{}u_{4}=0$ with
$4.~{}x*x*y=u_{1},5.~{}x*y*y=u_{2},6.~{}x*y=u_{3},7.~{}y*y*y=u_{4}\\}.$ These
equations with new constant symbols are now pure with $1,2,3$ purely in $+$
and $4,5,6,7$ purely in $*$. Let a total ordering on all constants be:
$u_{1}\gg u_{2}\gg u_{4}\gg u_{3}\gg y\gg x$. Extend it using the degree-
lexicographic ordering on $+$-monomials as well as $*$-monomials, Orienting
$*$ equations: $R_{*}=\\{4.~{}x*x*y\rightarrow u_{1},5.~{}x*y*y\rightarrow
u_{2},6.~{}x*y\rightarrow u_{3},7.~{}y*y*y\rightarrow u_{4}\\}$. Orienting $+$
equations, $R_{+}=\\{1.~{}7~{}u_{1}\rightarrow 3~{}x,2.~{}4~{}u_{2}\rightarrow
u_{3},3.~{}3~{}u_{4}\rightarrow 0\\}$.

SingleACCompletion on the ground equations on $*$ terms generates a reduced
canonical rewrite system: $R_{*}=\\{6.~{}x*y\rightarrow
u_{3},7.~{}y*y*y\rightarrow u_{4},8.~{}u_{3}*y\rightarrow
u_{2},9.~{}u_{3}*x\rightarrow u_{1},10.~{}u_{2}*x\rightarrow
u_{3}*u_{3},11.~{}u_{1}*y\rightarrow u_{3}*u_{3},12.~{}u_{1}*u_{4}\rightarrow
u_{2}*u_{2},13.~{}u_{4}*x*x\rightarrow
u_{2}*u_{3},14.~{}u_{3}*u_{3}*u_{3}\rightarrow u_{1}*u_{2}\\}$. This system
captures relationships among all product monomials appearing in the input.
This computation can be utilized in other examples specified using the same
monomial set and thus does not have to be repeated.

SingleACCompletion generating the congruence closure of the ground equations
on $+$ terms produces $R_{+}$ as a reduced canonical rewrite system:
$\\{1.~{}7~{}u_{1}\rightarrow 3~{}x,2.~{}4~{}u_{2}\rightarrow
u_{3},3.~{}3~{}u_{4}\rightarrow 0\\}$.

There are no constant rewrite rules to be propagated from one rewrite system
to another.

Due to the distributivity property relating $*$ and $+$, however, there are
interactions between a rule on $*$ monomials and a rule on $+$ monomials. To
illustrate, a $+$ rule 3, $3~{}u_{4}\rightarrow 0$ and a $*$ rule 12 interact
$~{}u_{1}*u_{4}\rightarrow u_{2}*u_{2}$ because of the common constant
$u_{4}$; Using the distributivity property,
$u_{1}*(3~{}u_{4})=u_{1}*(u_{4}+u_{4}+u_{4})=u_{1}*u_{4}+u_{1}*u_{4}+u_{1}*u_{4}$.
Rule 3 applies on the pure $+$ monomial in the mixed term $u_{1}*(3*u_{4})$
giving another mixed term $u_{1}*0$ which rewrites to 0 using a property of
CR1; rule 12 repeatedly applies on the pure $*$ monomials in the mixed term
$u_{1}*u_{4}+u_{1}*u_{4}+u_{1}*u_{4}$, producing
$u_{2}*u_{2}+u_{2}*u_{2}+u_{2}*u_{2}$. Consequently, a new equality on mixed
terms $u_{2}*u_{2}+u_{2}*u_{2}+u_{2}*u_{2}=0$ is derived. It is purified
resulting in new relations on $*$-monomials as well as $+$-monomials. Using a
new constant, say $v_{1}$, for $u_{2}*u_{2}$,a new $+$ rule $3v_{1}\rightarrow
0$ is generated and added to $R_{+}$; $u_{2}*u_{2}\rightarrow v_{1}$ is added
to $R_{*}$.

The pair $(u_{1}*0,u_{2}*u_{2}+u_{2}*u_{2}+u_{2}*u_{2})$ is a new kind of
critical pair, generated from the superposition
$u_{1}*(3~{}u_{4})=u_{1}*(u_{4}+u_{4}+u_{4})=u_{1}*u_{4}+u_{1}*u_{4}+u_{1}*u_{4}$
due to the distributivity property.

In general, given a $+$ rule $c~{}u\rightarrow r\in R_{+}$ and a $*$ rule
$u*m\rightarrow r^{\prime}\in R_{*}$ with a common constant $u$, where
$c\in\mathbb{Z}-\\{0\\}$, $r$ is a $+$-monomial, $m$ is $*$-monomial, the
superposition is $(c~{}u)*m$ generating the critical pair
$(r*m,c~{}r^{\prime})$ on mixed terms. $r*m$ and $c~{}r^{\prime}$ are
normalized using distributivity and other rules. The resulting mixed terms are
purified by introducing new constants, leading to new equalities on $*$
monomials as well as $+$ monomials. They are, respectively, added to $R_{*}$
and $R_{+}$ after normalizing and orienting them using the monomial ordering.

It can be proved that only this superposition needs to be considered to check
local confluence of $R_{+}\cup R_{*}$ [Kap21].

For more examples of superpositions and critical pairs between rules in
$R_{+}$ and $R_{*}$, consider rules $3$ and $13$, which give a trivial
critical pair.

Considering all such superpositions, the method terminates, leading to the
canonical rewrite system $\\{3~{}x\rightarrow 0,u_{1}\rightarrow
0,u_{2}\rightarrow u_{3},3~{}u_{4}\rightarrow 0\\}$ among other rules. When
converted into the original signature, this is precisely the Gröbner basis
reported as generated using Kandri-Rody and Kapur’s algorithm:
$\\{3~{}x\rightarrow 0,x*x*y\rightarrow 0,x*y*y\rightarrow
x*y,3~{}y*y*y\rightarrow 0\\}$ as reported in [KK88].

A very simple example taken from [KK88] consisting of 2 ground equations
$2~{}x*x*y=y,~{}3~{}x*y*y=x$ is also discussed in [Mar96] illustrating how
normalized rewriting and completion can be used to simulate Gröbner basis
computations on polynomial equations over the integers. Due to lack of space,
a detailed comparison could not be included; an interested reader is invited
to contrast the proposed approach with the one there. Comparing with [Mar96],
it is reported there that a completion procedure using normalized rewriting
generated 90 critical pairs in contrast to AC completion procedure [PS81]
which computed 1990 critical pairs. In the proposed approach, much fewer
critical pairs are generated in contrast, without needing to use any AC
unification algorithm or extension rules.

The above illustrates the power and elegance of the proposed combination
framework. The proposed approach can also be used to compute Gröbner basis of
polynomial ideals with coefficients over finite fields such as
$\mathbb{Z}_{p}$ for a prime number $p$, as well as domain with zero divisor
such as $\mathbb{Z}_{4}$ [KC09]. For example, in case of $\mathbb{Z}_{5}$,
another rule is added: $1+1+1+1+1\rightarrow 0$ added to the input basis.

##  12\. Conclusion

A modular algorithm for computing the congruence closure of ground equations
expressed using multiple AC function symbols and uninterpreted symbols is
presented; an AC symbol could also have additional semantic properties
including idempotency, nilpotency, identity, cancelativity as well as being an
Abelian group.

The algorithms are derived by generalizing the framework first presented in
[Kap97] for generating the congruence closure of ground equations over
uninterpreted symbols. The key insight from [Kap97]–flattening of subterms by
introducing new constants and congruence closure on constants, is generalized
by flattening mixed AC terms and purifying them by introducing new constants
to stand for pure AC terms in every (single) AC symbol. The result of this
transformation on a set of equations on mixed ground terms is a set of
constant equations, a set of flat equations relating a nonconstant term in an
uninterpreted symbol to a constant, and a set of equations on pure AC terms in
a single AC symbol possibly with additional semantic properties. Such
decomposition and factoring enable using congruence closure algorithms for
each of the pure subproblems on a single AC symbol and/or uninterpreted
symbols independently, which propagate equalities on shared constants. Once
the propagation of constant equalities stabilizes (reaches a fixed point), the
result is (i) a unique reduced canonical rewrite system for each subproblem
and finally, (ii) a unique reduced canonical rewrite system for the congruence
closure of a finite set of ground equations over multiple AC symbols and
uninterpreted symbols. The algorithms extend easily when AC symbols have
additional properties such as idempotency, identity, nilpotency, cancelativity
as well as their combination, as well as an AC symbol satisfying the
properties of a group.

Unlike previous proposals based on the Knuth-Bendix completion procedure and
its generalizations to AC theories and other equational theories, the proposed
algorithms do not need to use extensions rules, AC/E unification or AC/E
compatible termination ordering on arbitrary terms. Instead, the proposed
framework provides immense flexibility in using different termination
orderings for terms in different AC symbols and uninterpreted symbols.

The modularity of the algorithms leads to easier and simpler correctness and
termination proofs in contrast to those in [BTV03, Mar96]. The complexity of
the procedure is governed by the complexity of generating a canonical rewrite
system for AC ground equations on constants.

The proposed algorithm is a direct generalization of Kapur’s algorithm for the
uninterpreted case, which has been shown to be efficiently integrated into SMT
solvers including BarcelogicTools [NO07]. We believe that the AC congruence
closure can also be effectively integrated into SMT solvers. Unlike other
proposals, the proposed algorithm neither uses specialized AC compatible
orderings on nonground terms nor extension rules often needed in AC/E
completion algorithms and AC/E-unification, thus avoiding explosion of
possible critical pairs for consideration.

By-products of this research are new insights into combination algorithms for
ground canonical rewrite systems– (i) making orderings on ground terms more
flexible and general in which constants are allowed to be greater than
nonconstant terms as in orderings for Gröbner basis algorithms, as well as
(ii) a new kind of superpositions to incorporate semantic properties such as
cancelativity, which can be viewed as a superposition between the left side of
a rule with the right side of another rule. An instantiation of the proposed
combination framework to an AC symbol with identity and another AC symbol
enriched with the properties of an Abelian group, leads to a new way to view a
Gröbner basis algorithm for polynomial ideals over integers, as a combination
congruence closure algorithm. Canonical bases of polynomial ideals over the
integers (and more general rings with zero divisors) on extended signature can
be computed for which rewriting (simplification) of a polynomial by another
polynomial, when viewed on the original signature, need not terminate, opening
possible new research directions.

Integration of additional semantic properties of AC symbols is done on a case
by case basis by modifying critical pair constructions but avoiding extension
rules or complex AC termination orderings compatible with these semantic
properties. Using the proposed framework, other types of combinations of
congruence closure algorithms in the presence of other semantic properties
relating different AC symbols may be possible as well. A general approach
based on the properties of the semantic properties that avoids extension rules
and sophisticated machinery due to the AC properties, needs further
investigation.

Acknowledgments: Thanks to the FSCD conference and LMCS referees for their
reports which substantially improved the presentation. Thanks to Christophe
Ringeissen for detailed discussions which clarified the contributions of the
proposed approach in contrast to other approaches, particularly [BT97,
ABRS09]. Thanks also to Jose Abel Castellanos Joo for help in formulating some
of the algorithms.

## References

  * [ABRS09] Alessandro Armando, Maria Paola Bonacina, Silvio Ranise, and Stephan Schultz.  New results on rewrite-based satisfiability procedures.  ACM Transactions on Computational Logic, 10(1):4:1–4:51, 2009.
  * [BD88] Leo Bachmair and Nachum Dershowitz.  Critical pair criteria for completion.  J. Symb. Comput., 6(1):1–18, 1988.  doi:10.1016/S0747-7171(88)80018-X.
  * [BK20] Franz Baader and Deepak Kapur.  Deciding the word problem for ground identities with commutative and extensional symbols.  In Nicolas Peltier and Viorica Sofronie-Stokkermans, editors, Automated Reasoning - 10th International Joint Conference, IJCAR 2020, Paris, France, July 1-4, 2020, Proceedings, Part I, volume 12166 of Lecture Notes in Computer Science, pages 163–180. Springer, 2020.
  * [BL81] A. Michael Ballantyne and Dallas Lankford.  New decision algorithms for finitely presented commutative semigroups.  Computers and Mathematics Applications, 7:159–165, 1981.
  * [BN98] Franz Baader and Tobias Nipkow.  Term rewriting and all that.  Cambridge University Press, 1998.
  * [BRTV00] Leo Bachmair, I. V. Ramakrishnan, Ashish Tiwari, and Laurent Vigneron.  Congruence closure modulo associativity and commutativity.  In Hélène Kirchner and Christophe Ringeissen, editors, Frontiers of Combining Systems, Third International Workshop, FroCoS 2000, Nancy, France, March 22-24, 2000, Proceedings, volume 1794 of Lecture Notes in Computer Science, pages 245–259. Springer, 2000.  doi:10.1007/10720084\\_16.
  * [BT97] Franz Baader and Cesare Tinelli.  A new approach for combining decision procedure for the word problem, and its connection to the Nelson-Oppen combination method.  In William McCune, editor, Automated Deduction - CADE-14, 14th International Conference on Automated Deduction, Townsville, North Queensland, Australia, July 13-17, 1997, Proceedings, volume 1249 of Lecture Notes in Computer Science, pages 19–33. Springer, 1997.  doi:10.1007/3-540-63104-6\\_3.
  * [BTV03] Leo Bachmair, Ashish Tiwari, and Laurent Vigneron.  Abstract congruence closure.  J. Autom. Reason., 31(2):129–168, 2003.  doi:10.1023/B:JARS.0000009518.26415.49.
  * [BWK93] Thomas Becker, Volker Weispfenning, and Heinz Kredel.  Gröbner bases - a computational approach to commutative algebra, volume 141 of Graduate texts in mathematics.  Springer, 1993.
  * [DST80] Peter J. Downey, Ravi Sethi, and Robert Endre Tarjan.  Variations on the common subexpression problem.  J. ACM, 27(4):758–771, 1980.  doi:10.1145/322217.322228.
  * [Kap97] Deepak Kapur.  Shostak’s congruence closure as completion.  In Hubert Comon, editor, Rewriting Techniques and Applications, 8th International Conference, RTA-97, Sitges, Spain, June 2-5, 1997, Proceedings, volume 1232 of Lecture Notes in Computer Science, pages 23–37. Springer, 1997.  doi:10.1007/3-540-62950-5\\_59.
  * [Kap19] Deepak Kapur.  Conditional congruence closure over uninterpreted and interpreted symbols.  J. Syst. Sci. Complex., 32(1):317–355, 2019.  doi:10.1007/s11424-019-8377-8.
  * [Kap21] Deepak Kapur.  Weird Gröbner Bases: An application of associative-commutative congruence closure algorithm.  Tech report under preparation, Department of Computer Science, University of New Mexico, May 2021.
  * [KB70] Donald Knuth and Peter Bendix.  Simple word problems in universal algebras.  In Leech, editor, Computational Problems in Abstract Algebra, pages 263–297. Pergamon Press, 1970.
  * [KC09] Deepak Kapur and Yongyang Cai.  An algorithm for computing a Gröbner basis of a polynomial ideal over a ring with zero divisors.  Math. Comput. Sci., 2(4):601–634, 2009.  doi:10.1007/s11786-009-0072-z.
  * [KK83] Abdelilah Kandri-Rody and Deepak Kapur.  On relationship between Buchberger’s Gröbner basis algorithm and the Knuth-Bendix completion procedure.  Technical report no. ge-83crd286, General Electric Corporate Research and Development, Schenectady, NY, November 1983.
  * [KK88] Abdelilah Kandri-Rody and Deepak Kapur.  Computing a Gröbner basis of a polynomial ideal over a Euclidean domain.  J. Symb. Comput., 6(1):37–57, 1988.  doi:10.1016/S0747-7171(88)80020-8.
  * [KKN85] Abdelilah Kandri-Rody, Deepak Kapur, and Paliath Narendran.  An ideal-theoretic approach to work problems and unification problems over finitely presented commutative algebras.  In Jean-Pierre Jouannaud, editor, Rewriting Techniques and Applications, First International Conference, RTA-85, Dijon, France, May 20-22, 1985, Proceedings, volume 202 of Lecture Notes in Computer Science, pages 345–364. Springer, 1985.  doi:10.1007/3-540-15976-2\\_17.
  * [KKW89] Abdelilah Kandri-Rody, Deepak Kapur, and Franz Winkler.  Knuth-Bendix procedure and Buchberger algorithm: A synthesis.  In Gaston H. Gonnet, editor, Proceedings of the ACM-SIGSAM 1989 International Symposium on Symbolic and Algebraic Computation, ISSAC ’89, Portland, Oregon, USA, July 17-19, 1989, pages 55–67. ACM, 1989.  doi:10.1145/74540.74548.
  * [KMN88] Deepak Kapur, David R. Musser, and Paliath Narendran.  Only prime superpositions need be considered in the Knuth-Bendix completion procedure.  J. Symb. Comput., 6(1):19–36, 1988.  doi:10.1016/S0747-7171(88)80019-1.
  * [KN92] Deepak Kapur and Paliath Narendran.  Double-exponential complexity of computing a complete set of AC-unifiers.  In Proceedings of the Seventh Annual Symposium on Logic in Computer Science (LICS ’92), Santa Cruz, California, USA, June 22-25, 1992, pages 11–21. IEEE Computer Society, 1992.  doi:10.1109/LICS.1992.185515.
  * [KZ95] Deepak Kapur and Hantao Zhang.  An overview of Rewrite Rule Laboratory (RRL).  Computers and Mathematics Applications, 29:91–114, 1995.
  * [LeC83] Philippe LeChenadec.  Canonical forms in the finitely presented algebras.  Ph.D. Dissertation, U. of Paris 11, 1983.
  * [Mar91] Claude Marché.  On ground AC-completion.  In Ronald V. Book, editor, Rewriting Techniques and Applications, 4th International Conference, RTA-91, Como, Italy, April 10-12, 1991, Proceedings, volume 488 of Lecture Notes in Computer Science, pages 411–422. Springer, 1991.  doi:10.1007/3-540-53904-2\\_114.
  * [Mar96] Claude Marché.  Normalized rewriting: An alternative to rewriting modulo a set of equations.  J. Symb. Comput., 21(3):253–288, 1996.  doi:10.1006/jsco.1996.0011.
  * [MM82] Ernst W. Mayr and Albert Meyer.  The complexity of the word problems for commutative semigroups and polynomial ideals.  Advances in Mathematics, 46:305–439, 1982.
  * [NO80] Greg Nelson and Derek C. Oppen.  Fast decision procedures based on congruence closure.  J. ACM, 27(2):356–364, 1980.  doi:10.1145/322186.322198.
  * [NO07] Robert Nieuwenhuis and Albert Oliveras.  Fast congruence closure and extensions.  Inf. Comput., 205(4):557–580, 2007.  doi:10.1016/j.ic.2006.08.009.
  * [NR91] Paliath Narendran and Michaël Rusinowitch.  Any ground associative-commutative theory has a finite canonical system.  In Ronald V. Book, editor, Rewriting Techniques and Applications, 4th International Conference, RTA-91, Como, Italy, April 10-12, 1991, Proceedings, volume 488 of Lecture Notes in Computer Science, pages 423–434. Springer, 1991.  doi:10.1007/3-540-53904-2\\_115.
  * [PS81] Gerald E. Peterson and Mark E. Stickel.  Complete sets of reductions for some equational theories.  J. ACM, 28(2):233–264, 1981.  doi:10.1145/322248.322251.
  * [Sho78] Robert E. Shostak.  An algorithm for reasoning about equality.  Commun. ACM, 21(7):583–585, 1978.  doi:10.1145/359545.359570.
  * [Win84] Franz Winkler.  The Church-Rosser Property in Computer Algebra and Special Theorem Proving: An Investigation of Critical-Pair/Completion Algorithms.  Ph.d. dissertation, University of Linz, 1984.
  * [WNW+21] Max Willsey, Chandrakana Nandi, Yisu Remy Wang, Oliver Flatt, Zachary Tatlock, and Pavel Panchekha.  egg: Fast and extensible equality saturation.  Proc. ACM Program. Lang., 5(POPL):1–29, 2021.  doi:10.1145/3434304.
  * [Zha92] Hantao Zhang.  Implementing contextual rewriting.  In Michaël Rusinowitch and Jean-Luc Rémy, editors, Conditional Term Rewriting Systems, Third International Workshop, CTRS-92, Pont-à-Mousson, France, July 8-10, 1992, Proceedings, volume 656 of Lecture Notes in Computer Science, pages 363–377. Springer, 1992.  doi:10.1007/3-540-56393-8\\_28.


# Modelling COVID-19-III: endemic spread in India

Madhuchhanda Bhattacharjee School of Mathematics and Statistics University of
Hyderabad Gachibowli Hyderabad 500046 India email<EMAIL_ADDRESS>Arup
Bose111Research supported by J.C. Bose National Fellowship, Department of
Science and Technology, Government of India. Stat-Math Unit, Kolkata Indian
Statistical Institute 203 B.T. Road Kolkata 700108 India email:
<EMAIL_ADDRESS>
###### Abstract

A disease in a given population is termed endemic when it exhibits a steady
prevalence. We address the pertinent question as to what extent COVID-19 has
turned endemic in India.

There are several existing models for studying endemic behaviour, such as the
extensions of the traditional temporal SIR model or the spatio-temporal
endemic-epidemic model of Held et al. (2005) and its extensions.

We propose a “spatio-temporal Gravity model” in a state of the art generalised
linear model set up that can be deployed at various spatial resolutions. In
absence of routine and quality covariates in the context of COVID-19 at finer
spatial scales, we make use of extraneous covariates like air-traffic
passenger count that enables us to capture the local mobility and social
interactions effectively. This makes the proposed model different from the
existing models. The proposed gravity model not only produces consistent
estimators, but also outperforms the other models when applied to Indian
COVID-19 data.

Keywords. Air-passenger traffic data, COVID-19, distance decay, endemic-
epidemic model, endemic phase, gravity model, Indian sub-continent,
seasonality, spatial auto-regression, spatio-temporal model.

AMS 2010 Subject Classification. Primary 62P10; Secondary 92D30

##  1 Introduction

COVID-19 hit the world in December 2019, and there has been a large volume of
research articles devoted to its study from numerous angles. For a brief
background on the history and spread of COVID-19 in India, the reader is
referred to Bhattacharjee and Bose (2021), and the references therein. Data on
COVID-19 in India is primarily available in the form of three time series of
Cases ($C$), Recoveries ($R$), and Deaths ($D$).

The most commonly used model for this data is some form of the temporal SIR
model–with added compartments depending either on any additional information
which may be available, or on the required policy prescriptions. When fitted
dynamically, these models are capable of very good short term predictions.
Most models on the spread of COVID-19 in India are built from the temporal
view point.

India is often referred to as a sub-continent for its size and diverse nature
of the regions. The progression of COVID-19 cases in the sub-continent over
time has shown a wide variation across its 36 states and union territories,
henceforth collectively referred to as regions. This is evident from the plots
and discussions in Bhattacharjee (2020), Ranjan (2020) and Bhattacharjee et
al. (2021), as well as other sources.

It is naturally expected that the propagation of COVID-19 to neighboring
locations would lead to a spatial pattern. Areas which are connected, either
geographically, or via movement of people by railways, road, air and
waterways, would influence each others’ caseload. Moreover, due to many
factors, including the natural topography of India, regions of India have
varying degrees of accessibility and connectivity. As is well known, temporal
SIR and related models cannot unearth the spatial dependence that may be
present in a time series data.

On the other hand, it is widely accepted that taking spatial dependence into
account can significantly impact statistical analysis of epidemics. See for
example Clif and Ord (1981). Kirby et al. (2017) provides an excellent review
of spatio-temporal methods in epidemiological data. Incidentally, spatio-
temporal analysis of US and European COVID-19 data have been made by several
authors. See for instance Sannigrahi et al. (2020) and Mollalo et al. (2020).
For Indian COVID-19 data, it is natural to consider spatio-temporal models at
the regional level, with districts (sub-regions) being considered as the
primary units in each region.

To develop these models for India is a challenge for several reasons, not the
least of which is the paucity and unreliability of data on the relevant
covariates, and on the auxiliary variables. The extent and accuracy of this
information varies across the regions, and hence not every model is
meaningfully implementable in every region. We are aware of a few attempts
such as Bhattacharjee et al. (2021) and Ganesan and Subramani (2020, 2021).

When a disease within a given geographic region exhibits a steady prevalence
in the population, then it is termed endemic. An endemic behaviour is
connected to it being seasonal, so that the disease incidence sees spikes
around the same time every year–an obvious example is the common flu in
different parts of the world. There is a strong opinion (see for instance
Phillips (2021)) that COVID-19 will eventually turn endemic. A pertinent
question is whether we are witnessing the onset of its endemic behaviour. Some
researchers, for instance Bjórnstad et al. (2020), have estimated the time to
reach an endemic equilibrium for COVID-19 as four years.

The traditional endemic models are based on modified versions of the SIR
models. The endemic-epidemic (EE) models are more sophisticated spatio-
temporal models. They are time series models for multivariate surveillance
counts, and were introduced by Held et al. (2005). Subsequently, these have
been extended by several authors. For these models, the counts for the
different variables are postulated to progress over time in a Markovian
manner. While the parameters across the regions are in general different, some
of them may turn out to be identical. However, the regional variations and
apparent disharmonious occurrence of COVID-19 in different parts of India
cause difficulties with the direct applications of such models to this sub-
continent.

COVID-19 has been observed to have a strong epidemic component (see for
example Giuliani et al. (2020)), with noticeable spatial variation. A natural
approach would be to employ a spatio-temporal model at a sub-regional level.
This would allow for possible forecasts on various levels of spatial
aggregation, namely sub-regional, regional and also perhaps national. Thus it
would be worthwhile to derive a meaningful modelling framework at a finer
geographic resolution that would provide consistent and meaningful inference
on some of the essential aspects of this pandemic.

We have succeeded in developing spatio-temporal models for the Indian district
level data within each region and have investigated their performance with two
specific objectives. The first is to connect public mobility, social
connectivity, and local behaviour with the extent of occurrence of COVID-19.
The second objective is to assess how far it has reached an endemic state at
the individual regional levels as well as collectively.

We consider the publicly available daily COVID-19 incidence data at the
district level. Unfortunately, there are multiple issues with the available
data, and in particular, it is affected by many factors that are extraneous to
the epidemic. Thus we are forced to fall back on data on surrogate covariates
that are affected by, as well as are influencing the epidemic. Data on these
surrogates could possibly be collated by independent resources. While
theoretically there are many possibly useful covariates, in practice there is
a stark lack of data on such covariates at required space and time resolution.

The air-traffic passenger data is publicly available, and we consider this as
an important surrogate. It stands to reason that the observed spatial
fluctuations in the counts for cases are often influenced by the local social
events, such as festivals. At the same time, these very factors also influence
the corresponding “air-traffic” covariate. Hence, in absence of any other
reliable local spatial predictors, this surrogate covariate could play a
significant role. However, air-traffic data is available only on a monthly
basis, and there is a lag of a few weeks in its availability. To use this data
meaningfully, we aggregated the COVID-19 incidence data at the monthly level.

For the second aspect, since the extent of endemicity is somewhat subjective,
we use different models to assess this, and compare our estimates across the
different novel models that we employ.

The remaining part of this article is organised as follows. In Section 2 we
review the available literature on spatio-temporal models for epidemiological
data, especially those with endemic components. Section 3 contains our
proposed gravity models. Along with these novel models we also consider
modifications of the existing models. In Section 4, we present the findings
from these models, based of COVID-19 data from the Indian sub-continent, and
demonstrate the superiority of their performance The final Section 5 contains
discussion along with possible extensions of the current work.

##  2 Modeling endemic behaviour

There have been the two major paths explored for modelling endemic patterns.
The first of these include the mathematical and statistical models that extend
the well-known SIR model. We give an example of such an effort by Bjórnstad et
al. (2020). The second is the wholly probabilistic modelling pioneered by Held
et al. (2005), which we shall discuss in details.

Extension of SIR models

As mentioned earlier, the SIR models are among the most popular epidemic
models. The original SIR models can also be extended by taking into account
the number of exposed individuals (see Brauer et al. (2019)), say $E$. Since
there is a latent period between being infected and becoming infectious, the
trajectory of counts for $C$ will be in a lag with that for $E$ by this
latency period. With the reproduction number $R_{0}>1$, both $E$ and $I$ are
predicted to eventually stabilize at an endemic equilibrium (see Bjórnstad et
al. (2020) for an illustration).

Endemic-Epidemic (EE) models

A study by Colón-González et al. (2018) indicates that use of multiple data
streams arising from surveillance activities can be a useful approach to
disease detection. They opine that, syndromic surveillance complements
traditional public health surveillance, by collecting and analysing health
indicators in near real time. It is imperative that appropriate statistical
techniques are used to analyse such data.

The endemic-epidemic (EE) model was introduced by Held et al. (2005). These
were subsequently extended in Paul et al. (2008), Held and Paul (2012), Meyer
and Held (2014) and Meyer et al. (2017). These models have also been
effectively applied to several other epidemiological problems. See Dunbar and
Held (2020) for a review.

The EE model is based on the motivation that, while the epidemic component
enables the capture of occasional outbreaks, the endemic component should
explain the baseline rate of cases that is persistent with a stable temporal
pattern (see Adegboye and Adegboye (2017)). Suppose we have weekly data on
counts for different variables in different regions. Then, we assume the
Markovian structure, so that the counts for the $t$th week depend only on
those for the $(t-1)$th week. This entails assuming implicitly that the time
between appearance of the symptoms in successive generations is the same, and
equals the observation interval of a week. However, in reality, this serial
interval may vary randomly across infection events and may be longer than a
single observation interval (see Becker (2015), p.156). Other factors can also
introduce dependencies on time points $t-2,\ldots,t-D$ and $D$ can be quite
high. This limitation was addressed by Bracher and Held (2020) who introduced
flexible weighting schemes for past incidences. In particular, they suggested
the shifted Poisson, triangular and geometric weights based on the different
empirical distributions of the serial intervals.

Consider a specific region with $P$ districts. Let the number of cases in
district $i$ for week $t$ be denoted by $Y_{i,t}$, $i=1,\ldots,P$ and let
$\textbf{Y}_{t}$ be the vector of these values. A simple example of an EE
model is the following: first, conditionally on the past, $Y_{i,t}$ and
$Y_{j,t}$ from different districts (units) $i$ and $j$ at time $t$ are assumed
to be independent. Second, each $Y_{i,t}$, conditionally on the past, is
assumed to follow a negative binomial distribution. That is,

$Y_{i,t}|\textbf{Y}_{t-1},\textbf{Y}_{t-2},\ldots\sim\text{NB}(\lambda_{it},\psi_{i}),\
1\leq i\leq P,\ 1\leq t\leq T,$ (2.1)

with conditional mean $\lambda_{it}$ and over dispersion parameter $\psi_{i}$.
The $\\{\lambda_{it}\\}$ is further modeled as

$\lambda_{it}=\nu_{it}+\phi_{it}\sum_{j=1}^{P}\lfloor
w_{ji}\rfloor\textbf{Y}_{j,t-1},\ 1\leq i\leq P,\ 1\leq t\leq T.$ (2.2)

Here the parameter $\nu_{it}$ is the endemic component, and captures the
number of infections that are not directly linked to the observed cases from
the previous week. The remaining autoregressive term in (2.2) is the epidemic
component and describes how the incidence in district $i$ is linked to
previous cases $Y_{j,t-1}$ in districts $j=1,\ldots,P$.

The conditional variance is thus $\lambda_{it}+\psi_{i}\lambda_{it}^{2}$ where
$\\{\psi_{i}\\}$ are the over dispersion parameters. When $\psi_{i}=0$, the
conditional variance reduces to $\lambda_{it}$. A common simplification is to
assume that all the over dispersion parameters $\psi_{i}$ are identical, say,
equal to $\psi$.

The parameters $\\{\nu_{it}\\}$ and $\\{\phi_{it}\\}$ are constrained to be
non-negative and modeled in a log-linear fashion, for instance with sine-
cosine terms to account for seasonality (Held and Paul (2012)). Long-term
temporal trends and covariates such as meteorological conditions (Cheng et al.
(2016), Bauer and Wakefield (2018)) or vaccination coverage (Herzog et al.
(2011)) could also be included. The coupling between the districts is achieved
by the $P^{2}$ weights $w_{ji}$ which enter in (2.2) after normalization as:

$\lfloor w_{ji}\rfloor=\dfrac{w_{ji}}{\sum_{h=1}^{P}w_{jh}}.$ (2.3)

Unrestricted estimation of all $P^{2}$ weights is usually unstable in
practice, so more parsimonious and epidemiologically meaningful
parameterizations have been introduced. Specifically, the weights can be based
on social contact data for spread across age groups (Meyer and Held (2017)),
or on the geographical distance between cases to describe spatio-temporal
spread (Meyer and Held (2014)). In the latter case the weights can be
specified through a power law

$w_{ji}=(o_{ji}+1)^{\rho},$ (2.4)

where $o_{ji}$ is the path distance between the regions $j$ and $i$
($o_{ii}=0$ for all $i$, $o_{ji}=1$ when $i$ and $j$ are direct neighbors, and
so on), and $\rho$ is a decay parameter to be estimated from the data. The
power law formulation is motivated by human movement behaviour (see Brockmann
et al. (2006)), and has been found to be an efficient way of capturing spatial
dependence.

Modified EE model

The EE model of Berlamann and Haustein (2020) is more advanced. It allows for
possible over dispersion due to under-reporting or unobserved covariates that
affect the disease incidence. As before, suppose that the conditional
distribution of $Y_{i,t}$, given the past history up to time $t-1$ which is
symbolized as the $\sigma$-field $F_{t-1}$, is negative binomial with
parameters $\mu_{i,t}$ and $\psi_{i}$. We write

$Y_{i,t}|F_{t-1}\sim NB(\mu_{i,t},\psi_{i}),\ 1\leq i\leq P,\ 1\leq t\leq T.$

The parameter $\psi_{i}>0$ is the over dispersion parameter, so that the
conditional variance of $Y_{i,t}$ is given by

$\text{Var}[Y_{i,t}|F_{t-1}]=\mu_{i,t}\left(1+\psi_{i}\mu_{i,t}\right).$

The mean $\mu_{i,t}$ is given by

$\mu_{i,t}=\lambda_{i,t}\sum_{d=1}^{D}u_{d}Y_{i,t-d}+\phi_{i,t}\sum_{d=1}^{D}\sum_{j\neq
i}u_{d}\omega_{j,i}Y_{j,t-d}+\nu_{t},\ 1\leq i\leq P,\ 1\leq t\leq T.$ (2.5)

The three summands on the right side of (2.5) correspond to the epidemic-
within districts, the epidemic-between districts and the endemic components.

The within component first summand is modeled as an autoregression of the
number of cases in district $i$, on a weighted sum of the past number of cases
up to day $t-D$. The autoregression parameter is $\lambda_{i,t}$.

The between component second summand attempts to capture the spread of a
disease across different locations. It is taken as a regression of the number
of cases in district $i$ on a weighted sum of the past number of cases up to
time $t-D$ in other districts $j\neq i$. The autoregressive parameter
$\phi_{i,t}$ is district specific, and may depend on covariates. The component
is completed by a system of weights that have two factors; while $\\{u_{d}\\}$
are the same as for the within component, the parameters $\\{\omega_{j,i}\\}$
account for the spatial distance between districts.

The last term $\nu_{t}$ is the endemic component that describes the
seasonality of the difference between regions, and typically it is modeled by
capturing the seasonal variations of the disease through a harmonic wave:

$\log\nu^{END}_{t}=\alpha_{0}+\eta t+\gamma\sin(\omega t)+\delta\cos(\omega
t),$ (2.6)

with $\alpha_{0}$ being a constant, $\eta t$ being a time trend and $\gamma
sin(\omega t)+\delta cos(\omega t)$ capturing possible seasonal variation of
the endemic component as is typical for many viral diseases. The parameter
$\alpha_{0}$ could also be equal to $\alpha_{i}$ so that it varies across
districts. The parameter $w$ could then be $2\pi/365$ as Fourier frequencies
corresponding to daily data.

The EE modelling approach has been used in Grimée et al. (2022). More
recently, Celani and Giudici (2021) applied this model to data from Italy,
taking several covariates, such as the population density, a Stringency Index,
a Testing Policy index, and an indicator for weekend days. Two of their key
covariates were available at only the national level.

##  3 Proposed models

As mentioned earlier, the objectives here are twofold. The first is to be able
to connect public mobility, social connectivity, and local behaviour with the
extent of occurrence of COVID-19. In a noticeable contrast to theoretical
possibilities of such candidate covariates, in reality there is a stark lack
of data at required space and time resolution for the relevant covariates.

The second objective is surveillance of this epidemic to assess whether it has
reached a manageable state. Since the Indian sub-continent is extremely
heterogeneous in its socio-demographic profile, this would be key in assessing
the status of the epidemic at the individual regional levels as well as
collectively.

As mentioned earlier we will model the COVID-19 incidence data using surrogate
covariates like air-traffic passenger data which is publicly available and
collated by independent resources. It stands to reason that the observed
spatial fluctuations in the counts for cases are often influenced by the local
social events, such as festivals. At the same time, these same factors also
influence the corresponding “air-traffic” covariates. Hence, in absence of any
other reliable local spatial predictors, this data could play a significant
role. Unfortunately, the air-traffic data is available only on a monthly
basis, and in addition there is a lag of a few weeks in its availability. To
use this data gainfully, we aggregated the COVID-19 incidence data at the
monthly level.

Gravity model

As is commonly done, we apply a negative binomial ‘observation-driven’ model
for count data. The proposed model for (log-)mean has effectively three
components. The first component captures the purely time effects, the second
captures the purely spatial effect, and the third is a combination of spatio-
temporal effects reflecting the local variations.

We also borrow ideas from the gravity models. These models can be traced to
Ravenstein (1889) and Zipf (1946), and they have become popular due to the
subsequent work of Tinbergen (1962). See Anderson (2011) for a review of these
models. In a regression context, we may frame this model as follows: consider
locations $i$ and $j$ between which some form of transfer of population, say
$Y_{ij}$ occurs. Assume that information on the total movement from $i$ is
available as $U_{i}$ and total arrival at $j$ is available as $V_{j}$.
Additionally let $d_{ij}$ be the “distance” between the two locations. Then we
postulate that

$Y_{i,j}=\theta_{0}U_{i}^{\theta_{1}}V_{j}^{\theta_{2}}{d_{ij}}^{\theta_{3}}\eta_{ij},$
(3.1)

where $\theta$’s are unknown parameters and $\eta_{ij}$ is an error term.
Since $Y_{ij}$ are non-negative, a modified version of the above model
stipulates that the mean $\mu_{ij}$ of $Y_{ij}$ is given by:

$\log\mu_{ij}=\theta_{0}+{\theta_{1}}\times U_{i}+{\theta_{2}}\times
V_{j}+{\theta_{3}}\times{d_{ij}}.$ (3.2)

Even though the above model assumes a power decay $d_{ij}^{\theta_{3}}$
between the districts $i$ and $j$, more generally, it may be modeled by a
decreasing distance decay function $f$. Some specific choices of $f$ that have
been used in the literature are the (i) power decay
$f(d_{ij})=d_{ij}^{-\alpha}$, (ii) exponential-normal decay
$f(d_{ij})=e^{-\alpha d_{ij}^{2}}$, (iii) exponential-square-root decay
$f(d_{ij})=e^{-\alpha\sqrt{d_{ij}}}$, and (iv) exponential decay
$f(d_{ij})=e^{-\alpha d_{ij}}.$ For other choice of distances, see Berlamann
and Haustein (2020).

For count data regression when an offset variable is used, then the
corresponding regression coefficient is prefixed to be 1. Such a model then
represents rates in comparison to counts. In the literature it has been
suggested to use the population size as an offset term (see Cox (1981)).
However due to various shortcomings in the data available, we do not use the
population size as an offset, and instead assign a set of geographic region
specific (e.g. district level, regional level) parameters, which would
effectively serve the same purpose. An endemic component is also included in
the model as suggested earlier.

We employ the gravity models in two variations. Suppose there are $K_{i}$
airports within the $i$-th region. In variation 1, at the $t$-th time point,
for the $j$-th district of the $i$-th region, the log-mean of incidence is
given by

$\lambda_{ijt}=\nu_{it}+\phi_{ij;}+\sum_{k=1}^{K_{i}}\theta_{ik}D_{ijk}*X_{ikt},\
1\leq j\leq R_{i},\ 1\leq i\leq P,\ 1\leq t\leq T,$ (3.3)

where $D$ contains the distances with exponential decay, and $X$ contains the
air-passenger-traffic data.

In variation 2, we have a composite gravity model, where,

$\lambda_{ijt}=\nu_{t}+\phi_{ij}+\sum_{k=1}^{K_{i}}\theta_{ik}D_{ijk}*X_{ikt},\
1\leq j\leq R_{i},\ 1\leq i\leq P,\ 1\leq t\leq T,$ (3.4)

and $D$ and $X$ are as above.

Both models (3.3) and (3.4) can also be implemented with an intercept term,
while adjusting the $\nu_{t}$ and $\phi_{ij}$ parameters accordingly.

Endemic-epidemic type models

For local COVID-19 infection counts, following the ideas of Held et al.
(2005), we develop a class of spatio-temporal regression models that have both
endemic and epidemic components.

The scarcity of useful covariates led us to use air-traffic data for the
Gravity models described earlier, and this data is available only at a monthly
frequency. Note however that the EE-type models are (spatial-)autoregressive
in nature. Since COVID-19 is a rapidly evolving epidemic, such an
autoregressive model at the monthly time scale is not expected to capture its
evolution. Therefore the relevant parameters in these models are then
estimated using the daily occurrence data from various regions.

In our models, the expected conditional mean spread of the infectious disease
is decomposed into endemic and epidemic components. However, deviating from
the existing methods, we implemented these models at various levels of
aggregation of the parameters.

For example, for each individual region $i$, the expected conditional mean
$\mu_{it}$ at a time $t$ for that region is postulated as:

$\log\mu_{it}=\nu^{EPI}_{it}+\nu^{END}_{it},$

with an over-dispersion parameter $\psi>0$.

In contrast, in a separate model, we assume that the endemic component
$\nu^{END}_{it}$ at the time $t$, is common to all the regions, regardless of
the infection history of district $i$ and its neighbors. Let the common value
be $\nu^{END}_{t}$. In that case $\log\mu_{it}$ reduces to

$\log\mu_{it}=\nu^{EPI}_{it}+\nu^{END}_{t}.$

Interestingly Berlamann and Haustein (2020) further multiply $\nu^{END}_{t}$
by the size of the local population. While this would be a preferred
structure, shortcomings in the data prevents us from using this modification.

As described earlier, most commonly used structure for the epidemic part of
the regression are primarily autoregressive, where the spatial influences are
also modeled as a lagged regression (see Held et al. (2005), Bracher and Held
(2020)). This enables the capture of path dependencies and of self-exciting
behaviour that are known to be common with infectious diseases. Thus for the
epidemic part of our model, we adhere to autoregressive components similar to
those proposed in Held et al. (2005). For the spatial autoregression
component, we use the entries of the spatial lag-$1$ adjacency matrix as
weights.

Depending on the context, additional covariates such as age and gender-
specific effects, contagion in nearby districts in the geographic and social
space, as well as latent heterogeneity between the districts are taken into
account. See Fritz and Kauermann (2020)) for an example. They also provide an
extensive covariate based model for the endemic component. In our models we
have used region specific parameters.

Berlamann and Haustein (2020) use a region specific random effects component
in the epidemic part of their model, in conjunction with the spatial
autoregressive part. In contrast we have employed a fixed effect parameter for
the effects of regions, and the varying population sizes of the individual
regions are subsumed in this parameter.

##  4 Results

The data sources used for implementation of the models are
https://data.incovid19.org for COVID-19 incidence data, and
https://www.aai.aero for air traffic passenger data. We have compared various
alternative model specifications on these data. For model comparison and model
fit to assessment, we have used both, the Akaike Information Criterion (AIC)
and the Nagelkerke pseudo-$R^{2}$. Due to widely different nature of the
models implemented in this data, we have presented the pseudo-$R^{2}$ values.
The fact that these values are naturally limited to a common range of values,
eases the comparison across various models. It is to be noted that the gravity
models have been implemented with monthly aggregated data whereas the spatial-
autoregressive-type EE model uses daily incidence data.

The main commonality of the proposed Gravity and EE-type models is that they
all have an endemic component, at regional or national level, denoted by
$\nu_{.}$. The remaining components of these models pertain to capturing the
epidemic behaviour. We have applied various alternative specifications of the
epidemic components of the model, which serve to check the robustness of our
conclusions with regard to the estimation of the endemic pattern.

To assess the presence of a spatially-consistent endemic component, we
implemented the gravity model at the regional level with individual region
specific endemic parameters $\\{\nu_{it}\\}$, as well as an overall model
based on all regions with shared endemic parameters $\\{\nu_{t}\\}$. We
compared our findings with the competing EE-type model. Further, to assure us
against spuriousness, we also used multiple data sources and checked the
consistency of the endemic-parameter estimates, in terms of both, the pattern
as well as scale of the estimated values.

Gravity model: offset parameter estimates. As mentioned earlier, we did not
include the (log)-population in the model, and instead used a separate
parameter to capture the location. We have validated our choice of structure
for the various parts of our model (see expression (3.4)). For example, based
on the proposed model (3.3), for each region, we obtained estimates for
locations as sum of the estimated intercept parameter(s), and the average of
district level parameters from the respective region. From Figure 1 we observe
that, while there is a reasonable relationship between these estimated
location parameters and the corresponding log-population values, there are
some departures. This justifies our choice of the model structure in this
respect. We have also assessed the relationship with the population density,
and the observed pattern is similar to that in Figure 1.

Figure 1: Gravity model for individual regions: log population (in horizontal
axis) against estimated location parameter (in vertical axis).

Gravity model: spatial parameter estimates. In the epidemic part of the
proposed model, the very first components are the spatial/regional effect
parameters, $\phi$’s. In Figure 2 we have presented the estimated parameters
for various districts in India under the Gravity model implemented for the 17
regions where districts level COVID-19 data is available, in a combined
manner. The interpretation of the parameters (and their estimates) are on the
log-scale.

Figure 2: Gravity model at the regional level: estimated district level
parameters for the 17 regions. For the districts from other regions, the
overall average has been used. Figure 3: Gravity model at the regional level:
estimates of quantiles of the location adjusted district parameters for each
of the 17 regions.

Gravity model: distributional commonality of the spatial parameters. For each
region, once we take away the location effect as described above, the centered
estimates of the district level parameters appear to have similar
distribution. In Figure 3 we have provided the estimated quantiles of the
centered district level parameters for each region. This demonstrates the
commonality of the underlying behaviour of COVID-19 occurrences across
regions.

Gravity model: robustness of spatial pattern against implementation level.
These models are flexible enough to be implementable either in an individual
region or as a combined model over multiple regions. The distributional
behaviour of the district level estimated parameters appear to be extremely
similar in nature for the individual and the combined models. These results
are not shown to avoid repetition. Since the findings appear to be not
sensitive to the level of implementation, it provides the assurance of the
robustness of the pattern captured by the Gravity model.

Gravity model: effect of gravity components in model performance. As described
earlier, the gravity component is expected to capture the difference in the
local behaviour of the social and mobility aspects of the populations. Air-
traffic has been considered as a proxy variable, since data on other possible
covariates are not available publicly/systematically for modeling purposes.
One may raise the issue that since life has returned to near-normalcy, whether
these variables can still be useful for modeling COVID-19. In Table 1, we have
summarised the performance of a selected set of models. We can safely conclude
from these values that, these variables definitely contribute significantly in
the various model combinations applied on this data.

Table 1: Effect of gravity components in model fit, assessed by pseudo-$R^{2}$. | Gravity component
---|---
Other covariates in the model | Without | With
Regional level effects | 0.182 | 0.290
District level effects | 0.363 | 0.497
District level effects and seasonal components | 0.457 | 0.556
District level effects and endemic components | 0.826 | 0.997

Gravity model: gravity component with complete connection of sources and
destinations. Another relevant aspect in modeling the epidemic component is
whether it was possible to travel across the regional boundaries to a
neighboring region for air travel. We implemented a separate model which
allows for the inter-regional travel, leading to a gravity model with
completely connected source-destination network dampened only by distance
decay. We found the model estimates to be similar but the overall fit was
poorer. This is not surprising since for a considerable period of the
pandemic, various travel restrictions in each region had made cross-region
travel to airports quite infeasible for most people.

Gravity model: robustness of estimates against data source. The availability
of COVID-19 data for India has been volunteer driven. Thus there could be
skepticism whether our conclusions are data source specific. We have tried our
models on the popular/regular sources of COVID-19 data. We also extracted and
cleaned data from the MS-Bing website https://www.bing.com/covid/local/india.
Irrespective of the data source, the estimated parameters and patterns have
been extremely similar. Unfortunately since October 2021, data in the MS-Bing
resource has become erratic, and hence full comparison was not possible. The
partial results have not been shown here.

EE model: lag structure of regional level history. Before we arrived at an
estimate of $\nu_{it}$ or $\nu_{t}$ based on the EE model proposed by Celani
and Giudici (2021), it was essential to assess the lag effects of the time
series of regional COVID-19 counts on each other. For this purpose, we used
the daily COVID-19 occurrence data from each region, and implemented various
lag structures for the autoregressive and the spatial effects parts of the
models. We measured the effects of such changes by estimating the
pseudo-$R^{2}$ for each of these models. Figure 4 shows that there is little
impact of higher order lags on the EE model, when measured in terms of the
relative gain in Nagelkerke pseudo-$R^{2}$. Thus for subsequent
implementations, we used a lag of order 1 for both local and spatial
autoregression.

Figure 4: Nagelkerke pseudo-$R^{2}$ for regional endemic-epidemic type models
with daily count data, for different autoregression lag choices.

Endemic parameters: EE model with multi-region common parameters. In the same
EE model framework, instead of implementing the model for each region with
their own endemic parameters $\\{\nu_{it}\\}$, we implemented a combined model
with common endemic parameters $\\{\nu_{t}$}.

Endemic parameters: gravity model with robustness against implementation
level. Figure 5 shows the estimated $\\{\nu_{it}\\}$ parameters, based on
Gravity models applied to each region individually, and also the estimated
$\\{\nu_{t}\\}$ parameters under the combined model. We observe that these
parameter estimates show a strongly similar pattern over time, irrespective of
the wide difference in the behaviour of the pandemic in different regions of
India.

Figure 5: Endemic parameter estimates under Gravity model for individual
regions (labeled by region names), and all regions combined (labeled as
’Combined model’).

Similar comparison between the common $\\{\nu_{t}\\}$ parameters estimated
under the two models, namely the EE model of Held et al. (2005) (see Adegboye
and Adegboye (2017)) and the Gravity models, are quite similar in nature (see
Figure 6).

Endemic parameters: consistency under different model choices. In Figure 6 we
present a comparative plot of the estimated endemic parameters, assumed to be
common across the regions of India, estimated using the EE model as well as
the proposed model (3.4). It is to be noted that the EE model has been
implemented on the daily series, whereas the Gravity model utilized the
monthly data. Thus the unit of time $t$ for $\nu_{t}$ are different for the
two implementations. After a thorough comparison, it was felt that the mid-
month estimates of $\nu_{t}$ based on EE-model applied to daily data is best
for comparison with the estimates of the same from the Gravity model applied
to monthly data. The estimates from the EE model and the Gravity model are
shown in Figure 6, and it can be seen that they concur and show noticeable
similarities. Thus, the estimates of the endemic parameters are persistently
similar under different models/implementation and data sources. This enhances
our confidence in the model and the estimation procedure.

Figure 6: Endemic parameter estimates under EE model and Gravity model based
on multi region combined data.

Endemic parameters: assessment of seasonality. One of the natural indicators
of whether the epidemic has reached an endemic phase would be a strong
seasonal component. To assess this, we modified the models (3.3) and (3.4),
where instead of a time point specific $\nu_{t}$ parameter, we implement a
month specific set of parameters. If the endemic phase has indeed been
reached, then this modified model would be able to achieve the same modeling
precision with lesser number of parameters. In Figure 7 we have presented the
estimated effects of months on COVID-19 log-mean counts, based on (i) the
Gravity model for monthly data with seasonal parameters, (ii) the Gravity
model based on monthly data with time point specific parameters, and (iii) the
EE model based on daily data with daily parameters. We also applied a moving
average of order 2 on the seasonal estimates, which provides a reasonable fit.
The estimates of months from the Gravity and EE models are arrived at by
combining the estimates from the respective months. From this figure we notice
a consistency of pattern under the three modeling frameworks.

Figure 7: Seasonal endemic parameter estimates under various models.

The pseudo-$R^{2}$ values are listed in Table 2. These values show that while
all three modeling attempts arrive at a similar pattern, the seasonal model
performs far worse than the models (3.3) and (3.4). This indicates that
reduction in the number of parameters within the endemic part of the model is
not yet possible. In other words, the endemic phase of COVID-19 has not been
reached in India yet.

Table 2: Model fit assessment for selected models using Nagelkerke pseudo-$R^{2}$. Gravity model | EE model
---|---
Seasonal | Non-seasonal | Non-seasonal
0.556 | 0.997 | 0.878

##  5 Conclusions

The contributions of the Gravity type models proposed here are twofold. First,
we used a novel set of covariates that account for temporal dynamic, latent
effects and other covariates. This enhances our understanding of the spread of
COVID-19 at a local level. Second, we have built this model within a state-of-
the-art regression framework that allows ease of implementation for other
applications.

The extraneous covariates yielded an unexpected utility. Questions about the
nature of collected and reported COVID-19 data have persisted over time. Many
news articles have questioned what appeared to be the systematic under-
reporting of COVID-19 cases from certain regions. Since the models presented
here use externally collected covariates, these departures/lack-of-
compatibilities within the data from temporal and/or spatial perspective
become clearly visible. A detailed reporting of such findings are beyond the
scope of this article and hence have been omitted.

The underlying spatial pattern of COVID-19 is expected to be reasonably
smooth, although local mobility and social behaviour can strongly influence
the outcome. For most of the regions, this hypothesis is supported by plots of
the estimated coefficients for the district level effects on the map. There is
an evident degree of similarity in the underlying distributions of the
district level parameters. Additionally, these district level parameters
appear to be spatially smooth. In this regard the parameterization in our
model overlaps with that proposed by Alipour et al. (2021), although their
focus was to estimate the effect of “working from home”.

This Gravity model is in spirit similar to the surveillance model introduced
by Held et al. (2005). Their extension of the generalized linear models to
analyze surveillance data from epidemic outbreaks was further expanded to
handle multivariate surveillance data by Paul et al. (2008). Useful
modifications to account for seasonality and spatial heterogeneity was
proposed in Held and Paul (2012), and neighborhood information from social
contact data was included in Meyer and Held (2017). It is agreed that these
models have strong mathematical foundation.

However, the proposed Gravity type model has novel explanatory features in
addition to being mathematically rigorous. Our gravity model uses a part of
the socio-economic behaviour in explaining the occurrence of cases, instead of
using the evolution of disease process solely. This would be critical in
situations where disease data collection mechanism is erratic, as is the case
for India. In the end, the gravity component of the model successfully
captures the local fluctuations, and the final fit to the data is extremely
well. Berlamann and Haustein (2020) had concluded that their disaggregated
data showed that there is quite some variety in the relative importance of the
endemic, the autoregressive epidemic, and the spatial endemic components. In
our case, the endemic component is seen to be systematically similar across
various regions of India, and the main variations are observed to come from
local sources.

Further modifications to our models could be explored but would require
additional investigation. In one of our earlier works, we had employed novel
clustering algorithms that explored and captured similarities between
districts that are not necessarily spatially adjacent. Some possible
explanations of such behaviour would be the shared latent non-spatial
variables like climate, transport, and local culture. This could have been
exploited in our current model too. In the literature, for example in D’Urso
et al. (2019) and D’Urso and Vitale (2020), such techniques have been used on
different proximity dimensions, such as the geographical and social space, to
identify similar districts while taking into account spatial dependencies.
Spatial dependencies have also been directly incorporated in the correlation
structure in the context of spatial econometric models. See for example LeSage
and Pace (2009). This could also be a possible avenue to explore in extending
these models.

## References

  * Adegboye and Adegboye (2017) O. Adegboye and M. Adegboye.  Spatially correlated time series and ecological niche analysis of cutaneous leishmaniasis in Afghanistan.  _Int. J. Env. Res. Public Health_ , 14(3):309, 2017.  doi: 10.3390/ijerph14030309.
  * Alipour et al. (2021) J. Alipour, H. Fadinger, and J. Schymik.  My home is my castle: the benefits of working from home during a pandemic crisis. Evidence from Germany.  _Jour. Pub. Econ._ , 196:104373, 2021.  doi: 10.1016/j.jpubeco.2021.104373.
  * Anderson (2011) J. E. Anderson.  The gravity model.  _Annual Rev. in Econ._ , 3:133–160, 2011.  doi: 10.1146/annurev-economics-111809-125114.
  * Bauer and Wakefield (2018) C. Bauer and J. Wakefield.  Stratified space-time infectious disease modelling, with an application to hand, foot and mouth disease in china.  _Jour. Royal Stat. Soc. Ser. C, Applied Statistics_ , 67(5):1379–1398, 2018.  doi: 10.1111/rssc.12284.
  * Becker (2015) N. Becker.  _Modeling to Inform Infectious Disease Control_.  CRC Press, Boca Raton, 2015.
  * Berlamann and Haustein (2020) M. Berlamann and E. Haustein.  Right and yet wrong: a spatio-temporal evaluation of Germany’s covid-19 containment policy.  _CESifo Working Paper, 8446_ , 2020.  URL https://www.cesifo.org/DocDL/cesifo1_wp8446.pdf.
  * Bhattacharjee (2020) M. Bhattacharjee.  COVID-19 Website (plots, graphs and other resources).  _http://mathstat.uohyd.ac.in/people/profile/m-bhattacharjee_ , 2020\.
  * Bhattacharjee and Bose (2021) M. Bhattacharjee and A. Bose.  Modelling COVID-19 data–I: a dynamic SIR(D) model and application to Indian data.  _J. Indian Stat. Assoc._ , 58(2):113–133, 2021\.
  * Bhattacharjee et al. (2021) M. Bhattacharjee, K. Divya, and A. Bose.  Modelling COVID-19 data–II: spatio-temporal models with application to Kerala data.  _J. Indian Stat. Assoc._ , 1111(2222):3333, 2021\.
  * Bjórnstad et al. (2020) O. Bjórnstad, K. Shea, M. Krzywinski, and N. Altman.  The SEIRS model for infectious disease dynamics.  _Nature Methods_ , 17:557–558, 2020.  doi: 10.1038/s41592-020-0856-2.  Erratum: 18, 321, 2021.
  * Bracher and Held (2020) J. Bracher and L. Held.  Endemic-epidemic models with discrete-time serial interval distributions for infectious disease prediction.  _Int. J. Forecast._ , 38(3):1221–1233, 2020.  doi: 10.1016/j.ijforecast.2020.07.002.
  * Brauer et al. (2019) F. Brauer, C. Castillo-Chavez, and Z. Feng.  Endemic disease models.  _Mathematical Models in Epidemiology_ , 69:63–116, 2019\.  doi: 10.1007/978-1-4939-9828-9˙3.
  * Brockmann et al. (2006) D. Brockmann, L. Hufnagel, and T. Geisel.  The scaling laws of human travel.  _Nature_ , 439(7075):462–465, 2006.  doi: 10.1038/nature04292.
  * Celani and Giudici (2021) A. Celani and P. Giudici.  Endemic-epidemic models to understand COVID-19 spatio-temporal evolution.  _Spatial Stat._ , 2021.  doi: 10.1016/j.spasta.2021.100528.
  * Cheng et al. (2016) Q. Cheng, X. Lu, J. Wu, Z. Liu, and J. Huang.  Analysis of heterogeneous dengue transmission in Guangdong in 2014 with multivariate time series model.  _Scientific Rep._ , 6:33755, 2016.  doi: 10.1038/srep33755.
  * Clif and Ord (1981) A. Clif and J. Ord.  _Spatial processes: Models and applications_.  Pion Limited, London, 1981.
  * Colón-González et al. (2018) F. Colón-González, I. Lake, R. Morbey, and et al.  A methodological framework for the evaluation of syndromic surveillance systems: a case study of England.  _BMC Public Health_ , 18:544, 2018.  doi: 10.1186/s12889-018-5422-9.
  * Cox (1981) D. R. Cox.  Statistical analysis of time series: Some recent developments.  _Scand. Jour. of Stat._ , 8(2):93–115, 1981.
  * Dunbar and Held (2020) M. Dunbar and L. Held.  Endemic-epidemic framework used in covid-19 modelling (discussion on the paper by nunes, caetano, antunes and dias).  _Revstat. Stat. J._ , 18:565–574, 2020.  URL https://www.ine.pt/revstat/pdf/REVSTAT_v18-n5-02.pdf.
  * D’Urso and Vitale (2020) P. D’Urso and V. Vitale.  A robust hierarchical clustering for georeferenced data.  _Spatial Statistics_ , 35:100407, 2020.  doi: 10.1016/j.spasta.2020.100407.
  * D’Urso et al. (2019) P. D’Urso, L. De Giovanni, M. Disegna, and R. Massari.  Fuzzy clustering with spatial–temporal information.  _Spat. Stat._ , 30:71–102, 2019.  doi: 10.1016/j.spasta.2019.03.002.
  * Fritz and Kauermann (2020) C. Fritz and G. Kauermann.  On the interplay of regional mobility, social connectedness, and the spread of covid-19 in Germany.  _arXiv, https://arxiv.org/abs/2008.03013_ , 2020.
  * Ganesan and Subramani (2020) S. Ganesan and D. Subramani.  IISc-model website.  _https://cmg.cds.iisc.ac.in/covid/_ , 2020.
  * Ganesan and Subramani (2021) S. Ganesan and D. Subramani.  Spatio-temporal predictive modeling framework for infectious disease spread.  _Nature Scientific Rep._ , 11(6741), 2021.  doi: 10.1038/s41598-021-86084-7.
  * Giuliani et al. (2020) D. Giuliani, M. Dickson, G. Espa, and F. Santi.  Modelling and predicting the spatio-temporal spread of coronavirus disease 2019 (COVID-19) in Italy.  _BMC Infectious Diseases_ , 20:700, 2020.  doi: 10.1186/s12879-020-05415-7.
  * Grimée et al. (2022) M. Grimée, M. Bekker-Nielsen Dunbar, F. Hofmann, and et al.  Modelling the effect of a border closure between switzerland and italy on the spatiotemporal spread of covid-19 in switzerland.  _Spat. Stat._ , 49(100552), 2022.  doi: 10.1016/j.spasta.2021.100552.
  * Held and Paul (2012) L. Held and M. Paul.  Modeling seasonality in space-time infectious disease surveillance data.  _Biometrical Journal._ , 54(6):824–843, 2012\.  doi: 10.1002/bimj.201200037.
  * Held et al. (2005) L. Held, M. Hòhle, and M. Hofmann.  A statistical framework for the analysis of multivariate infectious disease surveillance counts.  _Statistical Modelling_ , 5(3):187–199, 2005\.
  * Herzog et al. (2011) S. Herzog, M. Paul, and L. Held.  Heterogeneity in vaccination coverage explains the size and occurrence of measles epidemics in German surveillance data.  _Epidemiology and Infection_ , 139(4):505–515, 2011.  doi: 10.1017/S0950268810001664.
  * Kirby et al. (2017) R. S. Kirby, E. Delmelle, and J. M. Eberth.  Advances in spatial epidemiology and geographic information systems.  _Annals of Epidemiology_ , 279(1):1–9, 2017.  doi: 10.1016/j.annepidem.2016.12.001.
  * LeSage and Pace (2009) J. LeSage and K. Pace.  _Introduction to Spatial Econometrics_.  Chapman and Hall/CRC, London, 2009.
  * Meyer and Held (2014) S. Meyer and L. Held.  Power-law models for infectious disease spread.  _Ann. Appl. Stat._ , 8(23):1612–1639, 2014.  doi: 10.1214/14-AOAS743.
  * Meyer and Held (2017) S. Meyer and L. Held.  Incorporating social contact data in spatio-temporal models for infectious disease spread.  _Biostatistics_ , 18:338–351, 2017.
  * Meyer et al. (2017) S. Meyer, L. Held, and M. Hòhle.  Spatio-temporal analysis od epidemic phenomena using the R package surveillance.  _J. of Stat. Software_ , 77(11):1–55, 2017.  doi: 10.18637/jss.v077.i11.
  * Mollalo et al. (2020) A. Mollalo, B. Vahedi, and K. Rivera.  GIS-based spatial modeling of covid-19 incidence rate in the continental united states.  _Science of the Total Environment_ , 728:138884, 2020.  doi: 10.1016/j.scitotenv.2020.138884.
  * Paul et al. (2008) M. Paul, L. Held, and A. Toschke.  Multivariate modelling of infectious disease surveillance data.  _Stat. in Med._ , 27(29):6250–6267, 2008.  doi: 10.1002/sim.3440.
  * Phillips (2021) N. Phillips.  The coronavirus will become endemic.  _Nature_ , 590:382–384, 2021.  doi: 10.1038/d41586-021-00396-2.
  * Ranjan (2020) R. Ranjan.  COVID-19 spread in India: dynamics, modeling, and future projections.  _J. Indian Stat. Assoc._ , 58(1):47–65, 2020\.
  * Ravenstein (1889) E. G. Ravenstein.  The laws of migration.  _J. Royal Stat. Soc._ , 52(2):241–305, 1889.  doi: 10.2307/2979333.
  * Sannigrahi et al. (2020) S. Sannigrahi, F. Pilla, B. Basu, A. Sarkar Basu, and A. Molter.  The overall mortality caused by covid-19 in the European region is highly associated with demographic composition: A spatial regression-based approach.  _Sustainable Cities and Society_ , 62:102418, 2020.  doi: 10.1016/j.scs.2020.102418.
  * Tinbergen (1962) J. Tinbergen.  _Shaping the World Economy: suggestions for an International Economic Policy_.  The Twentieth Century Fund, New York, 1962.
  * Zipf (1946) G. K. Zipf.  The P1 P2/D hypothesis: on the intercity movement of persons.  _American Sociological Review_ , 11(6):677–686, 1946.  doi: 10.2307/2087063.