# Bernoulli factories and duality in Wright–Fisher and Allen–Cahn models of
population genetics

Jere Koskela*, Krzysztof Łatuszyński, and Dario Spanò
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Department of Statistics, University of Warwick, Coventry CV4 7AL, UK
*Corresponding author

###### Abstract

Mathematical models of genetic evolution often come in pairs, connected by a
so-called duality relation. The most seminal example are the Wright–Fisher
diffusion and the Kingman coalescent, where the former describes the
stochastic evolution of neutral allele frequencies in a large population
forward in time, and the latter describes the genetic ancestry of randomly
sampled individuals from the population backward in time. As well as providing
a richer description than either model in isolation, duality often yields
equations satisfied by unknown quantities of interest. We employ the so-called
Bernoulli factory—a celebrated tool in simulation-based computing—to derive
duality relations for broad classes of genetics models. As concrete examples,
we present Wright–Fisher diffusions with general drift functions, and
Allen–Cahn equations with general, nonlinear forcing terms. The drift and
forcing functions can be interpreted as the action of frequency-dependent
selection. To our knowledge, this work is the first time a connection has been
drawn between Bernoulli factories and duality in models of population
genetics.

Keywords: Allen–Cahn equation, Bernoulli factory, duality, frequency-dependent
selection, Wright–Fisher diffusion

2020 MSC: 35C99, 60J70, 60J90, 92D10

##  1 Introduction

The Bernoulli factory problem is to construct a realisation of a
Bernoulli$(f(p))$ random variable (or an $f(p)$-coin) using an almost surely
finite number of independent $p$-coins, where $f:[0,1]\mapsto[0,1]$ is a known
function but $p\in[0,1]$ is unknown. The special case $f(p)=1/2$ was
formulated and solved by John von Neumann [von51]. Later, Keane and O’Brien
provided a necessary and sufficient condition for a given function $f$ to have
a Bernoulli factory [KO94]. In brief, $f$ has a Bernoulli factory if and only
if it is continuous and _polynomially bounded_ :

$\min\\{f(p),1-f(p)\\}\geq\min\\{p,1-p\\}^{n}$ (1)

for all $p\in[0,1]$ and some $n\geq\mathbb{N}$, or identically equal to zero
or one. However, the proof of Keane and O’Brien is only partly constructive:
it relies on a recursively defined sequence whose explicit solution is
intractable. Constructions of algorithms have relied of approximating $f$ by
Bernstein polynomials, which are naturally associated with $p$-coins [HNP11,
ŁKPR11, NP05], or with other series expansions of $f$ with non-negative
coefficients [Men19].

Many seminal models of population genetics rely on the random propagation of
alleles from one generation to the next. A prototypical example is the
Wright–Fisher model, in which a population of fixed size $N\in\mathbb{N}$
evolves in discrete generations. Individuals carry one of two alleles, $a$ or
$A$, and each individual inherits the allele of a parent which it samples
independently and uniformly from the previous generation. The mechanism of
sampling alleles by sampling parents ensures that the model carries
information of the forward-in-time evolution of allele frequencies, as well as
the backward-in-time genealogies of samples of individuals. Frequently, these
two modelling perspectives satisfy a duality relation which renders both
models more tractable than they would be in isolation. We direct interested
readers to e.g. [Dur08] for an introduction to Wright–Fisher and genealogical
models in population genetics.

In the absence of mutation, inheriting an allele from a uniformly sampled
parent models _neutral_ evolution, where the conditional mean allele frequency
in a generation equals that in the previous generation. Non-neutral models in
which the mean is not constant can be obtained by sampling offspring alleles
as $f(p)$-coins when the allele frequency in the previous generation is $p$,
and $f$ models so-called _frequency-dependent selection_. In order to retain
the aforementioned backward-in-time genealogical picture and its associated
duality, it is desirable to generate the $f(p)$-coins in a given generation by
sampling parental alleles from the previous generation. Since the allele
frequency in the parental generation is $p$, parental alleles can be thought
of as $p$-coins, motivating a connection to Bernoulli factories. Our
contribution is to use a Bernoulli factory to extend two standard models of
population genetics to more general settings than has been done previously.
They are (i) the non-neutral Wright–Fisher diffusion and its ancestral
selection graph dual [CHS22, EGT10, GS18, KN97], for which the function $f$
models frequency-dependent selection as described above, and (ii) the
Allen–Cahn PDE which models stationary allele frequencies in a spatial
continuum, in which $f$ appears as external forcing [EGL22].

The remainder of the manuscript is organised as follows. In Section 2, we
review the Bernoulli factory of [KO94] which turns out to be convenient for
our purposes. In Sections 3 and 4 we introduce the Wright–Fisher diffusion
with frequency–dependent selection and the Allen–Cahn equation. In each case,
we also demonstrate how Bernoulli factories facilitate the construction of
very large classes of these models. Section 5 concludes with a discussion on
connections to earlier results, as well as some potential extensions.

##  2 The Keane–O’Brien factory

Let $f:[0,1]\mapsto[0,1]$ be continuous and polynomially bounded as in (1).
Let $\bar{X}_{n}(p)$ be the sample mean of $n$ independent $p$-coins. Define
sequences of functions $f_{k}$ and integers $\eta(f,k)$ as follows: set
$f_{1}(p):=f(p)$, and

$f_{k+1}(p):=\frac{4}{3}\Bigg{(}f_{k}(p)-\frac{1}{4}\mathbb{P}\Big{(}f_{k}(\bar{X}_{\eta(f,k)}(p))\geq\frac{1}{2}\Big{)}\Bigg{)},$

where each $\eta(f,k)$ is finite, independent of $p,$ and large enough, so
that

$f_{k}(p)-\frac{1}{4}\mathbb{P}\Big{(}f_{k}(\bar{X}_{\eta(f,k)}(p))\geq\frac{1}{2}\Big{)}\in[0,3/4].$

Such $\eta(f,k)$ exists by [KO94, page 218]. Let
$L\sim\operatorname{Geo}(1/4)$ be independent of all $p$-coins. Then

$\mathds{1}\Big{\\{}f_{L}(\bar{X}_{\eta(f,L)}(p))\geq\frac{1}{2}\Big{\\}}\sim\operatorname{Ber}(f(p))$

is a Bernoulli factory for $f$ [KO94, pages 217–219]. It is based on the
series expansion

$\displaystyle f(p)$
$\displaystyle=\sum_{k=1}^{\infty}\Big{(}\frac{3}{4}\Big{)}^{k-1}\frac{1}{4}\mathbb{P}\Big{(}f_{k}(\bar{X}_{\eta(f,k)}(p))\geq\frac{1}{2}\Big{)}$
$\displaystyle=\sum_{k=1}^{\infty}\Big{(}\frac{3}{4}\Big{)}^{k-1}\frac{1}{4}\sum_{j=0}^{\eta(f,k)}\binom{\eta(f,k)}{j}p^{j}(1-p)^{\eta(f,k)-j}\mathds{1}\Big{\\{}f_{k}\Big{(}\frac{j}{\eta(f,k)}\Big{)}\geq\frac{1}{2}\Big{\\}},$
(2)

which converges uniformly in $p$. Pseudocode for this Keane–O’Brien factory is
shown in Algorithm 1.

Algorithm 1 Bernoulli factory for continuous $f:[0,1]\mapsto[0,1]$ satisfying
(1).

1:Sequences $\\{f_{k}\\}_{k\geq 1}$, $\\{\eta(f,k)\\}_{k\geq 1}$.

2:Sample $L\sim\operatorname{Geo}(1/4)$.

3:Sample $\bar{X}_{\eta(f,L)}(p)\sim\operatorname{Bin}(\eta(f,L),p)$.

4:if $f_{L}(\bar{X}_{\eta(f,L)}(p)/L)\geq 1/2$ then

5: Return 1.

6:else

7: Return 0.

8:end if

It is worth highlighting two features of the Keane–O’Brien factory which will
turn out to be essential.

######  Remark 1.

The Keane–O’Brien factory is defined for all $p\in[0,1]$ and $f(p)\in[0,1]$.
Factories based on Bernstein polynomial approximation typically have to
constrain the range (and sometimes the domain) of $f$ to (subsets of) $(0,1)$
to avoid degenerate distributions [HNP11, ŁKPR11, NP05]. In population
genetics applications, it is essential to allow the whole range $p\in[0,1]$,
and desirable to allow $f(0)=0$ and $f(1)=1$ to model fixation.

######  Remark 2.

The random variable $L$, and hence the number of $p$-coins $\eta(f,L)$ needed
to determine the allele of a child, is independent of the realisations of the
coins. This will facilitate the construction of ancestral graphs containing
all possible ancestors before the alleles carried by any of those ancestors
are known.

##  3 The frequency-dependent Wright–Fisher diffusion

Consider a population of $N$ individuals $\\{\mathbf{Z}_{k}^{(N)}\\}_{k\geq
0}$ evolving in discrete time, where the state of the system at time $k$ is
$\mathbf{Z}_{k}^{(N)}:=(Z_{1}^{(N)}(k),\ldots,Z_{N}^{(N)}(k))$ with
$Z_{k}^{(N)}\in\\{a,A\\}$. Let

$\tilde{Y}_{k}^{(N)}:=\frac{1}{N}\sum_{i=1}^{N}\mathds{1}_{\\{A\\}}(Z_{i}^{(N)}(k))$

