# Nonsense associations in Markov random fields with pairwise dependence

Sohom Bhattacharya  S. Bhattacharya Department of Statistics
University of Florida
Gainesville, FL 32611, USA<EMAIL_ADDRESS>,  Rajarshi Mukherjee  R.
Mukherjee Departments of Biostatistics
Harvard T.H. Chan School of Public Health
Boston, MA 02115, USA<EMAIL_ADDRESS>and  Elizabeth Ogburn  E.
Ogburn Departments of Biostatistics
Johns Hopkins Bloomberg School of Public Health
Baltimore, MD, 21205, USA<EMAIL_ADDRESS>
###### Abstract.

Yule (1926) identified the issue of “nonsense correlations” in time series
data, where dependence within each of two random vectors causes overdispersion
– i.e. variance inflation – for measures of dependence between the two. During
the near century since then, much has been written about nonsense correlations
– but nearly all of it confined to the time series literature. In this paper
we provide the first, to our knowledge, rigorous study of this phenomenon for
more general forms of (positive) dependence, specifically for Markov random
fields on lattices and graphs. We consider both binary and continuous random
vectors and three different measures of association: correlation, covariance,
and the ordinary least squares coefficient that results from projecting one
random vector onto the other. In some settings we find variance inflation
consistent with Yule’s nonsense correlation. However, surprisingly, we also
find variance deflation in some settings, and in others the variance is
unchanged under dependence. Perhaps most notably, we find general conditions
under which OLS inference that ignores dependence is valid despite positive
dependence in the regression errors, contradicting the presentation of OLS in
countless textbooks and courses.

##  1\. Introduction

It is well known that two independent, non-stationary time series often appear
to be correlated with one another and that tests of independence have inflated
type-I error [49, 26]; this phenomenon is known variously as nonsense,
spurious, or volatile correlations and has been the subject of extensive
research. Less well-known is that related phenomena can occur in many other
dependent data settings, e.g. dependence due to spatial, batch, population,
genetic, or network structure, and that similar phenomena can occur for
different measures of association. We will use the term association to be more
general than correlation and to encompass correlation coefficients, estimated
covariances, and regression coefficients as measures of linear non-
independence.

The problem of nonsense correlations in time series was first articulated by
[60], building upon earlier works by [32, 61, 59]. Much work has been on pre-
whitening data, e.g. extracting independent increments, before attempting to
learn about associations between time series. [26] derived the exact
asymptotic distribution of the estimated correlation of two independent Weiner
processes, showing overdispersion that does not vanish asymptotically and
“solving” the problem first put forth by [60].

However, outside of the time series econometrics literature, this idea seems
to be largely unknown. We have been able to find a smattering of references in
the broader literature: of [37, Section 12.2.3] illustrates nonsense
associations due to temporal autocorrelation in linear regression. [31]
discussed the importance of dealing with nonsense correlations in the study of
neurophysiological processes in neuroscience. [52] noted that natural
population processes (e.g. population growth) are susceptible to nonsense
correlations. Moving beyond time series data, [61] first noted that, in
analogy to temporal processes, spatial processes are also susceptible to more
highly variable estimates of correlation, and that the problem there can
similarly be solved by extracting independent increments. [13] and [23]
derived the effective sample size for tests of correlation between two spatial
processes under strong parametric assumptions. Section 12.4 of 24 includes a
brief discussion of ephemeral predictions, which is a manifestation of
nonsense associations in the prediction setting. Recently [41] made a case for
the importance of nonsense associations in social network data. [5] studied
the power and type I error of a family of tests of association commonly used
in genome-wide association studies (GWAS), where the independent (e.g.
phenotype) and dependent (e.g. SNP) variables both exhibit dependence due to
family or population structure; they showed that standard tests of association
can have type I error bounded away from the nominal level or even converging
to 1.

In addition to statistical papers, a handful of papers from the social
sciences have pointed out phenomena that appear to be instances of nonsense
associations, though none of these papers make the connection explicitly:
“Galton’s problem,” a well-known concept in cross-cultural studies and
anthropology, refers to the fact that within-culture dependence makes it
difficult to perform inference using cross-cultural data, and at least one
group of authors have connected this to network dependence [20]. [17]
describes how incidental clustering of beliefs or behaviors in social networks
can look like systematic associations. [38] shows that many findings in the
economics literature on “persistence” or “deep origins,” purporting to find
connections between characteristics of a place in modern day and
characteristics of the same place in the distant past, can be explained by
spatial autocorrelation rather than any meaningful association between past
and present characteristics.

MRFs have been used crucially for modeling spatially correlated data sets,
which appear naturally in population genetics [3, 56], image segmentation [9,
12], climate data modeling [30, 50]. In this research in conditional
independence and structure testing has gained popularity recently [45, 53], it
is worth noting that the literature on testing for the independence of two
MRFs typically assumes multiple independent observations from the MRF
distribution, with the number of observations going to infinity; in this paper
we are concerned with the depenedent data setting where there is only one draw
from each of two MRF distributions with the dimension going to infinity.

In this paper, we derive the first, to our knowledge, characterization of
nonsense associations for a class of network-dependent data. We consider the
asymptotic distributions of estimates and tests of association between two
random vectors, $\mathbf{X}$ and $\mathbf{Y}$, where the coordinates of
$\mathbf{X}$ and of $\mathbf{Y}$ exhibit dependence but $\mathbf{X}$ is
independent of $\mathbf{Y}$. We focus on dependence that is governed by graph
or network structure and generated from a Markov random field (MRF) model with
pairwise positive dependence. We will demonstrate that under some MRF models,
the asymptotic variance of estimates of correlation or covariance is inflated
relative to the i.i.d. setting; this is analogous to Yule’s nonsense
correlations and would result in inflated type I error were i.i.d. tests to be
used without correction. However, in other settings dependence results in
deflated asymptotic variance; this is in stark contrast to Yule’s finding and
contradicts the heuristic that positive dependence results in smaller
effective sample sizes and therefore higher variances.

In a particularly surprising result, we challenge the common wisdom that
positive dependence in the errors of a simple linear regression model should
result in standard i.i.d. inference for the regression coefficient being
anticonservative. In fact, we show that in order for i.i.d. inferential
procedures to be invalid, the variance-covariance matrices of the dependent
and the independent variables must share some structure. Otherwise, e.g. when
the independent variable is itself i.i.d., i.i.d. inferential procedures
achieve nominal type-I error rates. This contradicts the claim, presented in
countless regression classes and textbooks, that correlated regression errors
should be handled with modified methods in order for inference to be valid
[55, 29, 10].

##  2\. Setting

Consider two independent random vectors observed on $n$ nodes in a network,
$\mathbf{X}=(X_{1},\ldots,X_{n})$ and $\mathbf{Y}=(Y_{1},\ldots,Y_{n})$, where
the nodes exhibit positive pairwise dependence whenever they are connected by
a network tie. More precisely, for $\mathcal{Z}\subseteq\mathbb{R}$, we assume
$\mathbf{X},\mathbf{Y}\in\mathcal{Z}^{n}$ are independently distributed with
density of the form

$\mathbb{P}_{\mathbf{Q}_{n}}(z)\propto\exp{\left(\mathbf{z}^{\top}\mathbf{Q}_{n}\mathbf{z}\right)},\quad
z\in\mathcal{Z}^{n}$ (2.1)

where $\mathbf{Q}_{n}$ is some $n$ by $n$ matrix related to the adjacency
matrix of the underlying network, which can be different for
$\mathbf{X},\mathbf{Y}$. We drop the subscript $n$ when context allows it.
Such families are typical examples of Markov Random Fields (MRFs) which are a
natural class of models for dependent collections of random variables and have
been used extensively across an array of disciplines to explore statistical
inference in such dependent systems, including to model social networks (see,
e.g., 57, 19, 1, 39, 18, 14, 15, 36, 47). The family (2.1) includes two
classical examples which we will study in this paper: Gaussian graphical
models corresponding to $\mathcal{X}=\mathbb{R}$ and Ising models
corresponding to $\mathcal{X}=\pm 1$. We will extensively focus on these two
cases since the distribution of the statistics we aim to study here remains
asymptotically standard normal under no dependence within the coordinates of
$\mathbf{X}$ or $\mathbf{Y}$.

Throughout we consider a single pair $(\mathbf{X},\mathbf{Y})$ with $n$ going
to infinity and inference on the distribution of the $n$ pairs of observations
$(X_{i},Y_{i})$ with dependence across units. This distinguishes our framework
from the more typical framework in which a draw from an MRF represents one
single, high-dimensional observation and inference would typically require
multiple such observations. We consider three quantities of interest. The
sample covariance $T_{n}$ is defined by

$T_{n}:=\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\overline{\mathbf{X}}_{n})(Y_{i}-\overline{\mathbf{Y}}_{n}),$
(2.2)

where $\overline{\mathbf{X}}_{n}:=\frac{1}{n}\sum_{i=1}^{n}X_{i}$,
$\overline{\mathbf{Y}}_{n}:=\frac{1}{n}\sum_{i=1}^{n}Y_{i}$ are the sample
means. We will also consider the sample correlation

$\rho_{n}:=\frac{\sum_{i=1}^{n}(X_{i}-\overline{\mathbf{X}}_{n})(Y_{i}-\overline{\mathbf{Y}}_{n})}{\sqrt{\sum_{i=1}^{n}(X_{i}-\overline{\mathbf{X}}_{n})^{2}}\sqrt{\sum_{i=1}^{n}(Y_{i}-\overline{\mathbf{Y}}_{n})^{2}}}.$
(2.3)

Although tests of independence between $\mathbf{X}$ and $\mathbf{Y}$ can be
done using the sample covariance, it is of interest to note when the
asymptotic behavior of $\rho_{n}$ differs from that of $T_{n}$. Finally, for
continuous data we also consider $\widehat{\beta}$, the estimated coefficient
for $X$ in a simple linear regression of $Y$ on $X$.

In all of the settings we consider below, it is easy to show that the
asymptotic null distribution of both $T_{n}$ and $\rho_{n}$ under an i.i.d.
data-generating process is $N(0,1)$. Similarly, when both $\mathbf{X}$ and
$\mathbf{Y}$ are i.i.d., $\widehat{\beta}$ has a known normal asymptotic
distribution. Unlike the time series setting, where the null distribution of
$\rho_{n}$ under dependence generated by Weiner processes is not generally
Gaussian [26], we find that in many settings with MRF-generated dependence,
Gaussianity is preserved for various measures of association. Even in these
settings, we find that the asymptotic variance of $\rho_{n}$ may be inflated
or deflated relative to the nominal value of $1$. When variance is inflated,
the type I error of any procedure that ignores dependence will be inflated,
which is another hallmark of nonsense associations.

Proofs of all new results are in Section 6.

##  3\. Results for binary data

Let $\mathbf{Z}=(z_{1},\ldots,z_{n})^{\top}\in\\{\pm 1\\}^{n}$ be a random
vector with the joint distribution of $\mathbf{Z}$ given by an Ising model
defined as:

$\displaystyle\mathbb{P}_{\beta,\mathbf{Q}_{n}}(\mathbf{Z}=\mathbf{z}):=\frac{1}{Z(\beta,\mathbf{Q}_{n})}\exp{\left(\frac{\beta}{2}\mathbf{z}^{\top}\mathbf{Q}_{n}\mathbf{z}\right)},\qquad\forall\mathbf{z}\in\\{\pm
1\\}^{n},$ (3.1)

where $\mathbf{Q}_{n}$ is an $n\times n$ symmetric and hollow matrix,
$\beta\geq 0$ is a non-negative real number usually referred to as the
“inverse temperature,” and $Z(\beta,\mathbf{Q}_{n})$ is a normalizing
constant. The pair $(\beta,\mathbf{Q}_{n})$ characterizes the dependence among
the coordinates of $\mathbf{X}$, and $X_{i}$’s are independent if and only if
$\beta\mathbf{Q}_{n}=\mathbf{0}_{n\times n}$. The matrix $\mathbf{Q}_{n}$ will
be associated with a certain labeled graph: Let $\mathbb{G}_{n}=(V_{n},E_{n})$
with vertex set $V_{n}=\\{1,\dots,n\\}$ and undirected edge set
$E_{n}\subseteq V_{n}\times V_{n}$, then the corresponding
$\mathbf{Q}_{n}=nG_{n}/2|E_{n}|$, where $G_{n}$ is the adjacency matrix of
$\mathbb{G}_{n}$.

First, we present results for sparse graphs, specifically for lattices, where
we see that (in high temperature regimes) both the sample correlation and
sample covariance exhibit variance inflation. This result generalizes nonsense
associations from time series, where dependence is informed by distance in
$\mathbb{R}^{1}$, to higher dimensional dependence such as would be found in
spatial or network data, and we believe to be the first such generalization
beyond one-dimensional dependence. We suspect these results will hold for
other connected, bounded-neighborhood graphs, but we leave the proof to future
work.

We next present results for complete and dense (i.e., $E_{n}=\Theta(n^{2})$)
regular graphs. This includes some specific examples of mean-field type
models, see e.g. [4, 34]. In a surprising contrast to other settings,
dependence in these dense graphs sometimes results in variance deflation for
the sample covariance, while the sample correlation has the same asymptotic
distribution as under independence.

