# Ancestral reproductive bias in continuous time branching trees under various
sampling schemes

Jan Lukas Igelbrink  Jan Lukas Igelbrink
Institut für Mathematik, Johannes Gutenberg-Universität Mainz and Goethe-
Universität Frankfurt, Germany<EMAIL_ADDRESS>and  Jasper
Ischebeck  Jasper Ischebeck, Institut für Mathematik, Goethe-Universität
Frankfurt, Germany<EMAIL_ADDRESS>
###### Abstract.

Cheek and Johnston [CJ23] consider a continuous-time Bienaymé-Galton-Watson
tree conditioned on being alive at time $T$. They study the reproduction
events along the ancestral lineage of an individual randomly sampled from all
those alive at time $T$. We give a short proof of an extension of their main
results [CJ23, Theorems 2.3 and 2.4] to the more general case of Bellman-
Harris processes. Our proof also sheds light onto the probabilistic structure
of the rate of the reproduction events. A similar method will be applied to
explain (i) the different ancestral reproduction bias appearing in work by
Geiger [Gei99] and (ii) the fact that the sampling rule considered by Chauvin,
Rouault and Wakolbinger in [CRW91, Theorem 1] leads to a time homogeneous
process along the ancestral lineage.

###### Key words and phrases:

branching processes, spines, reproductive bias, inspection paradox,sampling
schemes

###### 2020 Mathematics Subject Classification:

Primary 60J80; secondary 60K05, 92D10

We thank Anton Wakolbinger for bringing the work [CJ23] to our attention. We
are grateful to him and also to Matthias Birkner, Götz Kersting and Marius
Schmidt for stimulating discussions and valuable hints. A substantial part of
this work was done during the 2023 seminar week of the Frankfurt probability
group in Haus Bergkranz.

##  1\. Introduction

Consider a continuous-time branching process with $N_{\\_}t$ individuals alive
at time $t$, started with one individual at time $0$. At the end of its
lifetime, an individual is replaced by a random number of offspring with
distribution $\left(p_{\\_}k\right)_{\\_}{k\geq 0}$. Cheek and Johnston [CJ23]
study the process of reproduction times and family sizes along the ancestral
lineage of an individual sampled from all those alive at a given time $T>0$,
conditioned on the event $\\{N_{\\_}T>0\\}$. While [CJ23] treat the (Bienaymé-
Galton-Watson-) case of i.i.d. exponentially distributed lifetimes, we
consider here the more general case of Bellman-Harris processes in which the
lifetimes are i.i.d. with an arbitrary distribution $\mu$ on
$\mathbb{R}_{\\_}+$, and give a short and conceptual probabilistic proof of
the main results of [CJ23]. The core idea of this proof is as follows:

On the event $\left\\{N_{\\_}T>0\right\\}$, assign independent, uniformly on
$[0,1]$ distributed markers to the individuals alive at time $T$. Then the
individual whose marker is largest constitutes a uniformly distributed random
pick from all the individuals alive at time $T$. As we will see, the argument
$s$ of the generating functions that appear in the analytic arguments of
[CJ23] corresponds to the realisation of the largest marker. Sections 2 to 4
will be devoted to formulating and proving Theorem 2.1.

Relating to work of Chauvin, Rouault and Wakolbinger [CRW91], in Section 5 we
will consider the case of potentially dependent but identically and atomless
distributed markers and conditioning on one marker taking the prescribed value
$s$. In contrast to the above, in this case one does not observe a time-
inhomogeneity along the sampled ancestral lineage.

In Section 6 we will consider a planar embedding of the Bellman-Harris tree
conditioned to survive up to time $T$, and analyse the leftmost ancestral
lineage among those surviving until time $T$. Here we follow Geiger [Gei99],
who gave a representation of discrete-time Galton-Watson processes conditioned
to survive up to a given number of generations. With this sampling rule we
observe a time-inhomogeneity of ancestral reproduction events that is
different from the one in [CJ23].

In Section 7 we briefly resume the discussion from [CJ23] on a possible
relation between the ancestral rate bias and the rate of mutations per cell
division in embryogenesis, and illustrate the various sampling schemes from a
more biological perspective.

##  2\. Sampling an ancestral line at random

On the event $\\{N_{\\_}T>0\\}$, let the individual $V$ be sampled as
described in the Introduction, and let $S$ be its mark. We define the process
$(N_{\\_}t)_{\\_}{t\geq 0}$ to be right continuous with left limits. As a
consequence, if $T_{\\_}1$ is the lifetime of the root individual, then
$N_{\\_}{T_{\\_}1}$ has distribution $\left(p_{\\_}k\right)_{\\_}{k\geq 0}$.
Let $J$ be the random number of reproduction events and
$0<T_{\\_}1<T_{\\_}2<\cdots<T_{\\_}J\leq T$ be the random times of
reproduction events along the ancestral lineage of $V$. Let
$L_{\\_}1,\ldots,L_{\\_}J$ be the offspring sizes in these reproduction events
and let $0<\tau_{\\_}1<\tau_{\\_}2<\cdots$ be the random arrival times in a
renewal process with interarrival time distribution $\mu$. Denote by
$\operatorname{\mathbf{P}}$ and $\operatorname{\mathbf{E}}$ the probability
measure and expectation for $N_{\\_}0=1$.

######  Theorem 2.1.

For $j\geq 0$, $0<t_{\\_}1<\ldots<t_{\\_}j\leq T\in\mathbb{R}$ and
$\ell_{\\_}1,\ldots,\ell_{\\_}j\in\mathbb{N}$ we have

(2.1) $\displaystyle\begin{aligned}
&\operatorname{\mathbf{P}}\left(N_{\\_}T>0,J=j,\,T_{\\_}1\in\operatorname{d\\!}t_{\\_}1,\ldots
T_{\\_}j\in\operatorname{d\\!}t_{\\_}j,\,L_{\\_}1=\ell_{\\_}1,\ldots,L_{\\_}j=\ell_{\\_}j,\,S\in\operatorname{d\\!}s\right)\\\
&=\operatorname{\mathbf{P}}\left(\tau_{\\_}1\in\operatorname{d\\!}t_{\\_}1,\ldots,\tau_{\\_}j\in\operatorname{d\\!}t_{\\_}j,\tau_{\\_}{j+1}>T\right)\prod_{\\_}{i=1}^{j}\left(\ell_{\\_}ip_{\\_}{\ell_{\\_}i}\operatorname{\mathbf{E}}\left[s^{N_{\\_}{T-t_{\\_}i}}\right]^{\ell_{\\_}i-1}\right)\operatorname{d\\!}{s}.\end{aligned}$

######  Corollary 2.2.

When integrated over $s\in(0,1)$, (LABEL:fullformula) reveals that the process
$(T_{\\_}1,L_{\\_}1),\ldots,(T_{\\_}J,L_{\\_}J)$ of reproduction times and
offspring sizes along the ancestral lineage of the uniformly chosen individual
(conditioned on $\\{N_{\\_}T>0\\}$) is a mixture of (what could be called)
“biased compound renewal processes”.

######  Remark 2.3.

  1. a)

When the lifetime distribution $\mu$ is the exponential distribution with
parameter $r$, then $\tau_{\\_}1,\tau_{\\_}2,\ldots$ are the points of a rate
$r$ Poisson point process. In this case Corollary 2.2 together with
(LABEL:fullformula) becomes a reformulation of the statements of [CJ23,
Theorems 2.3 and 2.4], and at the same time reveals the probabilistic role of
the mixing parameter $s$ in the mixture of biased compound Poisson processes
that appear in the “Cox process representation” of [CJ23].

Let us write (as in [CJ23]) $F_{\\_}t(s)\mathrel{\mathop{\mathchar
58\relax}}=\mathbf{E}[s^{N_{\\_}t}]$, and abbreviate

(2.2) $B(t,T,\ell)\mathrel{\mathop{\mathchar
58\relax}}=\frac{1}{1-F_{\\_}T(0)}\int_{\\_}0^{1}F_{\\_}{T-t}(s)^{\ell-1}F^{\prime}_{\\_}T(s)\operatorname{d\\!}s.$

[CJ23, Theorem 2.4] (as well as Theorem 2.1) says that the rate of size $\ell$
reproduction along the uniform ancestral lineage at time $t$ is $r\ell
p_{\\_}\ell\,B(t,T,\ell)$. In this sense the factor $B(t,T,\ell)$ can be
interpreted as an (ancestral) rate bias, on top of the classical term $r\ell
p_{\\_}\ell$. Indeed, the factor $B(t,T,\ell)$ is absent in trees that are
biased with respect to their size at time $T$. Galton-Watson trees of this
kind have been investigated (also in the multitype case) by Georgii and Baake
[GB03, Section 4]; they are continuous-time analogues of the size-biased trees
analysed by Lyons et al. [LPP95] and Kurtz et al. [KLPP97].

In the critical and supercritical case one can check that, for all fixed $u<T$
and $\ell\in\mathbb{N}$ one has the convergence $B(T-u,T,\ell)\to 1$ as
$T\to\infty$. In the supercritical case this stabilisation along the sampled
ancestral lineage corresponds to the “retrospective viewpoint” that has been
taken in [GB03] and, in the more general situation of Crump-Mode-Jagers
processes, by Jagers and Nerman [JN96].

  2. b)

The choice $\mu=\delta_{\\_}1$ renders the case of discrete time Galton-Watson
processes, starting with one individual at time $0$ and with reproduction
events at times $1,2,\ldots$. Then, with $T=n\in\mathbb{N}$, and
$L_{\\_}1,\ldots,L_{\\_}{n}$ being the family sizes along the ancestral
lineage of the sampled individual $V$, the formula (LABEL:fullformula)
specialises to

(2.3)
$\displaystyle\operatorname{\mathbf{P}}\left(N_{\\_}n>0,\,L_{\\_}1=\ell_{\\_}1,\ldots,L_{\\_}{n}=\ell_{\\_}{n},\,S\in\operatorname{d\\!}s\right)=\left(\prod_{\\_}{i=1}^{n}\ell_{\\_}ip_{\\_}{\ell_{\\_}i}\operatorname{\mathbf{E}}\left[s^{N_{\\_}{n-i}}\right]^{\ell_{\\_}i-1}\right)\operatorname{d\\!}{s}.$

##  3\. Maxima of i.i.d. random markers

As a preparation for the short probabilistic proof of Theorem 2.1 given in the
next section, we recall the following well-know fact:

For $\ell\in\mathbb{N}$, let $\widetilde{S}$ be the maximum of $\ell$
independent Unif$[0,1]$-distributed random variables
$U_{\\_}1,\ldots,U_{\\_}\ell$. Then the density of $\widetilde{S}$ is

(3.1)
$\operatorname{\mathbf{P}}\left(\widetilde{S}\in\operatorname{d\\!}s\right)=\ell
s^{\ell-1}\operatorname{d\\!}s,\quad 0\leq s\leq 1.$

Indeed, because of exchangeability,

$\operatorname{\mathbf{P}}\left(\widetilde{S}\in\operatorname{d\\!}s\right)=\ell\operatorname{\mathbf{P}}\left(U_{\\_}1\in\operatorname{d\\!}s\right)\operatorname{\mathbf{P}}\left(U_{\\_}2<s,\ldots,U_{\\_}\ell<s\right),$

which equals the r.h.s. of (3.1).

The following lemma specialises to (3.1) when putting $\tilde{N}\equiv 1$.

######  Lemma 3.1.

Let $\widetilde{N}$ be an $\mathbb{N}_{\\_}0$-valued random variable, and
$\widetilde{N}_{\\_}1,\widetilde{N}_{\\_}2,\ldots$ be i.i.d. copies of
$\tilde{N}$. Given $\widetilde{N}_{\\_}1,\widetilde{N}_{\\_}2,\ldots$ let
$U_{\\_}{1,1},\ldots U_{\\_}{1,\widetilde{N}_{\\_}1},U_{\\_}{2,1},\ldots
U_{\\_}{2,\widetilde{N}_{\\_}2},\ldots$ be independent Unif
$[0,1]$-distributed random variables, and write