be the frequency of the $A$ allele in generation $k$. Each individual in
generation $k+1$ samples its allele conditionally independently given
$\tilde{Y}_{k}^{(N)}$, where the probability of an $A$ allele is
$f^{(N)}(Y_{k}^{(N)})$ with

$f^{(N)}(p):=\Big{(}1-\frac{\sigma}{N}\Big{)}p+\frac{\sigma}{N}f(p)$

for a fixed constant $\sigma>0$ and continuous function $f:[0,1]\mapsto[0,1]$
satisfying (1).

Concretely, the realisation of each allele is implemented as follows. With
probability $1-\sigma/N$, a given individual samples one parent and inherits
its allele. With the complementary probability, the individual samples an
independent copy of $L$, followed by $\eta(f,L)$ parents chosen uniformly at
random from the population with replacement. Conditional on $L$ and
$\tilde{Y}_{k}^{(N)}$, the allele of the individual is $A$ if event
$\\{f_{L}(\bar{X}_{\eta(f,L)})\geq\tilde{Y}_{k}^{(N)}\\}$ happens. By (2), the
marginal probability of an $A$ allele is $f^{(N)}(\tilde{Y}_{k}^{(N)})$.

Our aim is to verify convergence of the allele frequency process
$\\{\tilde{Y}_{k}^{(N)}\\}_{k\geq 0}$ in the infinite population limit. To
that end, we define the continuous-time process $(Y_{t}^{(N)})_{t\geq 0}$ via

$Y_{t}^{(N)}:=\tilde{Y}_{\lfloor Nt\rfloor}^{(N)}.$

######  Theorem 1.

Let $f$ be polynomially bounded and Lipschitz continuous on $[0,1]$. Then, in
the Skorokhod topology, $(Y_{t}^{(N)})_{t\geq 0}\to(Y_{t})_{t\geq 0}$ weakly
as $N\to\infty$, where $Y_{t}$ solves

$\mathrm{d}Y_{t}=\sigma(f(Y_{t})-Y_{t})\mathrm{d}t+\sqrt{Y_{t}(1-Y_{t})}\mathrm{d}W_{t}$
(3)

subject to appropriate initial conditions, where $(W_{t})_{t\geq 0}$ is a
Brownian motion.

###### Proof.

Existence, uniqueness, and the Feller property of the putative limiting
process all follow from [EK86, Theorem 2.8] because the drift function
$\sigma(f(y)-y)$ is Lipschitz continuous.

The discrete-time process $\bar{Y}_{k}^{(N)}$ has generator

$\displaystyle
L^{N}h(y):=\sum_{k=0}^{N}\binom{N}{k}\Big{(}\frac{\sigma}{N}\Big{)}^{k}\Big{(}1-\frac{\sigma}{N}\Big{)}^{N-k}$
$\displaystyle\sum_{x=0}^{N-k}\binom{N-k}{x}y^{x}(1-y)^{N-k-x}$
$\displaystyle\times\sum_{z=0}^{k}\binom{k}{z}f(y)^{z}(1-f(y))^{k-z}h\Big{(}\frac{x+z}{N}\Big{)}-h(y).$

To show the claimed convergence, take $h$ to be a $C^{2}$ function with
bounded third derivatives, which is a convergence-determining class. Expanding
$h$ around $y$ on the right-hand side yields

$\displaystyle L^{N}h(y)={}$
$\displaystyle\sum_{k=0}^{N}\binom{N}{k}\Big{(}\frac{\sigma}{N}\Big{)}^{k}\Big{(}1-\frac{\sigma}{N}\Big{)}^{N-k}\sum_{x=0}^{N-k}\binom{N-k}{x}y^{x}(1-y)^{N-k-x}$
$\displaystyle\times\sum_{z=0}^{k}\binom{k}{z}f(y)^{z}(1-f(y))^{k-z}\Bigg{[}\Big{(}\frac{x+z}{N}-y\Big{)}h^{\prime}(y)+\frac{1}{2}\Big{(}\frac{x+z}{N}-y\Big{)}^{2}h^{\prime\prime}(y)\Bigg{]},$

up to terms which are of lower order since $h^{\prime\prime\prime}$ is
bounded. Evaluating the binomial expectations yields

$NL^{N}h(y)=\sigma(f(y)-y)h^{\prime}(y)+\frac{1}{2}y(1-y)h^{\prime\prime}(y)+o(1),$

as $N\to\infty$. The proof of weak convergence is completed by [Kal02, Theorem
19.28] because the limiting diffusion is Feller. ∎

######  Remark 3.

The requirement of a Lipschitz drift could be slightly relaxed by using the
more cumbersome conditions for a Feller semigroup given in [XYZ19, Theorem
3.3]. They cover drifts satisfying a condition akin to

$|f(y)-f(z)|\leq-C|y-z|\log(|y-z|),$

or small variations thereof, where $C>0$ is a constant. See equations (14) and
(15), as well as Remark 2.3, of [XYZ19] for a precise class of non-Lipschitz
functions which can be handled.

Next we consider a sample of $n\in\mathbb{N}$ individuals from a given
generation (which we say lived at time 0) in the pre-limiting Wright–Fisher
model $\\{\mathbf{Z}_{k}^{(N)}\\}_{k\geq 0}$. We define the _ancestral
process_ $A_{k}^{(N)}$ as the number of lineages which are ancestral to the
sample $k$ generations in the past. The number of lineages can decrease by one
when two lineages find a common ancestor, or increase by $\eta(f,L)-1$
whenever a lineage samples $\eta(f,L)$ ancestors. In the pre-limiting particle
system, any number of these events can co-occur in one generation,
particularly when $n\geq 3$. But transitions other than isolated binary
mergers and single multifurcations turn out to vanish in a suitably rescaled
infinite population limit. To that end, we define the continuous time Markov
jump process

$A_{t}:=\lim_{N\to\infty}A_{\lfloor Nt\rfloor}^{(N)}$ (4)

whose existence we prove next.

######  Theorem 2.

The limit in (4) exists, and $(A_{t})_{t\geq 0}$ has generator

$Gh(n)=\binom{n}{2}[h(n-1)-h(n)]+\sigma
n\sum_{k=1}^{\infty}\Big{(}\frac{3}{4}\Big{)}^{k-1}\frac{1}{4}[h(n+\eta(f,k)-1)-h(n)].$

###### Proof.

The pre-limiting, discrete-time process $\\{A_{k}^{(N)}\\}_{k\geq 0}$
undergoes a myriad of transitions involving subsets of individuals finding
common ancestors, as well as individuals branching into many potential
ancestors, potentially in the same generation. However, only transitions with
a per-generation probability $\Theta(1/N)$ will contribute to the time-
rescaled limit with a finite rate. Events with probability $o(1/N)$ will not
occur in the limit at all, while events with probability $\omega(1/N)$ will
appear at a dense set of times. However, it will turn out that the latter only
result in identity transitions in $A_{t}$, and hence do not affect the limit.

The probability of two lineages originating from a common ancestor one
generation earlier, with neither being involved in a branching event, is

$\Big{(}1-\frac{\sigma}{N}\Big{)}^{2}\frac{1}{N}=\frac{1}{N}+o(1/N),$ (5)

implying that two lineages will merge to a common ancestor at rate 1 in the
limit. A triple merger, or more than one simultaneous merger, has probability
at most

$\Big{(}1-\frac{\sigma}{N}\Big{)}^{3}\frac{1}{N^{2}}=o(1/N),$

and hence will not appear in the limit.

A single individual branches into $\eta(f,k)$ ancestral lineages with
probability

$\frac{\sigma}{N}\Big{(}\frac{3}{4}\Big{)}^{k-1}\frac{1}{4},$ (6)

while any event involving more than two lineages branching in one generation
has probability at most

$\Big{(}\frac{\sigma}{N}\Big{)}^{2}=o(1/N).$

Hence, only isolated branching events appear in the limit. It is also clear
from (5) and (6) that the probability of at least one merger and branching
event in one generation is $o(1/N)$.

All other transitions involve no mergers or branching events, and hence do not
affect the limiting ancestral process. Noting that there are $\binom{n}{2}$
pairs of individuals to merge with probability (5) and $n$ individuals to
branch with probabilities (6) yields the claimed generator $G$. ∎

Duality between the Wright–Fisher diffusion (3) and the ancestral process (4)
is a relation

$\mathbb{E}_{y}[h(Y_{t},n)|Y_{0}=y]=\mathbb{E}_{n}[h(y,A_{t})|A_{0}=n],$

where $y\in[0,1]$ and $n\in\mathbb{N}$ are respective initial conditions, and
$h$ is a _duality function_ , the specification of which will require some
exposition.

We follow [CHS22] and define the random function $P_{t}(y)$ as the conditional
probability that all $n$ leaves at time zero in the ancestral process carry
allele $A$, given that the $A$ allele frequency at time $t$ in the past is
$y\in[0,1]$, and given the realisation of the ancestral process
$A_{0:t}:=(A_{s})_{s\in[0,t]}$ started from the $n$ lineages. For example, in
the absence of branching events we have $P_{t}(y)=y^{A_{t}}$, while when
$n=1$, a single branching event into $\eta(f,k)$ ancestors and no mergers in
the history $A_{0:t}$ yields