###  3.1. Lattice models

Let the data points $i=[n]$ be vertices of $d$-dimensional hypercubic lattice
$\Lambda_{n,d}$ and the underlying graph $\mathbb{G}_{n}$ be the nearest
neighbor (in sense of Euclidean distance) graph on these vertices. This
nearest neighbor graph corresponds to the matrix
$\mathbf{Q}_{ij}=\mathbf{1}(\|i-j\|_{1}=1)$.

The lattice Ising model typically undergoes a “thermodynamic” phase transition
at a value of the temperature parameter that is determined by the dimension
$d$ of the lattice: there exists $\beta_{c}:=\beta_{c}(d)$ which governs the
macroscopic behavior of the system of observations $X_{i},i\in\Lambda_{n,d}$
as $n\rightarrow\infty$. In particular, the sample average
$\overline{\mathbf{X}}_{n}=n^{-1}\sum X_{i}$ converges to $0$ in probability
for $\beta<\beta_{c}(d)$ and to an equal mixture of two delta-dirac random
variables $m_{+}(\beta)$ and $m_{-}(\beta)=-m_{+}(\beta)$, for
$\beta>\beta_{c}(d)$ ([40]). The value of $\beta_{c}(d)$ (which is known to be
strictly positive for any fixed $d\geq 1$) is hard to derive precisely. For
$d=1$, it is known from the first work in this area [33] that
$\beta_{c}(1)=+\infty$ and consequently the Ising model in 1-dimension is said
to have no phase transitions. The seminal work of [48] derived
$\beta_{c}(2)=2\log(1+\sqrt{2})$. Obtaining an analytical formula for
$\beta_{c}(d)$ for $d\geq 3$ remains open. We refer the interested reader to
the excellent expositions in [27, 21] for more details on this topic.

Now, we state the main result of this section.

######  Theorem 1.

Suppose $\mathbf{X}\sim\mathbb{P}_{\beta_{1},\mathbf{Q}(\Lambda_{n,d})}$,
$\mathbf{Y}\sim\mathbb{P}_{\beta_{2},\mathbf{Q}(\Lambda_{n,d})}$ for some
$\beta_{1},\beta_{2}<\beta_{c}(d)$ with $\mathbf{X}$ independent of
$\mathbf{Y}$. Then, there exists $v=v(\beta_{1},\beta_{2})\in\mathbb{R}^{+}$
such that $\sqrt{n}T_{n}\xrightarrow{d}N(0,v^{2})$ and the asymptotic variance
satisfies $v^{2}>1$. The sample correlation has the same asymptotic
distribution as the sample covariance:
$\sqrt{n}\rho_{n}\xrightarrow{d}N(0,v^{2})$. Moreover,
$v(\beta_{1},\beta_{2})$ is a non-decreasing function in either coordinate.

Obtaining the precise limiting constant $v^{2}$ in Theorem 1 is an important
question and we leave it as future research direction. Since $v^{2}>1$, our
result demonstrates that both covariance and correlation exhibit variance
inflation. Moreover, in Section 5 we present evidence from simulations
suggesting that $v^{2}$ is strictly increasing in $\beta_{1},\beta_{2}$ for
$0<\beta_{1},\beta_{2}<\beta_{c}$. This result indicates that a researcher
ignoring the dependence in $\mathbf{X}$ and $\mathbf{Y}$ would perform
anticonservative inference, with type I error increasing as a function of
dependence; this is a direct analogue to nonsense correlations in time series.

###  3.2. Complete and dense graphs

The Ising model on the complete graph, where $E_{n}=V_{n}\times V_{n}$, is
known as the Curie-Weiss model [35, 46]. This model corresponds to having
$\mathbf{Q}_{n}:=\mathbf{Q}^{\rm
CW}=\frac{1}{n}(\mathbf{1}\mathbf{1}^{\top}-I)$. Similar to the lattice Ising
model, the Curie-Weiss model also undergoes a phase transition at
$\beta_{c}=1$. Namely, $\overline{\mathbf{X}}_{n}\xrightarrow[]{p}0$ for
$\beta\leq 1$ and
$\overline{\mathbf{X}}_{n}\xrightarrow{d}\frac{1}{2}(\delta_{m}+\delta_{-m})$
if $\beta>1$, i.e. an equal mixture of two Dirac delta random variables with
mass on $\pm m$ where $m=m(\beta)$ is the unique positive root of the equation
$m=\tanh(\beta m)$. Under this model, we precisely characterize the asymptotic
distribution of sample covariance $T_{n}$ and sample correlation $\rho_{n}$.

######  Theorem 2.

Suppose $\mathbf{X}\sim\mathbb{P}_{\beta_{1},\mathbf{Q}^{\text{CW}}}$,
$\mathbf{Y}\sim\mathbb{P}_{\beta_{2},\mathbf{Q}^{\text{CW}}}$ and $\mathbf{X}$
independent of $\mathbf{Y}$. Then

$\sqrt{n}T_{n}\xrightarrow{d}N(0,(1-m^{2}_{1})(1-m^{2}_{2})),$ (3.2)

where if $\beta_{i}\leq 1$, then $m_{i}=0$ and otherwise $m_{i}$ is the unique
positive root of the equation $m=\tanh(\beta_{i}m)$, $i=1,2$. Noting that
$m_{i}$ must lie in $[-1,1]$, this implies
$\lim\limits_{n\rightarrow\infty}\text{Var}(\sqrt{n}T_{n})=1$ iff
$\beta_{1},\beta_{2}\leq 1$, otherwise
$\lim\limits_{n\rightarrow\infty}\text{Var}(\sqrt{n}T_{n})<1$.

The implication of the above theorem is two-fold: first, it shows that the
distibution of (scaled) sample covariance is asymptotically normal
irrespective of the choice of $\beta_{1},\beta_{2}$. Second, it characterizes
asymptotic variance as a function of $\beta$ and shows a phase transition at
$\beta_{i}=1$, $i=1,2$. In fact, there is a variance deflation unless
$\max\\{\beta_{1},\beta_{2}\\}\leq 1$. To illustrate our findings, we plot the
asymptotic variance in Figure 1 when $\beta_{1}=\beta_{2}$.

Figure 1. Asymptotic variance of $\sqrt{n}T_{n}$ when
$\beta_{1}=\beta_{2}=\beta$.

The proof of Theorem 2 shows how the result can be generalized to a large
family of graphs as long as one of the variables is sampled from a Curie-Weiss
model. More precisely, if $\mathbf{X}\sim\mathbb{P}_{\beta_{1},\mathbf{Q}}$
satisfies $\overline{\mathbf{X}}^{2}_{n}\xrightarrow{\mathbb{P}}m^{2}_{1}$ and
$\mathbf{Y}\sim\mathbb{P}_{\beta_{2},\mathbf{Q}^{\text{CW}}}$ independent of
$\mathbf{X}$, then
$\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\overline{\mathbf{X}}_{n})^{2}\xrightarrow{\mathbb{P}}1-m^{2}_{1}$
and the conclusion of Theorem 2 still holds. Example of such interaction
matrices $\mathbf{Q}$ are scaled adjacency matrices of dense mean field graphs
(cf. [16, Theorem 1.1]) which includes Erdős-Rényi graphs, balanced stochastic
block models, sparse regular graphons, block spin Ising models and sparse
graphs like lattice Ising models.

A consequence of Theorem 2 is that the asymptotic distribution of sample
correlation $\rho_{n}$ is independent of the choice of $\beta_{1},\beta_{2}$.
The result is stated via the following corollary:

######  Corollary 3.

Suppose $\mathbf{X}\sim\mathbb{P}_{\beta_{1},\mathbf{Q}^{\text{CW}}}$,
$\mathbf{Y}\sim\mathbb{P}_{\beta_{2},\mathbf{Q}^{\text{CW}}}$ and
$\mathbf{X},\mathbf{Y}$ are independent. Then

$\sqrt{n}\rho_{n}\xrightarrow{d}N(0,1),$ (3.3)

for any $\beta_{1},\beta_{2}\geq 0$.

The asymptotic variance of sample correlation remains equal to $1$ for any
values of $\beta_{1},\beta_{2}$. It is notable that in this setting, unlike
the lattice setting, the asymptotic variances of the sample covariance and
correlation differ in nature. A practitioner unaware of dependence in their
data could still perform valid inference about the association between $X$ and
$Y$ using the sample correlation, but not necessarily using the sample
covariance.

These results for the Curie-Weiss model can be extended to mean-field type
models (cf. [4, 34]) for dense regular graphs. We make the following
assumptions on the dependency matrix $\mathbf{Q}_{n}$:

######  Assumption 1.

$\mathbf{Q}_{n}$ is a regular matrix, i.e., for all $i=1,2,\ldots,n$, we have
$\sum_{j=1}^{n}Q_{n,ij}=1$.

######  Assumption 2.

There exists a constant $\kappa>0$ such that

$\max_{1\leq i,j\leq n}Q_{n,ij}<\frac{\kappa}{n}.$

######  Assumption 3.

There exists a constant $\gamma\in\mathbb{R}$ such that

$\|\mathbf{Q}_{n}\|^{2}_{F}=\sum_{i,j=1}^{n}Q_{n,ij}\rightarrow\gamma.$

We will need to assume the sequence of matrices $\\{\mathbf{Q}_{n}\\}$
converge in cut-metric. To this end, we have to introduce the basics of
graphons and cut-distance convergence and refer the curious readers to [7, 8,
42].

A graphon is a symmetric bounded measurable function from $[0,1]^{2}$ to
$[0,1]$. Let $\mathcal{W}$ be the set of all graphons. One can equip the space
of graphons $\mathcal{W}$ by the cut-distance, given by

$d_{\square}(W_{1},W_{2}):=\sup_{S,T\subseteq[0,1]}\Big{|}\int_{S\times
T}(W_{1}(x,y)-W_{2}(x,y))dxdy\Big{|},$

where the supremum is taken over all measurable subsets $S,T$ of $[0,1]$.
Define the weak cut distance by

$\delta_{\square}(W_{1},W_{2}):=\inf_{\sigma}(W^{\sigma}_{1},W_{2})=\inf_{\sigma}(W_{1},W^{\sigma}_{2})$

where $\sigma$ ranges from all measure preserving bijections
$[0,1]\rightarrow[0,1]$ and $W^{\sigma}(x,y)=W(\sigma(x),\sigma(y))$.

Given an $n\times n$ matrix $A_{n}$, we can define the corresponding graphon
$W^{A_{n}}$ as a piecewise constant function given by

$W^{A_{n}}(x,y)=A_{n,\lceil nx\rceil,\lceil ny\rceil}.$

We call a set of matrices $A_{n}$, $n\in{N}$ converge in cut distance to a
graphon $W\in\mathcal{W}$ if $\delta_{\square}(W^{A_{n}},W)\rightarrow 0$.

Finally, given $W\in\mathcal{W}$, define its Hilbert-Schmidt operator $T_{W}$
from $[0,1]$ to $[0,1]$ given by

$T_{W}(f)(x)=\int_{[0,1]}W(x,y)f(y)dy,\quad x\in[0,1].$

Since the operator $T_{W}$ is a compact operator, it is has countably many
eigenvalues, which we denote by $\lambda_{1},\lambda_{2},\ldots$ sorted by
decreasing order in absolute values
$|\lambda_{1}|\geq|\lambda_{2}|\geq\ldots$.

Now, we are in a position to state our final assumption on the sequence of
matrices $\\{\mathbf{Q}_{n}\\}$:

######  Assumption 4.

The sequence of matrices $\\{n\mathbf{Q}_{n}\\}$ converges in cut distance to
a graphon $W$ (i.e., $\delta_{\square}(W^{n\mathbf{Q}_{n}},W)\rightarrow 0$).
Moreover, if $T_{W}$ has eigenvalues $\lambda_{1},\lambda_{2},\ldots$, we have
$\lambda_{1}>\sup_{i\geq 2}\lambda_{i}$.

Let us briefly discuss the spectral gap condition we have assumed here. By our
regularity Assumption 1, the largest eigenvalue of $\mathbf{Q}_{n}$ is exactly
$1$. This, along with [42, Theorem 11.54] yields that $\lambda_{1}=1$. Hence,
we are assuming $\sup_{i\geq 2}\lambda_{i}<1$. Since $|\lambda_{i}|\leq 1$, we
allow the possibility that $\lambda_{2}=-1$, which can happen for regular
connected bipartite graphs. One can expect any universality result regarding
the asymptotic distribution of sample covariance might fail if one does not
assume such a spectral gap condition.

Before we mention our main result of this section, we briefly provide some
examples of matrices $\mathbf{Q}_{n}$ which satisfy our Assumptions 1-4. The
examples arise as the scaled adjacency matrices of regular graphs. First, note
that the previously discussed Curie-Weiss model,
$\mathbf{Q}_{n}=\frac{1}{n}(\mathbf{1}\mathbf{1}^{\top}-I)$ satisfies all the
assumptions. This is because the corresponding graphons $f^{\mathbf{Q}_{n}}$
converges in cut distance to the constant graphon $W\equiv 1$, yielding
$\lambda_{1}=1$ and $\lambda_{i}=0$ for $i\geq 2$. Second, the scaled
adjacency matrix of the complete bipartite graph also satisfies the
Assumptions 1-4. To see this, note that,