$\displaystyle S_{\\_}k$ $\displaystyle\mathrel{\mathop{\mathchar 58\relax}}=$
$\displaystyle\max\left\\{U_{\\_}{k,1},\ldots,U_{\\_}{k,\widetilde{N}_{\\_}k}\right\\},\quad
k=1,2,\ldots$ $\displaystyle S^{(\ell)}$
$\displaystyle\mathrel{\mathop{\mathchar 58\relax}}=$
$\displaystyle\max\left\\{S_{\\_}1,\ldots,S_{\\_}\ell\right\\},\quad\ell\in\mathbb{N}$

where we put $\max(\emptyset)\mathrel{\mathop{\mathchar 58\relax}}=-\infty$.
Then, for all $\ell\in\mathbb{N}$, the density of $S^{(\ell)}$ is

(3.2)
$\operatorname{\mathbf{P}}\left(\widetilde{N}_{\\_}1+\ldots+\widetilde{N}_{\\_}\ell>0,\,S^{(\ell)}\in\operatorname{d\\!}s\right)=\ell\,\operatorname{\mathbf{E}}\left[s^{\widetilde{N}}\right]^{\ell-1}\operatorname{\mathbf{P}}\left(\widetilde{N}_{\\_}1>0,\,S_{\\_}1\in\operatorname{d\\!}s\right),\quad
0\leq s\leq 1.$

###### Proof.

Again because of exchangeability, the l.h.s. of (3.2) equals

(3.3)
$\ell\operatorname{\mathbf{P}}\left(\widetilde{N}_{\\_}1>0,\,S_{\\_}1\in\operatorname{d\\!}s\right)\operatorname{\mathbf{P}}\left(S_{\\_}2<s,\ldots,S_{\\_}\ell<s\right).$

Since by assumption the $S_{\\_}k$ are i.i.d. copies of $S_{\\_}1$, the
rightmost factor in (3.3) equals

$\operatorname{\mathbf{P}}\left(S_{\\_}1<s\right)^{\ell-1}=\operatorname{\mathbf{E}}\left[\operatorname{\mathbf{P}}\left(S_{\\_}1<s\mid\widetilde{N}_{\\_}1\right)\right]^{\ell-1}=\operatorname{\mathbf{E}}\left[s^{\widetilde{N}_{\\_}1}\right]^{\ell-1}=\operatorname{\mathbf{E}}\left[s^{\widetilde{N}}\right]^{\ell-1}.$

Hence, (3.3) equals the r.h.s. of (3.2), completing the proof of the lemma. ∎

The following corollary is immediate.

######  Corollary 3.2.

Let $L$ be an $\mathbb{N}_{\\_}0$-valued random variable that is independent
of all the random variables appearing in Lemma 3.1, with
$\operatorname{\mathbf{P}}(L=\ell)=p_{\\_}\ell$, $\ell\in\mathbb{N}_{\\_}0$.
Then we have for all $\ell\in\mathbb{N}_{\\_}0$,

(3.4)
$\operatorname{\mathbf{P}}\left(L=\ell,\,\widetilde{N}_{\\_}1+\ldots+\widetilde{N}_{\\_}\ell>0,\,S^{(\ell)}\in\operatorname{d\\!}s\right)=\ell
p_{\\_}\ell\,\operatorname{\mathbf{E}}\left[s^{\widetilde{N}}\right]^{\ell-1}\operatorname{\mathbf{P}}\left(\widetilde{N}_{\\_}1>0,\,S_{\\_}1\in\operatorname{d\\!}s\right),\quad
0\leq s\leq 1.$

##  4\. Proof of Theorem 2.1

For $j=0$, both sides of (2.1) are equal to
$\mu((T,\infty))\,\operatorname{d\\!}s$. For $j\geq 1$, a decomposition at the
first reproduction event of the branching process (which on the event
$\\{N_{\\_}T>0\\}$ necessarily is also the first reproduction event along the
sampled ancestral lineage) leads us directly to the situation of Corollary
3.2. Here, the random variable $L$ in Corollary 3.2 takes the role of the
$L_{\\_}1$ in (2.1), and the random variable $N_{\\_}{T-t_{\\_}1}$ from (2.1)
becomes the $\tilde{N}$ in Corollary 3.2. We thus obtain from Corollary 3.2
for $0<t_{\\_}1\leq T$, $0\leq s\leq 1$, $\ell_{\\_}1\in\mathbb{N}$

(4.1)
$\displaystyle\begin{split}&\operatorname{\mathbf{P}}\left(N_{\\_}T>0,J\geq
1,\,T_{\\_}1\in\operatorname{d\\!}t_{\\_}1,L_{\\_}1=\ell_{\\_}1,S\in\operatorname{d\\!}s\right)\\\
&=\operatorname{\mathbf{P}}\left(\tau_{\\_}1\in\operatorname{d\\!}t_{\\_}1\right)\ell_{\\_}1p_{\\_}{\ell_{\\_}1}\operatorname{\mathbf{E}}\left[s^{N_{\\_}{T-t_{\\_}1}}\right]^{\ell_{\\_}1-1}\operatorname{\mathbf{P}}\left(N_{\\_}{T-t_{\\_}1}>0,\,S_{\\_}1\in\operatorname{d\\!}s\right),\end{split}$

where, on the event $\\{N_{\\_}{T-t_{\\_}1}>0\\}$, $S_{\\_}1$ is the largest
among the marks assigned in an i.i.d. Unif$[0,1]$ manner to the
$N_{\\_}{T-t_{\\_}1}$ many individuals. Thanks to the independence and self-
similarity properties inherent in the branching processes, this can be
iterated, leading directly to (2.1).

##  5\. Conditioning on a marker value

Chauvin, Rouault and Wakolbinger [CRW91] consider a Markov process with an
atomless transition probability indexed by a continuous-time Galton-Watson-
tree and condition on an individual at time $T$ to be at a given location.

To relate this to the framework described in the Introduction, we assume that
each individual alive at time $T$ in the Bellmann-Harris tree carries a mark
in some standard Borel space $E$ and these random marks have the following
properties:

  1. (M1)

Their marginal distributions (denoted by $\nu$) are identical and do not
depend on the reproduction events

  2. (M2)

a.s. no pair of marks is equal.

Think for example of branching Brownian motion: The positions of the particles
depend on each other via the genealogy, but the movements after a branching
event are independent. At every fixed point in time a.s. no pair of particles
will be at the same position.

We now condition on $\left\\{N_{\\_}T>0\right\\}$ and, for given $s\in E$, on
one of the $N_{\\_}T$ individuals having marker value $s$. Denote by $V$ the
individual having marker $s$. Let $J$ be the random number of reproduction
events along the ancestral lineage of $V$ and
$0<T_{\\_}1<T_{\\_}2<\cdots<T_{\\_}J<T$ be the random times of these
reproduction events. Let $L_{\\_}1,\ldots,L_{\\_}J$ be the offspring sizes in
these reproduction events and let $0<\tau_{\\_}1<\tau_{\\_}2<\cdots$ be the
random arrival times in a renewal process with interarrival time distribution
$\mu$.

######  Theorem 5.1.

For $j\geq 0$, $0<t_{\\_}1<\ldots<t_{\\_}j<T$ and
$\ell_{\\_}1,\ldots,\ell_{\\_}j\in\mathbb{N}$ we have for $\nu$-almost all $s$

(5.1)
$\displaystyle\begin{split}&\operatorname{\mathbf{P}}\left(\left.J=j,\,T_{\\_}1\in\operatorname{d\\!}t_{\\_}1,\ldots
T_{\\_}j\in\operatorname{d\\!}t_{\\_}j,\,L_{\\_}1=\ell_{\\_}1,\ldots,L_{\\_}j=\ell_{\\_}j\right|N_{\\_}T>0,\exists\,\mathrm{mark}\in\operatorname{d\\!}s\right)\\\
&=\frac{1}{\operatorname{\mathbf{E}}[N_{\\_}T]}\operatorname{\mathbf{P}}\left(\tau_{\\_}1\in\operatorname{d\\!}t_{\\_}1,\ldots,\tau_{\\_}j\in\operatorname{d\\!}t_{\\_}j,\tau_{\\_}{j+1}\geq
T\right)\prod_{\\_}{i=1}^{j}\ell_{\\_}{i}p_{\\_}{\ell_{\\_}i}.\end{split}$

###### Proof.

Because of properties (M1), (M2) we have

$\operatorname{\mathbf{P}}(N_{\\_}T>0,\exists\,\mathrm{mark}\in\operatorname{d\\!}s)=\operatorname{\mathbf{E}}[N_{\\_}T]\nu(\operatorname{d\\!}s),\quad
s\in E.$

Hence (LABEL:fullformula2cw) is equivalent to

(5.2)
$\displaystyle\begin{split}&\operatorname{\mathbf{P}}\left(J=j,\,T_{\\_}1\in\operatorname{d\\!}t_{\\_}1,\ldots
T_{\\_}j\in\operatorname{d\\!}t_{\\_}j,\,L_{\\_}1=\ell_{\\_}1,\ldots,L_{\\_}j=\ell_{\\_}j,N_{\\_}T>0,\exists\,\mathrm{mark}\in\operatorname{d\\!}s\right)\\\
&=\operatorname{\mathbf{P}}\left(\tau_{\\_}1\in\operatorname{d\\!}t_{\\_}1,\ldots,\tau_{\\_}j\in\operatorname{d\\!}t_{\\_}j,\tau_{\\_}{j+1}\geq
T\right)\prod_{\\_}{i=1}^{j}\ell_{\\_}{i}p_{\\_}{\ell_{\\_}i}\,\nu(\operatorname{d\\!}s).\end{split}$

We prove the statement (LABEL:fullformula3cw) by induction over $j$,
simultaneously over all time horizons $T>0$. We write $\mathbf{P}^{T}$ for the
probability referring to time horizon $T$; this well be helpful in the
induction step where will encounter two different time horizons. For $j=0$ the
statement is true, since

$\operatorname{\mathbf{P}}^{T}(J=0,N_{\\_}T>0,\exists\,\mathrm{mark}\in\operatorname{d\\!}s)=\operatorname{\mathbf{P}}\left(\tau_{\\_}1\leq
T\right)\,\nu(ds).$

Assume we have proved (LABEL:fullformula3cw) for all time horizons
$T^{\prime}$ with $j-1$ (in place of $j$), for all times
$t_{\\_}1^{\prime},\dots,t_{\\_}{j-1}^{\prime}\leq T^{\prime}$, sizes
$\ell_{\\_}1^{\prime},\dots,\ell_{\\_}{j-1}^{\prime}\in\mathbb{N}$ and marker
distributions with the same marginal $\nu$ that satisfy conditions (M1), (M2).
Turning to (LABEL:fullformula3cw) as it stands, we note that on
$\\{T_{\\_}1=t_{\\_}1,L_{\\_}1=\ell_{\\_}1\\}$, the descendants of the
$\ell_{\\_}1$ siblings in the first branching event form $\ell_{\\_}1$
independent and identically distributed trees on the time interval
$[t_{\\_}1,T]$. Let $\mathcal{U}_{\\_}k,\,k=1,\dots,\ell_{\\_}1$, be the set
of markers of the individuals at time $T$ that descend from the $k$-th
sibling. By randomly permuting these $\ell_{\\_}1$ siblings, we can assume
that the set-valued random variables
$\mathcal{U}_{\\_}k,\,k=1,\ldots,\ell_{\\_}1$, are exchangeable. Note that the
markers in each $\mathcal{U}_{\\_}k$ satisfy conditions (M1), (M2). Because
the markers are a.s. pairwise different by assumption, the mark $s$ belongs to
at most one of those $\mathcal{U}_{\\_}k$, so

(5.3)
$\mathbf{1}_{\\_}{\\{\exists\,\mathrm{mark}\in\operatorname{d\\!}s\\}}=\sum_{\\_}{k=1}^{\ell_{\\_}1}\mathbf{1}_{\\_}{\\{\mathcal{U}_{\\_}k\cap\operatorname{d\\!}s\neq\emptyset\\}}\,\mbox{
a.s.}$