$P_{t}(y):=\sum_{j=0}^{\eta(f,k)}\binom{\eta(f,k)}{j}y^{j}(1-y)^{\eta(f,k)-j}\mathds{1}\Big{\\{}f_{k}\Big{(}\frac{j}{\eta(f,k)}\Big{)}\geq\frac{1}{2}\Big{\\}},$

where $f_{k}$ is as specified in Section 2. Other patterns of merger and
branching events will result in a more complicated Bernstein polynomial of
degree $A_{t}$,

$P_{t}(y):=\sum_{j=0}^{A_{t}}V_{t}(j)\binom{A_{t}}{j}y^{j}(1-y)^{A_{t}-j},$

where $V_{t}(j)$ is the random Bernstein coefficient which equals the
probability that the $n$ leaves all carry allele $A$, given that $j$ of
$A_{t}$ roots do. As in [CHS22, Definition 2.12 and Proposition 2.13], the
vector $(\bm{V}_{t})_{t\geq 0}:=(V_{t}(0),\ldots,V_{t}(n))_{t\geq 0}$ is also
a Markov jump process with transitions

$\displaystyle\bm{v}\mapsto\Bigg{(}\sum_{j=0}^{i\wedge\eta(f,k)}\frac{\binom{i}{j}\binom{n+\eta(f,k)-1-j}{\eta(f,k)-j}}{\binom{n+\eta(f,k)-1}{i}}[$
$\displaystyle\mathds{1}\\{f_{h}(j/\eta(f,k))\geq 1/2\\}v_{i+1-j}$
$\displaystyle+\mathds{1}\\{f_{k}(j/\eta(f,k))<1/2\\}v_{i-j}]\Bigg{)}_{i=0}^{n+\eta(f,k)-1}$
(7)

at rate $n\sigma\mathbb{P}(L=k)$, and

$\bm{v}\mapsto\Big{(}\frac{i}{n-1}v_{i+1}+\frac{n-1-i}{n-1}v_{i}\Big{)}_{i=0}^{n-1}$
(8)

at rate $\binom{n}{2}$. Also, let

$\bm{B}_{n}(x):=\Bigg{(}\binom{n}{i}x^{i}(1-x)^{n-i}\Bigg{)}_{i=1}^{n}$

be the vector of order $n$ Bernstein polynomials, and for vectors $\bm{u}$ and
$\bm{v}$ of common length $n+1$, let

$\langle\bm{u},\bm{v}\rangle:=\sum_{i=0}^{n}u_{i}v_{i}$

be the usual inner product.

######  Theorem 3.

The processes $(Y_{t})_{t\geq 0}$ and $(V_{t})_{t\geq 0}$ are dual, in that
for all $t\geq 0$,

$\mathbb{E}[\langle\bm{B}_{n}(Y_{t}),\bm{v}\rangle|Y_{0}=y]=\mathbb{E}[\langle\bm{B}_{\text{dim}(\bm{V}_{t})-1}(y),\bm{V}_{t}\rangle|\bm{V}_{0}=\bm{v}],$
(9)

for each $y\in[0,1]$, every $n\in\mathbb{N}$, and each
$\bm{v}\in\mathbb{R}^{n+1}$.

###### Proof.

The proof is a small adaptation of that of [CHS22, Theorem 2.14] to our more
general setting. Duality between the coalescing mechanism of the ancestral
process (8) and the diffusion coefficient of the Wright–Fisher diffusion (3)
is standard, and we omit it to focus on establishing the same relation for the
branching mechanism (7) and the drift term in (3).

Following [CHS22, Section 4.3], let $Y_{k}^{y}\sim\text{Bin}(k,y)$ and
$K_{k,i}^{n}\sim\text{Hyp}(n+k-i,k,i)$ be independent random variables, where
$\text{Hyp}(a,b,c)$ denotes the hypergeometric distribution with $b$ draws
without replacement from a population of size $a$ containing $c$ successes.
Let $H:(y,\bm{v})\mapsto\langle\bm{B}_{\text{dim}(\bm{v})-1}(y),\bm{v}\rangle$
be the putative duality function. Then, for $\bm{v}\in\mathbb{R}^{n+1}$ we
have