$Q_{n,ij}=\begin{cases}\frac{2}{n}\text{ if }1\leq i\leq\frac{n}{2}\text{ and
}\frac{n}{2}+1\leq j\leq n\\\ \frac{2}{n}\text{ if }1\leq
j\leq\frac{n}{2}\text{ and }\frac{n}{2}+1\leq j\leq n\\\ 0\text{
otherwise.}\end{cases}$

This is the (scaled) adjacency matrix of a regular graph with degree
$\frac{n}{2}$. Also, $n\mathbf{Q}_{n}$ converges, in cut norm to the limiting
graphon $W$ given by

$W(x,y)=\begin{cases}2\text{ if }0<x<\frac{1}{2}\text{ and
}\frac{1}{2}<y<1,\\\ 2\text{ if }\frac{1}{2}<x<1\text{ and
}0<y<\frac{1}{2},\\\ 0\text{ otherwise}.\end{cases}$

This limiting graphon $W$ has two non-zero eigenvalues
$\lambda_{1}=1,\lambda_{2}=-1$ satisying Assumption 4. In fact, our
assumptions continue to hold true for general complete multipartite graphs.
For our final example, we consider random $d_{n}$ regular graphs, where
$\frac{d_{n}}{n}\rightarrow d\in(0,1)$ and
$Q_{n,ij}=\frac{1}{d_{n}}\mathbf{1}_{(i,j)\in E_{n}}$. By definition,
$\mathbf{Q}_{n}$ is regular, with $n\mathbf{Q}_{n}$ converging in cut-norm to
the constant graphon $W\equiv 1$ [11, Theorem 1.1], hence the only non-zero
eigenvalue is $\lambda_{1}=1$.

Note that the Ising model with interaction matrix $\mathbf{Q}_{n}$ satisfying
Assumptions 1-4 satisfies the exact same phase transition as the Curie-Weiss
model $\mathbb{P}_{\beta,\mathbf{Q}^{\rm CW}}$ at $\beta_{c}=1$. We now state
our universality result.

######  Theorem 4.

Let $\\{\mathbf{Q}_{n}\\}$ be a sequence of matrices satisfying Assumptions
1-4. Suppose $\mathbf{X}\sim\mathbb{P}_{\beta_{1},\mathbf{Q}_{n}}$ and
$\mathbf{Y}\sim\mathbb{P}_{\beta_{2},\mathbf{Q}_{n}}$, independent of
$\mathbf{X}$. Then
$\sqrt{n}T_{n}\xrightarrow{d}N(0,(1-m^{2}_{1})(1-m^{2}_{2}))$, where $m_{i}=0$
if $\beta_{i}\leq 1$, and otherwise $m_{i}$ is the unique positive root of the
equation $m=\tanh(\beta_{i}m)$, $i=1,2$. Moreover, the sample correlation has
asymptotic distribution $\sqrt{n}\rho_{n}\xrightarrow{d}N(0,1)$.

The above result shows that the asymptotic distribution of (scaled) sample
covariance is the same for all dense regular graphs (under Assumptions 1-4).
Therefore, the phenomenon we observed in the previous section, i.e., variance
deflation for sample covariance, is not limited to complete graphs. It is very
natural to expect such a result to continue to hold for other mean-field
graphs like dense Erdős-Rényi graphs. However, the current proof crucially
relies on the fact that $\mathbf{Q}_{n}$ is exactly regular and we leave any
such extension for future research.

These results illustrate an important and surprising dichotomy between dense
regular graphs and lattice Ising models. While on lattices we see the familiar
phenomenon of variance inflation for both sample correlation and sample
covariance, analogous to nonsense correlations in time series, the dense
regime exhibits a variance deflation for the sample covariance and no change
in distribution for sample correlation.

##  4\. Results for continuous data

We now turn our attention to Gaussian graphical models
$\mathbf{X},\mathbf{Y}\sim N(0,\Sigma)$, corresponding to
$\mathcal{X}=\mathbb{R}$ and $\mathbf{Q}_{n}=-\frac{1}{2}\Sigma^{-1}$ in model
(2.1).

In addition to studying the asymptotic behavior of $T_{n}$ and $\rho_{n}$
under this model, we also study the behavior of the coefficient $\beta$ for
$X$ in a simple linear regression of $Y$ on $X$, allowing $\mathbf{X}$ and
$\mathbf{Y}$ to have different covariance matrices: $\mathbf{Y}\sim
N(0,\Sigma_{Y})$ and $\mathbf{X}\sim N(0,\Sigma_{X})$. We find that,
surprisingly and contrary to dozens of presentations on OLS, inference about
$\beta$ need not be anti-conservative when the regression errors exhibit
positive dependence. If $\Sigma_{X}$ is diagonal, or if a particular condition
holds on the product of the eigenvalues of $\Sigma_{X}$ and $\Sigma_{Y}$, then
inference under OLS is exact, i.e. it achieves the nominal type I error rate.

We find conditions under which $T_{n}$ and $\rho_{n}$ are asymptotically
normal, and, in line with the phenomenon of nonsense correlations, we find
that $\rho_{n}$ exhibits variance inflation. However, surprisingly, under
these same conditions $T_{n}$ may exhibit either variance inflation or
deflation.

###  4.1. Correlation and Covariance

Let $\mathbf{X},\mathbf{Y}\sim N(0,\Sigma)$ and define
$J:=I_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{\top}$, where $\mathbf{1}$ is the
$n$ length vector of all $1$’s. The asymptotic distribution of $T_{n}$ and
$\rho_{n}$ will depend on the eigenvalues of
$\widetilde{\Sigma}:=\Sigma^{1/2}J\Sigma^{1/2}$. Let
$\widetilde{\lambda}_{1}\geq\widetilde{\lambda}_{2}\ldots\geq\widetilde{\lambda}_{n}\geq
0$ be the eigenvalues of $\widetilde{\Sigma}$ in non-increasing order. The
asymptotic distributions of $T_{n}$ and $\rho_{n}$ depend on whether the
eigenvalues are all of similar order. We formalize this notion in the
following Theorem.

######  Theorem 5.

Suppose $\mathbf{X},\mathbf{Y}\sim N(0,\Sigma)$ are independent. Define
$a_{n}:=\frac{\sum_{i=1}^{n}\widetilde{\lambda}_{i}}{\sqrt{n\sum_{i=1}^{n}\widetilde{\lambda}^{2}_{i}}}$.

  1. (i)

If $\widetilde{\lambda}^{2}_{1}\ll\sum_{i=1}^{n}\widetilde{\lambda}^{2}_{i}$,
then

$\frac{nT_{n}}{\sqrt{\sum_{i=1}^{n}\widetilde{\lambda}^{2}_{i}}}\xrightarrow{d}N(0,1)$

and

$\sqrt{n}\rho_{n}a_{n}\xrightarrow{d}N(0,1);$

hence, $\liminf_{n\rightarrow\infty}\text{Var}(\sqrt{n}\rho_{n})\geq 1$.

  2. (ii)

On the other hand, if $\widetilde{\lambda}_{1}\gg n\widetilde{\lambda}_{2}$,
then

$\frac{nT_{n}}{\widetilde{\lambda}_{1}}\xrightarrow{d}AB,$

where $A\sim N(0,1)$, $B\sim\sqrt{\chi^{2}_{1}}$ and $A,B$ are independent.
Moreover, $\rho_{n}$ converges weakly to a Rademacher distribution, i.e.,

$\rho_{n}\xrightarrow{d}\frac{1}{2}(\delta_{1}+\delta_{-1}).$

Since $a_{n}\leq 1$, the Theorem implies variance inflation for the sample
correlation. In fact, for any $\sigma^{2}\geq 1$, one can construct a sequence
of covariance matrix $\\{\Sigma_{n}\\}$ such that
$\sqrt{n}\rho_{n}\xrightarrow{d}N(0,\sigma^{2})$. To see this, set
$l_{0}=\lceil\frac{n}{\sigma^{2}}\rceil$ and set
$\widetilde{\lambda}_{1}=\ldots=\widetilde{\lambda}_{l}=1$ and
$\widetilde{\lambda}_{i}=0$ if $i>l_{0}$. Then
$a_{n}\rightarrow\frac{1}{\sigma}$. Hence, if
$\frac{1}{\sqrt{n}}\mathbf{1},v_{1},\ldots,v_{n-1}$ forms an orthonormal
basis, choosing
$\Sigma_{n}:=\frac{1}{n}\mathbf{1}\mathbf{1}^{\top}+\sum_{i=2}^{n}\widetilde{\lambda}_{i-1}v_{i}v^{\top}_{i}$
would yield $\sqrt{n}\rho_{n}\xrightarrow{d}N(0,\sigma^{2})$.

Although the sample correlation exhibits variance inflation, the sample
covariance can exhibit either variance inflation or deflation. Corollary 6
below provides an example of variance deflation; to see example of variance
inflation choose $\widetilde{\lambda}_{1}=\sqrt{n}\log n$ and
$\widetilde{\lambda}^{2}_{1}\ll\sum_{i=1}^{n}\widetilde{\lambda}^{2}_{i}$.
This will yield
$\lim\limits_{n\rightarrow\infty}\text{Var}(\sqrt{n}T_{n})\rightarrow\infty$.
In future work it will be of interest to characterize the settings in which
$T_{n}$ and $\rho_{n}$ show divergent asymptotic behavior, with the latter
exhibiting variance inflation while the former exhibits deflation.

Theorem 5 part (ii) shows more pathological behavior when one eigenvalue is of
larger order than the others. In this case the limiting distribution of
$T_{n}$ is not normal and the variance will typically be larger than $1$, and
the sample correlation converges to a bimodal Rademacher distribution centered
at, but concentrated away from, the true value of 0. These phenomena are
analogous to nonsense associations in other settings, where measures of
association have been shown to follow asymptotic distributions that are
overdispersed relative to the Normal distribution [26] or to be concentrated
away from the truth [41].

For the special case of equicorrelation we have the following Corollary to
Theorem 4:

######  Corollary 6.

Suppose $\mathbf{X},\mathbf{Y}\sim N(0,\Sigma)$ are independent, where
$\Sigma$ is an equicorrelation matrix, i.e.,
$\Sigma=(1-\rho)I+\rho\mathbf{1}\mathbf{1}^{\top}$, $0<\rho<1$.

  1. (i)

The sample covariance satisfies
$\sqrt{n}T_{n}\xrightarrow{d}N(0,(1-\rho)^{2})$.

  2. (ii)

$\sqrt{n}\rho_{n}\xrightarrow{d}N(0,1)$.

The above result exhibits an interesting phenomenon: the asymptotic variance
of $\sqrt{n}T_{n}$ is decreasing in $\rho$ but the asymptotic variance of
$\sqrt{n}\rho_{n}$ is independent of the choice of $\rho$.

###  4.2. Ordinary Least squares111The authors are grateful to Youjin Lee,
Brian Gilbert, and Abhirup Datta for their work on simulations and informal
results that preceded the formal results below.

Let $\mathbf{X}$ and $\mathbf{Y}$ be multivariate normal with simultaneously
diagonalizable variance matrices $\Sigma_{Y}=V\Lambda_{Y}V^{T}$ and
$\Sigma_{X}=V\Lambda_{X}V^{T}$. Then, letting $\boldsymbol{\varepsilon}$ be
the error after subtracting the least squares projection of $\mathbf{Y}$ onto
$\mathbf{X}$, $\boldsymbol{\varepsilon}\perp\mathbf{X}$ is also multivariate
normal and with simultaneously diagonalizable with variance matrix
$\Sigma_{\varepsilon}=V\Lambda_{\varepsilon}V^{T}$ and

$\displaystyle\mathbf{Y}=\mathbf{X}\beta+\boldsymbol{\varepsilon}.$

In order to derive interpretable results, we assume $\lambda_{X,i}=f(i/n)$ and
$\lambda_{\varepsilon,i}=g(i/n)$ for integrable functions
$f,g:[0,1]\rightarrow\mathbb{R}_{+}$.

We will assess when inference based on standard i.i.d. inference is
nevertheless valid, via an asymptotic comparison of the naive to the true
variance of $\widehat{\beta}_{\rm ols}$, where the estimated (naive) standard
error is given by the standard i.i.d. formula
$\widehat{\text{Var}_{\text{naive}}}=\frac{\boldsymbol{\varepsilon}^{T}(\mathbb{I}-P_{\mathbf{X}})\boldsymbol{\varepsilon}}{n\|\mathbf{X}\|^{2}}=\frac{\boldsymbol{\varepsilon}^{T}(\|\mathbf{X}\|^{2}-\mathbf{X}\mathbf{X}^{\top})\boldsymbol{\varepsilon}}{n\|\mathbf{X}\|^{4}}$,

######  Theorem 7.

Assume $\mathbf{X}$ and $\mathbf{Y}$ are simultaneously diagonalizable with
eigenvalues as described above. Then naive i.i.d. inference for $\beta$ is
valid, achieving asymptotic Type I error less than or equal to a nominal
level, if and only if the following condition holds:

$\displaystyle\int f(x)g(x)dx\leq\int f(x)dx\int g(x)dx.$ (4.1)