Putting $t_{\\_}1^{\prime}\mathrel{\mathop{\mathchar
58\relax}}=t_{\\_}2-t_{\\_}1,\ldots,t_{\\_}{j-1}^{\prime}\mathrel{\mathop{\mathchar
58\relax}}=t_{\\_}j-t_{\\_}1$ we thus infer, using the branching property of
the Bellman-Harris tree, that the left hand side of (LABEL:fullformula3cw)
equals

(5.4)
$\operatorname{\mathbf{P}}(\tau_{\\_}1\in\operatorname{d\\!}t_{\\_}1)p_{\\_}{\ell_{\\_}1}\ell_{\\_}1\operatorname{\mathbf{P}}^{T-t_{\\_}1}\left(J=j-1,\ldots,T_{\\_}1\in\operatorname{d\\!}t_{\\_}1^{\prime},\ldots,T_{\\_}{j-1}\in\operatorname{d\\!}t_{\\_}{j-1}^{\prime},N_{\\_}{T-t_{\\_}1}>0,\exists\,{\rm
mark}\in\operatorname{d\\!}s\right)$

By the induction assumption this is equal to

(5.5)
$\operatorname{\mathbf{P}}(\tau_{\\_}1\in\operatorname{d\\!}t_{\\_}1)p_{\\_}{\ell_{\\_}1}\ell_{\\_}1\operatorname{\mathbf{P}}\left(\tau_{\\_}1^{\prime}\in\operatorname{d\\!}t_{\\_}1^{\prime},\ldots,\tau_{\\_}{j-1}^{\prime}\in\operatorname{d\\!}t_{\\_}{j-1}^{\prime},\tau_{\\_}{j}^{\prime}\geq
T-t_{\\_}1\right)\prod_{\\_}{i=2}^{j}\ell_{\\_}{i}p_{\\_}{\ell_{\\_}i}\,\nu(\operatorname{d\\!}s),$

where $(\tau_{\\_}1^{\prime},\tau_{\\_}2^{\prime},\ldots)$ have the same
distribution as $(\tau_{\\_}1,\tau_{\\_}2,\ldots)$. Obviously (5.5) equals the
r.h.s. of (LABEL:fullformula3cw), which completes the induction step and
concludes the proof. ∎

######  Remark 5.2.

If $\mu$ is the exponential distribution with parameter $r$, then
$\tau_{\\_}1,\tau_{\\_}2,\ldots$ are again the points of a rate $r$ Poisson
point process and (LABEL:fullformula2cw) implies that reproduction events
along the ancestral lineage of $V$ happen according to a time-homogeneous
Poisson process with rate $r\sum_{\\_}{\ell}\ell p_{\\_}{\ell}$. This
corresponds to the description of the events along the ancestral line of $V$
given in [CRW91, Theorem 1].

##  6\. Sampling the left-most ancestral lineage

We now aim to obtain results about what Geiger [Gei99] calls the leftmost
surviving ancestral lineage in a planar embedding of the tree: At any
reproduction event we assign independent uniformly on $[0,1]$ distributed
markers to all children. An individual can now be uniquely determined by the
markers along its ancestral lineage. On the event $\\{N_{\\_}T>0\\}$, let $V$
be the individual whose markers along the entire ancestral lineage comes first
in the lexicographic ordering. Let $J$ be the random number of reproduction
events and $0<T_{\\_}1<T_{\\_}2<\cdots<T_{\\_}J\leq T$ be the random times of
reproduction events along the ancestral lineage of $V$. Let
$L_{\\_}1,\ldots,L_{\\_}J$ be the offspring sizes in these reproduction events
and let $0<\tau_{\\_}1<\tau_{\\_}2<\cdots$ be the random arrival times in a
renewal process with interarrival time distribution $\mu$. Denote by
$K_{\\_}i$ the number of siblings born at reproduction event number $i$ along
the ancestral lineage of $V$ which have a lower lexicographic order than $V$
and whose descendants hence die out before time $T$.

######  Theorem 6.1.

For $j\geq 0$,
$0<t_{\\_}1<\ldots<t_{\\_}j<T,\,\ell_{\\_}1,\ldots,\ell_{\\_}j\in\mathbb{N}$
and $k_{\\_}i\in\left\\{1,\ldots,\ell_{\\_}i-1\right\\}$ we have

(6.1)
$\begin{split}&\operatorname{\mathbf{P}}\left(N_{\\_}T>0,J=j,\,T_{\\_}1\in\operatorname{d\\!}t_{\\_}1,\ldots
T_{\\_}j\in\operatorname{d\\!}t_{\\_}j,\,L_{\\_}1=\ell_{\\_}1,\ldots,L_{\\_}j=\ell_{\\_}j,K_{\\_}1=k_{\\_}1,\ldots,K_{\\_}j=k_{\\_}j\right)\\\
&=\operatorname{\mathbf{P}}\left(\tau_{\\_}1\in\operatorname{d\\!}t_{\\_}1,\ldots,\tau_{\\_}j\in\operatorname{d\\!}t_{\\_}j,\tau_{\\_}{j+1}\geq
T\right)\prod_{\\_}{i=1}^{j}\left(p_{\\_}{\ell_{\\_}i}\operatorname{\mathbf{P}}\left(N_{\\_}{T-t_{\\_}i}=0\right)^{k_{\\_}i}\right).\end{split}$

###### Proof.

The proof of the theorem works in analogy to the one of Theorem 2.1, but using
following analogue of Lemma 3.1. ∎

######  Lemma 6.2.

Let $\widetilde{N}$ be an $\mathbb{N}_{\\_}0$-valued random variable, and
$\widetilde{N}_{\\_}1,\widetilde{N}_{\\_}2,\ldots$ be i.i.d. copies of
$\tilde{N}$. Given $\widetilde{N}_{\\_}1,\widetilde{N}_{\\_}2,\ldots$ let
$U_{\\_}{1},U_{\\_}{2},\ldots$ be independent Unif $[0,1]$-distributed random
variables, and write

$\displaystyle S^{(\ell)}$ $\displaystyle\mathrel{\mathop{\mathchar
58\relax}}=$ $\displaystyle\max\left\\{U_{\\_}{k}\mid N_{\\_}k\geq
1,k=1,\ldots,\ell\right\\},$ $\displaystyle K^{(\ell)}$
$\displaystyle\mathrel{\mathop{\mathchar 58\relax}}=$
$\displaystyle\left|\left\\{U_{\\_}k\mid
U_{\\_}k>S^{(\ell)},k=1,\ldots,\ell\right\\}\right|$

where we put $\max(\emptyset)\mathrel{\mathop{\mathchar 58\relax}}=-\infty$.
Then, for all $k<\ell\in\mathbb{N}$ we have

(6.2)
$\operatorname{\mathbf{P}}\left(\widetilde{N}_{\\_}1+\ldots+\widetilde{N}_{\\_}\ell>0,\,K^{(\ell)}=k\right)=\operatorname{\mathbf{P}}\left(\widetilde{N}=0\right)^{k}\operatorname{\mathbf{P}}\left(\widetilde{N}>0\right).$

###### Proof.

Because $S^{(\ell)}$ and $K^{(\ell)}$ are symmetric in
$U_{\\_}1,\dots,U_{\\_}\ell$, we can use exchangeability to assume that
$U_{\\_}1>U_{\\_}2>\dots>U_{\\_}\ell$. For $K^{(\ell)}$ to be $k$,
$S^{(\ell)}$ has then to be $U_{\\_}{k+1}$. This is exactly the case if
$\widetilde{N}_{\\_}1,\dots,\widetilde{N}_{\\_}k=0$ and
$\widetilde{N}_{\\_}{k+1}>0$. ∎

##  7\. Biological perspectives

Cheek and Johnston [CJ23, Section 5] discuss recent studies ([PMK+21],
[CTAS+19]) which suggest that certain mutation rates are elevated for the
earliest cell divisions in embryogenesis. Under the assumptions that (1) cell
division times vary and (2) mutations arise not only at but also between cell
divisions, Cheek and Johnston argue that this early rate elevation might be
parsimoniously explained by their finding that in the supercritical case with
no deaths the rate of branching events along a uniformly chosen ancestral
lineage is increasing in $t\in[0,T]$ (which is a corollary to their Theorem
2.4).

The two-stage sampling rule first sample a random tree (“an adult”) that
survives up to time $T$, then sample an individual from this tree (“a cell
from this adult”) at time $T$ seems adequate for the situation discussed in
Cheek and Johnston [CJ23, Section 5]. In other modeling situations, again with
a large collection of i.i.d. Galton-Watson trees, one may think of a different
sampling rule: choose individuals at time $T$ uniformly at random from the
union of all the trees. This makes it more probable that the sampled
individuals belong to larger trees, and in fact corresponds to the size-
biasing of the random trees at time $T$ ([GB03, Section 4]) which we discussed
in Remark 2.3 a). As mentioned there, the rate bias (2.2) is not present in
this sampling scheme.

As can be seen from [CRW91, Theorem 1] (and Theorem 5.1), the rate bias (2.2)
is also absent along the ancestral lineage of an individual whose marker has a
prescribed value $s$, if one considers a situation in which a neutral marker
evolves along the trees in small (continuous) mutation steps, and if one
takes, for the prescribed value $s$, the collection of trees so large that one
individual at time $T$ has a marker value close to (ideally: precisely at)
$s$.

The sampling rule that appears in [Gei99] (and Theorem 6.1) leads to a rate
(and reproduction size) bias along the ancestral lineage that is different
from the ones we just discussed. This sampling rule can be defined via i.i.d.
real-valued valued neutral markers that are created at each birth and passed
to the offspring. The individual sampled at time $T$ (from the tree
conditioned to survive up to time $T$) is the one whose marker sequence is the
largest in lexicographic order among the individuals that live in the tree at
time $T$. This interpretation appears of less biological relevance, except in
the pure birth (or cell division) case, where one might think of one single
marker that is passed on in each generation to a randomly chosen daughter
cell.

## References

  * [CJ23] D. Cheek and S. G. G. Johnston.  Ancestral reproductive bias in branching processes.  Journal of Mathematical Biology, 86(5):70, 2023.  doi:10.1007/s00285-023-01907-7.
  * [CRW91] B. Chauvin, A. Rouault, and A. Wakolbinger.  Growing conditioned trees.  Stochastic Processes and their Applications, 39(1):117–130, 1991\.  doi:10.1016/0304-4149(91)90036-C.
  * [CTAS+19] T. H. Coorens, T. D. Treger, R. Al-Saadi, L. Moore, M. G. Tran, T. J. Mitchell, S. Tugnait, C. Thevanesan, M. D. Young, T. R. Oliver, et al.  Embryonal precursors of Wilms tumor.  Science, 366(6470):1247–1251, 2019.  doi:10.1126/science.aax1323.
  * [GB03] H.-O. Georgii and E. Baake.  Supercritical multitype branching processes: the ancestral types of typical individuals.  Advances in Applied Probability, 35(4):1090–1110, 2003.  doi:10.1239/aap/1067436336.
  * [Gei99] J. Geiger.  Elementary new proofs of classical limit theorems for Galton-Watson processes.  Journal of Applied Probability, 36(2):301–309, 1999.  doi:10.1239/jap/1032374454.
  * [JN96] P. Jagers and O. Nerman.  The asymptotic composition of supercritical, multi-type branching populations.  Séminaire de probabilités de Strasbourg, 30:40–54, 1996\.  doi:10.1007/BFb0094640.
  * [KLPP97] T. Kurtz, R. Lyons, R. Pemantle, and Y. Peres.  A conceptual proof of the Kesten-Stigum theorem for multi-type branching processes, pages 181–185.  Springer New York, New York, NY, 1997.  doi:10.1007/978-1-4612-1862-314.
  * [LPP95] R. Lyons, R. Pemantle, and Y. Peres.  Conceptual proofs of L log L criteria for mean behavior of branching processes.  The Annals of Probability, 23(3):1125–1138, 1995.  doi:10.1214/aop/1176988176.
  * [PMK+21] S. Park, N. M. Mali, R. Kim, J.-W. Choi, J. Lee, J. Lim, J. M. Park, J. W. Park, D. Kim, T. Kim, et al.  Clonal dynamics in early human embryogenesis inferred from somatic mutation.  Nature, 597(7876):393–397, 2021.  doi:10.1038/s41586-021-03786-8.