$\displaystyle\partial_{y}H(y,\bm{v})$
$\displaystyle=n\mathbb{E}[v_{Y_{n-1}^{y}+1}-v_{Y_{n-1}^{y}}],$ $\displaystyle
H(y,\bm{v})$
$\displaystyle=\mathbb{E}[v_{Y_{n}^{y}}]=\mathbb{E}[(1-y)v_{Y_{n-1}^{y}}+yv_{Y_{n-1}^{y}+1}],$
$\displaystyle f(y)$
$\displaystyle=\mathbb{E}[\mathds{1}\\{f_{L}(\bar{X}_{\eta(f,L)}\geq 1/2\\}],$

where $L\sim\text{Geo}(1/4)$ and the third equality is due to (2). We also
have

$(\bar{X}_{\eta(f,L)},Y_{n-1}^{y})\stackrel{{\scriptstyle
d}}{{=}}(K_{\eta(f,k),Y_{n+\eta(f,L)-1}^{y}}^{n},Y_{n+\eta(f,L)-1}^{y}-K_{\eta(f,k),Y_{n+\eta(f,L)-1}^{y}}^{n}).$

Thus the drift of the Wright–Fisher diffusion (3) satisfies

$\displaystyle\sigma(f(y)-y)\partial_{y}H(y,\bm{v})$
$\displaystyle=\sigma\mathbb{E}[\mathds{1}\\{f_{L}(\bar{X}_{\eta(f,L)}(y))\geq
1/2\\}-y]n\mathbb{E}[v_{Y_{n-1}^{y}+1}-v_{Y_{n-1}^{y}}]$
$\displaystyle=n\sigma\mathbb{E}[\mathds{1}\\{f_{L}(\bar{X}_{\eta(f,L)}(y))\geq
1/2\\}v_{Y_{n-1}^{y}+1}-\mathds{1}\\{f_{L}(\bar{X}_{\eta(f,L)}(y))\geq
1/2\\}v_{Y_{n-1}^{y}}-yv_{Y_{n-1}^{y}+1}$
$\displaystyle\phantom{=n\sigma\mathbb{E}[}+yv_{Y_{n-1}^{y}}+v_{Y_{n-1}^{y}}-v_{Y_{n-1}^{y}}]$
$\displaystyle=n\sigma(\mathbb{E}[\mathds{1}\\{f_{L}(\bar{X}_{\eta(f,L)}(y))\geq
1/2\\}v_{Y_{n-1}^{y}+1}+\mathds{1}\\{f_{L}(\bar{X}_{\eta(f,L)}(y))<1/2\\})v_{Y_{n-1}^{y}}]-H(y,\bm{v}))$
$\displaystyle=n\sigma\Bigg{(}\sum_{k=1}^{\infty}\Big{(}\frac{3}{4}\Big{)}^{k-1}\frac{1}{4}\sum_{j=0}^{n+\eta(f,k)-1}\mathbb{E}\Bigg{[}\mathds{1}\Bigg{\\{}f_{k}\Bigg{(}\frac{K_{\eta(f,k),j}^{n}}{\eta(f,k)}\Bigg{)}\geq
1/2\Bigg{\\}}v_{j-K_{\eta(f,k),j}^{n}+1}$
$\displaystyle\phantom{\ldots\ldots\ldots\ldots}+\mathds{1}\Bigg{\\{}f_{k}\Bigg{(}\frac{K_{\eta(f,k),j}^{n}}{\eta(f,k)}\Bigg{)}<1/2\Bigg{\\}}v_{j-K_{\eta(f,k),j}^{n}}\Bigg{]}\binom{n+\eta(f,k)-1}{j}y^{j}(1-y)^{n+\eta(f,k)-1-j}$
$\displaystyle\phantom{=n\sigma\mathbb{E}[}-H(y,\bm{v})\Bigg{)},$

which is precisely (7) applied to the $\bm{v}$-argument of
$H(y,\bm{v})=\langle\bm{B}_{\text{dim}(\bm{v})-1}(y),\bm{v}\rangle$. ∎

##  4 The Allen–Cahn model of spatial genetics

The Allen–Cahn equation on a Lipschitz domain $\Omega\subseteq\mathbb{R}^{d}$
is

$\partial_{t}u-\Delta u=\lambda f(u),$ (10)

for a given $\lambda>0$ and $f:[0,1]\mapsto[0,1]$, and subject to suitable
initial and boundary conditions. In [EGL22], the authors consider a model of
the spatially structured frequency $u(x,t)$ of an allele $A$ governed by

$\displaystyle\partial_{t}u-\Delta u$
$\displaystyle=\frac{1}{\varepsilon^{2}}u(1-u)(2u-1+\nu\varepsilon)\text{ for
}x\in\Omega,t>0,$ (11) $\displaystyle\partial_{n}u$ $\displaystyle=0\text{ for
}x\in\partial\Omega,t>0,$ $\displaystyle u(x,0)$ $\displaystyle=u_{0}(x)\text{
for }x\in\Omega,$

where $\varepsilon>0$, $\nu>0$, $\partial_{n}$ denotes the normal derivative
at the boundary, and $u_{0}:\Omega\to[0,1]$ (see also [Goo18]). They construct
a solution to (11) by using a particle system in which a single particle
started at a location $x\in\Omega$ undergoes a Brownian motion with speed 2,
branches into three independent copies at rate
$(1+\varepsilon\nu)/\varepsilon^{2}$, and particles are reflected from the
boundary. At a given end time $t>0$, a leaf particle at location $z\in\Omega$
samples one of two alleles, $a$ and $A$, with respective probabilities
$(1-u_{0}(z),u_{0}(z))$. All particles sample their alleles independently.
Then, particles propagate their alleles rootwards along the realisation of the
Brownian tree. The allele of an internal branch is decided by a majority vote
among its three children unless exactly one child carries allele $A$, in which
case the parent branch carries $A$ with probability
$2\nu\varepsilon/(3+3\nu\varepsilon)$, and $a$ otherwise. The probability that
the root particle $x\in\Omega$ carries allele $A$ under these dynamics solves
(11) [EGL22, Proposition 2.4].

The ingredients of the particle system can be read off from (11). The
Laplacian $\Delta$ is the generator of Brownian motion run at speed 2, the
branching rate $(1+\varepsilon\nu)/\varepsilon^{2}$ is an upper bound for the
right-hand side, and

$\frac{1}{\varepsilon^{2}}u(1-u)(2u-1+\nu\varepsilon)=\frac{1+\nu\varepsilon}{\varepsilon^{2}}\Bigg{(}u^{3}+3u^{2}(1-u)+\frac{2\nu\varepsilon}{3+3\nu\varepsilon}3u(1-u)^{2}-u\Bigg{)},$

which demonstrates that the nonlinearity in (11) has an interpretation as the
voting system described above. It would be straightforward to adapt the proof
of [EGL22, Proposition 2.4] to any other polynomial right-hand side of the
form

$\lambda\sum_{j=0}^{m}p_{j}\binom{m}{j}u^{j}(1-u)^{m-j}-u,$

for some $m\in\mathbb{N}$, $\lambda>0$, and coefficients $p_{j}\in[0,1]$ by
setting $\lambda$ as the branching rate into $m$ particles, and suitably
adapting the voting scheme. Our contribution is to prove an analogous result
for (10), subject to the same initial and boundary conditions as (11), when
$f$ is merely continuous and polynomially bounded.

Consider a branching Brownian motion reflected off the boundary
$\partial\Omega$, with branching rate $\lambda$, started from a single
particle at $x\in\Omega$ at time $0$. At a branching event, the number of
offspring is given by $\eta(f,L)$, where $L\sim\text{Geo}(1/4)$. At a terminal
time $t>0$, a leaf particle at $z\in\Omega$ samples allele $a$ (resp. $A$)
with probability $1-u_{0}(z)$ (resp. $u_{0}(z)$), and each leaf carries out
this choice independently. Alleles are propagated rootwards along the tree: a
branch with $\eta(f,k)$ offspring, $j$ of whom carry allele $A$, carries
allele $A$ if $f_{k}(j/\eta(f,k))\geq 1/2$, and carries allele $a$ otherwise.
Let $F(x,t)$ be the allele carried by the root particle at position
$x\in\Omega$ when the branching process is run until time $t>0$.

######  Theorem 4.

Suppose $f$ is continuous and polynomially bounded as in (1). Viewed as a
function of $x\in\Omega$, the probability $\mathbb{P}(F(x,t)=A)$ solves (10)
subject to the initial and boundary conditions in (11).

###### Proof.

The proof is an adaptation of that of [EGL22, Proposition 2.4] to our more
general setting. To verify that $q(x,t):=\mathbb{P}(F(x,t)=A)$ solves (10) on
the interior of $\Omega$, let $S$ denote the first branching time of the
initial particle and $(W_{t})_{t\geq 0}$ be a Brownian motion. For small
$h>0$,

$\displaystyle q(x,t+h)$
$\displaystyle=\mathbb{P}(F(x,t+h)=A|S>h)\mathbb{P}(S>h)+\mathbb{P}(F(x,t+h)=A|S\leq
h)\mathbb{P}(S\leq h)$
$\displaystyle=\mathbb{E}_{x}[q(W_{h},t)|S>h]\mathbb{P}(S>h)$
$\displaystyle\phantom{=}+\mathbb{P}(S\leq
h)\sum_{k=1}^{\infty}\Big{(}\frac{3}{4}\Big{)}^{k-1}\frac{1}{4}\sum_{j=0}^{\eta(f,k)}\binom{\eta(f,k)}{j}\mathbb{E}_{x}[q(W_{S},t)^{j}(1-q(W_{S},t))^{\eta(f,k)-j}]$
$\displaystyle\phantom{=+\mathbb{P}(S\leq
h)\sum_{k=1}^{\infty}\Big{(}\frac{3}{4}\Big{)}^{k-1}\frac{1}{4}\sum_{j=0}^{\eta(f,k)}}\times\mathds{1}\Bigg{\\{}f_{k}\Big{(}\frac{j}{\eta(f,k)}\Big{)}\geq\frac{1}{2}\Bigg{\\}}$
$\displaystyle=\mathbb{E}_{x}[q(W_{h},t)|S>h]e^{-\lambda
h}+\mathbb{E}_{x}[f(q(W_{S},t))|S\leq h](1-e^{-\lambda h}),$

where the subscript in $\mathbb{E}_{x}$ denotes the starting point of
$(W_{t})_{t\geq 0}$ and the last equality follows via (2). Since $f$ is
continuous, regularity of the heat semigroup and the tower law (viewing
$q(x,t)$ as the expectation of an indicator function) yield

$\mathbb{E}_{x}[f(q(W_{S},t))|S\leq
h]=f(\mathbb{E}_{x}[q(W_{s},t)])+O(h)=f(q(x,t))+O(h).$

Hence,

$\displaystyle\partial_{t}\mathbb{P}(F(x,t)=A)$ $\displaystyle=\lim_{h\to
0}\frac{\mathbb{E}_{x}[\mathbb{P}(F(W_{h},t)=A)|S>h]-\mathbb{P}(F(x,t)=A)}{h}e^{-\lambda
h}$ $\displaystyle\phantom{=}+\lim_{h\to
0}\frac{f(\mathbb{P}(F(x,t)=A))-\mathbb{P}(F(x,t)=A)}{h}(1-e^{-\lambda h})$
$\displaystyle=\Delta\mathbb{P}(F(x,t)=A)+\lambda[f(\mathbb{P}(F(x,t)=A))-\mathbb{P}(F(x,t)=A)].$

The boundary condition is inherited from reflecting Brownian motion [BH91]
since branching events occur at a finite rate and hence will not take place on
the boundary. ∎

##  5 Discussion

In [CHS22], the authors prove analogues of Theorems 1–3 for the case

$f(x)-x:=\sum_{\ell=2}^{m}\beta_{\ell}\sum_{i=0}^{\ell}\binom{\ell}{i}x^{i}(1-x)^{\ell-i}\Big{(}p_{i,\ell}-\frac{i}{\ell}\Big{)},$
(12)

for a fixed constant $m\in\mathbb{N}$, positive coefficients
$\\{\beta_{\ell}\\}_{\ell=2}^{m}$, and a sequence of $[0,1]$-valued
coefficients $\\{\\{p_{i,\ell}\\}_{i=0}^{\ell}\\}_{\ell=2}^{m}$. Earlier work
by González Casanova and Spanò also covered the case

$f(x):=\sum_{j=0}^{\infty}\pi_{j}x^{j},$

where $\\{\pi_{j}\\}_{j=0}^{\infty}$ is a probability mass function [GS18].
Our results cover all Lipschitz continuous, polynomially bounded functions
$f:[0,1]\mapsto[0,1]$, which includes both of these classes as special cases.
Furthermore, [CHS22, Section 2.10] mentions that their approach should extend
to the $m=\infty$ case. Our Bernoulli factory approach demonstrates that that
is true, and also that there is no further difficulty in handling our more
general class of drift functions.

The works of [GS18] and [CHS22] are more general than our results in two ways:
they tackle sequences of drift functions $f_{N}\to f$ as $N\to\infty$, and
they incorporate jumps in the limiting forward-in-time Wright–Fisher
diffusion, along with multiple mergers into the reverse-time ancestral
process. Both of these generalisations could be incorporated into our model at
the cost of increased technicality. We have chosen to omit them to focus on
the class of drift functions, which is our main interest.

The link our work establishes between individual-based models and diffusive
scaling limits provides rigorous justification for a range of diffusive
approximations which have been obtained for non-neutral finite-population
models [TN06, LL07]. In addition, we provide an associated genealogical
description and a formal duality relation between the two processes.
Genealogical processes of Wright–Fisher diffusions with frequency-dependent
selection are also considered in [CG04], though their approach is to condition
on the allele frequency trajectory and hence avoid the need for branching
events in the ancestral process. Their setting is formulated for generic drift
functions of the form $\beta(x)x(1-x)$, though in practice they focus on
Bernstein polynomials of degree no more than two.

Approaches based on moment duality, with duality function $h(x,n)=x^{n}$, have
been used to obtain series expansions of Wright–Fisher diffusion transition
functions [BEG00], and [CHS22] use Bernstein duality in (9) to study fixation
when the drift is of the form (12) using the Bernstein coefficient process
$(\mathbf{V}_{t})_{t\geq 0}$. It is unclear whether similar results can be
usefully obtained in our setting in practice because the sequences
$\\{f_{k}\\}_{k\geq 1}$ and $\\{\eta(f,k)\\}_{k\geq 1}$ are difficult to
compute. Hence, so are the Bernstein coefficients.

The construction of the solution to the Allen–Cahn equation (10) from a
branching Brownian motion is also an example of duality between these two
processes. This result could also be generalised by replacing the Laplacian
with a more general second order differential operator, provided that it has
suitable regularity and generates a diffusion with tractable reflecting
behaviour at boundaries.

The key advantage of the Keane–O’Brien factory is that it covers the whole
class of continuous, polynomially bounded functions $f$. Other factories could
be used to obtain different ancestral processes, Wright–Fisher diffusions, and
constructions of the solution to the Allen–Cahn equation. However, it is
essential that the number of coins needed by the factory is independent of the
realisations of those coins, or at least that there is a almost surely finite
upper bound on the number of coins regardless of realisations. Otherwise, the
trick of keeping track of all resulting branches to fill in alleles (or votes)
later cannot work, because the number of branches cannot be determined until
earlier alleles have been resolved. As far as we are aware, the only other
somewhat general Bernoulli factory with this independence property is that of
[Men19], which is identical to the branching mechanism used in [GS18] and
applies to exactly the same functions. Seminal Bernoulli factories, such as
that of Nacu and Peres [NP05] based on Bernstein polynomial approximation of
$f$, inherently link the decision to keep flipping coins to the outcomes of
earlier coins.

Finally, there are multivariate analogues of Bernoulli factories, in which
independent, $m$-sided dice with probability mass function
$(p_{1},\ldots,p_{m})$ are used to construct a $v$-sided die with mass
function $f(p_{1},\ldots,p_{m})$. To date, attention has focused on domains
which exclude boundaries, and where the coordinates of $f$ are rational
functions [Mor21, MŁNW22], or where $v=1$ [PLS]. We believe a suitable
extension of such _dice enterprises_ to cases including boundaries could be
used to obtain analogues of the convergence and duality results presented here
in the case with more than two non-neutral alleles.

## Funding and data sharing statements

JK was supported by EPSRC research grant EP/V049208/1. Data sharing is not
applicable to this article as no new data were generated or analysed.

## References

  * [BEG00] A. D. Barbour, S. N. Ethier, and R. C. Griffiths.  A transition function expansion for a diffusion model with selection.  Ann. Appl. Probab., 10:123–162, 2000.
  * [BH91] R. F. Bass and P. Hsu.  Some potential theory for reflecting brownian motion in Holder and Lipschitz domains.  Ann. Probab., 19(2):486–508, 1991.
  * [CG04] G. Coop and R. C. Griffiths.  Ancestral inference on gene trees under selection.  Theor. Popul. Biol., 66:219–232, 2004.
  * [CHS22] F. Cordero, S. Hummel, and E. Schertzer.  General selection models: Bernstein duality and minimal ancestral structures.  Ann. Appl. Probab., 32(3):1499–1556, 2022.
  * [Dur08] R. Durrett.  Probability models for DNA sequence evolution.  Springer, 2008.
  * [EGL22] A. M. Etheridge, M. D. Gooding, and I. Letter.  On the effects of a wide opening in the domain of the (stochastic) Allen–Cahn equation and the motion of hybrid zones.  Electron. J. Probab., 27:Article 1, 2022.
  * [EGT10] A. M. Etheridge, R. C. Griffiths, and J. E. Taylor.  A coalescent dual process in a moran model with genic selection, and the Lambda coalescent limit.  Theor. Popln Biol., 78(2):77–92, 2010.
  * [EK86] S. N. Ethier and T. G. Kurtz.  Markov Processes: Characterization and Convergence.  John Wiley & Sons, Inc., 1986.
  * [Goo18] M. Gooding.  Long term behaviour of spatial population models with heterozygous or asymmetric homozygous selection.  PhD thesis, Oxford University, 2018.
  * [GS18] A. Gonzalez Casanova and D. Spanò.  Duality and fixation in $\Xi$-Wright–Fisher processes in frequency-dependent selection.  Ann. Appl. Probab., 28(1):250–284, 2018.
  * [HNP11] O. Holtz, F. Nazarov, and Y. Peres.  New coins from old, smoothly.  Constr. Approx., 33:331–363, 2011.
  * [Kal02] O. Kallenberg.  Foundations of Modern Probability.  Springer New York, NY, 2nd edition, 2002.
  * [KN97] S. M. Krone and C. Neuhauser.  Ancestral processes with selection.  Theor. Popln Biol., 51(3):210–237, 1997.
  * [KO94] M. S. Keane and G. L. O’Brien.  A Bernoulli factory.  AMC Trans. Modeling and Computer Simulation, 4:213–219, 1994.
  * [ŁKPR11] K. Łatuszyński, I. Kosmidis, O. Papaspiliopoulos, and G. O. Roberts.  Simulating events of unknown probabilities via reverse-time martingales.  Random Struct. Algorithms, 38(4):441–452, 2011.
  * [LL07] S. Lessard and V. Ladret.  The probabilty of fixation of a single mutant in an exchangeable selection model.  J. Math. Biol., 54:721–744, 2007.
  * [Men19] L. Mendo.  An asymptotically optimal Bernoulli factory for certain functions that can be expressed as power series.  Stoch. Proc. Appl., 129:4366–4384, 2019.
  * [MŁNW22] G. Morina, K. Łatuszyński, P. Nayar, and A. Wendland.  From the Bernoulli factory to a dice enterprise via perfect sampling of Markov chains.  Ann. Appl. Probab., 32(1):327–359, 2022.
  * [Mor21] G. Morina.  Extending the Bernoulli factory to a dice enterprise.  PhD thesis, University of Warwick, 2021.
  * [NP05] Ş. Nacu and Y. Peres.  Fast simulation of new coins from old.  Ann. Appl. Probab., 15:93–115, 2005.
  * [PLS] R. Paes Leme and J. Schneider.  Multiparameter Bernoulli factories.  Preprint, arXiv:2202.07216.
  * [TN06] C. Taylor and M. A. Nowak.  Evolutionary game dynamics with non-uniform interaction rates.  Theor. Popul. Biol., 69:243–252, 2006.
  * [von51] J. von Neumann.  Various techniques used in connection with random digits.  Appl. Math. Ser., 12:36–38, 1951.  Reprinted in _von Neumann’s Collected Works 5_ , Oxford University Press, 1963, pp. 768–770.
  * [XYZ19] F. Xi, G. Yin, and C. Zhu.  Regime-switching jump diffusions with non-Lipschitz coefficients and countably many switching states: existence and uniqueness, Feller, and strong Feller properties.  In G. Yin and Q. Zhang, editors, Modelling, Stochastic Control, Optimization, and Applications. Springer, Cham., 2019.


1]GESIS-Leibniz Institute for Social Sciences, Germany 2]Nokia Bell Labs,
Cambridge, UK 3]CUSP, King’s College, London, UK 4]University College London,
UK 5]King’s College Business School, London, UK<EMAIL_ADDRESS>
# Insider stories: Analyzing Internal Sustainability Efforts of major US
companies from online reviews