Notably, this condition is met whenever $\mathbf{X}$ is i.i.d., even if
$\mathbf{Y}$ and $\varepsilon$ are not. In this case ratio of the true to the
naive variance is asymptotically equal to $1$, meaning that inference using
the naive variance estimator will asymptotically exactly achieve the nominal
type I error rate.

If $f,g$ are both increasing or decreasing, (4.1) is violated (by the FKG
inequality), yielding anticonservative inference. However, if one is
increasing and the other is decreasing, then the condition is met and niave
i.i.d. inference is valid.

Note that even when (4.1) holds and naive inference is valid or exact, this
inference is not typically efficient (unless $\Sigma_{\varepsilon}$ is the
identity matrix). The efficient estimator in these cases is generalized least
squares (GLS). Under some conditions the estimated regression coefficient from
generalized least squares (GLS) will coincide with that from ordinary least
squares (OLS) [51]; in future work it will be interesting to explore the
relationship between the conditions under which OLS inference is valid and the
conditions under which OLS estimates are efficient.

##  5\. Numerical Experiments

Figures 2 and 3 illustrate the claims in Theorem 3. We fixed $G_{n}$ to be a
$1000$ by $1000$ lattice and used the Potts package in R to simulate the Ising
model on $G$, varying $\beta$ between $0$ and $\beta_{c}=1.76$.

Figure 4 illustrates the Rademacher distribution of Theorem 4 (ii).

Figure 5 shows confidence interval of $\beta$ under $X,Y$ multivariate normal
with non-diagonal covariance matrix. The columns correspond to $\Sigma_{X}$
having eigenvalues $(i/n)^{2}$, $i=1,2,\ldots,n$ and $\Sigma_{Y}$ have
eigenvalues (a) $e^{i/n^{0.85}}$, (b) $e^{i/n}$, (c) $e^{-i/n}$, (d)
$e^{-i/n^{0.85}}$, (e) $\Sigma_{X}=I$. In (a) and (b) we have anticonservative
inference, in (c) and (d) we have conservative inference, and in (e) inference
is exact (though not efficient). This illustrates the theoretical findings in
Theorem 7.

Figure 2. Empirical distribution of $\sqrt{n}\rho_{n}$ over 1000 simulation
replicates with $n=10,000$. The inverse temperature parameter $\beta$ ranges
from $0$ (i.i.d. data) to the critical value. The standard deviation of the
distribution increases monotonically with $\beta$, as does the type I error
rate $\alpha$ for a naive null hypothesis test of $\rho=0$. Figure 3.
Empirical distribution of $\sqrt{n}T_{n}$ over 1000 simulation replicates with
$n=10,000$. The inverse temperature parameter $\beta$ ranges from $0$ (i.i.d.
data) to the critical value. The standard deviation of the distribution
increases monotonically with $\beta$, as does the type I error rate $\alpha$
for a naive null hypothesis test of $T=0$. Figure 4. Empirical distribution of
$\sqrt{n}\rho_{n}$ over 1000 simulation replicates with $n=200$ with
$\mathbf{X},\mathbf{Y}\sim N(0,\Sigma)$ with $\widetilde{\lambda_{1}}\gg
n\widetilde{\lambda_{2}.}$

(a) Anticonservative inference

(b) Anticonservative inference

(c) Conservative inference

(d) Conservative inference

(e) Exact inference

Figure 5. Each column corresponds to $95\%$ naive OLS confidence intervals for
$\beta$ in the OLS regression of $Y$ onto $X$. In (a) and (b) the eigenvalues
of $\Sigma_{X}$ and $\Sigma_{Y}$ are both increasing; in (c) and (d) the
eigenvalues for $\Sigma_{X}$ are decreasing and $\Sigma_{Y}$ increasing. In
(e) $\Sigma_{X}$ is the identity matrix. Each column represents $500$
simulated replications of sample size $n=200$.

##  6\. Proofs

###  6.1. Proof of Theorem 1

For notational simplicity, we will assume throughout the proof that
$\beta_{1}=\beta_{2}=\beta$. The proof for general $\beta_{1},\beta_{2}$ goes
through with minor modification.

Recall that $T_{n}$ is the sample covariance of $\mathbf{X}$ and $\mathbf{Y}$
and $\mathbf{Q}_{n}$ is the adjaceny matrix on $\Lambda_{n,d}$. To obtain the
asymptotic distribution of $\sqrt{n}T_{n}$, we want to invoke [43, Theorem 6]
which provides a sufficient condition for asymptotic normality.

Following the notation of [43, Theorem 6], we set $A=\\{i\\}$ with $i\in[n]$
and $X^{(n)}_{A}=\frac{1}{\sqrt{n}}(X_{i}-\overline{\mathbf{X}}_{n})$, and
$\sigma_{A}=Y_{i}$. Then, we have
$\sqrt{n}T_{n}=\sum_{i=1}^{n}X^{n}_{\\{i\\}}\sigma_{\\{i\\}}$. The conditions
of the theorem are easily satisfied since,

  1. (1)

$\sum_{i=1}^{n}(X^{n}_{\\{i\\}})^{2}=\sum_{i=1}^{n}\frac{1}{n}(X_{i}-\overline{\mathbf{X}}_{n})^{2}\leq
1$, and

  2. (2)

$\max_{i\in n}|X^{n}_{\\{i\\}}|\leq\frac{2}{\sqrt{n}}\rightarrow 0$.

Therefore, taking $h=0$, we have, by [43, Theorem 6],
$\sqrt{n}T_{n}\xrightarrow{d}N(0,v)$ once we can show

$\lim\limits_{n\rightarrow\infty}v_{n}:=\frac{1}{n}\mathbb{E}\left(\sum_{i=1}^{n}(X_{i}-\overline{\mathbf{X}}_{n})Y_{i}\right)^{2}\rightarrow
v.$ (6.1)

To show (6.1), we have,

$\displaystyle v_{n}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}(X_{i}-\overline{\mathbf{X}}_{n})^{2}+\frac{1}{n}\sum_{i\neq
j}\mathbb{E}(X_{i}-\overline{\mathbf{X}}_{n})(X_{j}-\overline{\mathbf{X}}_{n})\mathbb{E}(Y_{i}Y_{j})$
$\displaystyle=1-\mathbb{E}(\overline{\mathbf{X}}_{n}^{2})+\frac{1}{n}\sum_{i\neq
j}\mathbb{E}(X_{i}-\overline{\mathbf{X}}_{n})(X_{j}-\overline{\mathbf{X}}_{n})\mathbb{E}(Y_{i}Y_{j}).$
(6.2)

Since $\mathbf{X},\mathbf{Y}\sim\mathbb{P}_{\beta,\mathbf{Q}}$, we have
$\mathbb{E}(Y_{i}Y_{j})=\mathbb{E}(X_{i}X_{j})$. Hence, upon expanding, the
second summand of (6.1) becomes,

$\frac{1}{n}\Big{(}\sum_{i\neq
j}\text{Cov}^{2}(X_{i},X_{j})+\mathbb{E}(\overline{\mathbf{X}}_{n}^{2})\sum_{i\neq
j}\text{Cov}(X_{i},X_{j})-2\sum_{i\neq
j}\text{Cov}(X_{i},X_{j})\mathbb{E}(X_{i}\overline{\mathbf{X}}_{n})\Big{)}.$
(6.3)

Now we want to show the third summand of the above display is $o(1)$. To see
this, note that,

$\displaystyle\frac{1}{n}\sum_{i\neq
j}\text{Cov}(X_{i},X_{j})\mathbb{E}(X_{i}\overline{\mathbf{X}}_{n})=\frac{1}{n^{2}}\sum_{i\neq
j}\text{Cov}(X_{i},X_{j})\times\sum_{l=1}^{n}\mathbb{E}(X_{i}X_{l})$
$\displaystyle=\frac{1}{n^{2}}\Big{(}\sum_{i\neq
j}\text{Cov}(X_{i},X_{j})\Big{(}\sum_{l\neq
i,j}\mathbb{E}(X_{i}X_{l})+1+\mathbb{E}(X_{i},X_{j})\Big{)}$
$\displaystyle=\frac{1}{n^{2}}\Big{(}\sum_{i\neq j\neq
l}\text{Cov}(X_{i},X_{j})\text{Cov}(X_{i},X_{l})+\sum_{i\neq
j}\text{Cov}(X_{i},X_{j})+\sum_{i\neq j}\text{Cov}^{2}(X_{i},X_{j})\Big{)}$
(6.4)

To bound the second and third summands of (6.4), we begin by stating the
following fact: there exists a constant $c=c(\beta,d)$ such that for all
$j\in\mathbb{Z}^{d}$,

$0\leq\mathrm{Cov}(X_{i},X_{j})\leq e^{-c\|i-j\|}$ (6.5)

if $\beta<\beta_{c}$ and $\|\cdot\|$ is the $L^{1}$ norm in ${Z}^{d}$. The
first inequality is known as the GKS inequality (cf [27]) and the second
inequality is proved in [2, 6, 22] for $\beta<\beta_{c}$. Using (6.5), one can
show that [6, Theorem 5]

$\sum_{i\neq j}\text{Cov}(X_{i},X_{j})=O(n).$ (6.6)

Now, using (6.6)and $\text{Cov}(X_{i},X_{j})\geq 0$ by (6.5), we further have
$\sum_{i\neq j}\text{Cov}^{2}(X_{i},X_{j})=O(n)$, yielding the second and
third summand of (6.4) to be $o(1)$. Finally, the first summand of (6.4) can
be simplified as

$\displaystyle\frac{1}{n^{2}}\sum_{i\neq j\neq
l}\text{Cov}(X_{i},X_{j})\text{Cov}(X_{i},X_{l})$
$\displaystyle=O\left(\frac{1}{n^{2}}\sum_{i\neq
j}\text{Cov}(X_{i},X_{j})\right)$
$\displaystyle=O\left(\frac{1}{n}\right)=o(1).$

This yields the third summand of (6.3) is $o(1)$. Now, we focus on the second
summand of (6.3). Since $\overline{\mathbf{X}}\xrightarrow{\mathbb{P}}0$ as
$\beta<\beta_{c}$, by the dominated convergence theorem, we have
$\mathbb{E}(\overline{\mathbf{X}}_{n}^{2})\rightarrow 0$. Again, invoking
(6.6), we have $\sum_{i\neq j}\text{Cov}(X_{i},X_{j})=O(n)$, yielding the
second summand to be $o(1)$. Plugging this into (6.1), we have

$v_{n}=1+\frac{1}{n}\sum_{i\neq j}\text{Cov}^{2}(X_{i},X_{j})+o(1)$ (6.7)

Now we claim the following:
$\lim\limits_{n\rightarrow\infty}\frac{1}{n}\sum_{i\neq
j}\text{Cov}^{2}(X_{i},X_{j})$ exists and is strictly positive. This will show
asymptotic normality and the fact that the asymptotic variance is strictly
larger than $1$, concluding our proof.

To see the first part of the claim, note that
$\lim\limits_{n\rightarrow\infty}\frac{1}{n}\sum_{i\neq
j}\text{Cov}(X_{i},X_{j})$ exists using [6, Theorem 5]. Since
$0\leq\mathrm{Cov}^{2}(X_{0},X_{j})\leq e^{-2c\|j\|}$, we have the same proof
for existence of limit.

To prove the second part of the claim, we plan to show, there exists a
constant $C:=C(\beta,d)>0$ such that

$\inf_{i\sim j}\text{Cov}(X_{i},X_{j})\geq C,$ (6.8)

where $i\sim j$ means $i$ and $j$-th vertices are neighbors in
$\Lambda_{n,d}$. This, in turn will show that $\frac{1}{n}\sum_{i\neq
j}\text{Cov}^{2}(X_{i},X_{j})\geq\frac{C}{n}\sum_{i\neq j}\mathbf{1}_{i\sim
j}=\Theta(1)$, since the number of neighbors on $\Lambda_{n,d}$ is
$\Theta(1)$, yielding the asymptotic distribution of $\sqrt{n}T_{n}$.

We now turn our attention to proving (6.8). To this end, we invoke the
Edwards-Sokal coupling between the measures $\mathbb{P}_{\beta,\mathbf{Q}}$
and random-cluster model. The random cluster model, on a finite graph
$G=(V,E)$ is a probability measure on $\\{0,1\\}^{E}$ defined by:

$\phi_{G,p,q}(\omega)=p^{|\\{e:\omega_{e}=1\\}|}(1-p)^{|\\{e:\omega_{e}=0\\}|}q^{k(\omega)}\qquad\omega\in\\{0,1\\}^{E}.$
(6.9)

where $k(\omega)$ denotes the number of connected components of $\omega$
(where we think of $\omega$ as the graph $(V,\\{e:\omega_{e}=1\\})$). The
Edwards-Sokal coupling (cf. [28, Theorem 4.91]) states that

$\text{Cov}(X_{i},X_{j})=\phi_{G,p,2}(i\leftrightarrow j),$ (6.10)

where the notation $i\leftrightarrow j$ implies that the vertices $i$, $j$ are
connected in the random-cluster model via open edges and $p:=1-e^{2\beta}$.
Now, we invoke [28, (3.23)] with $p_{1}=p$, $p_{2}=p/2$, $q_{1}=2$ and
$q_{2}=1$. Hence

$\phi_{G,p,2}(i\leftrightarrow j)\geq\phi_{G,p/2,1}(i\leftrightarrow j)=p/2,$
(6.11)

if $i\sim j$ in $G$ since the random-cluster model $\phi_{G,p,q}$ with $q=1$
corresponds to Bernoulli percolation with each edge being open independently
with probability $p$. Hence, if $i\sim j$, we have, by (6.10),
$\text{Cov}(X_{i},X_{j})\geq\frac{1}{2}(1-e^{2\beta})$, proving (6.8).

Moreover, since $\overline{\mathbf{X}}_{n}\xrightarrow{\mathbb{P}}0$ as
$\beta<\beta_{c}$, we have
$\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\overline{\mathbf{X}}_{n})^{2}=1-\overline{\mathbf{X}}^{2}_{n}\xrightarrow{\mathbb{P}}1$.
Therefore, by Slutsky’s theorem, we have
$\sqrt{n}\rho_{n}\xrightarrow{d}N(0,v^{2})$.
Finally, from (6.7), we know for $\beta(p_{1})<\beta(p_{2})<\beta_{c}$, where
$\beta{(p_{j})}$ solves $p_{j}=\frac{1}{2}(1-e^{2\beta(p_{j})})$, $j=1,2$,