$\begin{array}[]{lllllllllllll}\left|\mathbb{P}(\bm{S}_{\mathscr{N}}\in\mathscr{A})-\Phi(\bm{Z}\in\mathscr{A})\right|&\leq&\dfrac{83}{\xi_{\epsilon^{\star}}^{3/2}}\,\sqrt{\dfrac{p^{7/2}}{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}}+\dfrac{4}{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}+\dfrac{8\,\left[D_{g}\right]^{+}}{\left(\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}\right)^{2}}.\end{array}$

$\hfill\blacksquare$

Note that the asymptotic multivariate normality can be established provided

$\begin{array}[]{lllllllllllll}\lim\limits_{N\to\infty}\,\left[\dfrac{83}{\xi_{\epsilon^{\star}}^{3/2}}\,\sqrt{\dfrac{p^{7/2}}{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}}+\dfrac{4}{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}+\dfrac{8\,\left[D_{g}\right]^{+}}{\left(\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}\right)^{2}}\right]&=&0,\end{array}$

resulting in following the asymptotic convergence in distribution:

$\begin{array}[]{lllllllllllll}\bm{S}_{\mathscr{N}}&\overset{D}{\longrightarrow}&\bm{Z}&\sim&\text{MvtNorm}\left(\textbf{0},\,\bm{I}_{p}\right).\end{array}$

##  Appendix F Proof of Theorem 2

Proof of Theorem 2. In order to demonstrate the feasibility of the normal
approximation for maximum likelihood estimators $\widehat{\bm{\theta}}$ of
$\bm{\theta}^{\star}$, we first start with a standard Taylor expansion of the
negative score equation:

$\begin{array}[]{lllllllllllllllll}-\nabla_{\bm{\theta}}\,\ell(\widehat{\bm{\theta}};\bm{x},\bm{y})&=&-\nabla_{\bm{\theta}}\,\ell(\bm{\theta}^{\star};\bm{x},\bm{y})-\nabla_{\bm{\theta}}^{2}\,\ell(\bm{\theta}^{\star};\bm{x},\bm{y})\,(\widehat{\bm{\theta}}-\bm{\theta}^{\star})-\bm{R},\end{array}$
(15)

where $\bm{R}\in\mathbb{R}^{p}$ is the vector of remainders given in the
Lagrange form. Denoting by $R_{i}$,
$(\widehat{\bm{\theta}}-\bm{\theta}^{\star})_{i}$, and
$(\nabla_{\bm{\theta}}\,\ell(\bm{\theta};\bm{x},\bm{y}))_{i}$ the
$i^{\text{th}}$ component of $\bm{R}$,
$(\widehat{\bm{\theta}}-\bm{\theta}^{\star})$, and the score function
$\nabla_{\bm{\theta}}\,\ell(\bm{\theta};\bm{x},\bm{y})$, respectively. The
remainder term $R_{i}$ ($i=1,\ldots,p$) is given by

$\begin{array}[]{lllllllllllll}R_{i}&=&\displaystyle\sum\limits_{j=1}^{p}\,\dfrac{1}{2}\,\dfrac{\partial^{2}\,(\nabla_{\bm{\theta}}\,\ell(\dot{\bm{\theta}}_{i};\bm{x},\bm{y}))_{i}}{\partial\,\theta_{j}^{2}}(\widehat{\bm{\theta}}-\bm{\theta}^{\star})_{j}^{2}+\;\displaystyle\sum\limits_{1\leq
j<k\leq
p}\,\dfrac{\partial^{2}\,(\nabla_{\bm{\theta}}\ell(\dot{\bm{\theta}}_{i};\bm{x},\bm{y}))_{i}}{\partial\,\theta_{j}\,\partial\,\theta_{k}}\,(\widehat{\bm{\theta}}-\bm{\theta}^{\star})_{j}\,(\widehat{\bm{\theta}}-\bm{\theta}^{\star})_{k},\end{array}$

where
$\dot{\bm{\theta}}_{i}=t_{i}\,\widehat{\bm{\theta}}+(1-t_{i})\,\bm{\theta}^{\star}$
(for some $t_{i}\in[0,1]$). By Proposition 2.1,

$\begin{array}[]{lllllllllllll}\ell(\bm{\theta};\bm{x},\bm{y})&=&\log\,\mathbb{P}_{\bm{\theta}}(\bm{X}=\bm{x}\,|\,\bm{Y}=\bm{y})+\log\,g(\bm{y}).\end{array}$

By Lemma 4, the probability mass function
$\mathbb{P}_{\bm{\theta}}(\bm{X}=\bm{x}\,|\,\bm{Y}=\bm{y})$ belongs to a
minimal exponential family with sufficient statistic vector
$s(\bm{x})=(s_{1}(\bm{x}),\ldots,s_{p}(\bm{x}))$ given by equation (11) in
Lemma 4. We then have,

$\begin{array}[]{lllllllllllll}&-\nabla_{\bm{\theta}}\,\ell(\bm{\theta};\bm{x},\bm{y})&=&-(s(\bm{x})-\mathbb{E}_{\bm{\theta}}^{\bm{y}}\,s(\bm{X}))\vspace{0.25cm}\\\
&-\nabla_{\bm{\theta}}^{2}\,\ell(\bm{\theta};\bm{x},\bm{y})&=&\mathbb{V}_{\bm{\theta}}^{\bm{y}}\,s(\bm{X})\;=\;I(\bm{\theta}^{\star})\,|\\!|\bm{y}|\\!|_{1},\end{array}$

where $\mathbb{E}_{\bm{\theta}}^{\bm{y}}$ and
$\mathbb{V}_{\bm{\theta}}^{\bm{y}}$ are the conditional expectation and
variance operators, respectively, of the conditional distribution of $\bm{X}$
given $\bm{Y}=\bm{y}$, and by using standard formulas for exponential families
[e.g., Proposition 3.8, pp. 29, Sundberg, 2019] and the results of Lemma 1.
Note $\nabla_{\bm{\theta}}\,\ell(\widehat{\bm{\theta}};\bm{x},\bm{y})=0$, as
the maximum likelihood estimator $\widehat{\bm{\theta}}$ solves the score
equation by definition. We re-arrange (15) and multiply both sides by
$(I(\bm{\theta}^{\star})\,|\\!|\bm{Y}|\\!|_{1})^{-1/2}$ to obtain

$\begin{array}[]{lllllllllllllllll}&&(I(\bm{\theta}^{\star})\,|\\!|\bm{Y}|\\!|_{1})^{1/2}\,(\widehat{\bm{\theta}}-\bm{\theta}^{\star})-\;(I(\bm{\theta}^{\star})\,|\\!|\bm{Y}|\\!|_{1})^{-1/2}\,\bm{R}\vspace{0.25cm}\\\
&=&(I(\bm{\theta}^{\star})\,|\\!|\bm{Y}|\\!|_{1})^{-1/2}\,(s(\bm{X})-\mathbb{E}^{\bm{Y}}\,s(\bm{X})).\end{array}$
(16)

Define
$\tilde{\bm{R}}\coloneqq(I(\bm{\theta}^{\star})\,|\\!|\bm{Y}|\\!|_{1})^{-1/2}\,\bm{R}$.
Let $\mathscr{A}\subset\mathbb{R}^{p}$ be any measurable convex subset of
$\mathbb{R}^{p}$ and $\bm{Z}\sim\text{MvtNorm}(\bm{0}_{p},\,\bm{I}_{p})$. We
are interested in bounding the quantity

$\begin{array}[]{lllllllllllll}\left|\mathbb{P}((I(\bm{\theta}^{\star})\,|\\!|\bm{Y}|\\!|_{1})^{1/2}\,(\widehat{\bm{\theta}}-\bm{\theta}^{\star})-\tilde{\bm{R}}\in\mathscr{A})-\Phi(\bm{Z}\in\mathscr{A})\right|.\end{array}$

Then from (LABEL:eq:bridge),

$\begin{array}[]{lllllllllllll}\mathbb{P}\left((I(\bm{\theta}^{\star})\,|\\!|\bm{Y}|\\!|_{1})^{1/2}\,(\widehat{\bm{\theta}}-\bm{\theta}^{\star})-\tilde{\bm{R}}\in\mathscr{A}\right)\,=\,\mathbb{P}\big{(}(I(\bm{\theta}^{\star})\,|\\!|\bm{Y}|\\!|_{1})^{-1/2}\,(s(\bm{X})-\mathbb{E}^{\bm{Y}}\,s(\bm{X}))\in\mathscr{A}\big{)}.\end{array}$

Applying Proposition E, for all measurable convex sets
$\mathscr{A}\subseteq\mathbb{R}^{p}$,

$\begin{array}[]{lllllllllllll}&&\left|\mathbb{P}\left((I(\bm{\theta}^{\star})\,|\\!|\bm{Y}|\\!|_{1})^{-1/2}\,(s(\bm{X})-\mathbb{E}^{\bm{Y}}\,s(\bm{X}))\in\mathscr{A}\right)-\Phi(\bm{Z}\in\mathscr{A})\right|\vspace{0.25cm}\\\
&\leq&\dfrac{83}{\xi_{\epsilon^{\star}}^{3/2}}\,\sqrt{\dfrac{p^{7/2}}{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}}+\dfrac{4}{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}+\dfrac{8\,\left[D_{g}\right]^{+}}{\left(\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}\right)^{2}}.\end{array}$

Hence,

$\begin{array}[]{lllllllllllll}&&\left|\mathbb{P}((I(\bm{\theta}^{\star})\,|\\!|\bm{Y}|\\!|_{1})^{1/2}\,(\widehat{\bm{\theta}}-\bm{\theta}^{\star})-\tilde{\bm{R}}\in\mathscr{A})-\Phi(\bm{Z}\in\mathscr{A})\right|\vspace{0.25cm}\\\
&\leq&\dfrac{83}{\xi_{\epsilon^{\star}}^{3/2}}\,\sqrt{\dfrac{p^{7/2}}{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}}+\dfrac{4}{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}+\dfrac{8\,\left[D_{g}\right]^{+}}{\left(\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}\right)^{2}}.\vspace{0.25cm}\end{array}$

We lastly handle the term $\tilde{\bm{R}}$ by showing that
$|\\!|\tilde{\bm{R}}|\\!|_{2}$ is small with high probability. We first use
standard vector/matrix norm inequalities to bound

$\begin{array}[]{lllllllllllll}|\\!|\widetilde{\bm{R}}|\\!|_{2}&=&|\\!|(I(\bm{\theta}^{\star})\,|\\!|\bm{Y}|\\!|_{1})^{-1/2}\,\bm{R}|\\!|_{2}&\leq&\dfrac{|\\!|\\!|I(\bm{\theta}^{\star})^{-1/2}|\\!|\\!|_{2}}{\sqrt{|\\!|\bm{Y}|\\!|_{1}}}\;|\\!|\bm{R}|\\!|_{2}&\leq&\dfrac{|\\!|\bm{R}|\\!|_{2}}{\sqrt{\xi_{\epsilon^{\star}}|\\!|\bm{Y}|\\!|_{1}}},\end{array}$

noting that the spectral norm
$|\\!|\\!|I(\bm{\theta}^{\star})^{-1/2}|\\!|\\!|_{2}$ is equal to the largest
eigenvalue of $I(\bm{\theta}^{\star})^{-1/2}$ which will be the reciprocal of
the smallest eigenvalue of $I(\bm{\theta}^{\star})^{1/2}$, which is bounded
below by $\sqrt{\xi_{\epsilon^{\star}}}$. Using a standard result from the
Taylor theorem for functions with multiple variables, if for each
$i=1,\ldots,p$, there exists constants $M_{i}>0$ such that

$\begin{array}[]{lllllllllllll}\sup\limits_{\bm{\theta}\in\mathbb{R}^{p}\,:\,|\\!|\bm{\theta}-\bm{\theta}^{\star}|\\!|_{1}\leq|\\!|\widehat{\bm{\theta}}-\bm{\theta}^{\star}|\\!|_{1}}\,\left|\,\dfrac{\partial^{2}(\nabla_{\bm{\theta}}\,\ell(\bm{\theta};\bm{x},\bm{y}))_{i}}{\partial\,\theta_{j}\,\partial\,\theta_{k}}\,\right|&\leq&M_{i},&&1\leq
j\leq k\leq p,\end{array}$