Indira Sen  Daniele Quercia  Licia Capra  Matteo Montecchi  Sanja Šćepanović
[  [  [  [  [  [

###### Abstract

It is hard to establish whether a company supports internal sustainability
efforts (ISEs) like gender equality, diversity, and general staff welfare, not
least because of lack of methodologies operationalizing these internal
sustainability practices, and of data honestly documenting such efforts. We
developed and validated a six-dimension framework reflecting Internal
Sustainability Efforts (ISEs), gathered more than 350K employee reviews of 104
major companies across the whole US for the (2008-2020) years, and developed a
deep-learning framework scoring these reviews in terms of the six ISEs.
Commitment to ISEs manifested itself at micro-level — companies scoring high
in ISEs enjoyed high stock growth. This new conceptualization of ISEs offers
both theoretical implications for the literature in corporate sustainability,
and practical implications for companies and policy makers. To further explore
these implications, researchers need to add potentially missing ISEs, to do so
for more companies, and establish the causal relationship between company
success and ISEs.

## Introduction

Figure 1: The Wheel of Internal Sustainability Efforts (ISEs). The wheel
includes the two macro categories (financial benefits _vs._ staff welfare)
under which the six ISEs are classified. The wheel’s outer layer reports
keywords representative of each ISE.

Investments in sustainability are becoming paramount as many companies are
under constant pressure to reduce the social and environmental impact of their
operations [bai2020supply] and increase accountability towards stakeholders
and the wider society [de2021reimagining, serafeim2020social]. While corporate
sustainability efforts tend to focus primarily on external stakeholders (e.g.,
customers, supply-chain partners, governmental organizations)
[doi:10.1177/0022242921992052], internal stakeholders (e.g., employees)
represent a critical, and sometimes overlooked, target group to ensure
effective corporate engagement with the sustainability agenda
[martin2021exploring, paine2014sustainability].

Internal sustainability efforts (ISEs) encompass a wide range of corporate
policies directed towards internal stakeholders, including, for example,
promoting a healthy employee work-life balance [kelliher2019all], investing in
gender equality and diversity [NADEEM2017874], and ensuring an harassment-free
working environment [cassino2019race]. These ISEs can reduce staff turnover
[giauque2019stress] and improve market competitiveness [wang2012explaining].
However, although many companies openly advertise their commitment to internal
sustainability, employees often report contrasting accounts of their
experience of such efforts [peloza2011how], and the extent to which ISEs
successfully propagate throughout the organization remains unclear.

We partly tackle those issues by running a large-scale assessment of
organizational practices aligned with ISEs. Since there exists no agreed-upon
definition of corporate ISEs, to inform our research we started from the
United Nations (UN) World Commission on Environment and Development (WCED)
definition of sustainability as a strategy oriented towards “meeting the needs
of the present without compromising the ability of future generations to meet
their own needs” [wced1987world]. This definition is operationalized in the UN
17 Sustainable Development Goals (SDGs), which represent both a framework and
a call-to-action for organizations to invest in addressing critical societal
issues such as “good health and well-being”, “decent work and economic
growth”, and “peace, justice and strong institutions”
[nations2015transforming]. Not all 17 UN SDGs are relevant to a company’s
internal stakeholders (this is the case, for example, for SDG “life under
water”). To identify the relevant ones and sharpen their definitions in the
internal corporate context, we developed and validated a mixed-method approach
that ended up paraphrasing the broad UN SDGs into six corporate-relevant ISEs.
These ISEs concerned health, education, diversity, monetary benefits,
infrastructure, and atmosphere (Figure 1). Core to the approach is a state-of-
the-art Natural Language Processing (NLP) framework that processed more than
350K geo-referenced reviews about 104 S&P 500 companies.

## Data

Our aim was to understand and capture the microfoundations of ISEs; we did so
in a bottom-up fashion, starting from the perspectives of employees. More
specifically, we collected data from a popular company reviewing site, where
current and, more likely, former employees write reviews about their own
corporate experiences, ranging from job interviews to salaries to workplace
culture. These reviews have been recently used in studies exploring corporate
culture at scale [das2020modeling]. As of 2021, there are 50M monthly visitors
on the platform, and 70M reviews of 1.3M companies. To ensure quality reviews,
the site: a) performs both automatic and manual content moderation; b) allows
for full access to content only to users who register on the site and write at
least one review each (encouraging neutral and unbiased reviews); and c)
allows for posting maximum one review per employee per year. Our dataset
consisted of reviews published over twelve years, from 2008 to 2020.

Each review consists of a title; a ‘pro’ portion (i.e., positive aspects of
the company); a ‘con’ portion (i.e., its negative aspects); a set of four
ratings on a [0,5] scale scoring the company’s _balance_ , _career_ ,
_culture_ , and _management_ ; and a final _overall_ rating of the company.
Since reviewers have the option to include their location, we were able to
identify the states for part of the reviews. To ensure the robustness of our
text processing method, we retained companies that had at least 1,000 reviews
and were present in at least 10 states, leaving us with a dataset of 358,527
reviews of 104 US-based companies (which represented 88.7% of the original
dataset); 80% of these are S&P 500. As detailed in Supplementary Information,
these 104 companies offer the same level of representativeness as the S&P 500
companies, in terms of the distribution of industry sectors and the geographic
distribution across states. In addition to the reviews, we collected yearly
stock growth values of the 104 companies from the Yahoo Finance Portal[WinNT].

## Methods

### The three-step mixed-method approach for defining ISEs

We developed a mixed-method approach to operationalize ISEs. This approach
unfolded in three main steps, which are detailed in Supplementary Information
and summarized here as follows (Figure 2):

Figure 2: Summary of the three-step mixed-method approach for defining ISEs.
Starting with the 17 UN SDGs, three annotators unanimously discarded those
that did not apply to the corporate context (step 1, pre-selection): 13 SDGs
were left. From these, the subset of SDGs accurately captured by our NLP deep-
learning framework was identified (step 2, unsupervised discovery): 8 SDGs
were selected. Finally, three annotators merged the goals that, in the context
of company reviews, ended up being paraphrased with very similar meanings
(step 3, consolidation): this final step resulted in the identification of the
six ISEs.

Step 1 - Pre-selection of goals:

Using a deductive content analysis [elo2008qualitative], three independent
annotators assessed each of the UN seventeen goals’ definitions and decided
whether they applied to the corporate context or not. We took a conservative
approach and discarded the goals that the annotators _unanimously_ discarded,
which ended up being four, leaving us with thirteen potentially relevant goals
(after step 1 in Figure 2).

Step 2 - Unsupervised discovery of goals:

An unsupervised deep-learning framework based on the sentence-level BERT
algorithm [reimers2019sentence] was developed (its technical architecture is
discussed in Supplementary Information). This framework scored each employee’s
review against the 13 goals found in the previous step. The framework
identified the five reviews most relevant to each goal, and three other
independent annotators then manually assessed the relevance of these reviews.
To conservatively retain only the goals that were accurately identified by the
framework, we discarded any goal for which the majority of the annotators
marked less than 4 of the goal’s 5 reviews as relevant (overall, the agreement
among the annotators was high, i.e., Fleiss $K=0.83$). As a result, five goals
were dropped; these had more to do with environmental sustainability (e.g.,
“clean water”, “climate change”) than with internal corporate affairs. This
left us with eight goals (after step 2 in Figure 2).

Step 3 - Consolidation of goals:

Finally, the three annotators assessed if any of the eight goals ended up
acquiring very similar meanings in company reviews. Two pairs were merged,
ultimately leaving us with six ISEs (after step 3 in Figure 2). Table 1
reports the names of these ISEs (first column), corresponding original UN SDGs
(second column), and related excerpts of real reviews (third column).

### Metrics

We studied the six ISEs at the company-level $u$ to test whether commitment to
ISEs manifests itself at a micro-level (e.g., in a company’s growth). To that
end, we computed the score $s(u,i)$ of the $i^{th}$ ISE for company $u$ as the
fraction of $u$’s reviews that mentioned $i$:

$s(u,i)=\frac{\sum_{p\in R(u)}sim_{t}(v_{p},v_{i})}{|R(u)|}$ (1)

where $R(u)$ is the set of $u$’s reviews, $v_{i}$ is the SBERT (Sentence-BERT)
vector of ISE $i$ (the six vectors/phrases for the ISEs are in Supplementary
Information in Table 4), and $sim_{t}(v_{p},v_{i})$ is the _thresholded_ SBERT
similarity score [reimers2019sentence] between the SBERT vector of review $p$
and the SBERT vector of ISE $i$. More precisely, $sim_{t}(v_{p},v_{i})$ is
defined as:

$sim_{t}(v_{p},v_{i})=\begin{cases}sim(v_{p},v_{i}),&\text{if
}sim(v_{p},v_{i})>0.31\textrm{ AND }sim(v_{p},v_{i})\textrm{>}95\%(i)\\\
0,&\text{otherwise}\end{cases}$ (2)

We chose the threshold of $0.31$ by computing the mean SBERT similarity for
each of the 8 goals left after stage 2 of our 3-step ISE selection procedure
as we had established that the NLP method worked well for these 8 goals. We
then took the average value of the eight means (which was $0.31$). Based on
further validation, we also established that the SBERT values for all ISEs
were not equally distributed and, as such, the fixed generalized threshold of
$0.31$ had to be paired with an ISE-specific threshold: based on our
experiments reported in Supplementary Information, this latter threshold value
(denoted as $95\%(i)$) was the 95% percentile of the ISE’s distribution, which
is the very same threshold found in previous studies [choi20ten]. We finally
ranked companies by their score $s(u,i)$ for each $i^{th}$ ISE. Note that, by
review, we mean the pro portion of the review. That is because we were mostly
interested in positive initiatives (pros) rather than shortcomings (cons). In
Supplementary Information, we indeed show that, if we were to instead take
cons (or combine cons with pros together), our deep-learning framework would
perform worse in the two validation steps of our mixed-method approach (steps
2 and 3).

Internal Sustainability Efforts (ISEs) | UN Goal | Example of Review Sentence
---|---|---
| Monetary
---
| decent work and
---
economic growth
| “Professional growth, training, co-workers, mutuality,
---
income, entrepreneurship”.
Health | good health and wellbeing |  | “Excellent work-life balance. Great information
---
offered to improve health and equality of life”.
| Education
---
quality education |  | “Encourage continual education and offer multiple
---
learning opportunities”..
| Diversity
---
gender equality | “Respect for gender equality”.
| Infrastructure
---
| industry, innovation,
---
and infrastructure
“Good infrastructure to support the work environment”.
| Atmosphere
---
| peace, justice, and
---
strong instituions
| “Collaborative environment, excellent benefits,
---
opportunity for growth and development”.
Table 1: The six internal sustainability efforts resulting from the three-step
mixed-method approach for defining ISEs.

## Results

Figure 3: Overview of the research questions investigated in this work and the
sources of data we used for them.

We identified each ISE’s keywords from all reviews associated with it (e.g.,
the keyword ‘salary’ for the ‘monetary’ ISE), and ascertained through a
principled linguistic validation that the keywords are semantically related to
the ISE (RQ1). After establishing that our ISE scoring is valid, we scored the
companies, and studied the relationship between a company’s ISE scores and its
success in the forms of company ratings and stock growth (RQ2), and uncovered
ISE scores variability across industry sectors (RQ3). Figure 3 summarizes our
analyses and the data used for them.

### RQ1: Does our machine learning method capture Internal Sustainability
Efforts?

We validated our deep learning method for detecting ISEs based on a
triangulation approach [denzin2012triangulation], during which we first
established its face validity by inspecting the language used in reviews, and
subsequently examined our results with respect to external reports. We discuss
the former next, while the latter is detailed in Supplementary Information.

To establish the face validity of the proposed ISE detection method, we took
the linguistic approach explored by [das2020modeling]. First, for each of the
six ISEs, we obtained the most frequent keywords — 1, 2, 3, and 4-grams from
the reviews deemed relevant by our method. We then computed the TF-IDF scores
for such n-grams, where each document was comprised of all shortlisted reviews
for each ISE. Finally, we ranked keywords for each ISE based on their TF-IDF
score. This allowed us to find the keywords judged to be important for a
certain ISE by our embedding-based method. The top ranked keywords for the six
ISEs are visualized as a heatmap in Figure 4.

Figure 4: Top n-grams in sentences expressing ISEs. Darker colors (higher
normalized TF-IDF score) indicate greater relative relevance to a particular
ISE.

We observed many keywords to be highly discriminative of the ISE they
associated with: for example, keywords ‘pay good’ and ‘salary’ were
(correctly) ranked highly for ISE ‘monetary’ only; ‘health’, ‘health
benefits’, and ‘take care’ were ranked highly for ISE ‘health’ instead.
Keywords ‘opportunities learn’, ‘experience’, ‘good train’, and ‘program’ were
uniquely strongly associated with ISE ‘education’. Keyword ‘flexibility’ was
highly discriminative of ISE ‘diversity’; ‘industry’ and ‘technology’ were
strongly associated with ISE ‘infrastructure’; lastly, n-grams like ‘positive
work environment’ and ‘friendly work environment’ strongly associated with ISE
‘atmosphere’. Other keywords ranked highly in more than one ISE instead: this
was the case, for example, for keywords ‘benefit worklife balance’ and
‘family’, which were highly associated with both the ‘health’ ISE (as one
might expect), and to the ‘diversity’ ISE. Health-enhancing factors like work-
life balance and flexible working conditions options have been shown to
facilitate gender equality and improve the diversity of employees
[chung2018flexible, lyonette2015part], therefore it was promising that our
method was capable of picking up these semantically-related concepts too.

Indeed the 6 ISEs we identified were not mutually exclusive concerns (and
neither are the UN SDGs), and one may wonder to what extent they are
semantically related. To shed light onto this question, we conducted a
Principal Component Analysis (PCA) on s(u, i) at a company level to assess how
much of the variance in the data could be explained by different principal
components, and how those components related to the 6 ISEs. We found that, at
company level, just two components explained 88% of the variance —
specifically, the first component explained 73% and the second component
explained 15%. We report the correlation between the first two PCA components
and the six ISEs in the last two columns of Table 2.

We observed that all ISEs with the exception of ‘monetary’ were strongly
correlated with the first component and weakly negatively correlated with the
second component; on the other hand, ‘monetary’ was moderately correlated with
both the first and second principal components. These two findings suggested
that the ‘monetary’ ISE was orthogonal to the other five, and that these other
five were strongly interconnected with one another. Indeed, one may expect
that improving work-life balance has a positive impact on both the ‘health’
ISE and the ‘diversity’ ISE; on the other hand, improving monetary conditions
may not directly affect other aspects of corporate internal sustainability.
Overall, we thus found two main facets of employee-centred sustainability — a
_staff welfare_ -related one ($PC_{1}$) and a _financial benefits_ -related
one ($PC_{2}$). To avoid multicollinearity, we used these two main facets of
ISEs (rather than the six individual ones) to answer the following research
questions.

| Monetary | Health | Education | Diversity | Infrastructure | Atmosphere |  | Staff Welfare
---
(PC1)
| Financial Benefits
---
(PC2)
Monetary | 1.00 | 0.74 | 0.39 | 0.43 | 0.41 | 0.55 | 0.67 | 0.68
Health | 0.74 | 1.00 | 0.55 | 0.76 | 0.74 | 0.85 | 0.90 | 0.32
Education | 0.39 | 0.55 | 1.00 | 0.72 | 0.77 | 0.73 | 0.82 | -0.39
Diversity | 0.43 | 0.76 | 0.72 | 1.00 | 0.82 | 0.92 | 0.90 | -0.19
Infrastructure | 0.41 | 0.74 | 0.77 | 0.82 | 1.00 | 0.92 | 0.90 | -0.25
Atmosphere | 0.55 | 0.85 | 0.73 | 0.92 | 0.92 | 1.00 | 0.96 | -0.07
Staff Welfare (PC1) | 0.67 | 0.90 | 0.82 | 0.90 | 0.90 | 0.96 | 1.00 | 0.00
Financial Benefits (PC2) | 0.68 | 0.32 | -0.39 | -0.19 | -0.25 | -0.07 | 0.00 | 1.00
Table 2: Cross-correlation between the 6 ISE scores and the 2 principle
components obtained via PCA at a company level.

### RQ2: Is sustainability associated with company success?

There are several ways to measure a company’s success. We considered two
complementary ones: the online ratings it received from its employees
(available from the company reviewing site), and its financial position
(measured as stock growth).

Sustainability and company online ratings. Employees have the option to rate
the company they are reviewing based on four different facets — balance,
career, culture, management, plus a fifth company’s overall one. We thus
investigated to what extent a company’s success across these five facets could
be predicted based on the company’s commitments to the ISEs. We did so by
first aggregating ISE scores and ratings at company level, then conducting an
OLS regression using our two main sustainability facets as predictors (‘staff
welfare’ and ‘financial benefits’), while also controlling for a company’s
total number of reviews. As reported in Table 3, we found that these two
sustainability facets could explain up to 64% of the variance in a company’s
ratings; particularly noteworthy was that the staff welfare facet of corporate
internal sustainability was strongly positively correlated with all aspects of
a company’s success, including balance and culture, in line with previous
research findings [rao2017work, isensee2020relationship].

| Balance | Career | Culture | Management | Overall
---|---|---|---|---|---
Const |  16.465∗∗ |  29.995∗∗∗ |  30.165∗∗∗ |  24.863∗∗∗ | 5.651
| (6.286) | (4.723) | (5.990) | (5.614) | (4.740)
Staff Welfare (PC1) |  0.732∗∗∗ |  0.681∗∗∗ |  0.721∗∗∗ |  0.615∗∗∗ |  0.766∗∗∗
| (0.095) | (0.071) | (0.090) | (0.085) | (0.067)
Financial Benefits (PC2) |  0.183∗ |  -0.164∗∗ | -0.160 |  -0.201∗∗ |  0.313∗∗∗
| (0.109) | (0.082) | (0.104) | (0.097) | (0.076)
Total Reviews |  |  |  |  | 0.102
|  |  |  |  | (0.063)
Observations | 84 | 84 | 84 | 84 | 84
$R^{2}$ | 0.443 | 0.534 | 0.444 | 0.405 | 0.658
Adjusted $R^{2}$ | 0.429 | 0.523 | 0.430 | 0.390 | 0.645
Residual Std. Error | 16.283(df = 81) | 12.234(df = 81) | 15.515(df = 81) | 14.541(df = 81) | 11.355(df = 80)
F Statistic |  32.174∗∗∗ |  46.461∗∗∗ |  32.289∗∗∗ |  27.518∗∗∗ |  51.319∗∗∗
| (df = 2.0; 81.0) | (df = 2.0; 81.0) | (df = 2.0; 81.0) | (df = 2.0; 81.0) | (df = 3.0; 80.0)
Note: |  ∗p$<$0.1; ∗∗p$<$0.05; ∗∗∗p$<$0.01
Table 3: Predicting company online ratings from the two main facets of
sustainability (staff welfare and financial benefits) using a stepAIC analysis
on an OLS regression.

Sustainability and company stock growth. We obtained stock data for 84 of the
104 companies in our dataset, from 2009 to 2019, using the Yahoo Finance
portal [WinNT]. For each company, we calculated the geometric mean of its
stock growth during such period; we used the geometric mean since the
distribution of stock growth values across companies was heavy-tailed (as
reported in Supplementary Information). To inspect whether a company’s
financial success (measured as stock growth) was associated with its
sustainability efforts, we then plotted in Figure 5 the geometric mean of its
stock growth ($y$ axis) against its ranking in terms of the staff welfare
facet of sustainability and the financial benefits facet of sustainability
($x$ axis). We also included in the figure the total number of reviews, to
check whether stock growth was merely associated with the company’s popularity
rather than its internal sustainability practices.

As showed in Figure 5, companies that focused on both staff welfare and
financial benefits sustainability tended to have high stock growth; between
the two facets, it was staff welfare that most strongly associated with high
stock growth, in line with previous research [bcg]. Notably, companies with
high stock growth did not invest as heavily in financial sustainability only,
bolstering previous work which noted that focusing on staff welfare
sustainability could lead to greater stakeholder engagement even without high
pay [ziegler2007effect]. Overall, our results suggest that a company’s
financial success is associated with its investment in internal sustainability
practices, but only if they focus on a holistic approach to sustainability
that tackles both staff welfare and financial benfits.

Figure 5: Geometric mean of stock growth values for increasing ISE score
ranking.

### RQ3: Is sustainability associated with specific industry sectors?

Figure 6: Sustainability and industry sector. Boxplots showing the
distribution of the staff welfare and financial benefits ISE scores across
different industry sectors.

To examine whether certain industry sectors were leading the corporate
sustainability agenda, we plotted in Figure 6 the distribution of the two
facets of sustainability for each industry sector. We further conducted a
MANOVA analysis and found the differences in sustainability scores to be
significant across sectors. In terms of staff welfare sustainability efforts,
we found Industrials and IT to lead, possibly due to recent investment in this
type of sustainability initiatives [higon2017ict]. The Financial sector
followed, while the health care one exhibited very high variability. This
could be explained by health care professionals often sacrificing personal
well-being and work-life balance due to the highly demanding nature of their
work [schwartz2019work, shanafelt2015changes]. We found consumer staples and
consumer discretionary to lag significantly behind. This was also the case
when looking at financial benefits, although differences between sectors were
smaller along this facet of internal sustainability efforts.

Figure 7: Scatterplot of the scores of each company’s staff welfare _vs._
financial benefits. The size of a company’s dot represents its stock growth.
We highlighted in blue some of the companies to assess them qualitatively.
Consumer staples and discretionary companies like Kmart, Macy’s, and Kohl’s
scored low for both types of sustainability. Traditional IT companies like
Infosys, IBM, and Accenture scored high for staff welfare sustainability but
not for financial benefits sustainability.

Figure 6 offered an overview of engagement with sustainability efforts at
industry sector level. To reveal more nuanced variations within the same
sector, we plotted individual companies’ engagement with each of the two main
sustainability facets in Figure 7. Notable variations emerged: our sector-
based analysis revealed low sustainability scores for consumer discretionary
and staples companies overall; upon closer inspection, we found some companies
(e.g., Dollar General, K-mart) to indeed score low on both facets of ISEs,
while others (e.g., Costco) to score low on staff welfare but high in the
financial benefits facet of sustainability, a phenomenon noted in previous
work too [cascio2006decency]. Variations emerged also within the IT sector,
previously showed to be leading sustainability efforts on both dimensions: a
more nuanced investigation revealed high sustainability scores on both
financial benefits and staff welfare ISEs for companies like Microsoft,
Google, and Apple; however, more traditional IT companies like Infosys, IBM,
and Cognizant scored high on staff welfare sustainability only.

One must be mindful that comparisons among (same-sector) companies were
further affected by the type of employees reviewing their employer. In our
study, this was apparent for companies like Amazon, that enjoyed high stock
growth but surprisingly scored low for both types of sustainability. Despite
the company employing a large number of software engineers as well as
warehouse workers, upon close inspection, we found the most common roles of
Amazon employees in our reviews to be ‘warehouse associate’ and ‘warehouse
worker’. Previous research did find logistic workers at Amazon to face poor
working conditions [amzn, chan2015examining], corroborating the low
sustainability scores that our method computed for this company. Furthermore,
previous literature noted that Amazon’s lack of focus on sustainability
practices has yet to hurt its profitability [amzn, chan2015examining].

## Discussion

Many companies are under constant pressure to invest in a wide range of
internal sustainability practices designed to enhance working conditions
[barko2022shareholder, jakob2022like]. However, the benefits of such
investments for both the company and its stakeholders are often difficult to
assess. By examining how employees form perceptions of their company’s
engagement with ISEs, this research spells out the microfoundations of
internal sustainability and provides evidence of the strategic importance of
investing in business practices and policies geared towards ISEs
[de2021reimagining, zhao2022influence].

By examining how the wider UN SDGs agenda can be translated into diverse
internal corporate efforts directed towards employees, our work offers
substantive methodological, conceptual, and empirical contributions to
internal sustainability research and managerial practice. More specifically,
it offers two main theoretical contributions. The first has to do with the
conceptualization of ISEs. We have shown how the sustainability agenda brought
forward by the introduction of the UN SDGs informs and shapes six
sustainability efforts within a company. Efforts to do with health, education,
diversity, monetary benefits, supporting infrastructure, and a supportive
atmosphere. While the existing literature often presents sustainability as a
monolithic construct [chen2020corporate, liu2020too], our two-factor
conceptualization of ISEs delineated the two core strategic aspects that
companies should carefully balance when implementing ISEs: one aspect had to
do with traditional financial benefits (e.g., salary, bonuses), and the other
had to do with broader aspects of staff welfare (e.g., diversity, atmosphere).
The second theoretical implication enhances the understanding of what makes
companies economically successful and how internal sustainability practices
differ by sector, especially in emerging sectors like IT.

This work also offers practical implications, and it does so for three main
stakeholders. The first stakeholder consists of scholars. Our method is
grounded in the UN SDGs and performed consistently well across several rounds
of external validation. By providing a robust framework for examining mentions
of ISEs through automated text analysis, new textual datasets could be
academically studied in the future.

The second stakeholder consists of policy makers. We showed that high levels
of ISEs engagement (not only for financial aspects but also for general staff
welfare) were associated with high economic growth. This result supports
policies in recent years that have fostered a corporate culture that goes
beyond financial rewards and are oriented towards equality and well-being
[triana2019perceived]. Beyond company efforts, policy makers themselves would
be able to strategically decide which ISEs to incentivize with taxation
schemes or set out a legislation agenda that would attract workers who care
about specific ISEs. To inform more targeted interventions, we also showed
that the impact of engaging with ISEs varies across sectors: companies in the
IT and business-to-business industrial goods sectors outperformed companies
that produce and commercialize consumer goods. This finding is noteworthy as
previous research shows that sustainability signals tend to be stronger in
business-to-consumer than in business-to-business market contexts (c.f.,
[hoejmose2012green])

The third stakeholder consists of company managers. By reflecting employees’
perceptions, our analytical framework represents an invaluable tool to
operationalize the microfoundations of internal sustainability, assess how
corporate efforts in this area directly impact employees, and quantify and
qualify the extent to which corporate engagement with ISEs becomes visible to
employees across different organizational levels.

Our work comes with five main limitations though. The first is that our list
of ISEs may not be accurate or exhaustive. While the corporate sustainability
literature has focused on initiatives that are external to a company and have
an impact on the wider world’s sustainability, the practices that are internal
to a company and have an impact on employees received less attention. As a
result we only found non-comprehensive frameworks for internal sustainability
practices suggested, such as those focusing on social aspects only
[baumgartner2010corporate]. To tackle that, we started from the well-grounded
definitions of the UN sustainability goals, used a principled mixed-method
approach to paraphrase those most relevant to the corporate context, and
validated the resulting list with both qualitative and quantitative
approaches. These approaches are generalizable, in that they could be used to
study other constructs appearing in reviews in the future (e.g., how employees
in a company deal with stress).

The second limitation is that, since the reviewing site was founded in 2008,
key financial events prior 2008 (e.g., the dot-com bubble in the late 1990s)
may have impacted our results but could not be accounted for because of lack
of data.

The third limitation is that the number of companies under study is invariably
limited. We were able to study 104 major companies, largely because the other
companies had a limited number of reviews that did not allow for automatic
processing. Future work should explore alternative mixed-method approaches
(likely qualitative ones) to study ISEs for these companies.

The fourth limitation is the lack of causal claims. Given our data, we could
not assess the causal direction between ISEs and socio-economic returns. More
specifically, we could not assess whether focusing on ISEs led to better
socio-economic returns (e.g., stock growth), whether better socio-economic
conditions created a breeding ground for fostering ISEs, or whether these two
causal relations were in a self-reinforcing cycle.

The fifth and final limitation has to do with the representativeness of our
data. Companies in certain sectors (e.g., IT) may have been reviewed more
often than those in other sectors (e.g., consumer discretionary). Despite that
possibility, in Supplementary Information, we show that our data was still
representative along three major dimensions: a) the distribution of industry
sectors of the S&P 500 companies, which our data matched without over-
representing any specific sector; b) official population in a state, which
scaled linearly with the number of employees in the state in our data; and c)
number of company headquarters in a state from official sources, which has a
nearly perfect correlation with the number of headquarters per state in our
data. Finally, despite the platform’s mechanisms to guarantee review quality,
as discussed in the section _Data_ , we acknowledge that potential self-
selection bias could cause our reviewers’ sample to be non-representative. To
reduce the impact of such a bias and ensure robustness, we restricted our
analyses to companies having at least 1000 reviews.

## Data availability

We made our code and data available in a readily usable format (http://social-
dynamics.net/InsiderStories/) for reproducibility. For each company, we shared
the following attributes: company name, #total reviews, #ISE reviews, average
rating, rating of work-life balance, rating of career prospects, rating of the
company, rating of the culture, rating of the management, ISE type, stock
values/growth for: 2009, 2012, 2014, 2019, and industry sector.

## Ethical approval and informed consent statement

This article does not contain any studies with human participants performed by
any of the authors.

## Contributions of the authors

All authors contributed to the research design and writing of the paper. IS
was mainly responsible for the linguistics and statistical analysis with
guidance from SS, LC, and DQ. SS cleaned, collected and pre-processed the
data. MM provided theoretical guidance. All authors wrote, read, and approved
the final manuscript.

## Competing interests

The authors declare no competing interests.

## Acknowledgements

We thank Dr. Edyta Bogucka for her help with some of the figures in this
manuscript.