$\displaystyle\lim\limits_{n\rightarrow\infty}v_{n}(\beta(p_{1}))$
$\displaystyle=1+\lim\limits_{n\rightarrow\infty}\frac{1}{n}\sum_{i\neq
j}\text{Cov}^{2}_{\beta(p_{1}),\mathbf{Q}}(X_{i},X_{j})$ $\displaystyle\leq
1+\lim\limits_{n\rightarrow\infty}\frac{1}{n}\sum_{i\neq
j}\text{Cov}^{2}_{\beta(p_{2}),\mathbf{Q}}(X_{i},X_{j})=\lim\limits_{n\rightarrow\infty}v_{n}(\beta(p_{2})),$

where the inequality is due to [6, Lemma 15]. This shows that the asymptotic
variance is a nondecreasing function of $\beta$, completing the proof.

###  6.2. Proof of Theorem 2

We want to show, for any $t\in\mathbb{R}$,

$\displaystyle\lim\limits_{n\rightarrow\infty}\mathbb{E}\Big{(}\exp(t\sqrt{n}T_{n})\Big{)}=\exp\left(\frac{t^{2}(1-m^{2}_{1})(1-m^{2}_{2})}{2}\right),$
(6.12)

where, recall that,
$T_{n}=\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\overline{\mathbf{X}}_{n})Y_{i}$ and
$m_{1},m_{2}$ are defined in the statement of Theorem 2. By Lemma 10, there
exists $Z_{2}$ such that $Z_{2}\sim
N\Big{(}\overline{\mathbf{Y}}_{n},\frac{1}{n\beta_{2}}\Big{)}$ and
$Y_{1},\ldots,Y_{n}$ are independent given $Z_{2}$, i.e.,
$\mathbb{P}(\mathbf{Y}=\mathbf{y}|Z_{2})\propto\prod_{i=1}^{n}e^{\beta_{2}Z_{2}y_{i}}$
where $\mathbf{y}=(y_{1},\ldots y_{n})$. Therefore, for any $i=1,2,\ldots,n$,
we have $\mathbb{E}(Y_{i}|Z_{2})=\tanh(\beta_{2}Z_{2})$ and
$\text{Var}(Y_{2}|Z_{2})=\operatorname{sech}^{2}(\beta_{2}Z_{2})$.

Moreover, by Lemma 10, $Z_{2}\xrightarrow{d}0$ if $\beta_{2}\leq 1$ and
$Z_{2}\xrightarrow{d}\frac{1}{2}(\delta_{m_{2}}+\delta_{-{m_{2}}})$ otherwise.
Invoking continuous mapping theorem, we have

$\operatorname{sech}^{2}(\beta_{2}Z_{2})\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}1-m^{2}_{2},$
(6.13)

where we used $\operatorname{sech}^{2}(\beta_{2}m_{2})=1-m^{2}_{2}$. Also,
since $\operatorname{sech}^{2}(x)\leq 1$, by the dominated convergence
theorem,
$\mathbb{E}\Big{(}\operatorname{sech}^{2}(\beta_{2}Z_{2})\Big{)}\rightarrow
1-m^{2}_{2}$.

Using Lemma 11, we obtain $\overline{\mathbf{X}}_{n}\xrightarrow{d}0$ if
$\beta_{1}\leq 1$ and
$\overline{\mathbf{X}}_{n}\xrightarrow{d}\frac{1}{2}(\delta_{m_{1}}+\delta_{-{m_{1}}})$
otherwise (using symmetry of the model around $0$). Therefore, by continuous
mapping theorem,

$\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\overline{\mathbf{X}}_{n})^{2}=1-\overline{\mathbf{X}}^{2}_{n}\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}1-m^{2}_{1}$
(6.14)

Since $\overline{\mathbf{X}}^{2}_{n}\leq 1$, by the dominated convergence
theorem,
$\mathbb{E}\Big{(}\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\overline{\mathbf{X}}_{n})^{2}\Big{)}\rightarrow
1-m^{2}_{1}$.

Now, conditioned on $\mathbf{X},Z_{2}$,

$\displaystyle\mathbb{E}\Big{(}\exp(tT_{n})|\mathbf{X},Z_{2}\Big{)}$
$\displaystyle=\exp\Big{(}\underbrace{\sum_{i=1}^{n}\Big{(}\log\cosh(\frac{t(X_{i}-\overline{\mathbf{X}}_{n})}{n^{1/2}}+\beta_{2}Z_{2})-\log\cosh(\beta_{2}Z_{2})}_{=:W_{n}}\Big{)}\Big{)}.$

By a Taylor series expansion,

$\displaystyle W_{n}-\frac{t^{2}}{2}(1-m^{2}_{1})(1-m^{2}_{2})$
$\displaystyle=\frac{t^{2}}{2}\left(\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\overline{\mathbf{X}}_{n})^{2}\operatorname{sech}^{2}(\beta_{2}Z_{2})-(1-m^{2}_{1})(1-m^{2}_{2})\right)+O\left(\frac{1}{\sqrt{n}}\right).$

Since
$\sum_{i=1}^{n}(X_{i}-\overline{\mathbf{X}}_{n})^{2}\xrightarrow{\mathbb{P}}1-m^{2}_{1}$
and $\operatorname{sech}^{2}(\beta_{2}Z_{2})\xrightarrow{d}1-m^{2}_{2}$ by
(6.14) and (6.13), we obtain

$\displaystyle|W_{n}-\frac{t^{2}}{2}(1-m^{2}_{1})(1-m^{2}_{2})|=o_{\mathbb{P}}(1).$
(6.15)

To show (6.12), using uniform integrability, it is enough to show

$\limsup\limits_{n\rightarrow\infty}\mathbb{E}(e^{2W_{n}-t^{2}(1-m^{2}_{1})(1-m^{2}_{2})})<\infty,$

which is immediate from the Taylor expansion above, since

$\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\overline{\mathbf{X}}_{n})^{2}\operatorname{sech}^{2}(\beta_{2}Z_{2})-(1-m^{2}_{1})(1-m^{2}_{2})\right|\leq
4.$

This completes the proof of the Theorem.

Proof of Corollary 3 When $\beta_{1}=\beta_{2}=0$, $X_{i}$, $Y_{i}$s are
independent and the proof is immediate. For $\beta_{i}>0$ case, the conclusion
follows combining Theorem 2 and (6.14) and using Slutsky’s theorem.

###  6.3. Proof of Theorem 4

Recall that $T_{n}$ is the sample covariance of $\mathbf{X}$ and $\mathbf{Y}$.
We need the following auxilliary lemma.

######  Lemma 8.

Suppose $\mathbf{Z}=(Z_{1},\ldots,Z_{n})$ are i.i.d. random variables with
mean $0$ and variance $\tau^{2}_{n}$ with converging to
$\tau^{2}\in(0,\infty)$. Assume the matrix $\mathbf{Q}_{n}$ satisfies
Assumptions 1-4. Set
$B_{n}=\mathbf{Q}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{\top}$. Suppose
$\mathbf{X}\sim\mathbb{P}_{\beta_{1},\mathbf{Q}_{n}}$ and $\mathbf{X}$ and
$\mathbf{Z}$ are independent.

Define
$\sqrt{n}T^{Z}_{n}=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(X_{i}-\overline{\mathbf{X}}_{n})Z_{i}$.
Then,

$(\sqrt{n}T^{Z}_{n},\mathbf{Z}^{\top}B_{n}\mathbf{Z})\xrightarrow{d}((1-m^{2}_{1})\tau_{0}W_{0},S_{0}),$
(6.16)

where $W_{0}\sim N(0,1)$, $S_{0}$ is some random variable, and $W_{0}$ and
$S_{0}$ independent, and $m_{1}$ is the unique positive solution of the
equation $m=\tanh(\beta_{1}m)$ if $\beta_{1}>1$ and $0$ otherwise.

We defer the proof of the above Lemma to Section Appendix. We now invoke the
Lemma 8 to prove the following convergence result for Curie-Weiss model.

######  Lemma 9.

Suppose $\mathbf{X}\sim\mathbb{P}_{\beta_{1},\mathbf{Q}_{n}}$,
$\mathbf{Y}\sim\mathbb{P}_{\beta_{2},\mathbf{Q}^{\text{CW}}}$,
$\mathbf{X},\mathbf{Y}$ are independent and $\mathbf{Q}_{n}$ satisfies
Assumptions 1-4. Then

$(\sqrt{n}T_{n},\mathbf{Y}^{\top}B_{n}\mathbf{Y})\xrightarrow{d}(W_{1},S_{0}),$

where $W_{1}\sim N(0,(1-m^{2}_{1})(1-m^{2}_{2}))$, $S_{0}$ defined in Lemma 8
and $W_{1}$ and $S_{0}$ are independent.

The proof of the Lemma is deferred to Section Appendix.

Now, we are in a position to prove the Theorem 4. To this end note the
following observations:

  1. (1)

The measures $P_{\beta_{1},\mathbf{Q}_{n}}\times P_{\beta_{2},\mathbf{Q}_{n}}$
and $P_{\beta_{1},\mathbf{Q}_{n}}\times P_{\beta_{2},\mathbf{Q}^{\text{CW}}}$
are mutually contiguous by [58, Lemma 4.4 part (b)].

  2. (2)

The log-likelihood ratio

$\log\frac{d(P_{\beta_{1},\mathbf{Q}_{n}}\times
P_{\beta_{2},\mathbf{Q}_{n}})}{d(P_{\beta_{1},\mathbf{Q}_{n}}\times
P_{\beta_{2},\mathbf{Q}^{\text{CW}}})}(\mathbf{X},\mathbf{Y})=\frac{\beta_{2}}{2}\mathbf{Y}^{\top}B_{n}\mathbf{Y}-Z_{n}(\beta_{2},\mathbf{Q}_{n})+Z_{n}(\beta_{2},\mathbf{Q}^{\text{CW}})$

is a function of $\mathbf{Y}^{\top}B_{n}\mathbf{Y}$.

  3. (3)

When
$(\mathbf{X},\mathbf{Y})\sim\mathbb{P}_{\beta_{1},\mathbf{Q}^{\text{CW}}}\times\mathbb{P}_{\beta_{2},\mathbf{Q}^{\text{CW}}}$,
$\sqrt{n}T_{n}$ has a non-degenerate asymptotic distribution independent of
$\mathbf{Y}^{\top}B_{n}\mathbf{Y}$.

Therefore, by Le Cam’s third Lemma [54], the asymptotic distribution of
$\sqrt{n}T_{n}$ is same under the measures $P_{\beta_{1},\mathbf{Q}_{n}}\times
P_{\beta_{2},\mathbf{Q}_{n}}$ and $P_{\beta_{1},\mathbf{Q}_{n}}\times
P_{\beta_{2},\mathbf{Q}^{\text{CW}}}$. This yields the asymptotic distribution
by using Lemma 9.

For the sample correlation, recall that
$\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\overline{\mathbf{X}}^{2}_{n})=1-\overline{\mathbf{X}}^{2}_{n}\xrightarrow{\mathbb{P}}1-m^{2}_{1}$
and
$\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-\overline{\mathbf{Y}}^{2}_{n})\xrightarrow{\mathbb{P}}1-m^{2}_{2}$.
This, along with Slutsky’s Theorem, completes the proof.

###  6.4. Proof of Theorem 5

(i) Here, we assume $\mathbf{X},\mathbf{Y}\sim N(0,\Sigma)$ are independent
and $T_{n}$ is their sample covariance. So,
$nT_{n}=\mathbf{X}^{\top}J\mathbf{Y}$ where
$J=I-\frac{1}{n}\mathbf{1}\mathbf{1}^{\top}$. Conditioned on $\mathbf{X}$, we
know $nT_{n}\sim N(0,\mathbf{X}^{\top}J\Sigma J\mathbf{X})$. Recall that
$\widetilde{\Sigma}:=\Sigma^{1/2}J\Sigma^{1/2}$ and let
$\widetilde{\Sigma}=P\widetilde{\Lambda}P^{\top}$ be the singular value
decomposition of $\widetilde{\Sigma}$. Therefore