then the Lagrange remainder is bounded above by

$\begin{array}[]{lllllllllllll}R_{i}&\leq&\dfrac{M_{i}}{2}\,|\\!|\widehat{\bm{\theta}}-\bm{\theta}^{\star}|\\!|_{1}^{2}\end{array}$

on the set
$\\{\bm{\theta}\in\mathbb{R}^{p}:|\\!|\bm{\theta}-\bm{\theta}^{\star}|\\!|_{1}\leq|\\!|\widehat{\bm{\theta}}-\bm{\theta}^{\star}|\\!|_{2}\\}$.
By Lemma 7, conditional on $\bm{Y}=\bm{y}$, we have, for all $i=1,\ldots,p$,
the bound $M_{i}\leq 2\,|\\!|\bm{y}|\\!|_{1}$. Hence,

$\begin{array}[]{lllllllllllllllll}|\\!|\tilde{\bm{R}}|\\!|_{2}&\leq&\dfrac{1}{\sqrt{\xi_{\epsilon^{\star}}\,|\\!|\bm{y}|\\!|_{1}}}\;\sqrt{\displaystyle\sum\limits_{i=1}^{p}\,R_{i}^{2}}&\leq&\dfrac{1}{\sqrt{\xi_{\epsilon^{\star}}\,|\\!|\bm{y}|\\!|_{1}}}\,\sqrt{\displaystyle\sum\limits_{i=1}^{p}\,|\\!|\bm{y}|\\!|_{1}^{2}\,|\\!|\widehat{\bm{\theta}}-\bm{\theta}^{\star}|\\!|_{1}^{4}}\vspace{0.25cm}\vspace{0.25cm}\\\
&\leq&\dfrac{1}{\sqrt{\xi_{\epsilon^{\star}}\,|\\!|\bm{y}|\\!|_{1}}}\,\sqrt{p\,|\\!|\bm{y}|\\!|_{1}^{2}\,|\\!|\widehat{\bm{\theta}}-\bm{\theta}^{\star}|\\!|_{1}^{4}}&\leq&\dfrac{\sqrt{p}\,|\\!|\bm{y}|\\!|_{1}\,|\\!|\widehat{\bm{\theta}}-\bm{\theta}^{\star}|\\!|_{1}^{2}}{\sqrt{\xi_{\epsilon^{\star}}\,|\\!|\bm{y}|\\!|_{1}}}\vspace{0.25cm}\vspace{0.25cm}\\\
&\leq&\dfrac{\sqrt{p}\,\sqrt{|\\!|\bm{y}|\\!|_{1}}\,p\,|\\!|\widehat{\bm{\theta}}-\bm{\theta}^{\star}|\\!|_{2}^{2}}{\sqrt{\xi_{\epsilon^{\star}}}}&\leq&\dfrac{p^{3/2}\,\sqrt{|\\!|\bm{y}|\\!|_{1}}\,|\\!|\widehat{\bm{\theta}}-\bm{\theta}^{\star}|\\!|_{2}^{2}}{\sqrt{\xi_{\epsilon^{\star}}}}.\end{array}$
(17)

By Chebyshev’s inequality—as in the proof of Lemma 2—we can establish that

$\begin{array}[]{lllllllllllllllll}\mathbb{P}\left(\left||\\!|\bm{Y}|\\!|_{1}-\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}\right|>\dfrac{1}{2}\,\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}\right)&\leq&\dfrac{4}{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}+\dfrac{8\,[D_{g}]^{+}}{(\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1})^{2}},\end{array}$
(18)

whereas Theorem 1 established there exists $N_{0}\geq 3$ such that, for all
$N\geq N_{0}$, the event

$\begin{array}[]{lllllllllllllllll}|\\!|\widehat{\bm{\theta}}-\bm{\theta}^{\star}|\\!|_{2}&\leq&\sqrt{\dfrac{3\,p\,\log
N}{\mathbb{E}|\\!|\bm{Y}|\\!|_{1}}}\;\dfrac{\sqrt{1+[D_{g}]^{+}}}{\xi_{\epsilon^{\star}}},\end{array}$
(19)

occurs with probability at least
$1-3\,(\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1})^{-1}$. Define $\mathscr{E}_{1}$ to be
the event

$\begin{array}[]{lllllllllllll}||\\!|\bm{Y}|\\!|_{1}-\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}|&\leq&\dfrac{1}{2}\,\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}\end{array}$

and $\mathscr{E}_{2}$ to be the event in (19), and define $\mathcal{R}$ to be
the corresponding values of $\tilde{\bm{R}}$ in the event
$(\bm{X},\bm{Y})\in\mathscr{E}_{1}\cap\mathscr{E}_{2}$, under which we have
the bound

$\begin{array}[]{lllllllllllllllll}|\\!|\tilde{\bm{R}}|\\!|_{2}&\leq&\dfrac{p^{3/2}\,\sqrt{|\\!|\bm{y}|\\!|_{1}}}{\sqrt{\xi_{\epsilon^{\star}}}}\,\dfrac{3\,p\,\log
N}{\mathbb{E}|\\!|\bm{Y}|\\!|_{1}}\;\dfrac{1+[D_{g}]^{+}}{\xi_{\epsilon^{\star}}^{2}}\vspace{0.25cm}\vspace{0.25cm}\\\
&=&\dfrac{3\,p^{5/2}\,\sqrt{2\,\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}\,\log\,N}{\mathbb{E}|\\!|\bm{Y}|\\!|_{1}}\,\dfrac{1+[D_{g}]^{+}}{\xi_{\epsilon^{\star}}^{2}}\vspace{0.25cm}\vspace{0.25cm}\\\
&=&\dfrac{3\,\sqrt{2}\,(1+[D_{g}]^{+})}{\xi_{\epsilon^{\star}}^{2}}\;\dfrac{p^{5/2}\,\log
N}{\sqrt{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}}.\end{array}$ (20)

combining the bounds in (17), (18), and (19) and using the fact that

$\begin{array}[]{lllllllllllll}|\\!|\bm{y}|\\!|_{1}&\leq&\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}+\dfrac{1}{2}\,\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}&\leq&2\,\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}\end{array}$

in the event $\bm{y}\in\mathscr{E}_{1}$. Moreover, a union bound shows that

$\begin{array}[]{lllllllllllll}\mathbb{P}(\tilde{\bm{R}}\not\in\mathcal{R})&\leq&\dfrac{7}{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}+\dfrac{8\,[D_{g}]^{+}}{(\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1})^{2}}.\end{array}$

Hence,

$\begin{array}[]{lllllllllllllllll}\mathbb{P}\left(|\\!|\tilde{\bm{R}}|\\!|_{2}\leq\dfrac{3\,\sqrt{2}\,(1+[D_{g}]^{+})}{\xi_{\epsilon^{\star}}^{2}}\;\dfrac{p^{5/2}\,\log
N}{\sqrt{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}}\right)&\geq&1-\dfrac{7}{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}-\dfrac{8\,[D_{g}]^{+}}{(\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1})^{2}}.\end{array}$
(21)

Taken together, we have shown under the assumptions of Theorem 1 that there
exists $N_{0}\geq 3$ such that, for all $N\geq N_{0}$, the error of the
multivariate normal approximation

$\begin{array}[]{lllllllllllll}\left|\mathbb{P}((I(\bm{\theta}^{\star})\,|\\!|\bm{Y}|\\!|_{1})^{1/2}\,(\widehat{\bm{\theta}}-\bm{\theta}^{\star})-\tilde{\bm{R}}\in\mathscr{A})-\Phi(\bm{Z}\in\mathscr{A})\right|\end{array}$

is bounded above by

$\begin{array}[]{lllllllllllll}\dfrac{83}{\xi_{\epsilon^{\star}}^{3/2}}\,\sqrt{\dfrac{p^{7/2}}{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}}+\dfrac{4}{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}+\dfrac{8\,\left[D_{g}\right]^{+}}{\left(\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}\right)^{2}}\end{array}$

where $\tilde{\bm{R}}$ satisfies

$\begin{array}[]{lllllllllllll}\mathbb{P}\left(|\\!|\tilde{\bm{R}}|\\!|_{2}\leq\dfrac{3\,\sqrt{2}\,(1+[D_{g}]^{+})}{\xi_{\epsilon^{\star}}^{2}}\;\dfrac{p^{5/2}\,\log
N}{\sqrt{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}}\right)&\geq&1-\dfrac{7}{\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1}}-\dfrac{8\,[D_{g}]^{+}}{(\mathbb{E}\,|\\!|\bm{Y}|\\!|_{1})^{2}}.\end{array}$

$\hfill\blacksquare$

###  F.1 Auxiliary results for proof of Theorem 2

###### Lemma 7

Consider multilayer networks satisfying (1) defined on a set of $N\geq 3$
nodes and $K\geq 1$ layers and denote by the log-likelihood function by
$\ell(\bm{\theta};\bm{x},\bm{y})$. Then, for each $i=1,\ldots,p$ and
$\epsilon>0$,

$\begin{array}[]{lllllllllllll}\sup\limits_{\bm{\theta}\in\mathbb{R}^{p}\,:\,|\\!|\bm{\theta}-\bm{\theta}^{\star}|\\!|_{2}\leq\epsilon}\quad\left|\dfrac{\partial^{2}\,(\nabla_{\bm{\theta}}\,\ell(\bm{\theta};\bm{x},\bm{y}))_{i}}{\partial\,\theta_{j}\,\partial\,\theta_{k}}\right|&\leq&2\,|\\!|\bm{y}|\\!|_{1},\end{array}$

where $(\nabla_{\bm{\theta}}\,\ell(\bm{\theta};\bm{x},\bm{y}))_{i}$ is the
$i^{th}$ component of the score function
$\nabla_{\bm{\theta}}\,\ell(\bm{\theta};\bm{x},\bm{y})$.

Proof of Lemma 7. By Proposition 2.1, given the observation $\bm{x}$ of
$\bm{X}$ (i.e., observation of the event $\bm{X}=\bm{x}$), $\bm{Y}$ is
predictable with unique value $\bm{y}\in\mathbb{Y}$ given by the formula in
Proposition 2.1, and $(\bm{x},\bm{y})$ is network concordant. Further, by
Proposition 2.1

$\begin{array}[]{lllllllllllll}\ell(\bm{\theta};\bm{x},\bm{y})&=&\log\,\mathbb{P}_{\bm{\theta}}(\bm{X}=\bm{x}\mid\bm{Y}=\bm{y})+\log\,g(\bm{y}),\end{array}$

where $\log\,\mathbb{P}_{\bm{\theta}}(\bm{X}=\bm{x}|\bm{Y}=\bm{y})$ is the
log-likelihood of a minimal, full, and regular exponential family. Thus, the
second order derivative of $\ell(\bm{\theta};\bm{x},\bm{y})$ with respect to
the $i^{\text{th}}$ and $j^{\text{th}}$ components of $\bm{\theta}$ correspond
to the variance (in the case $i=j$) or covariance (in the case of $i\neq j$)
of corresponding sufficient statistic(s) of the exponential family [e.g.,
Proposition 3.8, p. 29, Sundberg, 2019], with sufficient statistics given in
Lemma 4. For $\\{i,j\\}\subseteq\\{1,\ldots,p\\}$,

$\begin{array}[]{lllllllllllll}\dfrac{\partial\,(\nabla_{\bm{\theta}}\,\ell(\bm{\theta};\bm{x},\bm{y}))_{i}}{\partial\,\theta_{j}}&=&\dfrac{\partial^{2}\,\ell(\bm{\theta};\bm{x},\bm{y})}{\partial\,\theta_{i}\;\partial\,\theta_{j}}&=&\mathbb{C}_{\bm{\theta}}(s_{i}(\bm{X}),s_{j}(\bm{X})\,|\,\bm{Y}=\bm{y}),\end{array}$

and when $i=j\in\\{1,\ldots,p\\}$,

$\begin{array}[]{lllllllllllll}\dfrac{\partial\,(\nabla_{\bm{\theta}}\,\ell(\bm{\theta};\bm{x},\bm{y}))_{i}}{\partial\,\theta_{i}}&=&\dfrac{\partial^{2}\,\ell(\bm{\theta};\bm{x},\bm{y})}{\partial\,\theta_{i}^{2}}&=&\mathbb{V}_{\bm{\theta}}(s_{i}(\bm{X})\,|\,\bm{Y}=\bm{y}).\end{array}$

As a result, for $\\{i,j\\}\subseteq\\{1,\ldots,p\\}$ and
$k\in\\{1,\ldots,p\\}$,

$\begin{array}[]{lllllllllllll}\left|\,\dfrac{\partial^{2}\,(\nabla_{\bm{\theta}}\,\ell(\bm{\theta};\bm{x},\bm{y}))_{i}}{\partial\,\theta_{j}\,\partial\,\theta_{k}}\,\right|&=&\left|\,\dfrac{\partial\,\mathbb{C}_{\bm{\theta}}(s_{i}(\bm{X}),s_{j}(\bm{X})\,|\,\bm{Y}=\bm{y})}{\partial\,\theta_{k}}\,\right|,\end{array}$

and when $i=j\in\\{1,\ldots,p\\}$ and $k\in\\{1,\ldots,p\\}$,

$\begin{array}[]{lllllllllllll}\left|\,\dfrac{\partial^{2}\,(\nabla_{\bm{\theta}}\,\ell(\bm{\theta};\bm{x},\bm{y}))_{i}}{\partial\,\theta_{i}\,\partial\,\theta_{k}}\,\right|&=&\left|\,\dfrac{\partial\,\mathbb{V}_{\bm{\theta}}(s_{i}(\bm{X})\,|\,\bm{Y}=\bm{y})}{\partial\,\theta_{k}}\,\right|.\end{array}$

By Lemma 4 equation (11), conditional on $\bm{Y}=\bm{y}$, each sufficient
statistic $s_{i}(\bm{X})$ ($i\in\\{1,\ldots,p\\}$) can be decomposed into the
sum of conditionally independent statistics of each dyad $\bm{X}_{v,w}$, for
$\\{v,w\\}\subseteq\mathscr{N}$. We can then write

$\begin{array}[]{lllllllllllll}\mathbb{C}_{\bm{\theta}}(s_{i}(\bm{X}),s_{j}(\bm{X})\,|\,\bm{Y}=\bm{y})&=&\displaystyle\sum\limits_{\\{v,w\\}\subset\mathscr{N}}\,\mathbb{C}_{\bm{\theta}}(s_{i,v,w}(\bm{X}_{v,w}),\,s_{j,v,w}(\bm{X}_{v,w})\,|\,\bm{Y}=\bm{y}),\end{array}$

noting that by conditional independence
$\mathbb{C}_{\bm{\theta}}(s_{i,v,w}(\bm{X}_{v,w}),s_{j,r,t}(\bm{X}_{r,t})\,|\,\bm{Y}=\bm{y})=0$
whenever $\\{r,t\\}\neq\\{v,w\\}$, and when $i=j$, we can write

$\begin{array}[]{lllllllllllll}\mathbb{V}_{\bm{\theta}}(s_{i}(\bm{X})\,|\,\bm{Y}=\bm{y})&=&\displaystyle\sum\limits_{\\{v,w\\}\subset\mathscr{N}}\,\mathbb{V}_{\bm{\theta}}(s_{i,v,w}(\bm{X}_{v,w})\,|\,\bm{Y}=\bm{y}),\end{array}$

again appealing to the conditional independence given $\bm{Y}$ of the random
variables $s_{i,v,w}(\bm{X}_{v,w})$ ($\\{v,w\\}\subset\mathscr{N}$). As a
result, for $k\in\\{1,\ldots,p\\}$, it suffices to show that,

$\begin{array}[]{lllllllllllll}\left|\,\dfrac{\partial\,\mathbb{C}_{\bm{\theta}}(s_{i,v,w}(\bm{X}_{v,w}),\,s_{j,v,w}(\bm{X}_{v,w})\,|\,\bm{Y}=\bm{y})}{\partial\,\theta_{k}}\,\right|&\leq&2,\end{array}$

and

$\begin{array}[]{lllllllllllll}\left|\,\dfrac{\partial\,\mathbb{V}_{\bm{\theta}}(s_{i,v,w}(\bm{X}_{v,w})\,|\,\bm{Y}=\bm{y})}{\partial\,\theta_{k}}\,\right|&\leq&1.\end{array}$

Recall that the sufficient statistic $s_{i,v,w}(\bm{X})$ ($i=1,\ldots,p$) is
defined in Lemma 4 by

$\begin{array}[]{lllllllllllll}s_{i,v,w}(\bm{X}_{v,w})&=&\displaystyle\prod\limits_{t=1}^{h}\,X_{v,w}^{(k_{t})},&&\\{v,w\\}\subset\mathscr{N},\end{array}$

for some $h\in\\{1,\ldots,H\\}$ and
$\\{k_{1},\ldots,k_{h}\\}\subseteq\\{1,\ldots,K\\}$. Define the set
$S_{i,v,w}$ of components of the sufficient statistic vector
$\bm{s}_{v,w}(\bm{X})$ for $\\{v,w\\}\subset\mathscr{N}$ and $i=1,\ldots,p$ by

$\begin{array}[]{lllllllllllll}S_{i,v,w}\;\coloneqq\;\left\\{\displaystyle\prod\limits_{t=1}^{h^{\prime}}\,X_{v,w}^{(l_{t})}\;:\;\\{l_{1},\ldots,l_{h^{\prime}}\\}\subset\\{k_{1},\ldots,k_{h}\\},\;h^{\prime}<h\right\\},\end{array}$

where $h\in\\{1,\ldots,H\\}$ and
$\\{k_{1},\ldots,k_{h}\\}\subseteq\\{1,\ldots,K\\}$. The set $S_{i,v,w}$ is
the set of components of the sufficient statistic vector
$\bm{s}_{v,w}(\bm{X})$ of dyad $\\{v,w\\}\subset\mathscr{N}$ that have a value
of 1 when $s_{i,v,w}(\bm{X})=1$. For the ease of notation, let $I_{S_{i,v,w}}$
be the set of dimension indices whose corresponding components of the
sufficient statistic vector $\bm{s}_{v,w}(\bm{X})$ belong to the set
$S_{i,v,w}$:

$\begin{array}[]{lllllllllllll}I_{S_{i,v,w}}\;\coloneqq\;\\{j\in\\{1,\ldots,p\\}\;:\;s_{j,v,w}(\bm{X})\in
S_{i,v,w}\\}.\end{array}$

Define the conditional expectation of $s_{i,v,w}(\bm{X})$ given
$\bm{Y}=\bm{y}$ for any $i\in\\{1,\ldots,p\\}$ and
$\\{v,w\\}\subset\mathscr{N}$ by

$\begin{array}[]{lllllllllllll}P_{i,v,w}(\bm{\theta};\bm{X},\bm{y})\;\coloneqq\;\mathbb{P}_{\bm{\theta}}\,(s_{i,v,w}(\bm{X})=1\,|\,\bm{Y}=\bm{y}).\end{array}$

For further notation simplicity, denote by $L_{i}$ the set of layer indices
$\\{k_{1},\ldots,k_{h}\\}\subseteq\\{1,\ldots,K\\}$ that define the $i^{th}$
component $s_{i,v,w}(\bm{X}_{v,w})$ of the sufficient statistic vector
$\bm{s}_{v,w}(\bm{X})$ for any $\\{v,w\\}\subset\mathscr{N}$,
$j\in\\{1,\ldots,p\\}$, and some $h\in\\{1,\ldots,H\\}$. We then define

$\begin{array}[]{lllllllllllll}\bm{X}_{v,w}^{(L_{i})}&\coloneqq&\left\\{X_{v,w}^{(k_{1})},\ldots,X_{v,w}^{(k_{h})}\right\\},&&\bm{X}_{v,w}^{(-L_{i})}\;\coloneqq\;\bm{X}_{v,w}\setminus\bm{X}_{v,w}^{(L_{i})},\end{array}$

and the corresponding sample space

$\begin{array}[]{lllllllllllll}\mathbb{X}_{v,w}^{(L_{i})}&\coloneqq&\\{0,1\\}^{h},&&\mathbb{X}_{v,w}^{(-L_{i})}\;\coloneqq\;\\{0,1\\}^{H-h},\end{array}$

for some $h\in\\{1,\ldots,H\\}$. Then we can write

$\begin{array}[]{lllllllllllll}P_{i,v,w}(\bm{\theta};\bm{X},\bm{y})&=&\mathbb{P}_{\bm{\theta}}\left(\displaystyle\prod\limits_{l\in
L_{i}}\,X_{v,w}^{(l)}=1\,|\,\bm{Y}=\bm{y}\right)\vspace{0.25cm}\\\
&=&\dfrac{\displaystyle\sum\limits_{\mathbb{X}_{v,w}^{(-L_{i})}}\,\exp\,\left(\sum_{j\in
I_{S_{i,v,w}}}\,\theta_{j}+\sum_{j\in
I_{S_{i,v,w}}^{c}}\,\theta_{j}\,s_{j,v,w}(\bm{x})\right)}{\displaystyle\sum\limits_{\mathbb{X}_{v,w}}\,\exp\left(\sum_{j=1}^{p}\,\theta_{j}\,s_{j,v,w}(\bm{x})\right)}.\end{array}$

Let

$\begin{array}[]{lllllllllllll}Z(\bm{\theta})&\coloneqq&\displaystyle\sum\limits_{\mathbb{X}_{v,w}}\,\exp\left(\sum_{j=1}^{p}\,\theta_{j}\,s_{j,v,w}(\bm{x})\right),\end{array}$

and take the derivative of $P_{i,v,w}(\bm{\theta};\bm{x},\bm{y})$ with respect
to $\theta_{k}$ for $k=1,\ldots,p$. We have

$\begin{array}[]{lllllllllllllllll}&\dfrac{\partial\,P_{i,v,w}(\bm{\theta};\bm{X},\bm{y})}{\partial\,\theta_{k}}\vspace{0.25cm}\\\
\leq&\dfrac{\displaystyle\sum\limits_{\mathbb{X}_{v,w}^{(-L_{i})}}\,\exp\,\left(\sum_{j\in
I_{S_{i,v,w}}}\,\theta_{j}+\sum_{j\in
I_{S_{i,v,w}}^{c}}\,\theta_{j}\,s_{j,v,w}(\bm{x})\right)\,\left(Z(\bm{\theta})-\dfrac{\partial\,Z(\bm{\theta})}{\partial\,\theta_{k}}\right)}{Z(\bm{\theta})^{2}}\vspace{0.25cm}\\\
=&\dfrac{\displaystyle\sum\limits_{\mathbb{X}_{v,w}^{(-L_{i})}}\,\exp\,\left(\sum_{j\in
I_{S_{i,v,w}}}\,\theta_{j}+\sum_{j\in
I_{S_{i,v,w}}^{c}}\,\theta_{j}\,s_{j,v,w}(\bm{x})\right)\,\left(\displaystyle\sum\limits_{\mathbb{X}_{v,w}}\,\exp\left(\sum_{j=1}^{p}\,\theta_{j}\,s_{j,v,w}(\bm{x})\right)\,(1-s_{k,v,w}(\bm{x}))\right)}{Z(\bm{\theta})^{2}}\vspace{0.25cm}\\\
\leq&\dfrac{\displaystyle\sum\limits_{\mathbb{X}_{v,w}^{(-L_{i})}}\,\exp\,\left(\sum_{j\in
I_{S_{i,v,w}}}\,\theta_{j}+\sum_{j\in
I_{S_{i,v,w}}^{c}}\,\theta_{j}\,s_{j,v,w}(\bm{x})\right)\,Z(\bm{\theta})}{Z(\bm{\theta})^{2}}\vspace{0.25cm}\\\
\leq&1.\end{array}$ (22)

The first inequality is obtained because $s_{k,v,w}(\bm{x})\leq 1$, and the
last inequality is due to the fact that