$\mathbf{X}^{\top}J\Sigma
J\mathbf{X}={Z}^{\top}\widetilde{\Sigma}^{2}{Z}={Z}^{\top}P\widetilde{\Lambda}^{2}P^{\top}{Z}\stackrel{{\scriptstyle
D}}{{=}}{Z}^{\top}\widetilde{\Lambda}^{2}{Z},$ (6.17)

where ${Z}\sim N(0,I_{n})$. Now, we invoke Lemma 12 since,
$\widetilde{\lambda}^{2}_{1}\ll\sum_{i=1}^{n}\widetilde{\lambda}^{2}_{i}$, to
obtain
$\frac{{Z}^{\top}\widetilde{\Lambda}^{2}{Z}}{\mathbb{E}({Z}^{\top}\widetilde{\Lambda}^{2}{Z})}\xrightarrow{\mathbb{P}}1$.
Now,
$\mathbb{E}({Z}^{\top}\widetilde{\Lambda}^{2}{Z})=\text{Trace}(\widetilde{\Lambda}^{2})=\sum_{i=1}^{n}\widetilde{\lambda}^{2}_{i}$.
Hence, by Slutsky’s Theorem, we prove ((i)).

To show asymptotic distribution of sample correlation $\rho_{n}$, note that,
$\widetilde{\lambda}_{1}\ll\sum_{i=1}^{n}\widetilde{\lambda}_{i}$, hence by
Lemma 12, we have
$\frac{\mathbf{X}^{\top}J\mathbf{X}}{\mathbb{E}(\mathbf{X}^{\top}J\mathbf{X})}\xrightarrow{\mathbb{P}}1$.
Also,
$\mathbb{E}(\mathbf{X}^{\top}J\mathbf{X})=\sum_{i=1}^{n}\widetilde{\lambda}_{i}$.
Hence,

$\displaystyle\sqrt{n}\rho_{n}a_{n}$
$\displaystyle=\sqrt{n}\frac{nT_{n}a_{n}}{\sqrt{(\mathbf{X}^{\top}J\mathbf{X})({Y}^{\top}J{Y})}}=(1+o_{\mathbb{P}}(1))\sqrt{n}\frac{nT_{n}a_{n}}{\sum_{i=1}^{n}\widetilde{\lambda}_{i}}$
$\displaystyle=(1+o_{\mathbb{P}}(1)\frac{nT_{n}}{\sqrt{\sum_{i=1}^{n}\widetilde{\lambda}^{2}_{i}}}\xrightarrow{d}N(0,1).$

Also, by Cauchy-Schwarz inequality and the fact $\widetilde{\lambda}_{i}\geq
0$, we have $a_{n}\leq 1$. Therefore
$\liminf_{n\rightarrow\infty}\text{Var}(\sqrt{n}\rho_{n})\geq\lim\limits_{n\rightarrow\infty}\frac{1}{a_{n}}\geq
1$.

(ii) Now, we assume $\widetilde{\lambda}^{2}_{1}\gg
n\widetilde{\lambda}^{2}_{2}$. As the previous part, conditioned on
$\mathbf{X}$, $nT_{n}\sim N(0,\mathbf{X}^{\top}J\Sigma J\mathbf{X})$. Hence,
$\frac{nT_{n}}{\sqrt{\mathbf{X}^{\top}J\Sigma J\mathbf{X}}}\sim N(0,1)$ and it
is independent of $\mathbf{X}$. So, it is enough to show
$\zeta_{n}:=\frac{\mathbf{X}^{\top}J\Sigma
J\mathbf{X}}{\widetilde{\lambda}^{2}_{1}}\xrightarrow{d}\chi^{2}_{1}$. By
(6.17), we have

$\displaystyle\zeta_{n}\stackrel{{\scriptstyle
d}}{{=}}\frac{\sum_{i=1}^{n}\widetilde{\lambda}^{2}_{i}Z^{2}_{i}}{\widetilde{\lambda}^{2}_{1}}=Z^{2}_{1}+\sum_{i=2}^{n}\frac{\widetilde{\lambda}^{2}_{i}}{\widetilde{\lambda}^{2}_{1}}Z^{2}_{i}$
$\displaystyle\leq
Z^{2}_{1}+\frac{n\widetilde{\lambda}^{2}_{2}}{\widetilde{\lambda}^{2}_{1}}\Big{(}\frac{1}{n}\sum_{i=2}^{n}Z^{2}_{i}\Big{)}$
$\displaystyle=Z^{2}_{1}+o_{\mathbb{P}}(1),$

since $\widetilde{\lambda}^{2}_{1}\gg n\widetilde{\lambda}^{2}_{2}$ and
$\sum_{i=2}^{n}Z^{2}_{i}=O_{\mathbb{P}}(n)$. As $Z^{2}_{1}\sim\chi^{2}_{1}$,
we get the desired result.

To compute the asymptotic distribution of $\rho_{n}$, note that, following a
similar calculation as above, if $Z_{1},Z_{2}\sim N(0,I)$ are independent,
then

$\displaystyle\rho_{n}=\frac{\sum_{i=1}\widetilde{\lambda}_{i}Z_{1i}Z_{2i}}{\sqrt{\Big{(}\sum_{i=1}\widetilde{\lambda}_{i}Z^{2}_{1i}\Big{)}\Big{(}\sum_{i=1}\widetilde{\lambda}_{i}Z^{2}_{2i}\Big{)}}}=\frac{Z_{11}Z_{21}+o_{\mathbb{P}}(1)}{\sqrt{Z^{2}_{11}Z^{2}_{21}+o_{\mathbb{P}}(1)}}\xrightarrow{d}\frac{1}{2}(\delta_{1}+\delta_{-1}),$

concluding the proof of the Theorem.

Proof of Lemma 6 Here we assume the Contagion model, i.e.,
$\mathbf{X},\mathbf{Y}\sim N(0,\Sigma)$ are independent and
$\Sigma=(1-\rho)I+\rho\mathbf{1}\mathbf{1}^{\top}$, $0<\rho<1$.

(i) We want to compute the asymptotic distribution of $\sqrt{n}T_{n}$. To this
end, note that $J$ has eigenvalues $1$ with multiplicity $(n-1)$ and $0$ with
multiplicity $1$ and $J\Sigma=\Sigma J=(1-\rho)J$. Since the eigenvalues of
$(\Sigma^{1/2})^{\top}J\Sigma^{1/2}$ are same as the eigenvalues of $J\Sigma$,
we have, using the notation of Theorem 5,
$\widetilde{\lambda}_{1}=\ldots=\widetilde{\lambda}_{n-1}=(1-\rho)$ and
$\widetilde{\lambda}_{n}=0$. Therefore,
$\widetilde{\lambda}^{2}_{1}\ll\sum_{i=1}^{n}\widetilde{\lambda}^{2}_{i}$ and
we invoke Theorem 5 (i) to claim
$\frac{nT_{n}}{\sqrt{(n-1)(1-\rho)^{2}}}\xrightarrow{d}N(0,1)$, i.e.,
$\sqrt{n}T_{n}\xrightarrow{d}N(0,(1-\rho)^{2})$. Moreover,

$\displaystyle
a_{n}=\frac{(1-\rho)(n-1)}{\sqrt{n(n-1)(1-\rho)^{2}}}\rightarrow 1,$

yielding $\sqrt{n}\rho_{n}\xrightarrow{d}N(0,1)$.
(ii) We want to compute the asymptotic distribution of
$S_{n}=\mathbf{X}^{\top}{Y}$. Borrowing notation from Theorem 5, we have
$\lambda_{1}=1+(n-1)\rho$ and $\lambda_{2}=\ldots=\lambda_{n}=(1-\rho)$. Since
$\lambda^{2}_{1}\gg n\lambda^{2}_{2}$, we invoke Theorem 5 to obtain
$\frac{nS_{n}}{1+(n-1)\rho}=(1+o(1))\frac{S_{n}}{\rho}\xrightarrow{d}AB$,
where $A\sim N(0,1)$, $B\sim\sqrt{\chi^{2}_{1}}$ and $A,B$ are indendent. This
completes the proof.

###  6.5. Proof of Proposition 7

We need to compare the naive variance with the variance of
$\widehat{\beta}_{\rm ols}$. To this end, note the following conclusions:

$\displaystyle\frac{1}{n}\|\varepsilon\|^{2}\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}\int
g(x)dx,\qquad\frac{1}{n}\|\mathbf{X}\|^{2}\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}\int
f(x)dx,\qquad\frac{1}{n}\varepsilon^{\top}\mathbf{X}\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}0.$
(6.18)

The first convergence in (6.18) is true since
$\mathbb{E}(\frac{1}{n}\|\varepsilon\|^{2})=\frac{1}{n}\sum_{i=1}^{n}g(i/n)\rightarrow\int
g(x)dx$, and $\text{Var}(\frac{1}{n}\|\varepsilon\|^{2})\rightarrow 0$. The
other two convergences follow similarly. Moreover

$\frac{1}{n}\mathbf{X}^{T}\Sigma_{\varepsilon}\mathbf{X}\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}\int
f(x)g(x)dx,$ (6.19)

since
$\mathbb{E}(\frac{1}{n}\mathbf{X}^{T}\Sigma_{\varepsilon}\mathbf{X})=\frac{1}{n}\sum_{i=1}^{n}f(i/n)g(i/n){\rightarrow}\int
f(x)g(x)dx$, and
$\text{Var}(\frac{1}{n}\mathbf{X}^{T}\Sigma_{\varepsilon}\mathbf{X})\rightarrow
0$.

Hence,

$\displaystyle n\widehat{\text{Var}_{\text{naive}}}$
$\displaystyle=\frac{\boldsymbol{\varepsilon}^{T}(\|\mathbf{X}\|^{2}-\mathbf{X}\mathbf{X}^{\top})\boldsymbol{\varepsilon}}{\|\mathbf{X}\|^{4}}=\frac{\|\varepsilon\|^{2}}{\|\mathbf{X}\|^{2}}\left(1-\frac{(\varepsilon^{\top}\mathbf{X})^{2}}{\|\varepsilon\|^{2}\|\mathbf{X}\|^{2}}\right)\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}\frac{\int
g(x)dx}{\int f(x)dx},$

where the first convergence follows from (6.18). Further, we have,

$\displaystyle n\text{Var}(\widehat{\beta}_{\rm
ols})=\frac{n\mathbf{X}^{T}\Sigma_{\varepsilon}\mathbf{X}}{\|\mathbf{X}\|^{4}}\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}\frac{\int
f(x)g(x)dx}{(\int f(x)dx)^{2}}$

using (6.18) and (6.19). Naive inference is valid iff
$\frac{\text{Var}(\widehat{\beta}_{\rm
ols})}{\widehat{\text{Var}_{\text{naive}}}}<1$, which by the above two
displays equivalent to

$\int f(x)g(x)dx\leq\Big{(}\int f(x)dx\Big{)}\Big{(}\int g(x)dx\Big{)},$

concluding our proof.

##  7\. Conclusion

This paper is a first exploration of the conditions under which nonsense
associations – or the opposite phenomenon of valid / conservative inference –
occur for dependent data with dependence governed by MRF distributions. Our
results serve as a proof of concept that researchers ought to be concerned
about nonsense associations in settings beyond time series, but more work is
needed to fully characterize the effects of unit-level dependence on measures
of association. In particular, we report several surprising findings, and the
extent to which these can occur in other settings should be explored in future
work: (1) settings where the asymptotic behavior of the sample correlation is
qualitatively different from the asymptotic behavior of the sample covariance,
and (2) settings where unit-level positive dependence results in variance
deflation and conservative inference. The latter is especially surprising
because it contradicts the notion that positive dependence results in smaller
effective sample sizes for statistical inference. Finally, our results on OLS
upend the way that dependent residuals is typically presented in introductory
regression texts.

## Appendix

Here, we state the technical Lemmas along with proofs which we have used
throughout the paper. The following lemma follows from [44, Lemma 3].

######  Lemma 10.

Suppose $\mathbf{X}\sim\mathbb{P}_{\beta,\mathbf{Q}^{\rm CW}}$ with
$\mathbf{Q}^{\rm CW}_{i,j}=\mathbf{1}(i\neq j)/n$. Define a random variable
$Z_{n}$ which given $\mathbf{X}$ has a distribution
$N(\overline{\mathbf{X}}_{n},\frac{1}{n\beta})$. Then under
$\mathbb{P}_{\beta,\mathbf{Q}^{\rm CW}}$, given $Z_{n}=z_{n}$, each $X_{i}$’s
are i.i.d. with

$\displaystyle\mathbb{P}_{\beta,\mathbf{Q}^{\rm
CW}}(X_{i}=x_{i}|Z_{n}=z_{n})=\frac{e^{\beta z_{n}x_{i}}}{e^{\beta
z_{n}}+e^{-\beta z_{n}}}.$

Moreover, $Z\xrightarrow{d}0$ if $\beta\leq 1$ and
$Z\xrightarrow{d}\frac{1}{2}(\delta_{m}+\delta_{-{m}})$ otherwise, where
$m=m(\beta)$ is the unique positive root of $m=\tanh(\beta m)$.

The following lemma follows from [25, Theorems 1-3].

######  Lemma 11.