$\begin{array}[]{lllllllllllll}\displaystyle\sum\limits_{\mathbb{X}_{v,w}^{(-L_{i})}}\,\exp\,\left(\sum_{j\in
I_{S_{i,v,w}}}\,\theta_{j}+\sum_{j\in
I_{S_{i,v,w}}^{c}}\,\theta_{j}\,s_{j,v,w}(\bm{x})\right)&\leq&\displaystyle\sum\limits_{\mathbb{X}_{v,w}}\,\exp\left(\sum_{j=1}^{p}\,\theta_{j}\,s_{j,v,w}(\bm{x})\right).\end{array}$

Now we turn to show the derivative of the conditional variance and covariance
of the sufficient statistics of each dyad are bounded. Given $\bm{Y}=\bm{y}$,
for all $\\{i\\}\subset\\{1,\ldots,p\\}$, $s_{i,v,w}(\bm{X})$ are
conditionally independent across $\\{v,w\\}\subseteq\mathscr{N}$. Then we have

$\begin{array}[]{lllllllllllll}&\mathbb{C}_{\bm{\theta}}(s_{i,v,w}(\bm{X}_{v,w}),\,s_{j,v,w}(\bm{X}_{v,w})\,|\,\bm{Y}=\bm{y})\vspace{0.25cm}\\\
=&\mathbb{E}\,[s_{i,v,w}(\bm{X})\,s_{j,v,w}(\bm{X})\,|\,\bm{Y}=\bm{y}]-\mathbb{E}\,[s_{i,v,w}(\bm{X})\,|\,\bm{Y}=\bm{y}]\,\mathbb{E}\,[s_{j,v,w}(\bm{X})\,|\,\bm{Y}=\bm{y}]\vspace{0.25cm}\\\
=&\mathbb{P}_{\bm{\theta}}\,(s_{i,v,w}(\bm{X})=1,s_{j,v,w}(\bm{X})=1\,|\,\bm{Y}=\bm{y})-P_{i,v,w}(\bm{\theta};\bm{X},\bm{y})\,P_{j,v,w}(\bm{\theta};\bm{X},\bm{y})\vspace{0.25cm}\\\
=&\mathbb{P}_{\bm{\theta}}\left(\displaystyle\prod\limits_{l\in L_{i}\cup
L_{j}}\bm{X}_{v,w}^{(l)}=1\,|\,\bm{Y}=\bm{y}\right)-P_{i,v,w}(\bm{\theta};\bm{X},\bm{y})\,P_{j,v,w}(\bm{\theta};\bm{X},\bm{y}).\end{array}$

Using the inequality derived in (LABEL:cov_ineq) and suppressing the notation
of $\\{v,w\\}$ and $(\bm{X},\bm{y})$ in
$P_{i,v,w}(\bm{\theta};\bm{X},\bm{y})$, the derivative of the covariance with
respect to $\theta_{k}$, $k=1,\ldots,p$ is given by

$\begin{array}[]{lllllllllllll}&\left|\,\dfrac{\partial\,\mathbb{C}_{\bm{\theta}}(s_{i,v,w}(\bm{X}_{v,w}),\,s_{j,v,w}(\bm{X}_{v,w})\,|\,\bm{Y}=\bm{y})}{\partial\,\theta_{k}}\,\right|\vspace{0.25cm}\\\
=&\dfrac{\partial\,\mathbb{P}_{\bm{\theta}}\left(\displaystyle\prod\limits_{l\in
L_{i}\cup
L_{j}}\bm{X}_{v,w}^{(l)}=1\,|\,\bm{Y}=\bm{y}\right)}{\partial\,\theta_{k}}-\dfrac{\partial\,P_{i}(\bm{\theta})}{\partial\,\theta_{k}}\,P_{j}(\bm{\theta})-P_{i}(\bm{\theta})\,\dfrac{\partial\,P_{j}(\bm{\theta})}{\partial\,\theta_{k}}\vspace{0.25cm}\\\
\leq&2.\end{array}$

Using the same inequality and notation in (LABEL:cov_ineq), the derivative of
the variance of a Bernoulli random variable $s_{i,v,w}(\bm{X})$ is given by

$\begin{array}[]{lllllllllllll}\left|\,\dfrac{\partial\,\mathbb{V}_{\bm{\theta}}(s_{i,v,w}(\bm{X}_{v,w})\,|\,\bm{Y}=\bm{y})}{\partial\,\theta_{k}}\,\right|&=&\left|\,\left(1-2\,P_{i}(\bm{\theta})\right)\,\dfrac{\partial\,P_{i}(\bm{\theta})}{\partial\,\theta_{k}}\,\right|&\leq&1.\end{array}$

Finally, for $\\{i,j\\}\subseteq\\{1,\ldots,p\\}$ and $k\in\\{1,\ldots,p\\}$,
we obtain

$\begin{array}[]{lllllllllllll}\left|\,\dfrac{\partial^{2}\,(\nabla_{\bm{\theta}}\,\ell(\bm{\theta};\bm{x},\bm{y}))_{i}}{\partial\,\theta_{j}\,\partial\,\theta_{k}}\,\right|&=&\left|\,\dfrac{\partial\,\mathbb{C}_{\bm{\theta}}(s_{i}(\bm{X}),s_{j}(\bm{X})\,|\,\bm{Y}=\bm{y})}{\partial\,\theta_{k}}\,\right|\vspace{0.25cm}\\\
&\leq&\displaystyle\sum\limits_{\\{v,w\\}\subset\mathscr{N}}\,\left|\dfrac{\partial\,\mathbb{C}_{\bm{\theta}}(s_{i,v,w}(\bm{X}_{v,w}),\,s_{j,v,w}(\bm{X}_{v,w})\,|\,\bm{Y}=\bm{y})}{\partial\,\theta_{k}}\,\right|\vspace{0.25cm}\\\
&\leq&2\,|\\!|\bm{y}|\\!|_{1}\end{array}$

due to the fact that
$\mathbb{C}_{\bm{\theta}}(s_{i,v,w}(\bm{X}_{v,w}),\,s_{j,v,w}(\bm{X}_{v,w})\,|\,\bm{Y}=\bm{y})=0$
when $Y_{v,w}=0$ for $\\{v,w\\}\subset\mathscr{N}$. Similarly,
$\mathbb{V}_{\bm{\theta}}(s_{i,v,w}(\bm{X}_{i,v,w})\,|\,\bm{Y}=\bm{y})=0$ when
$Y_{v,w}=0$ for $\\{v,w\\}\subset\mathscr{N}$, and when
$i=j\in\\{1,\ldots,p\\}$ and $k\in\\{1,\ldots,p\\}$, we have

$\begin{array}[]{lllllllllllll}\left|\,\dfrac{\partial^{2}\,(\nabla_{\bm{\theta}}\,\ell(\bm{\theta};\bm{x},\bm{y}))_{i}}{\partial\,\theta_{i}\,\partial\,\theta_{k}}\,\right|&=&\left|\,\dfrac{\partial\,\mathbb{V}_{\bm{\theta}}(s_{i}(\bm{X})\,|\,\bm{Y}=\bm{y})}{\partial\,\theta_{k}}\,\right|\vspace{0.25cm}\\\
&\leq&|\\!|\bm{y}|\\!|_{1}.\end{array}$

$\hfill\blacksquare$

##  Appendix G Additional simulation results

We present additional simulation results that enhance those contained in
Section 5.

###  G.1 Normal approximation with different basis networks

The multivariate normality of $\widetilde{\bm{\theta}}$ is tested by Zhou-
Shao’s multivariate normal test [Zhou and Shao, 2014], and the p-values are
provided in tabel 4. Q-Q plots of $\widetilde{\bm{\theta}}$ estimated from 6
different model-generating parameters with a dense Bernoulli basis network, a
sparse Bernoulli basis network, a stochastic block model (SBM) generated basis
network, and a latent space model (LSM) generated basis network are shown in
Figure 6, 7, 8 and 9, respectively.

Table 4:  P-values of the Zhou-Shao’s test for multivariate normality of $\widetilde{\bm{\theta}}$ for 6 model-generating parameters ($\bm{\theta}^{\star}_{1}$, $\bm{\theta}^{\star}_{2}$, $\bm{\theta}^{\star}_{3}$, $\bm{\theta}^{\star}_{4}$, $\bm{\theta}^{\star}_{5}$, $\bm{\theta}^{\star}_{6}$) estimated from 250 network samples at size 1000 on four basis network structures. All p-values are larger than .1.  Basis network model | $\bm{\theta}^{\star}_{1}$ | $\bm{\theta}^{\star}_{2}$ | $\bm{\theta}^{\star}_{3}$ | $\bm{\theta}^{\star}_{4}$ | $\bm{\theta}^{\star}_{5}$ | $\bm{\theta}^{\star}_{6}$
---|---|---|---|---|---|---
Dense Bernoulli | .138 | .473 | .053 | .699 | .587 | .983
Sparse Bernoulli | .554 | .132 | .232 | .634 | .904 | .373
SBM | .65 | .891 | .982 | .975 | .871 | .674
LSM | .859 | .831 | .5 | .227 | .613 | .409

Figure 6: Q-Q plots and p-values of six components of
$\widetilde{\bm{\theta}}$ estimated from 250 multilayer network samples at
size 1000 on the dense Bernoulli basis network for 6 model-generating
parameters on each row.

Figure 7: Q-Q plots and p-values of six components of
$\widetilde{\bm{\theta}}$ estimated from 250 multilayer network samples at
size 1000 on the sparse Bernoulli basis network for 6 model-generating
parameters on each row.

Figure 8: Q-Q plots and p-values of six components of
$\widetilde{\bm{\theta}}$ estimated from 250 multilayer network samples at
size 1000 on the SBM generated basis network for 6 model-generating parameters
on each row.

Figure 9: Q-Q plots and p-values of six components of
$\widetilde{\bm{\theta}}$ estimated from 250 multilayer network samples at
size 1000 on the LSM generated basis network for 6 model-generating parameters
on each row.

###  G.2 Additional results on the false discovery rate

The false discovery rate (FDR) of the multiple testing correction procedures
of Bonferroni, Benjamini-Hochberg, Hochberg, and Holm to detect non-zero
components of $\bm{\theta}^{\star}$ at a family-wise significance level of
$\alpha=0.05$ with a sparse Bernoulli basis network, an SBM generated basis
network and an LSM generated basis network are provided in Table 5, 6 and 7,
respectively (recall that the third and the sixth component
$\theta_{1,3}^{\star}$ and $\theta_{3}^{\star}$ of $\bm{\theta}^{\star}$ are
set to 0).