Suppose $\mathbf{X}\sim\mathbb{P}_{\beta,\mathbf{Q}^{\rm CW}}$ with
$\mathbf{Q}^{\rm CW}_{i,j}=\mathbf{1}(i\neq j)/n$.

  1. (a)

If $\beta\in(0,1)$ then

$\sqrt{n}\overline{\mathbf{X}}_{n}\stackrel{{\scriptstyle
d}}{{\rightarrow}}N\left(0,\frac{1}{1-\beta}\right).$

  2. (b)

If $\beta>1$ then we have

$\sqrt{n}(\overline{\mathbf{X}}_{n}-m|\overline{\mathbf{X}}_{n}>0)\stackrel{{\scriptstyle
d}}{{\rightarrow}}N\left(0,\frac{1-m^{2}}{1-\beta(1-m^{2})}\right),$

where $m$ is the unique positive root of the equation $t=\tanh(\beta t)$.

  3. (c)

If $\beta=1$ then we have

$n^{1/4}\overline{\mathbf{X}}_{n}\stackrel{{\scriptstyle d}}{{\rightarrow}}Y,$

where $Y$ is a random variable with density proportional to $e^{-y^{4}/12}$.

Next, we prove a lemma regarding concentration inequalities of quadratic forms
of multivariate normal distributions.

######  Lemma 12.

Suppose $A$ is an $n\times n$ matrix with eigenvalues
$\lambda_{1}\geq\lambda_{2}\geq\ldots\geq 0$ and
$\frac{\lambda_{1}}{\sum_{i=1}^{n}\lambda_{i}}=o(1)$. If $Z\sim N(0,I_{n})$,
then
$\frac{Z^{\top}AZ}{\mathbb{E}(Z^{\top}AZ)}\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}1$.

###### Proof.

By Markov’s inequality, for any $t>0$,

$\displaystyle\mathbb{P}\Bigg{(}\Big{|}\frac{Z^{\top}AZ}{\mathbb{E}(Z^{\top}AZ)}-1\Big{|}\geq
t\Bigg{)}\leq\frac{\text{Var}(Z^{\top}AZ)}{t^{2}(\mathbb{E}(Z^{\top}AZ))^{2}}$
(7.1)

Note that, $\mathbb{E}(Z^{\top}AZ)=\text{Tr}(A)=\sum_{i=1}^{n}\lambda_{i}$ and
$\text{Var}(Z^{\top}AZ)=2\text{Tr}(A^{2})=2\sum_{i=1}^{n}\lambda^{2}_{i}$. So,
the RHS of (7.1) is

$\displaystyle\frac{2\sum_{i=1}^{n}\lambda^{2}_{i}}{t^{2}(\sum_{i=1}^{n}\lambda_{i})^{2}}\leq\frac{2\lambda_{1}}{t^{2}\sum_{i=1}^{n}\lambda_{i}}=o(1),$

by our assumption. This completes the proof. ∎

Next, we provide the proof of Lemmas omitted before.

Proof of Lemma 8 The proof is done in two steps.

  1. (1)

Suppose $G_{1},\ldots G_{n}$ are i.i.d. $N(0,\tau^{2}_{n})$ random variables.
Then the conclusion of the Lemma holds.

  2. (2)

For any $a,b\in\mathbb{N}$, we have

$\mathbb{E}(\sqrt{n}T^{Z}_{n})^{a}(\mathbf{Z}^{\top}B_{n}\mathbf{Z})^{b}-\mathbb{E}(\sqrt{n}T^{G}_{n})^{a}(\mathbf{G}^{\top}B_{n}\mathbf{G})^{b}\rightarrow
0,$

where
$\sqrt{n}T^{G}_{n}=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(X_{i}-\overline{\mathbf{X}}_{n})G_{i}$

To show $(1)$, we will throughout assume $\tau^{2}_{n}=1$ and the case of
general $\tau^{2}_{n}$ is immediate by change of scale. We begin by noting
that for $a,b\in\mathbb{R}$ with $b<1/2$, one has

$\mathbb{E}(\exp(aG_{1}+bG^{2}_{1}))=\frac{1}{\sqrt{1-2b}}\exp(\frac{a^{2}}{2(1-2b)}).$
(7.2)

Define $\upsilon_{n}=\frac{1}{\sqrt{n}}(X_{1}-\overline{\mathbf{X}}_{n},\ldots
X_{n}-\overline{\mathbf{X}}_{n})$ so that
$\sqrt{n}T^{G}_{n}=\upsilon^{T}_{n}\mathbf{G}$. Let
$\mathbf{Q}_{n}=P^{\top}\Lambda
P=\sum_{i=1}^{n}\lambda_{i,n}\mathbf{p}_{i}\mathbf{p}^{\top}_{i}$ be the
spectral decomposition of $\mathbf{Q}_{n}$, where the eigenvalues
$\lambda_{i,n}$, $i=1,2,\ldots,n$ be the eigenvalues of $\mathbf{Q}_{n}$
sorted in the decreasing order of absolute values. Since
$\sum_{j=1}^{n}Q_{n,ij}=1$ for all $i$, we have $\lambda_{1,n}=1$ and
$\mathbf{p}_{1}=\frac{1}{\sqrt{n}}\mathbf{1}$. Setting
$\widetilde{\mathbf{G}}=P\mathbf{G}$, where $\mathbf{G}=(G_{1},\ldots,G_{n})$,
we have,

$\Big{(}\upsilon^{T}_{n}\mathbf{G},\mathbf{G}^{\top}B_{n}\mathbf{G}\Big{)}\stackrel{{\scriptstyle
d}}{{=}}\Big{(}\upsilon^{T}_{n}P^{T}\widetilde{\mathbf{G}},\sum_{i=2}^{n}\lambda_{i,n}\widetilde{G}^{2}_{i}\Big{)}\stackrel{{\scriptstyle
d}}{{=}}\Big{(}\mathbf{a}^{\top}_{n}\mathbf{G},\sum_{i=2}^{n}\lambda_{i,n}G^{2}_{i}\Big{)},$

where we defined $\mathbf{a}_{n}=P\upsilon_{n}$ and used the fact
$P\mathbf{G}\stackrel{{\scriptstyle d}}{{=}}\mathbf{G}$. Now, invoking (7.2),
we have, for $s,t\in\mathbb{R}$ and $t<\frac{1}{2}$ and fixed $\mathbf{X}$,

$\displaystyle\log\mathbb{E}\Big{(}\exp(s\mathbf{a}^{\top}_{n}\mathbf{G}+t\sum_{i=2}^{n}\lambda_{i,n}G^{2}_{i})\Big{)}$
$\displaystyle=-\frac{1}{2}\sum_{i=2}^{n}\log(1-2t\lambda_{i,n})+\frac{s^{2}a^{2}_{1}}{2}$
(7.3) $\displaystyle+\sum_{i=2}^{n}\frac{s^{2}a^{2}_{i}}{2(1-2b_{i})}.$ (7.4)

Denote the first summand of RHS by $f_{n}(t)$. Following the proof of [58,
Lemma 4.1], we have $f_{n}(t)$ converges to the log of moment generating
function of a random variable $S_{0}$. Next we analyse the third summand of
the above display. To this end, we have

$\displaystyle\sum_{i=2}^{n}\frac{s^{2}a^{2}_{i}}{2(1-2t\lambda_{i,n})}=\sum_{i=2}^{n}\sum_{j=0}^{\infty}\frac{s^{2}a^{2}_{i}}{2}(2t\lambda_{i,n}))^{j}$
$\displaystyle=\sum_{i=2}^{n}\frac{s^{2}a^{2}_{i}}{2}+\underbrace{2ts^{2}\sum_{i=2}^{n}a^{2}_{i}\lambda_{i,n}+\sum_{j=2}^{\infty}\sum_{i=2}^{n}\frac{s^{2}a^{2}_{i}}{2}(2t\lambda_{i,n}))^{j}}_{R_{n}}$
(7.5)

To control the second summand above, $\mathbf{1}^{T}\upsilon_{n}=0$ yielding
$a_{1}=0$. Hence,

$\displaystyle\sum_{i=2}^{n}a^{2}_{i}\lambda_{i,n}=\sum_{i=1}^{n}a^{2}_{i}\lambda_{i,n}$
$\displaystyle=\upsilon^{\top}_{n}\mathbf{Q}_{n}\upsilon_{n}$
$\displaystyle=\frac{1}{n}\mathbf{X}^{\top}(I-\frac{1}{n}\mathbf{1}\mathbf{1}^{\top})\mathbf{Q}_{n}(I-\frac{1}{n}\mathbf{1}\mathbf{1}^{\top})\mathbf{X}$
$\displaystyle=\frac{1}{n}\mathbf{X}^{\top}B_{n}\mathbf{X}\xrightarrow{\mathbb{P}}0,$

since $\mathbf{Q}_{n}\mathbf{1}=\mathbf{1}$ and
$\mathbf{X}^{\top}B_{n}\mathbf{X}$ is tight by [58, Lemma 3.1] (note that,
[58, Lemma 3.1] states the result conditional on $\overline{\mathbf{X}}>0$,
but the same proof works for the distribution conditional on
$\overline{\mathbf{X}}\leq 0$). This yields that the second summand of
(Appendix) is $o_{\mathbb{P}}(1)$. Finally the third summand of (Appendix) by

$\displaystyle\sum_{j=2}^{\infty}\sum_{i=2}^{n}\frac{s^{2}a^{2}_{i}}{2}(2t\lambda_{i,n})^{j}$
$\displaystyle=\sum_{j=2}^{\infty}\frac{s^{2}(-2t)^{j}}{2}\sum_{i=2}^{n}a^{2}_{i}\lambda^{j}_{i,n}$
$\displaystyle\leq\sum_{j=2}^{\infty}\frac{s^{2}(2t)^{j}}{2}\sum_{i=2}^{n}a^{2}_{i}\lambda^{2}_{i,n}\leq\frac{s^{2}}{2(1-2t)}\sum_{i=2}^{n}a^{2}_{i}\lambda^{2}_{i,n},$

where the first inequality uses $|\lambda_{i,n}|\leq 1$. Note that,

$\sum_{i=2}^{n}a^{2}_{i}\lambda^{2}_{i,n}=\upsilon^{\top}Q^{2}_{n}\upsilon=\frac{1}{n}\mathbf{X}^{\top}B^{2}_{n}\mathbf{X}\xrightarrow{\mathbb{P}}0$

invoking [58, Lemma 3.1] as above. Therefore, we have, by (7.3),

$\log\mathbb{E}\Big{(}\exp(s\mathbf{a}^{\top}_{n}\mathbf{G}+t\sum_{i=2}^{n}\lambda_{i,n}G^{2}_{i})\Big{)}=f_{n}(t)+\frac{s^{2}}{2}\mathbf{a}^{\top}_{n}\mathbf{a}_{n}+R_{n},$

where $R_{n}=o_{\mathbb{P}}(1)$ and $|R_{n}|\leq
2|t|s^{2}+\frac{s^{2}}{2|1-2t|}$. Moreover,
$\mathbf{a}^{\top}_{n}\mathbf{a}_{n}=\upsilon^{T}_{n}\upsilon_{n}=\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\overline{\mathbf{X}}_{n})^{2}\rightarrow
1-m^{2}_{1}$ and $|\mathbf{a}^{\top}_{n}\mathbf{a}_{n}|\leq 4$. Therefore, by
the dominated convergence theorem, we have

$\displaystyle\lim\limits_{n\rightarrow\infty}\mathbb{E}\Big{(}\exp(s\mathbf{a}^{\top}_{n}\mathbf{G}+t\sum_{i=2}^{n}\lambda_{i,n}G^{2}_{i})\Big{)}=e^{\frac{s^{2}}{2}(1-m^{2}_{1})}\mathbb{E}e^{tS_{0}}.$

This shows,
$(\upsilon^{T}_{n}\mathbf{G},\mathbf{G}^{\top}B_{n}\mathbf{G})\xrightarrow{d}(S_{0},\sqrt{1-m^{2}_{1}}W_{0})$
where $W_{0}\sim N(0,1)$ and $S_{0},W_{0}$ are asymptotically independent
proving (1).

To show $(2)$, note that, $|B_{n}(i,j)|\leq\frac{4}{n}$, since
$\mathbf{Q}_{n}$ is the (scaled) adjacency matrix of $\mathcal{G}_{n,p}$ or
$\mathcal{G}^{d}_{n}$ and $|c_{i,n}|\leq\frac{\kappa}{\sqrt{n}}$ by our
assumption. Set $r=\frac{a}{2}+b$ and $S(l,2r)$ be the set of all positive
integer solutions to the equation $\sum_{i=1}^{l}\alpha_{i}=2r$. Then we have,

$\displaystyle|\mathbb{E}(\sqrt{n}T^{Z}_{n})^{a}(\mathbf{Z}^{\top}B_{n}\mathbf{Z})^{b}-\mathbb{E}(\sqrt{n}T^{G}_{n})^{a}(\mathbf{G}^{\top}B_{n}\mathbf{G})^{b}|$
$\displaystyle\leq
n^{-r}\kappa^{a}4^{b}\sum_{l=1}^{2r}n^{l}\sum_{\mathbb{\alpha}\in
S(l,2r)}|\mathbb{E}\prod\limits_{i=1}^{l}Z^{l}_{i}-\mathbb{E}\prod\limits_{i=1}^{l}G^{l}_{i}|,$