Table 5:  False discovery rates of four procedures for detecting non-zero effects of six model-generating parameters ($\bm{\theta}^{\star}_{1}$, $\bm{\theta}^{\star}_{2}$, $\bm{\theta}^{\star}_{3}$, $\bm{\theta}^{\star}_{4}$, $\bm{\theta}^{\star}_{5}$, $\bm{\theta}^{\star}_{6}$) estimated from 250 multilayer network samples at size 1000 on the sparse Bernoulli basis network. All FDRs are smaller than 0.05. Procedure | $\bm{\theta}^{\star}_{1}$ | $\bm{\theta}^{\star}_{2}$ | $\bm{\theta}^{\star}_{3}$ | $\bm{\theta}^{\star}_{4}$ | $\bm{\theta}^{\star}_{5}$ | $\bm{\theta}^{\star}_{6}$
---|---|---|---|---|---|---
Bonferroni | .002 | .003 | .003 | .003 | .003 | .011
Benjamini-Hochberg | .020 | .011 | .022 | .022 | .014 | .017
Hochberg’s | .009 | .008 | .012 | .010 | .010 | .014
Holm’s | .007 | .008 | .011 | .009 | .006 | .014
Table 6:  False discovery rates of four procedures for detecting non-zero effects of six model-generating parameters ($\bm{\theta}^{\star}_{1}$, $\bm{\theta}^{\star}_{2}$, $\bm{\theta}^{\star}_{3}$, $\bm{\theta}^{\star}_{4}$, $\bm{\theta}^{\star}_{5}$, $\bm{\theta}^{\star}_{6}$) estimated from 250 multilayer network samples at size 1000 on the SBM generated basis network. All FDRs are smaller than 0.05. Procedure | $\bm{\theta}^{\star}_{1}$ | $\bm{\theta}^{\star}_{2}$ | $\bm{\theta}^{\star}_{3}$ | $\bm{\theta}^{\star}_{4}$ | $\bm{\theta}^{\star}_{5}$ | $\bm{\theta}^{\star}_{6}$
---|---|---|---|---|---|---
Bonferroni | .002 | .002 | .003 | .001 | .001 | .004
Benjamini-Hochberg | .022 | .013 | .014 | .015 | .015 | .018
Hochberg’s | .009 | .014 | .01 | .008 | .011 | .014
Holm’s | .009 | .013 | .005 | .009 | .009 | .011
Table 7:  False discovery rates of four procedures for detecting non-zero effects of six model-generating parameters ($\bm{\theta}^{\star}_{1}$, $\bm{\theta}^{\star}_{2}$, $\bm{\theta}^{\star}_{3}$, $\bm{\theta}^{\star}_{4}$, $\bm{\theta}^{\star}_{5}$, $\bm{\theta}^{\star}_{6}$) estimated from 250 multilayer network samples at size 1000 on the LSM generated basis network. All FDRs are smaller than 0.05. Procedure | $\bm{\theta}^{\star}_{1}$ | $\bm{\theta}^{\star}_{2}$ | $\bm{\theta}^{\star}_{3}$ | $\bm{\theta}^{\star}_{4}$ | $\bm{\theta}^{\star}_{5}$ | $\bm{\theta}^{\star}_{6}$
---|---|---|---|---|---|---
Bonferroni | .004 | .006 | .000 | .005 | .003 | .004
Benjamini-Hochberg | .016 | .013 | .011 | .015 | .016 | .017
Hochberg’s | .009 | .014 | .009 | .011 | .010 | .011
Holm’s | .008 | .014 | .009 | .011 | .007 | .010

## References

  * Airoldi et al. [2008] E. Airoldi, D. Blei, S. Fienberg, and E. Xing.  Mixed membership stochastic blockmodels.  _Journal of Machine Learning Research_ , 9:1981–2014, 2008\.
  * Albert and Barabási [2002] R. Albert and A. L. Barabási.  Statistical mechanics of complex networks.  _Reviews of Modern Physics_ , 74(1):47, 2002.
  * Arroyo et al. [2021] J. A. Arroyo, A. Athreya, J. Cape, G. Chen, C. E. Priebe, and J. T. Vogelstein.  Inference for multiple heterogeneous networks with a common invariant subspace.  _The Journal of Machine Learning Research_ , 22(1):6303–6351, 2021.
  * Athreya et al. [2018] A. Athreya, D. E. Fishkind, M. Tang, C. E. Priebe, J. T. Vogelstein, K. Levin, V. Lyzinski, Y. Qin, D. L. Sussman, E. Fishkind, and Y. Park.  Statistical inference on random dot product graphs: a survey.  _Journal of Machine Learning Research_ , 18:1–92, 2018\.
  * Besag [1974] J. Besag.  Spatial interaction and the statistical analysis of lattice systems.  _Journal of the Royal Statistical Society, Series B_ , 36:192–225, 1974.
  * Bhamidi et al. [2011] S. Bhamidi, G. Bresler, and A. Sly.  Mixing time of exponential random graphs.  _The Annals of Applied Probability_ , 21:2146–2170, 2011\.
  * Block [2015] P. Block.  Reciprocity, transitivity, and the mysterious three-cycle.  _Social Networks_ , 40:163–173, 2015.
  * Butts [2020] C. T. Butts.  A dynamic process interpretation of the sparse ERGM reference model.  _Journal of Mathematical Sociology_ , 2020.
  * Cai et al. [2016] D. Cai, T. Campbell, and T. Broderick.  Edge-exchangeable graphs and sparsity.  _Advances in Neural Information Processing Systems_ , 29, 2016.
  * Caimo and Gollini [2020] A. Caimo and I. Gollini.  A multilayer exponential random graph modelling approach for weighted networks.  _Computational Statistics & Data Analysis_, 142:106825, 2020.
  * Caron and Fox [2017] F. Caron and E. B. Fox.  Sparse graphs using exchangeable random measures.  _Journal of the Royal Statistical Society, Series B (with discussion)_ , 79:1–44, 2017.
  * Chen et al. [2022] S. Chen, S. Liu, and Z. Ma.  Global and individualized community detection in inhomogeneous multilayer networks.  _The Annals of Statistics_ , 50(5):2664–2693, 2022.
  * Crane and Dempsey [2018] H. Crane and W. Dempsey.  Edge exchangeable models for interaction networks.  _Journal of the American Statistical Association_ , 113(523):1311–1326, 2018.
  * Frank [1980] O. Frank.  Transitivity in stochastic graphs and digraphs.  _Journal of Mathematical Sociology_ , 7:199–213, 1980.
  * Furi and Martelli [1991] M. Furi and M. Martelli.  On the mean value theorem, inequality, and inclusion.  _The American Mathematical Monthly_ , 98(9):840–846, 1991.
  * Geyer and Thompson [1992] C. J. Geyer and E. A. Thompson.  Constrained Monte Carlo maximum likelihood for dependent data.  _Journal of the Royal Statistical Society, Series B_ , 54:657–699, 1992.
  * Goodreau et al. [2009] S. M. Goodreau, J. A. Kitts, and M. Morris.  Birds of a feather, or friend of a friend? using exponential random graph models to investigate adolescent social networks.  _Demography_ , 46(1):103–125, 2009.
  * Hoff et al. [2002] P. D. Hoff, A. E. Raftery, and M. S. Handcock.  Latent space approaches to social network analysis.  _Journal of the American Statistical Association_ , 97:1090–1098, 2002.
  * Holland and Leinhardt [1972] P. W. Holland and S. Leinhardt.  Some evidence on the transitivity of positive interpersonal sentiment.  _American Journal of Sociology_ , 77:1205–1209, 1972.
  * Holland and Leinhardt [1981] P. W. Holland and S. Leinhardt.  An exponential family of probability distributions for directed graphs.  _Journal of the American Statistical Association_ , 76:33–65, 1981.
  * Holland et al. [1983] P. W. Holland, K. B. Laskey, and S. Leinhardt.  Stochastic block models: some first steps.  _Social Networks_ , 5:109–137, 1983.
  * Huang et al. [2022] S. Huang, H. Weng, and Y. Feng.  Spectral clustering via adaptive layer aggregation for multi-layer networks.  _Journal of Computational and Graphical Statistics_ , pages 1–15, 2022.
  * Hunter and Handcock [2006] D. R. Hunter and M. S. Handcock.  Inference in curved exponential family models for networks.  _Journal of Computational and Graphical Statistics_ , 15:565–583, 2006.
  * Hunter et al. [2008] D. R. Hunter, S. M. Goodreau, and M. S. Handcock.  Goodness of fit of social network models.  _Journal of the American Statistical Association_ , 103:248–258, 2008.
  * Johnson and Wichern [2002] R.A. Johnson and D.W. Wichern.  _Applied multivariate statistical analysis_.  Prentice hall, 2002.
  * Krivitsky and Kolaczyk [2015] P. N. Krivitsky and E. D. Kolaczyk.  On the question of effective sample size in network modeling: An asymptotic inquiry.  _Statistical Science_ , 30:184–198, 2015.
  * Krivitsky et al. [2009] P. N. Krivitsky, M. S. Handcock, A. E. Raftery, and P. D. Hoff.  Representing degree distributions, clustering, and homophily in social networks with latent cluster random effects models.  _Social Networks_ , 31:204–213, 2009.
  * Krivitsky et al. [2011] P. N. Krivitsky, M. S. Handcock, and M. Morris.  Adjusting for network size and composition effects in exponential-family random graph models.  _Statistical Methodology_ , 8:319–339, 2011.
  * Krivitsky et al. [2020] P. N. Krivitsky, L. M. Koehly, and C. S. Marcum.  Exponential-family random graph models for multi-layer networks.  _Psychometrika_ , 85(3):630–659, 2020.
  * Lazega [2001] E. Lazega.  _The collegial phenomenon: The social mechanisms of cooperation among peers in a corporate law partnership_.  Oxford University Press, 2001.
  * Lei et al. [2020] J. Lei, K. Chen, and B. Lynch.  Consistent community detection in multi-layer network data.  _Biometrika_ , 107(1):61–73, 2020.
  * Li et al. [2012] W. Li, Y. Xu, J Yang, and Z. Tang.  Finding structural patterns in complex networks.  In _2012 IEEE Fifth International Conference on Advanced Computational Intelligence_ , pages 23–27, 2012.
  * Lusher et al. [2013] D. Lusher, J. Koskinen, and G. Robins.  _Exponential Random Graph Models for Social Networks_.  Cambridge University Press, Cambridge, UK, 2013.
  * Maathuis et al. [2018] M. Maathuis, M. Drton, S. Lauritzen, and M. Wainwright.  _Handbook of graphical models_.  CRC Press, 2018.
  * McPherson et al. [2001] M. McPherson, L. Smith-Lovin, and J. M. Cook.  Birds of a feather: Homophily in social networks.  _Annual Review of Sociology_ , 27:415–444, 2001.
  * Ortega and Rheinboldt [2000] J. M. Ortega and W. C. Rheinboldt.  _Iterative solution of nonlinear equations in several variables_.  SIAM, 2000.
  * Raič [2019] M. Raič.  A multivariate Berry–Esseen theorem with explicit constants.  _Bernoulli_ , 25(4A):2824–2853, 2019.
  * Rohe et al. [2011] K. Rohe, S. Chatterjee, and B. Yu.  Spectral clustering and the high-dimensional stochastic block model.  _The Annals of Statistics_ , 39:1878–1915, 2011.
  * Schweinberger et al. [2020] M. Schweinberger, P. N. Krivitsky, C. T. Butts, and Jonathan Stewart.  Exponential-family models of random graphs: Inference in finite, super, and infinite population scenarios.  _Statistical Science_ , 35:627–662, 2020.
  * Sewell and Chen [2015] D. K. Sewell and Y. Chen.  Latent space models for dynamic networks.  _Journal of the American Statistical Association_ , 110:1646–1657, 2015.
  * Snijders et al. [2006] T. A. B. Snijders, P. E. Pattison, G. L. Robins, and M. S. Handcock.  New specifications for exponential random graph models.  _Sociological Methodology_ , 36:99–153, 2006.
  * Sosa and Betancourt [2022] J. Sosa and B. Betancourt.  A latent space model for multilayer network data.  _Computational Statistics & Data Analysis_, 169:107432, 2022.
  * Stewart et al. [2019] J. Stewart, M. Schweinberger, M. Bojanowski, and M. Morris.  Multilevel network data facilitate statistical inference for curved ERGMs with geometrically weighted terms.  _Social Networks_ , 59:98–119, 2019.
  * Stewart and Schweinberger [2021] J. R. Stewart and M. Schweinberger.  Pseudo-likelihood-based $M$-estimation of random graphs with dependent edges and parameter vectors of increasing dimension.  _arXiv preprint arXiv:2012.07167_ , 2021.
  * Strauss and Ikeda [1990] D. Strauss and M. Ikeda.  Pseudolikelihood estimation for social networks.  _Journal of the American Statistical Association_ , 85:204–212, 1990.
  * Sundberg [2019] R. Sundberg.  _Statistical modelling by exponential families_ , volume 12.  Cambridge University Press, 2019.
  * Sussman et al. [2014] D. Sussman, M. Tang, and C. Priebe.  Consistent latent position estimation and vertex classification for random dot product graphs.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 36:48–57, 2014.
  * Tang et al. [2013] M. Tang, D. L. Sussman, and C. E. Priebe.  Universally consistent vertex classification for latent positions graphs.  _The Annals of Statistics_ , 41:1406–1430, 2013.
  * Wainwright and Jordan [2008] M. J. Wainwright and M. I. Jordan.  Graphical models, exponential families, and variational inference.  _Foundations and Trends in Machine Learning_ , 1:1–305, 2008.
  * Zhou and Shao [2014] M. Zhou and Y. Shao.  A powerful test for multivariate normality.  _Journal of Applied Statistics_ , 41(2):351–363, 2014.