and the RHS in the above display can be bounded exactly as the proof of [58,
Lemma 4.1].

Proof of Lemma 9 We do the entire proof conditioning on $\mathbf{X}$. By Lemma
10, we have the existence of auxilliary variable $Z_{n,2}$ such that
$Y_{1},\ldots,Y_{n}$ are independent given $Z_{n,2}$ with

$\displaystyle\mathbb{E}(Y_{i}|Z_{n,2})=\tanh(\beta_{2}Z_{n,2})=:\mu,\quad\text{Var}(Y_{i}|Z_{n,2})=\operatorname{sech}^{2}(\beta_{2}Z_{n,2}),\text{
for all }i.$

Define $\boldsymbol{\mu}_{n}=\mu\mathbf{1}$ and recall from proof of Lemma 8,
we have
$\upsilon_{n}=\frac{1}{\sqrt{n}}(X_{1}-\overline{\mathbf{X}}_{n},\ldots
X_{n}-\overline{\mathbf{X}}_{n})$. Therefore, we have, using
$B_{n}\mathbf{1}=0$, given $Z_{n}$

$\displaystyle(\sqrt{n}T_{n},\mathbf{Y}^{\top}B_{n}\mathbf{Y})$
$\displaystyle=(\upsilon^{\top}_{n}\mathbf{Y},\mathbf{Y}^{\top}B_{n}\mathbf{Y})$
$\displaystyle=(\upsilon^{\top}_{n}(\mathbf{Y}-\boldsymbol{\mu}_{n}),(\mathbf{Y}-\boldsymbol{\mu}_{n})^{\top}B_{n}(\mathbf{Y}-\boldsymbol{\mu}_{n}))$
$\displaystyle\xrightarrow{d}(\sqrt{(1-m^{2}_{1})(1-m^{2}_{2})}W_{0},S_{0}),$

with $W_{0}\sim N(0,1)$, $S_{0}$ as in Lemma 8 and $W_{0},S_{0}$ being
independent by Lemma 8 since
$\operatorname{sech}^{2}(\beta_{2}Z_{n,2})\xrightarrow{\mathbb{P}}1-m^{2}_{2}$
by Lemma 10.

## References

  * [1] Amr Ahmed and Eric P Xing.  Recovering time-varying networks of dependencies in social and biological studies.  Proceedings of the National Academy of Sciences, 106(29):11878–11883, 2009.
  * [2] Michael Aizenman, David J Barsky, and Roberto Fernández.  The phase transition in a general class of ising-type models is sharp.  Journal of Statistical Physics, 47(3-4):343–374, 1987.
  * [3] Michael Banf and Seung Y Rhee.  Enhancing gene regulatory network inference through data integration with markov random fields.  Scientific reports, 7(1):41174, 2017.
  * [4] Anirban Basak and Sumit Mukherjee.  Universality of the mean-field for the potts model.  Probability Theory and Related Fields, pages 1–44, 2015.
  * [5] Sohom Bhattacharya, Rounak Dey, and Rajarshi Mukherjee.  Pc adjusted testing for low dimensional parameters.  arXiv preprint arXiv:2209.10774, 2022.
  * [6] Sohom Bhattacharya, Rajarshi Mukherjee, and Gourab Ray.  Sharp signal detection under ferromagnetic ising models.  arXiv preprint arXiv:2110.02949, 2021.
  * [7] Christian Borgs, Jennifer T Chayes, László Lovász, Vera T Sós, and Katalin Vesztergombi.  Convergent sequences of dense graphs i: Subgraph frequencies, metric properties and testing.  Advances in Mathematics, 219(6):1801–1851, 2008.
  * [8] Christian Borgs, Jennifer T Chayes, László Lovász, Vera T Sós, and Katalin Vesztergombi.  Convergent sequences of dense graphs ii. multiway cuts and statistical physics.  Annals of Mathematics, pages 151–219, 2012.
  * [9] Gilles Celeux, Florence Forbes, and Nathalie Peyrard.  Em procedures using mean field-like approximations for markov model-based image segmentation.  Pattern recognition, 36(1):131–144, 2003.
  * [10] Samprit Chatterjee and Jeffrey S Simonoff.  Handbook of regression analysis, volume 5.  John Wiley & Sons, 2013.
  * [11] Sourav Chatterjee, Persi Diaconis, and Allan Sly.  Random graphs with a given degree sequence.  The Annals of Applied Probability, pages 1400–1435, 2011.
  * [12] Sotirios P Chatzis and Theodora A Varvarigou.  A fuzzy clustering approach toward hidden markov random field models for enhanced spatially constrained image segmentation.  IEEE Transactions on Fuzzy Systems, 16(5):1351–1361, 2008.
  * [13] Peter Clifford, Sylvia Richardson, and Denis Hemon.  Assessing the significance of the correlation between two spatial processes.  Biometrics, pages 123–134, 1989.
  * [14] Constantinos Daskalakis, Nishanth Dikkala, and Ioannis Panageas.  Regression from dependent observations.  In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 881–889, 2019.
  * [15] Constantinos Daskalakis, Nishanth Dikkala, and Ioannis Panageas.  Logistic regression with peer-group effects via inference in higher-order ising models.  In International Conference on Artificial Intelligence and Statistics, pages 3653–3663. PMLR, 2020.
  * [16] Nabarun Deb and Sumit Mukherjee.  Fluctuations in mean-field ising models.  The Annals of Applied Probability, 33(3):1961–2003, 2023.
  * [17] Daniel DellaPosta, Yongren Shi, and Michael Macy.  Why do liberals drink lattes?  American Journal of Sociology, 120(5):1473–1511, 2015.
  * [18] Peter M DeMarzo, Dimitri Vayanos, and Jeffrey Zwiebel.  Persuasion bias, social influence, and unidimensional opinions.  The Quarterly journal of economics, 118(3):909–968, 2003.
  * [19] Pedro Domingos and Matt Richardson.  Mining the network value of customers.  In Proceedings of the seventh ACM SIGKDD international conference on Knowledge discovery and data mining, pages 57–66. ACM, 2001.
  * [20] Malcolm M Dow, Michael L Burton, Douglas R White, and Karl P Reitz.  Galton’s problem as network autocorrelation.  American Ethnologist, 11(4):754–770, 1984.
  * [21] Hugo Duminil-Copin.  Lectures on the ising and potts models on the hypercubic lattice.  arXiv preprint arXiv:1707.00520, 2017.
  * [22] Hugo Duminil-Copin and Vincent Tassion.  A new proof of the sharpness of the phase transition for bernoulli percolation and the ising model.  Communications in Mathematical Physics, 343(2):725–745, 2016.
  * [23] Pierre Dutilleul, Peter Clifford, Sylvia Richardson, and Denis Hemon.  Modifying the t test for assessing the correlation between two spatial processes.  Biometrics, pages 305–314, 1993.
  * [24] Bradley Efron and Trevor Hastie.  Computer age statistical inference, volume 5.  Cambridge University Press, 2016.
  * [25] Richard S Ellis and Charles M Newman.  The statistics of curie-weiss models.  Journal of Statistical Physics, 19(2):149–161, 1978.
  * [26] Philip A Ernst, Larry A Shepp, Abraham J Wyner, et al.  Yule’s “nonsense correlation” solved!  Annals of Statistics, 45(4):1789–1809, 2017.
  * [27] Sacha Friedli and Yvan Velenik.  Statistical mechanics of lattice systems: a concrete mathematical introduction.  Cambridge University Press, 2017.
  * [28] Geoffrey R Grimmett.  The random-cluster model, volume 333.  Springer Science & Business Media, 2006.
  * [29] Jürgen Groß.  Linear regression, volume 175.  Springer Science & Business Media, 2003.
  * [30] Dominique Guillot, Bala Rajaratnam, and Julien Emile-Geay.  Statistical paleoclimate reconstructions via Markov random fields.  The Annals of Applied Statistics, 9(1):324 – 352, 2015.
  * [31] Kenneth D Harris.  Nonsense correlations in neuroscience.  bioRxiv, pages 2020–11, 2020.
  * [32] RH Hooker.  On the correlation of successive observations.  Journal of the Royal Statistical Society, 68(4):696–703, 1905.
  * [33] Ernst Ising.  Beitrag zur theorie des ferromagnetismus.  Zeitschrift für Physik A Hadrons and Nuclei, 31(1):253–258, 1925.
  * [34] Vishesh Jain, Frederic Koehler, and Elchanan Mossel.  The mean-field approximation: Information inequalities, algorithms, and complexity.  In Conference On Learning Theory, pages 1326–1347. PMLR, 2018.
  * [35] Mark Kac.  Mathematical mechanisms of phase transitions.  Technical report, Rockefeller Univ., New York, 1969.
  * [36] Vardis Kandiros, Yuval Dagan, Nishanth Dikkala, Surbhi Goel, and Constantinos Daskalakis.  Statistical estimation from dependent data.  In International Conference on Machine Learning, pages 5269–5278. PMLR, 2021.
  * [37] Robert E Kass, Uri T Eden, and Emery N Brown.  Analysis of neural data, volume 491.  Springer, 2014.
  * [38] Morgan Kelly.  The standard errors of persistence.  2019\.
  * [39] Ross P Kindermann and J Laurie Snell.  On the relation between markov random fields and social networks.  Journal of Mathematical Sociology, 7(1):1–13, 1980.
  * [40] Joel L Lebowitz.  Coexistence of phases in ising ferromagnets.  Journal of Statistical Physics, 16(6):463–476, 1977.
  * [41] Youjin Lee and Elizabeth L Ogburn.  Network dependence can lead to spurious associations and invalid inference.  Journal of the American Statistical Association, 116(535):1060–1074, 2021.
  * [42] László Lovász.  Large networks and graph limits, volume 60.  American Mathematical Soc., 2012.
  * [43] Anders Martin-Löf.  Mixing properties, differentiability of the free energy and the central limit theorem for a pure phase in the ising model at low temperature.  Communications in Mathematical Physics, 32:75–92, 1973.
  * [44] Rajarshi Mukherjee, Sumit Mukherjee, and Ming Yuan.  Global testing against sparse alternatives under ising models.  The Annals of Statistics, 46(5):2062–2093, 2018.
  * [45] Matey Neykov and Han Liu.  Property testing in high-dimensional ising models.  The Annals of Statistics, 47(5):2472–2503, 2019.
  * [46] Hidetoshi Nishimori.  Statistical physics of spin glasses and information processing: an introduction.  Number 111. Clarendon Press, 2001.
  * [47] Elizabeth L Ogburn, Ilya Shpitser, and Youjin Lee.  Causal inference, social networks and chain graphs.  Journal of the Royal Statistical Society Series A: Statistics in Society, 183(4):1659–1676, 2020.
  * [48] Lars Onsager.  Crystal statistics. i. a two-dimensional model with an order-disorder transition.  Physical Review, 65(3-4):117, 1944.
  * [49] Peter CB Phillips.  Understanding spurious regressions in econometrics.  Journal of econometrics, 33(3):311–340, 1986.
  * [50] Marcos O. Prates, Douglas R. M. Azevedo, Ying C. MacNab, and Michael R. Willig.  Non-Separable Spatio-Temporal Models via Transformed Multivariate Gaussian Markov Random Fields.  Journal of the Royal Statistical Society Series C: Applied Statistics, 71(5):1116–1136, 05 2022.
  * [51] Simo Puntanen and George PH Styan.  The equality of the ordinary least squares estimator and the best linear unbiased estimator.  The American Statistician, 43(3):153–161, 1989.
  * [52] T Royama.  Moran effect on nonlinear population processes.  Ecological Monographs, 75(2):277–293, 2005.
  * [53] Helmut Schaeben.  A mathematical view of weights-of-evidence, conditional independence, and logistic regression in terms of markov random fields.  Mathematical Geosciences, 46(6):691–709, 2014.
  * [54] Aad W Van der Vaart.  Asymptotic statistics, volume 3.  Cambridge university press, 2000.
  * [55] Eric Vittinghoff, David V Glidden, Stephen C Shiboski, and Charles E McCulloch.  Regression methods in biostatistics: linear, logistic, survival, and repeated measures models.  2006\.
  * [56] Zhi Wei and Hongzhe Li.  A markov random field model for network-based analysis of genomic data.  Bioinformatics, 23(12):1537–1544, 2007.
  * [57] Robert West, Hristo S Paskov, Jure Leskovec, and Christopher Potts.  Exploiting social network structure for person-to-person sentiment analysis.  arXiv preprint arXiv:1409.2450, 2014.
  * [58] Yuanzhe Xu and Sumit Mukherjee.  Inference in Ising models on dense regular graphs.  The Annals of Statistics, 51(3):1183 – 1206, 2023.
  * [59] G Udny Yule.  On the time-correlation problem, with especial reference to the variate-difference correlation method.  Journal of the Royal Statistical Society, 84(4):497–537, 1921.
  * [60] G Udny Yule.  Why do we sometimes get nonsense-correlations between time-series?–a study in sampling and the nature of time-series.  Journal of the royal statistical society, 89(1):1–63, 1926.
  * [61] “STUDENT”.  Iv. the elimination of spurious correlation due to position in time or space.  Biometrika, 10(1):179–180, 1914.