# A note on efficient audit sample selection

Laura Boeschoten, Sander Scholtus & Arnout van Delden

###### Abstract

Auditing is a widely used method for quality improvement, and many guidelines
are available advising on how to draw samples for auditing. However,
researchers or auditors sometimes find themselves in situations that are not
straightforward and the standard sampling techniques are not sufficient, for
example when a selective sample has initially been audited and the auditor
desires to re-use as many cases as possible from this initial audit in a new
audit sample that is representative with respect to some background
characteristics. In this paper, we introduce a method that selects an audit
sample that re-uses initially audited cases by considering the selection of a
representative audit sample as a constrained minimization problem. In
addition, we evaluate the performance of this method by means of a simulation
study and we apply the method to draw an audit sample of establishments to
evaluate the quality of an establishment registry used to produce statistics
on energy consumption per type of economic activity.

Keywords: Audit sample, validation sample, gold standard, deviance,
minimization procedure

##  1 Introduction

Auditing is a widely used method for quality improvement. The goal of an audit
is to estimate the extent to which practice in a field of profession complies
with certain identified review criteria. The degree of compliance helps in
identifying areas where improvements can be made [9]. As practical reasons
often limit the ability to apply an audit to the total population, typically a
sample is drawn on which the audit is performed, the _audit sample_. The
procedure of performing an audit is known under different denominators. An
audit can be performed in the sense that a sample is thoroughly evaluated, or
a high quality test can be applied to the audit (which is then often known as
the _gold standard_). The study of comparing (potentially) error-prone
outcomes to true scores is also known as a _validation study_ and the set of
_true scores_ is then known as a _validation sample_. For supervised machine
learning a manually annotated set is created also referred to as a _labelled
set_. Audits are performed in different fields of profession.

In the field of clinical research, audits can for example be used to evaluate
the performance of diagnostics tests, as applied by [4]. In such cases, an
audit is used to evaluate the performance of a certain test that is used to
diagnose a specific type of disease. This typically means that a diagnostic
test is available and of high quality (‘gold standard’). However, it is often
expensive, difficult to apply or invasive and therefore has a high patient
burden. An alternative, cheaper, more easy to use, lower patient burden
diagnostic test is also available, although it has a lower quality in the
sense that more outcomes are either false positives or false negatives
compared to the gold standard test. To evaluate the quality of the new
diagnostic test, both tests are applied to the audit sample and the outcomes
of the new test are compared to outcomes of the ‘gold standard’ test.

In many organizations, financial audits are performed, sometimes because it is
required by law [17]. Here, the extent to which an organization’s financial
statements are presented fairly is evaluated [6]. Within an organization, all
financial transactions can be considered a population, from which a sample is
drawn, sometimes by means of a risk function, from which the total amount of
financial misstatements in the organization is inferred [7].

In official statistics, audits are also being used to estimate the quality of
administrative sources or of surveys. This can for example be done
simultaneously by using a structural equation model as shown by [15] including
an audit sample [16]. Another field where the use of an audit sample plays a
role is in supervised machine learning. There, a manually annotated set with
(error-free) labels is used to train and test a machine learning model. The
test set is used to estimate the performance of the model for the intended
population. Therefore, for the test set it is important that it is
representative for the target population. When one is interested in predicting
a class with a small proportion in the population, then using a simple random
sample to obtain an annotated set is not so efficient. One may for instance
use a set of key words to create a selective sample. The question then is how,
starting from a part of the selective sample, it can be supplemented to obtain
a representative test set. In the field of machine learning, this issue is
also known as ‘sample selection bias’ [20] and is for example manifested when
a road network is inferred from Weigh-in-Motion road sensors [11].

As an audit is often too expensive and/or invasive to be performed on the
complete population, typically a sample is selected on which the audit is
performed. Ideally, this sample is selected at random from the complete
population directly or from a Randomized Controlled Trial (RCT). In such
cases, it can be assumed that the audit sample is representative for the
population. In practice however, auditors often face situations that are more
complex. Examples are when an observational study is performed instead of an
RCT, when the auditor prefers to draw a judgement sample or when the register
that is available as a sampling frame is incomplete. In such cases, the sample
should be selected more carefully in order for it to be representative.
Representativity can for example be achieved by investigating background
characteristics of the sample units and by drawing a stratified sample.
Typically, guidelines are also available for drawing stratified samples, such
as provided by the British Medical Journal [12] or the American Institute of
Certified Public Accountants [2].

In practice however, auditors find themselves often in situations that are
more complex and for which the standard guidelines are not sufficient. A first
example of such a situation is when the auditor has determined a maximum
number of units for which auditing is feasible or when the auditor has
limitations in terms of budget. A second example of a complex situation is
when a small non-random sub-sample of the population has already been audited
which should ideally be included in the audit sample or when the auditor plans
to include particular units in the audit sample with certainty. In such
situations, using standard guidelines for auditing are not sufficient.

In this paper, we propose to approach the selection of an audit sample as a
constrained minimization problem. In the proposed framework, a representative
audit sample can be straightforwardly selected for the non-standard situations
described above. In Section 2, notation and the proposed framework are
introduced. In Section 3, it is illustrated how the framework can be used for
the complex audit situations. In Section 4, a simulation study is used to
investigate the performance of the method under different conditions. Section
5 illustrates how the proposed method can be applied in practice using ready-
made R scripts. Finally, section 6 concludes.

##  2 Background

###  2.1 A representative sample

When an audit is performed, the goal is to select a sample that is
representative with respect to the target population. We assume that this
audit sample is drawn from a subset of units in the target population for
which a possibly error-prone measure of the outcome variable of interest is
available. We will refer to this subset as the _observed sample_. In some
applications, the observed sample is in fact equal to the target population.
In general, the observed sample is a sample from the target population, and
standard sampling theory applies to both the audit sample with respect to the
observed sample and the observed sample with respect to the target population.

When an audit is performed for a diagnostic test in health care, the observed
sample refers to the sample on which an inferior test is applied. When an
algorithm is validated, the observed sample refers to the set on which the
prediction was performed. In official statistics, the observed sample can
refer to an administrative data set that is used to compile statistics; in
such cases the observed sample and target population are often assumed to be
the same. For all above scenarios it holds that both the observed sample and
the audit sample should be representative with respect to the population to
which inference is desired.

In what follows, all variables are assumed to be categorical. The observed
sample is denoted as $U$, with $U={1,2,\ldots,N}$. The true variable of
interest is denoted by $W$, with $W_{g}$ denoting the value of this variable
for element $g$ in $U$. The aim is to estimate certain parameters of the
distribution of $W$, for instance, its frequency distribution in the target
population. It is assumed that $W$ is observed _only_ for units that are
audited. For all units in the observed sample $U$, an error-prone version of
the variable of interest is measured, which is denoted by $X$. The values of
$X$ in the observed sample are denoted by $X_{1},X_{2},\ldots,X_{N}$. In
addition to estimating the distribution of $W$, a secondary aim may be to
estimate the association between $W$ and $X$, e.g., the error probabilities
$\text{Pr}(X_{g}=x\>|\>W_{g}=w)$. Knowledge of these error probabilities makes
it possible to estimate the distribution of $W$ based on the whole observed
sample, not just the audit sample, which is more efficient. In addition, these
error probabilities may be of interest in their own right, as quality measures
of the error-prone indicator in $X$.

In addition to $X$, we suppose that one or more background variables are
available in the observed sample, collectively denoted as $Y$, with values
$Y_{1},Y_{2},\ldots,Y_{N}$. We assume that the joint distribution of $Y$ is
known (or previously estimated) for the target population. The content of the
variables in $Y$ depends on the specific audit study of interest. In case of a
diagnostic test in healthcare, $Y$ can contain variables such as age, gender,
ethnic background and prior health status. In case of a financial audit, $Y$
can contain variables such as employer ID, monetary size or department in
firm.

Let $Z$ denote the selection indicator of the audit sample within $U$, where
$Z_{g}$ is the value for element $g$ with $g=1,2,\ldots,N$. $Z_{g}=1$ when the
element is included in the sample and $Z_{g}=0$ otherwise. The expected value
of $Z_{g}$ with respect to the sample design is denoted by

$\pi_{g}=E(Z_{g})=\text{Pr}(Z_{g}=1).$ (1)

The audit sample size, denoted by $n$, is equal to the sum of the values of
$Z$:

$n=\sum_{g=1}^{N}Z_{g}.$ (2)

In order to make appropriate inferences regarding the target population, it is
important that the audit sample is representative with respect to the target
population. If the observed sample was randomly drawn from that population and
is representative, the audit sample can also be randomly selected from the
observed sample by giving every unit an equal inclusion probability
($\pi_{g}=n/N$).

In practice, situations may occur in which the inclusion probabilities should
be allowed to depend on $Y$, for example because some type of elements are
more likely or are allowed to decline for audit, because a judgement sample
was selected, or because the observed sample was not representative with
respect to the target population, for example because some groups of elements
were so small that they were over-sampled in the observed sample. In such (in
practice very likely) cases a simple random sample without replacement for the
audit sample will not be sufficient to obtain representativity.

We will now give a more precise definition of representativity as it will be
used in this paper. As a theoretical starting point to investigate the
selectivity of an audit sample, the joint distribution in the observed sample
of the variables $(W,X,Y,Z)$ is of interest, where $W$ and $X$ represent the
true variable of interest and its error-prone observed version, $Z$ indicates
whether a unit is included in the audit sample and $Y$ represents a grouping
variable based on background characteristics. We make the following
simplifying assumption:

Assumption A. There is no _direct_ association between the true variable of
interest $W$ and the audit inclusion indicator $Z$, once $X$ and $Y$ are
accounted for.

For any audit sample that we draw ourselves, we can ensure that this
assumption is satisfied. Alternatively, if we are given a possibly non-random
sample of previously audited units, the assumption may require that the right
background variables are included in $Y$. The importance of this assumption
which will be seen below is that it allows the selectivity of the audit sample
to be investigated by analysing the joint distribution of $(X,Y,Z)$ instead of
$(W,X,Y,Z)$ (see the end of Section 2.2). In practice, the former distribution
is known, whereas inference about the latter distribution becomes possible
only once a representative audit sample has been obtained.

Selectivity in whether or not a unit is included in the audit sample (_audit
inclusion_), meaning that the distribution of $Z$ depends on $X$, is not
problematic as long as differences in inclusion probabilities are completely
related to the background characteristics in $Y$ for which the distribution in
the target population is known:

$\text{Pr}(Z_{g}=1\>|\>X_{g}=x,Y_{g}=y)=\text{Pr}(Z_{g}=1\>|\>Y_{g}=y).$ (3)

From a missing data perspective, it can be said that audit inclusion is not
problematic as long as exclusion is Missing At Random (MAR) [14]. However, if
the distribution of $Z$ depends on $X$ and this dependence is not explained by
$Y$, this can be problematic, in particular when the error probabilities
$\text{Pr}(X_{g}=x\>|\>W_{g}=w)$ are of interest. From a missing data
perspective, it can be said that exclusion is then Missing Not At Random
(MNAR) [14]. Therefore, our goal is to select an audit sample in such a way
that (3) holds; we refer to this as the _MAR selectivity requirement_.

###  2.2 Deviance as a criterion for representativity

To determine whether a given audit sample meets the MAR selectivity
requirement (3), we propose to analyse the joint distribution $(X,Y,Z)$ in the
observed sample using a non-saturated log-linear model $(XY)(YZ)$. Note that
the difference between the non-saturated model and a saturated model is that
the direct association $(XZ)$ and the three-way interaction $(XYZ)$ are
excluded. In the non-saturated model it is assumed that $X$ and $Z$ are
conditionally independent given $Y$, i.e. the observed outcome variable and
the audit inclusion probability are conditionally independent given the
background characteristics under consideration. Therefore, we refer to it as
the _independence model_.

The number of observations in cell $(X=i,Y=j,Z=k)$ is denoted as $n_{ijk}$ and
the number of observations estimated by the independence model is denoted as
$\hat{n}_{ijk}$. The fit of the independence model can be measured by the
likelihood ratio test statistic or deviance ($D$) comparing it to the
saturated model where $\hat{n}_{ijk}=n_{ijk}$ [1]:

$\displaystyle D=2\sum_{i,j,k}n_{ijk}\textrm{ log
}n_{ijk}-2\sum_{i,j,k}{n}_{ijk}\textrm{ log }\hat{n}_{ijk},$ (4)

where $D\geq 0$. A higher value for $D$ indicates a stronger conditional
dependence between $X$ and $Z$ given $Y$ in the observed data-set, which means
that the selectivity problem is more substantive. This suggests that an audit
sample should be selected for which $D$ is ‘sufficiently small’.

For the independence model, the estimated value for $\hat{n}_{ijk}$ can be
calculated directly, without the use of an iterative algorithm, as

$\hat{n}_{ijk}=\frac{n_{ij+}n_{+jk}}{n_{+j+}},$ (5)

see for example [1] or [3]. Here, $n_{ij+}=\sum_{k}n_{ijk}$,
$n_{+jk}=\sum_{i}n_{ijk}$, and $n_{+j+}=\sum_{i,k}n_{ijk}$.

The definition of $D$ can be expressed alternatively by incorporating (5):

$\displaystyle D$ $\displaystyle=2\sum_{i,j,k}n_{ijk}\textrm{ log
}n_{ijk}-2\sum_{i,j}n_{ij+}\textrm{ log }n_{ij+}$ (6)
$\displaystyle\quad-2\sum_{j,k}n_{+jk}\textrm{ log
}n_{+jk}+2\sum_{j}n_{+j+}\textrm{ log }n_{+j+}$
$\displaystyle=C+2\sum_{i,j,k}n_{ijk}\textrm{ log
}n_{ijk}-2\sum_{j,k}n_{+jk}\textrm{ log }n_{+jk},$

where

$C=2\sum_{j}n_{+j+}\textrm{ log }n_{+j+}-2\sum_{i,j}n_{ij+}\textrm{ log
}n_{ij+}$ (7)

is a constant term which depends only on the distribution of $(X,Y)$ and
therefore will be the same for every possible choice of audit sample. For the
second term in the first line of Equation (6) we used that
$\sum_{i,j,k}n_{ijk}\textrm{ log }n_{ij+}=\sum_{i,j}n_{ij+}\textrm{ log
}n_{ij+}$; for the third and fourth term alike.

Recall from Section 2.1 that the full (unobserved) distribution of interest is
$(W,X,Y,Z)$. We will now show that under Assumption A (there is no direct
association between $W$ and $Z$), the selectivity of the audit sample can be
analysed using the deviance defined in Equation (6), which is based on the
observed distribution $(X,Y,Z)$.

Let $n^{\prime}_{hijk}$ denote the number of observations in cell
$(W=h,X=i,Y=j,Z=k)$, with
$n_{ijk}=n^{\prime}_{+ijk}=\sum_{h}n^{\prime}_{hijk}$. Under Assumption A, the
maximal hierarchical log linear model for $n^{\prime}_{hijk}$ is $(WXY)(XYZ)$.
For this model, the predicted values $\hat{n}^{\prime}_{hijk}$ can also be
obtained directly:

$\hat{n}^{\prime}_{hijk}=\frac{n^{\prime}_{hij+}n^{\prime}_{+ijk}}{n^{\prime}_{+ij+}}=\frac{n^{\prime}_{hij+}n_{ijk}}{n_{ij+}},$

see [3].

Hence, the deviance that compares this model to the saturated model with
$\hat{n}^{\prime}_{hijk}=n^{\prime}_{hijk}$ can be written analogously to
Equation (6) as

$\displaystyle D_{W,max}$
$\displaystyle=2\sum_{h,i,j,k}n^{\prime}_{hijk}\textrm{ log
}n^{\prime}_{hijk}-2\sum_{h,i,j}n^{\prime}_{hij+}\textrm{ log
}n^{\prime}_{hij+}$ $\displaystyle\quad-2\sum_{i,j,k}n_{ijk}\textrm{ log
}n_{ijk}+2\sum_{i,j}n_{ij+}\textrm{ log }n_{ij+}.$

The maximal model contains a direct association between $X$ and $Z$, pointing
to selectivity in the audit sample. The analogue of the independence model for
$(W,X,Y,Z)$ is given by the log linear model $(WXY)(YZ)$. For this model, the
predicted values satisfy:

$\hat{n}^{\prime}_{hijk}=\frac{n^{\prime}_{hij+}n^{\prime}_{++jk}}{n^{\prime}_{++j+}}=\frac{n^{\prime}_{hij+}n_{+jk}}{n_{+j+}},$

see [3].

Hence, the deviance that compares this independence model to the saturated
model for $(W,X,Y,Z)$ can be written as

$\displaystyle D_{W}$ $\displaystyle=2\sum_{h,i,j,k}n^{\prime}_{hijk}\textrm{
log }n^{\prime}_{hijk}-2\sum_{h,i,j}n^{\prime}_{hij+}\textrm{ log
}n^{\prime}_{hij+}$ $\displaystyle\quad-2\sum_{j,k}n_{+jk}\textrm{ log
}n_{+jk}+2\sum_{j}n_{+j+}\textrm{ log }n_{+j+}.$

However, under Assumption A the saturated model is too large, and it makes
more sense to compare the fit of the independence model to that of the above
maximal model. This can be done using the conditional deviance for nested
models [3], which in this case is given by:

$\displaystyle D_{W}-D_{W,max}$ $\displaystyle=2\sum_{i,j,k}n_{ijk}\textrm{
log }n_{ijk}-2\sum_{i,j}n_{ij+}\textrm{ log }n_{ij+}$
$\displaystyle\quad-2\sum_{j,k}n_{+jk}\textrm{ log
}n_{+jk}+2\sum_{j}n_{+j+}\textrm{ log }n_{+j+}.$

Clearly, $D_{W}-D_{W,max}=D$ in Equation (6). Hence, under assumption A, to
analyse the fit of the independence model compared to the maximal model for
the full distribution $(W,X,Y,Z)$, it suffices to analyse the fit of the model
$(XY)(YZ)$ for the observed distribution $(X,Y,Z)$ .

##  3 Method

In this section, we introduce how a representative audit sample can be
selected by using the deviance defined in Equation (6) as a criterion for
representativity. The proposed method has two main advantages. First, the
auditor can determine a maximum number of units to be audited. Here, the
proposed method will provide an audit sample that is as representative as
possible given that maximum number. Second, if some units have already been
audited or the auditor is planning to include specific units in the audit
sample, the method attempts to include as many of these units into the final
audit sample as possible, given the constraint that this final sample should
be sufficiently representative.

###  3.1 Basic optimization problem

Let $n_{ijk}$ denote the number of observations in cell $(X=i,Y=j,Z=k)$ prior
to applying the method in this section. If previously audited units are
available, these have $Z=1$ and so it holds that $n_{ij1}>0$ for certain
units. If the selection of the audit sample starts with a ‘blank canvas’, then
$n_{ij1}=0$ and $n_{ij0}=n_{ij+}$ for all $i$ and $j$. We consider the general
situation where additional units may be selected for auditing (moved from
$Z=0$ to $Z=1$) and previously selected units may be excluded (moved from
$Z=1$ to $Z=0$). When applying the method, the user should specify a maximum
number of additional units to include in the audit sample ($M_{+}$) and a
maximum number of previously audited units to exclude ($M_{-}$). Special cases
are obtained by setting one of these bounds to zero. Moreover, in the special
case of no previously audited units, it automatically holds that $M_{-}=0$ and
$M_{+}$ indicates the maximal size of the audit sample.

After applying the method, the adjusted number of units in cell
$(X=i,Y=j,Z=k)$ is denoted as $m_{ijk}$. We write

$m_{ij1}=n_{ij1}+\delta^{+}_{ij}-\delta^{-}_{ij},$ (8)

where $\delta^{+}_{ij}$ and $\delta^{-}_{ij}$ indicate the number of
additional units to audit and the number of previously audited units to
exclude with $(X=i,Y=j)$, respectively. Note that during this procedure,
values in the marginal table $(X,Y)$ are not adjusted, only units are
transported from $Z=0$ to $Z=1$ and vice versa. Hence, $m_{ij+}=n_{ij+}$ for
all $i$ and $j$.

Now, the goal is to sample the (additional) units from different cells of
table $(X,Y,Z)$ in such a way that $D$ defined in Equation (6) is minimized.
This is a minimization problem for which the target function can be written
as:

$D(m)=C+2\sum_{i,j,k}m_{ijk}\textrm{ log }m_{ijk}-2\sum_{j,k}m_{+jk}\textrm{
log }m_{+jk},$ (9)

since $C$ remains equal to (7), as $m_{ij+}=n_{ij+}$ and $m_{+j+}=n_{+j+}$.
Therefore, $C$ can be ignored when minimizing $D(m)$.

When minimizing $D(m)$, a number of constraints apply. First, the user-
specified bounds on the number of additional units to include ($M_{+}$) and
the number of units to exclude ($M_{-}$) lead to the following constraints:

$\displaystyle\sum_{i,j}\delta^{+}_{ij}$ $\displaystyle\leq M_{+}$ (10)
$\displaystyle\sum_{i,j}\delta^{-}_{ij}$ $\displaystyle\leq M_{-}.$ (11)

Second, for each combination $(X=i,Y=j)$, the number of units in- and excluded
in the audit sample should add up to the known marginal value:

$\displaystyle m_{ij1}+m_{ij0}=n_{ij+},\qquad$ $\displaystyle i=1,\dots,I;$
(12) $\displaystyle j=1,\dots,J,$

where $I$ is the number of categories in $X$ and $J$ is the number of
categories in $Y$. Together with $m_{ij1}$ defined in Equation (8), this
restriction implies that

$m_{ij0}=n_{ij0}-\delta^{+}_{ij}+\delta^{-}_{ij}.$ (13)

Third, for each combination $(X=i,Y=j)$, there exist bounds on
$\delta_{ij}^{+}$ and $\delta_{ij}^{-}$ based on the initial counts $n_{ijk}$,
since no more units can be moved from $Z=0$ to $Z=1$ or vice versa than are
initially available:

$\displaystyle 0\leq\delta_{ij}^{+}\leq n_{ij0},\qquad$ $\displaystyle
i=1,\dots,I;$ $\displaystyle j=1,\dots,J.$ (14) $\displaystyle
0\leq\delta_{ij}^{-}\leq n_{ij1},\qquad$ $\displaystyle i=1,\dots,I;$
$\displaystyle j=1,\dots,J.$ (15)

In combination with the other restrictions, these bounds imply that
$m_{ijk}\geq 0$ for every value.

The minimization procedure can now be written as follows:

$\min\Big{\\{}C+2\sum_{i,j,k}m_{ijk}\textrm{ log
}m_{ijk}-2\sum_{j,k}m_{+jk}\textrm{ log }m_{+jk}\Big{\\}}$ (16)

under constraints

$\displaystyle m_{ij1}=n_{ij1}+\delta^{+}_{ij}-\delta^{-}_{ij};$
$\displaystyle m_{ij0}=n_{ij0}-\delta^{+}_{ij}+\delta^{-}_{ij};$
$\displaystyle\sum_{i,j}\delta^{+}_{ij}\leq M_{+};$
$\displaystyle\sum_{i,j}\delta^{-}_{ij}\leq M_{-};$ $\displaystyle
0\leq\delta_{ij}^{+}\leq n_{ij0};$ $\displaystyle 0\leq\delta_{ij}^{-}\leq
n_{ij1}.$

This is a non-linear minimization problem with linear constraints which can be
solved by standard software, for instance using the R function ‘constrOptim’
[13]. To find the optimal value efficiently, it is useful to have the gradient
of the objective function $D(m)$, i.e., the vector of partial derivatives with
respect to the free parameters $\delta_{ij}^{+}$ and $\delta_{ij}^{-}$. For
the elements of this gradient, we find [using (8) and (13) and some standard
calculus]:

$\displaystyle\frac{\partial D(m)}{\partial\delta_{ij}^{+}}$
$\displaystyle=2\bigg{(}\log m_{ij1}-\log m_{ij0}-\log m_{+j1}+\log
m_{+j0}\bigg{)};$ (17) $\displaystyle\frac{\partial
D(m)}{\partial\delta_{ij}^{-}}$ $\displaystyle=2\bigg{(}\log m_{ij0}-\log
m_{ij1}-\log m_{+j0}+\log m_{+j1}\bigg{)}.$ (18)

###  3.2 A procedure for optimizing the audit sample

In practice, it turns out that the objective function defined in Equation (16)
is difficult to minimize, because it has many local minima. To ensure that a
global minimum is found, multiple starting values have to be tried in the
optimization algorithm.

Algorithm 1 Optimize the audit sample

input: The number of observations $n_{ijk}$ ($i=1,\ldots,I;\ j=1,\ldots,J;\
k=0,1$).

1:Initialize $D_{\text{best}}$ (the best deviance found so far) as the value
of $D$ for the original counts $n_{ijk}$.

2:Determine the bounds $M_{+}$ and $M_{-}$.

3:Determine $N_{\text{attempts}}$, the maximal number of attempts to search
for the best solution.

4:for 1 : $N_{\text{attempts}}$ do

5: Select random starting values for the solution:

  * •

Draw a value $\widetilde{M}_{+}$ from a uniform distribution on
$[l_{+}M_{+},u_{+}M_{+}]$, where $l_{+}$ and $u_{+}$ are chosen constants with
$0<l_{+}<u_{+}<1$.

  * •

Assign random starting values $\delta_{ij}^{+}$ inside their feasible
intervals such that $\sum_{i,j}\delta_{ij}^{+}=\widetilde{M}_{+}$.

  * •

Draw a value $\widetilde{M}_{-}$ from a uniform distribution on
$[l_{-}M_{-},u_{-}M_{-}]$, where $l_{-}$ and $u_{-}$ are chosen constants with
$0<l_{-}<u_{-}<1$.

  * •

Assign random starting values $\delta_{ij}^{-}$ inside their feasible
intervals such that $\sum_{i,j}\delta_{ij}^{-}=\widetilde{M}_{-}$.

6: Run the optimization algorithm to solve (16) with these starting values.

7: Compute $D(m)$ for the current solution. If $D(m)<D_{\text{best}}$, then
set $D_{\text{best}}:=D(m)$ and store the current solution. Otherwise, discard
the current solution.

8:end for

9:If $D_{\text{best}}$ is still considered too large, return to step 3 and
change the bounds $M_{+}$ and/or $M_{-}$.

output: The number of additional units to audit $\delta_{ij}^{+}$ and remove
from the audit sample $\delta_{ij}^{-}$ ($i=1,\ldots,I;\ j=1,\ldots,J$) found
for the best solution, as well as the value of $D_{\mathrm{best}}$.

For solving Equation (16), we propose the practical procedure found under
Algorithm 1. To help decide when the best deviance value $D_{\text{best}}$
should be considered ‘too large’ in step 9, it may be noted that under the
independence model $(XY)(YZ)$ the deviance asymptotically follows a chi-square
distribution with $J\times(I-1)$ degrees of freedom [1]. Hence, the
$(1-\alpha)\times 100\%$ percentile $\chi_{J(I-1)}^{2}(1-\alpha)$ of this
distribution (e.g., with $\alpha=.05$) could be used as a cut-off point: we
accept the current solution in step 9 when
$D_{\text{best}}\leq\chi_{J(I-1)}^{2}(1-\alpha)$.

Once the optimal values $\delta_{ij}^{+}$ and $\delta_{ij}^{-}$ have been
obtained, a representative audit sample can be obtained by applying the
following two steps (in parallel) to each stratum $(X=i,Y=j)$:

  * •

If $\delta_{ij}^{+}>0$, draw a simple random sample without replacement of
size $\delta_{ij}^{+}$ from the $n_{ij0}$ units in this stratum with $Z=0$.
These units are selected for additional auditing, so moved to $Z=1$.

  * •

If $\delta_{ij}^{-}>0$, draw a simple random sample without replacement of
size $\delta_{ij}^{-}$ from the (original) $n_{ij1}$ units in this stratum
with $Z=1$. These previously audited units are removed from the audit sample,
so moved to $Z=0$.

Once these steps have been run, the final audit sample consists of all units
with $Z=1$.

###  3.3 Minimizing the deviance while maximizing the number of re-used cases

The approach introduced in Section 3.2 uses deviance as a primary criterion to
select an audit sample, where the sample is selected with the lowest value for
deviance. Recycling units that have been in the initial audit sample is
considered as a secondary criterion. However, in practice a deviance that is
as low as possible is not always essential, and recycling as many cases as
possible can be more important. Depending on sample size and number of
categories in $Y$ a critical value for deviance can be determined and all
selected samples with a deviance lower than the critical value can be
considered as being representative with respect to $Y$.

More formally stated, the specification of the optimization problem in (16)
allows that the optimal solution contains combinations $(i,j)$ where
$\delta_{ij}^{+}>0$ and $\delta_{ij}^{-}>0$ simultaneously. This corresponds
to a solution where in the same stratum some units are added to the audit
sample while other, previously audited units, are removed from the audit
sample. This seems inefficient from a practical point of view.

It is easy to correct this after the solution has been found, by moving to an
equivalent solution with

$\displaystyle\delta_{ij,alt}^{+}$
$\displaystyle=\max\left\\{0,\delta_{ij}^{+}-\delta_{ij}^{-}\right\\},$
$\displaystyle\delta_{ij,alt}^{-}$
$\displaystyle=\max\left\\{0,\delta_{ij}^{-}-\delta_{ij}^{+}\right\\}.$

Clearly,
$\delta_{ij,alt}^{+}-\delta_{ij,alt}^{-}=\delta_{ij}^{+}-\delta_{ij}^{-}$, so
this corresponds to the same solution in terms of $m_{ijk}$, but now it holds
that $\delta_{ij,alt}^{+}\delta_{ij,alt}^{-}=0$ in all cases, as desired.
However, it would be preferable to avoid this type of inefficient solution
altogether.

A possible solution is to replace the objective function $D(m)$ of the
minimization problem defined in Equation (16) by one of the following
alternatives:

$\displaystyle F_{1}(m)$
$\displaystyle=D(m)+\lambda\sum_{i,j}\big{(}\delta_{ij}^{+}+\delta_{ij}^{-}\big{)},$
(19) $\displaystyle F_{2}(m)$
$\displaystyle=D(m)+\exp\left\\{-\frac{D(m)}{\kappa}\right\\}\sum_{i,j}\big{(}\delta_{ij}^{+}+\delta_{ij}^{-}\big{)}.$
(20)

where $\lambda$ is a small positive constant (e.g., $\lambda=0.01$) and
$\kappa$ is a positive constant such that $D(m)/\kappa$ is larger than, say,
$10$ for any candidate solution with a deviance value for which the
independence model would be rejected (as discussed in Section 3.2).

Both alternative objective functions penalize candidate solutions with large
total numbers of additional units to audit and previously audited units to
remove (i.e., large values of
$\sum_{i,j}\big{(}\delta_{ij}^{+}+\delta_{ij}^{-}\big{)}$). This should make
any candidate solution with $\delta_{ij}^{+}\delta_{ij}^{-}>0$ unattractive,
because the value of $F_{1}(m)$ or $F_{2}(m)$ can be reduced by replacing
$(\delta_{ij}^{+},\delta_{ij}^{-})$ by
$(\delta_{ij,alt}^{+},\delta_{ij,alt}^{-})$ as above. The exponential term in
$F_{2}(m)$ is supposed to ensure that the penalty term becomes relevant only
for candidate solutions with small (i.e., acceptable) values of $D(m)$. For
large values of $D(m)$, it holds that $F_{2}(m)\approx D(m)$. The same can be
achieved with $F_{1}(m)$, provided that the constant $\lambda$ is chosen with
some care.

The problem of minimizing (16) with the target function replaced by (19) or
(20) can be solved by the same procedure as outlined in Section 3.2. To obtain
the gradient values for $F_{1}(m)$, we simply add a term $\lambda$ to (17) and
(18). For $F_{2}(m)$ we obtain:

$\frac{\partial F_{2}(m)}{\partial x}=\frac{\partial D(m)}{\partial
x}+\exp\left\\{-\frac{D(m)}{\kappa}\right\\}\left\\{1-\frac{1}{\kappa}\frac{\partial
D(m)}{\partial
x}\sum_{i,j}\big{(}\delta_{ij}^{+}+\delta_{ij}^{-}\big{)}\right\\},$

where $x=\delta_{ij}^{+}$ or $x=\delta_{ij}^{-}$.

##  4 Simulation study

###  4.1 Simulation approach for deviance and bias

To empirically evaluate the performance of the audit sample selection
procedure, we conducted a simulation study using R [13]. The code used for the
simulation study is available on https://github.com/lauraboeschoten/audit_2021
and a small illustrative example of the code applied to one situation is
available at
https://github.com/lauraboeschoten/audit_2021/tree/master/reproducible_example.

The main simulation study consists of three sets of conditions. For each set,
theoretical populations with varying relationships between $W$, $X$, $Y$ and
$Z$ are generated, where $X$ and $Y$ have three categories and $Z$ has two
categories by design. Note that in this simulation study, the target
population and observed sample are the same. The first set investigates how
the performance of the proposed procedure might be affected by the strength of
the relationship between $X$ and $Z$ before the start of the optimization
procedure, varying from no relationship to a strong relationship. The second
set investigates how the performance might be affected by the strength of the
relationship between $W$ and $X$, varying from a perfect relationship (meaning
that observed variable $X$ is a perfect measurement of the outcomes after
audit $W$) to an imperfect and asymmetrical relationship (meaning that
observed variable $X$ contains measurement error and that the probability of
an incorrect score differs for different scores of the audit variable $W$).
The third set investigates how the performance of the procedure might be
affected by the strength of the relationship between $W$ and $Y$, varying from
a strong relationship to a weak relationship. Note that for each set, the
conditions are ordered from most desired situation to least desired situation.
Each set contains four conditions resulting in a total of twelve simulation
conditions, which can all be found in Table 1. Note that when the different
relationships between $X$ and $Z$ are investigated, the first conditions
listed in Table 1 for $(W,X)$ and $(W,Y)$ are selected and similar approaches
are used when investigating $(W,X)$ and $(W,Y)$ to investigate the main
effects of these relationships. In the last part of the study, the
interactions between the selected relations are investigated by taking the
most desired and least desired condition for each relation, and investigating
their combinations in a full factorial design.

In all conditions, a joint probability density is generated by multiplying the
probabilities listed in Table 1. The probabilities generated here follow the
log-linear model $(WX)(WY)(XZ)$. This model is contained within the previously
discussed maximal model $(WXY)(XYZ)$ and it contains a term $(XZ)$ that cannot
be found in the independence model $(WXY)(YZ)$. Because of this term, it is
expected that this model results in a large value for the deviance and it
would therefore be beneficial to adjust the audit sample via the proposed
minimization procedure. For each of the described conditions, $1,000$ data-
sets of size $N=10,000$ are sampled from the generated joint probability
density. The size of the audit sample is $300$ in all conditions, as follows
also from the proportions described in Table 1.

Simulation condition number |  |  |  |  |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---
|  |  | 1 |  |  |  | 2 |  |  |  | 3 |  |  |  | 4 |
|  | X |  |  |  | X |  |  |  | X |  |  |  | X |  |
|  | 1 | 2 | 3 |  | 1 | 2 | 3 |  | 1 | 2 | 3 |  | 1 | 2 | 3
Z | 0 | .323 | .323 | .323 |  | .323 | .323 | .323 |  | .323 | .323 | .323 |  | .323 | .323 | .323
| 1 | .010 | .010 | .010 |  | .012 | .010 | .008 |  | .015 | .010 | .005 |  | .018 | .010 | .002
|  | W |  |  |  | W |  |  |  | W |  |  |  | W |  |
|  | 1 | 2 | 3 |  | 1 | 2 | 3 |  | 1 | 2 | 3 |  | 1 | 2 | 3
X | 1 | .333 | 0 | 0 |  | .267 | .033 | .033 |  | .333 | 0 | 0 |  | .300 | .017 | .017
| 2 | 0 | .333 | 0 |  | .033 | .267 | .033 |  | .017 | .300 | .017 |  | .033 | .267 | .033
| 3 | 0 | 0 | .333 |  | .033 | .033 | .267 |  | .033 | .033 | .267 |  | .050 | .050 | .233
Y | 1 | .267 | .033 | .033 |  | .200 | .067 | .067 |  | .300 | .017 | .017 |  | .267 | .033 | .033
| 2 | .033 | .267 | .033 |  | .067 | .200 | .067 |  | .033 | .267 | .033 |  | .067 | .200 | .067
| 3 | .033 | .033 | .267 |  | .067 | .067 | .267 |  | .050 | .050 | .233 |  | .100 | .100 | .133
Table 1: Overview of the bivariate relationships used to generate data for the
simulation study

Per condition, the procedure as described in Section 3.2 is applied on each
generated data-set. The bounds $M_{+}$ and $M_{-}$ that should be determined
at step 2 of Algorithm 1 are set at $100$ and $10$ respectively and the
maximal number of attempts for the optimization procedure to search for a
solution at step 3 is set at $N_{\text{attempts}}=200$. Note that step 9
(adjusting the bounds $M_{+}$ and $M_{-}$ if the obtained value for the
deviance is considered too large) is left out to appropriately compare the
results under different conditions. The settings described here remain
consistent for all simulation conditions.

We considered two types of target parameters to estimate in this simulation
study: the true proportion of units in the population with $W_{g}=w$ for each
category $w$,

$P^{W}_{w}=\frac{1}{N}\sum_{g=1}^{N}I(W_{g}=w),$ (21)

and the proportion of units in the population with true category $W_{g}=w$
that are observed in category $X_{g}=x$ (measurement error probabilities):

$P^{X|W}_{x|w}=\frac{\sum_{g=1}^{N}I(W_{g}=w,X_{g}=x)}{\sum_{g=1}^{N}I(W_{g}=w)}.$
(22)

Under the assumption that the independence model holds, $P^{W}_{w}$ defined in
(21) can be estimated without bias from the audit sample by

$\hat{P}^{W}_{w}=\sum_{y}P^{Y}_{y}\frac{\sum_{g=1}^{n}I(W_{g}=w,Y_{g}=y)}{\sum_{g=1}^{n}I(Y_{g}=y)}\equiv\sum_{y}P^{Y}_{y}p^{W|Y}_{w|y},$
(23)

where $P^{Y}_{y}$ is the (known) proportion of units with $Y_{g}=y$ in the
population and $p^{W|Y}_{w|y}$ denotes the (unweighted) observed proportion of
cases with $Y_{g}=y$ in the audit sample that also have $W_{g}=w$. (Here, for
convenience it is assumed that the audit sample consists of the first $n$
units.) Similarly, $P^{X|W}_{x|w}$ defined in (22) can be estimated
consistently from the audit sample by:

$\hat{P}^{X|W}_{x|w}=\frac{\sum_{y}P^{Y}_{y}\frac{\sum_{g=1}^{n}I(W_{g}=w,X_{g}=x,Y_{g}=y)}{\sum_{g=1}^{n}I(Y_{g}=y)}}{\sum_{y}P^{Y}_{y}\frac{\sum_{g=1}^{n}I(W_{g}=w,Y_{g}=y)}{\sum_{g=1}^{n}I(Y_{g}=y)}}\equiv\frac{\sum_{y}P^{Y}_{y}p^{WX|Y}_{w,x|y}}{\sum_{y}P^{Y}_{y}p^{W|Y}_{w|y}}.$
(24)

Note that this is a special case of the so-called combined ratio estimator
[5].

We investigate the bias of ${P}^{W}_{w}$,
$B({P}^{W}_{w})=E(\hat{P}^{W}_{w}-P^{W}_{w})$ and of $\hat{P}^{X|W}_{x|w}$,
$B(\hat{P}^{X|W}_{x|w})=E(\hat{P}^{X|W}_{x|w}-P^{X|W}_{x|w}$) while selecting
only $Z=1$ both before and after the procedure and compare these results.

###  4.2 Simulation approach for variance

In practice, it may be relevant to estimate the (design-based) variance of
$\hat{P}^{W}_{w}$ and $\hat{P}^{X|W}_{x|w}$ from the observed data. Under the
assumption that the independence model holds, we propose to estimate these
variances by treating the audit sample as a stratified simple random sample
with $Y$ as a stratifying variable; see the remark at the end of this
subsection for a justification. For $\hat{P}^{W}_{w}$ this yields the
following variance estimator, assuming for simplicity that the sampling
fraction in each stratum is small enough so that finite population corrections
can be neglected:

$\widehat{\mathrm{var}}(\hat{P}^{W}_{w})=\sum_{y}\frac{(P^{Y}_{y})^{2}}{n_{y}}p^{W|Y}_{w|y}(1-p^{W|Y}_{w|y}),$
(25)

where $n_{y}$ denotes the number of units with $Y_{g}=y$ in the audit sample.
Similarly, under the same assumptions the approximate variance of
$\hat{P}^{X|W}_{x|w}$ can be derived from expression (6.51) in [5]. Assuming
that $NP^{Y}_{y}\gg 1$ for all $y$, we obtain the following variance
estimator:

$\displaystyle\widehat{\mathrm{var}}(\hat{P}^{X|W}_{x|w})$
$\displaystyle=\frac{1}{(\hat{P}^{W}_{w})^{2}}\sum_{y}\frac{(P^{Y}_{y})^{2}}{n_{y}}\left\\{p^{WX|Y}_{w,x|y}(1-p^{WX|Y}_{w,x|y})+(\hat{P}^{X|W}_{x|w})^{2}p^{W|Y}_{w|y}(1-p^{W|Y}_{w|y})\right.$
$\displaystyle\phantom{=}\qquad\left.-2\hat{P}^{X|W}_{x|w}p^{WX|Y}_{w,x|y}(1-p^{W|Y}_{w|y})\right\\}.$
(26)

Note that these design-based variances treat the target population of size $N$
as fixed. To evaluate the performance of these variance estimators, a separate
set of simulations was run. Combining the least desired condition of $(W,X)$
in Table 1 (i.e., the right-most condition) with each of the four conditions
of $(W,Y)$, four fixed target populations of size $N=10,000$ were generated.
From each of these populations, $1,000$ initial audit samples were generated
according to the least desired condition of $(X,Z)$ in Table 1. Again, the
procedure from Section 3.2 was applied to each audit sample, with the same
settings as before. For each final audit sample, we computed $\hat{P}^{W}_{w}$
and $\hat{P}^{X|W}_{x|w}$ as well as the associated variance estimates from
(25) and (26). This allowed us to compare the estimated variances with the
empirical variances of $\hat{P}^{W}_{w}$ and $\hat{P}^{X|W}_{x|w}$ across
$1,000$ simulation rounds.

_Remark:_ To justify the assumption used to derive (25) and (26) — that the
final audit sample resembles a stratified simple random sample — we note that
the independence model $(WXY)(YZ)$ is equivalent to a logistic regression
model for $Z$ that, besides the constant term, uses only $Y$ as a predictor;
cf. [1, Section 9.5.1 ]. Hence, the distribution of the sample size $n_{y}$
per stratum does not depend on the target parameters $P^{W}_{w}$ and
$P^{X|W}_{x|w}$: the sample sizes $n_{y}$ are _ancillary statistics_ with
respect to these target parameters. As discussed by [10], it is preferable to
evaluate the variance of estimated parameters conditional on any ancillary
statistics. Hence, in our case the variances of $\hat{P}^{W}_{w}$ and
$\hat{P}^{X|W}_{x|w}$ should be evaluated conditional on the realised sample
sizes $n_{y}$. Moreover, if the independence model holds after the
optimization procedure, according to the above interpretation as a logistic
regression model, all units in the same stratum based on $Y$ have the same
(final) inclusion probability. Therefore, when the realised sample sizes
$n_{y}$ are treated as fixed, in the absence of more information it makes
sense to consider the design of the audit sample as equivalent to a stratified
simple random sample.

###  4.3 Results

####  4.3.1 Deviance

Figure 1: Results in terms of relative deviance (see text). Boxplots
demonstrate the spread of the 1,000 relative deviance values obtained per
simulation condition. The three panels represent the three sets of bivariate
relationships specified in the study, and the four columns per panel
illustrate the four alternative strengths of those relationships specified
(see Table 1).

Figure 1 illustrates the different values obtained for relative deviance under
the different simulation conditions. By relative deviance we mean the deviance
after applying the optimization procedure as a proportion of the deviance
before applying the procedure, so a value close to zero means that a new audit
sample is drawn that substantially improves the representativity of the audit.

In the left panel, the results are shown for the four different $WX$
relationships under which the initial sample was drawn. As representativity is
defined with respect to $Y$, not with respect to $W$, differences in relative
deviance between these $WX$ conditions are not expected and this is confirmed
here. In the middle panel, the results are shown for the four different $WY$
relationships under which the initial sample was drawn. The (lack of)
differences between these boxplots illustrates that the proposed method is
able to perform in situations where the relationship between background
variables and the variable of interest in the audit sample is strong and
possibly also unbalanced. In the right panel, results are shown for the four
different $XZ$ relationships under which the initial audit sample was drawn.
Here, a small increase in average relative deviance can be detected when
inclusion in the initial audit sample relates more strongly to scores on $X$,
while the spread of the results becomes drastically smaller. This is likely to
be caused by the fact that the deviance of the initial model was more
substantive.

Figure 2: Results in terms of relative deviance (see text). Boxplots
demonstrate the spread of the 1,000 relative deviance values obtained per
simulation condition for eight different combinations of different strengths
of relationships between the variables $(WX)$, $(WY)$ and $(YZ)$ (see Table
1).

Figure 2 illustrates the different values obtained for relative deviance under
interactions between the different simulation conditions. As indicated by the
boxplot labels, for each boxplot, the starting dataset is generated under a
different combination of conditions, for which the complete overview can be
found in Table 1. These boxplots illustrate that the scores for relative
deviance are stable over different data-generating conditions. It is
particularly noteworthy that a strong relationship for $(XZ)$ apparently has
the most substantial influence on scores for relative deviance.

####  4.3.2 Bias in W

Figure 3: Boxplots of the bias distribution of the estimated proportions of
$W$, based on 1,000 replicates per simulation condition. The three panels
represent the three sets of bivariate relationships specified in the study,
and the four columns per panel illustrate the four strengths of those
relationships (see Table 1).

Figure 3 illustrates the different values obtained for bias under the
different simulation conditions compared to when the initial audit sample
would have been used directly. Only the results for $W=1$ are shown as the
results for $W=2$ and $W=3$ behaved very similarly.

In the left panel, the results are shown for the four different $WX$
relationships under which the initial sample was drawn. In the first condition
$W=X$. As $X$ is the observed target variable of interest and $W$ its true
version, it is expected that no bias is present in results for $W$ if $W=X$.
In conditions 2-4, $W\neq X$ and small increases in spread can be detected as
a consequence. However, it can also be seen that if the procedure is not
applied in such cases, bias is induced in the results for $W$. In the middle
panel, the results are shown for the four different $WY$ conditions. Similar
to the results found for the deviance, it can be concluded that the
relationship between the variable of interest and background variables does
not affect the bias in $W$. However, bias is also not detected before the
procedure was applied, so a correction procedure is not required in such
situations. In the right panel, results are shown for the four different $XZ$
conditions. Here, a small increase in bias can be detected when inclusion in
the initial audit sample relates more strongly to scores on $X$, i.e. the
procedure finds it more difficult to obtain a sample that is unbiased with
respect to $W$ if inclusion in the initial sample becomes more unbalanced with
respect to $X$. However, not applying the procedure in such situation results
in very substantive amounts of bias in $W$.

Figure 4: Boxplots of the distribution of the bias of the estimated
proportions of $W$, based on 1,000 replicates per simulation condition for
eight combinations of different strengths of relationships between the
variables $(WX)$, $(WY)$ and $(YZ)$ (see Table 1).

Figure 4 illustrates the different values of bias in $W$ after applying the
method under interactions between the different simulation conditions. As
indicated by the boxplot labels, for each boxplot, the data is generated under
a different combination of conditions. These boxplots illustrate that although
the bias present in estimates of $W$ can be caused by both a weaker
relationship in $(WX)$ or a stronger relationship in $(XZ)$, the combination
of these conditions do not result in a multiplication in the effect in terms
of bias after applying the procedure. Furthermore, deviations in $(WX)$ result
in a wider spread while deviations from independence between $Z$ and $X$
result in more systematic bias. In addition, these boxplots again illustrate
that in situations of an imbalance in $XZ$, applying a procedure to obtain a
representative audit sample is essential.

####  4.3.3 Bias in XW

Figure 5: Boxplots of the distribution of the bias of the estimated
proportions of $X$ conditional on $W$, based on 1000 replicates per simulation
condition. The three panels represent the three sets of bivariate
relationships specified in the study, and the four columns per panel
illustrate the four alternative strengths of those relationships specified
(see Table 1).

Figure 5 illustrates the different values obtained for bias in $X\mid W$. Only
the results for $(X=1\mid W=1)$ are shown as the results for other proportions
were very similar. In the left panel, the results are shown for the four
different $WX$ relationships. Similarly as when evaluating the bias for $W$,
$X$ is the observed target variable of interest and $W$ its true version. It
is therefore expected that no bias is present in results for $X\mid W$ if
$W=X$. In conditions 2 and 4, $W\neq X$ and small increases in bias can be
detected as a consequence, while in condition 3 there is also no bias present
due to the fact that here $W=X$ for $X=1$. Similarly, no bias is present in
$X\mid W$ under the different simulation conditions for $WY$ and $XZ$, as
under these conditions $W=X$ as well. In cases of observed bias, the spread of
the bias is smaller when the procedure is applied compared to when the
procedure is not applied.

Figure 6: Boxplots of the distribution of the bias of the estimated
proportions of $X$ conditional on $W$ before and after applying the procedure,
based on 1000 replicates per simulation condition for eight combinations of
different strengths of relationships between the variables $(WX)$, $(WY)$ and
$(YZ)$ (see Table 1).

In Figure 6 it is illustrated that if $W\neq X$, generally more bias is
present, although the amount of bias remains limited if the audit sampling
procedure is applied. Furthermore, it is shown that combinations of
measurement error in $W$ ($W\neq X$), different strengths of $WY$ and
selectivity in the initial audit sample, $XZ$, do not cause for more
substantive bias in $X\mid W$ than already caused by $WX$. Again, in such
cases the bias is more substantive if the audit sample procedure has not been
applied.

####  4.3.4 Variance estimation

Table 2 and Table 3 summarize the results of the separate simulation study
mentioned in Section 4.2, to evaluate the performance of the variance
estimators defined in Equation (25) and Equation (26) for a fixed target
population. Each se-sd ratio in these tables represents the ratio of the
average standard error (se) and the empirical standard deviation (sd) across
1,000 simulations, for a particular condition and target parameter $P^{W}_{w}$
(Table 2) or $P^{X|W}_{x|w}$ (Table 3). Ideally, this se-sd ratio would be
equal to $1$.

For the estimated $P^{W}_{w}$ the empirical sd is also shown in Table 2. It is
seen that the sd values were similar for all four conditions where the initial
audit sample was selective. For the estimated $P^{X|W}_{x|w}$ the empirical sd
values are omitted here to save space; again, these values were similar across
all four conditions.

| True sd |  | Se-sd ratio
---|---|---|---
Condition | $w=1$ | $w=2$ | $w=3$ |  | $w=1$ | $w=2$ | $w=3$
WX4,WY1,XZ4 | 0.012 | 0.013 | 0.013 |  | 1.22 | 1.23 | 1.18
WX4,WY2,XZ4 | 0.012 | 0.013 | 0.015 |  | 1.50 | 1.38 | 1.29
WX4,WY3,XZ4 | 0.012 | 0.014 | 0.015 |  | 1.20 | 1.15 | 1.13
WX4,WY4,XZ4 | 0.012 | 0.014 | 0.014 |  | 1.40 | 1.35 | 1.38
WX4,WY1,XZ1* | 0.019 | 0.020 | 0.019 |  | 1.03 | 0.99 | 1.03
Table 2: Estimated true standard deviation (true sd) and se-sd ratio based on 1,000 simulations for the estimated proportions of $W$. Rows refer to different simulation conditions, columns to different categories of $W$. (* = no optimization.) | $w=1$ |  | $w=2$ |  | $w=3$
---|---|---|---|---|---
Condition | $x=1$ | $x=2$ | $x=3$ |  | $x=1$ | $x=2$ | $x=3$ |  | $x=1$ | $x=2$ | $x=3$
WX4,WY1,XZ4 | 1.18 | 1.11 | 1.19 |  | 1.20 | 1.27 | 1.27 |  | 1.27 | 1.41 | 1.45
WX4,WY2,XZ4 | 1.08 | 1.07 | 1.06 |  | 1.11 | 1.19 | 1.14 |  | 1.28 | 1.23 | 1.35
WX4,WY3,XZ4 | 1.21 | 1.18 | 1.17 |  | 1.19 | 1.27 | 1.27 |  | 1.22 | 1.33 | 1.35
WX4,WY4,XZ4 | 1.13 | 1.11 | 1.09 |  | 1.18 | 1.25 | 1.14 |  | 1.15 | 1.18 | 1.28
WX4,WY1,XZ1* | 0.99 | 0.98 | 0.95 |  | 1.02 | 1.02 | 0.97 |  | 1.00 | 1.01 | 1.02
Table 3: Results in terms of se-sd ratio based on 1,000 simulations for the
estimated proportions of $X$ conditional on $W$. Rows refer to different
simulation conditions, columns to different combinations of categories of $W$
and $X$. (* = no optimization.)

It is seen that for all conditions where the initial audit sample was
selective, the variance estimators as defined in Equation (25) and Equation
(26) tended to overestimate the true variance; the overestimation in terms of
standard error varied between $6\%$ and $50\%$, which in practice could be
considered a moderate bias. As a benchmark, we also included a condition where
$(X,Z)$ is given by the left-most matrix in Table 1, so that the initial audit
sample already satisfied the independence model. In this case we did not apply
the optimization procedure to alter the initial sample. Here it is seen that
both variance estimators performed much better (last line in Tables 2 and 3).
Thus, for a truly random sample from the independence model, our variance
estimators are approximately correct, and the overestimation of the variance
for the other conditions is due to the optimization procedure. Apparently, by
minimizing the deviance, the distribution of possible audit samples obtained
by this procedure is more restricted than a stratified simple random sampling
design, and this is not reflected by Equation (25) and Equation (26).

Within the selective conditions, it is seen that the degree of overestimation
in both tables differs between, on the one hand, the first and third condition
and, on the other hand, the second and fourth condition of $(W,Y)$. From Table
1, it is seen that this distinction corresponds to associations between $W$
and $Y$ that are relatively strong and relatively weak, respectively. For the
estimated proportions $\hat{P}^{W}_{w}$, the overestimation by variance
estimator (25) in Table 2 appears to be smaller when the association between
$W$ and $Y$ is stronger. Surprisingly, for the estimated error probabilities
$\hat{P}^{X|W}_{x|w}$ and variance estimator (26), the opposite effect is seen
in Table 3.

##  5 Application

Since the introduction of the sustainable development goals [19], the interest
in statistics regarding energy consumption has increased substantially,
including many statistics which have never been produced before. One statistic
of interest is the total energy consumption of establishments per economic
sector. To obtain these statistics, a statistical register is required
containing all companies, the economic sector they operate in and the address
of their location. To construct this database, multiple incomplete
administrative sources can be combined on unit level, where the unit in this
case is an establishment. Many of these different data sources contain
information about the economic activity of the establishment, and the observed
economic activity per establishment may vary over the sources. The economic
activity is classified according to the NACE rev 2. codes [8].

To investigate the quality of the constructed database, especially of the
final code for the economic activity, an audit should be performed on a
representative sample of the population of establishments. Some units (a non-
representative sample) have already been audited throughout the year, and to
efficiently perform an audit on a representative sample, as many of the
previously audited units should be included in this new audit sample.

The combined administrative data-set comprises $2,037,088$ establishments in
$2019$. Of these units, on $173,617$ either a previous audit has been
performed or the NACE codes were equal over three or more different data
sources. However, the distribution of audited units is not representative for
the distribution of total units with respect to background variable $Y$ and
target variable $X$. Here, the target variable $X$ concerns $21$ economic
sectors, which are first digit NACE codes, while background variable $Y$
contains six categories:

  * •

LV: Low voltage

  * •

HV: High voltage

  * •

SG: Small gas consumption

  * •

MG: Middle gas consumption

  * •

LG: Large gas consumption

  * •

OPC: other profiles

The upper heat map of Figure 7 illustrates the observed proportions of
economic sectors $X$ per profile category $Y$. The purpose of the intended
audit is to investigate the quality of these economic sector codes. In this
heatmap we for example see that in the ’other profiles’ category, most units
are in the economic sector “public administration and defense, compulsory
social security” (“O”). In the category of large gas consumers, most units are
in the agricultural sector (“A”) followed by industry (“C”). In the category
of small gas consumers, most units are in the sector wholesale and retail
(“G”).

Of all units shown in the upper heat map of Figure 7, a subset has previously
been audited, and the final economic sector codes for these units are known
and can be found in the middle heat map of Figure 7. By comparing the two heat
maps, it can be clearly seen that given certain profiles certain sectors have
been more thoroughly audited than others. For example, units in the economic
sector “human health and social work activities” (“Q”) have been audited
relatively thoroughly, while audits for sectors such as water and waste (“E”),
construction (“F”), finance (“K”) real estate (“L”), science (“M”) and
administration (“N”) are relatively scarce.

To obtain an audit sample that is more representative in terms of the economic
sectors and with respect to the background profiles, we applied the
methodology introduced in this manuscript. We first applied the methodology
discussed in Section 3. Here, we applied the method multiple times using audit
sample sizes $500$, $1000$, $1500$ and $2000$ ($M_{+}$), and for each sample
size we varied the number of cases to exclude from the initial sample
($M_{-}$). The size of $M_{-}$ was specified as a factor with respect to
$M_{+}$, and we used the values of $50$, $100$ and $120$ for this factor.
Here, we investigated the resulting deviance values and concluded that we
would probably be able to obtain a representative audit sample of a size
between $1000$ and $1500$ and we concluded that we would need to increase the
factor for removing cases. Therefore, we started with an audit sample of
$1500$ and were able to reduce the sample size to $1200$. We tried to reduce
the sample size even further to $1100$ but this did not result in a
representative sample by using our deviance criterion.

Now that we knew that we were able to obtain a representative sample by
sampling approximately $1200$ new cases, we applied the adjusted procedure for
the objective function $F_{2}$ as described in Section 3.3 to ensure that all
additionally sampled cases were indeed not in the initial audit. To assess the
deviance value, as suggested in Section 3.2, we compared it to a chi-square
distribution with $J\times(I-1)$ degrees of freedom, where $I$ and $J$ are the
number of categories of $X$ and $Y$, respectively, Here, with $I=21$ and
$J=6$, we used a chi-square distribution with $120$ degrees of freedom.
Choosing $\alpha=.05$, we found a value of 147 as an approximate cut-off
point. Thus, if we could draw an audit sample with a deviance lower than 147,
we would find this acceptable in terms of representativity.

Additionally, we tried to reduce the sample size even further, but this was
not successful. Finally, we selected a sample of $1200$ and we removed
$144000$ of the units from the initial audit sample (multiplication factor of
$120$), with a deviance of $117.27$.

The last heat map in Figure 7 illustrates the selected audit sample. This
audit sample is a combination of units from the initial audit (middle heat
map), units removed from this initial audit, and additionally selected units
for audit ($1200$ units). In this heat map, it can be seen that certain
economic sectors that were underrepresented in the initial audit are now more
prominently present, such as water and waste, construction, finance and
administration (“F”, “J”, “K” and “N”). Economic sectors real estate and
science (“L” and “M”) are still underrepresented. However, overall it can be
concluded that the distribution of the final audit sample is more similar to
the distribution in the combined administrative data set compared to the
initial audit sample, while the final audit comprises a smaller selection of
units.

Figure 7: Heatmap of the proportion of establishments per economic sector for
different profile categories: for the combined administrative data set (upper
panel), for the original audit sample (middle panel) and for the audit sample
adjusted with the proposed procedure (lower panel). Note that we excluded the
category ‘OPC’ as it highly affected the interpretability of the heatmap.

##  6 Conclusion

In this paper, we introduced a method that can be used to efficiently select a
representative audit sample if the situation is not straightforward. For
example, if a selective sample that is not representative with respect to some
background characteristics has already been audited, and a representative
audit sample is desired that re-uses the available sample as much as possible.

More specifically, our method uses the joint distribution between the observed
variable of interest, the background characteristics and an indicator for
initial audit inclusion. Furthermore, the method assumes that there is no
direct association between the true variable of interest and initial audit
inclusion conditional on the observed variable of interest and the background
characteristics. Here, we can analyse the previously mentioned joint
distribution using an independence model with respect to the background
characteristics. The fit in terms of deviance of this model can be compared to
that of the saturated model. Our method then numerically searches a new
solution that minimises the deviance by including new cases and excluding
already audited cases.

The simulation study illustrated that particularly when there is measurement
error in the observed target variable, this can cause biased estimates of the
target variable if the audit sample selection procedure is not applied to
obtain a more representative audit sample. In addition, the study illustrated
that without applying the procedure, bias in estimates of the target variable
is obtained in situations where inclusion in the initial audit is related to
the scores of the target variable of interest. This conclusion is particularly
relevant, as there is in practice often a relation between the scores on the
target variable and initial audit inclusion. With financial audits, the
largest companies are for example often audited by default. Furthermore, we
conclude from our simulation study that the variance when based on simple
random sampling is often overestimated compared to the true variance, due to
the fact that by minimizing the deviance, the distribution of possible audit
samples is more restricted than the distribution of simple random samples
would be.

The reason why we developed this procedure was that we wanted to perform an
audit on establishments to investigate whether classification errors occurred
in the economic sector that is used to for a publication on energy consumption
per economic sector. An audit sample was already available, but this was very
selective in the sense that not all economic sectors were evenly distributed.
Therefore we developed, investigated and applied the method introduced in this
paper. By applying this method, we were able to select an audit sample that
was more representative with respect to the economic sectors, that utilized
the cases that were already audited and we had control over the size of the
new audit sample.

We expect that the method is also interesting for other applications. For
instance, assume that we have trained and tested a supervised machine learning
model to derive a new variable. For instance, [18] describes that they were
interested to predict whether potential crimes in police records concern cyber
crime or not. Unfortunately, the annotated set that was used to train and test
the model was selective, because keywords were used to seek for cases likely
to be cyber crime rather than taking a randomised sample. Consequently, the
model prediction error in the test set might not be representative for the
true prediction error in the population. The police records also contain
background variables ($Y$) that relate to the target variable, such as crime
type and whether the case has been declared by a victim or not. One could then
use the predicted cyber crime in all police records by the originally trained
machine learning model as a proxy ($X$) for the true cyber variable ($W$)
which is available in the selective set. Next, one could apply the procedure
to select additional units to generate a test set that is more more
representative for the target population. The additionally selected units
could then be manually annotated to obtain $W$. Finally, one could use this
adjusted test set to obtain a more reliable estimate of the prediction error
of the model.

In this paper we illustrated that the method can be used in different ways.
For example, an audit sample can be selected that truly has the lowest
deviance, or a trade-off can be made between obtaining an acceptable value for
the deviance and maximizing the amount of cases from the initial audit to be
re-used. When applying the method in practice, the auditor should be aware of
the practical implications when selecting a certain amount of cases to leave
out or additional cases to include, particularly because including more cases
in the audit will improve the probability to meet the deviance criterion on
the one hand, but are also combined with increasing costs to perform the audit
in practice on the other hand.

Further research should give more insight in how the method performs when
different implementations of the method are chosen. In addition, more
practical applications should provide insight into what further improvements
are interesting for researchers. For example, are users of the method
typically interested in selecting a certain audit sample of a maximum size
that is within their budget? Or is this often more flexible and are users more
interested in an optimal combination of smallest sample size and appropriate
representativity? Finally, the results in terms of variance invite for a more
thorough investigation if the audit samples drawn using this method will also
be used to draw conclusions in terms of variance of estimates.

## References

  * [1] Alan Agresti  “Categorical data analysis”  John Wiley & Sons, 2013
  * [2] Timothy B Bell, Arnold M Bell and American Accounting Association  “Auditing practice, research, and education: A productive collaboration”  American Institute of Certified Public Accountants, 1995
  * [3] Yvonne M Bishop, Stephen E Fienberg and Paul W Holland  “Discrete Multivariate Analysis: Theory and Applications”  MIT Press, 1975
  * [4] J Chataway et al.  “Herpes simplex encephalitis: an audit of the use of laboratory diagnostic tests”  In _Qjm_ 97.6 Oxford University Press, 2004, pp. 325–330
  * [5] W G Cochran  “Sampling Techniques”  John Wiley & Sons, 1977
  * [6] Koen Derks, Jacques Swart, Eric-Jan Wagenmakers and Jan Wille  “JASP for Audit: Bayesian Tools for the Auditing Practice”  PsyArXiv, 2019
  * [7] Randal J Elder et al.  “Audit sampling research: A synthesis and implications for future research”  In _Auditing: A Journal of Practice & Theory_ 32.sp1 American Accounting Assocation, 2013, pp. 99–129
  * [8] Eurostat  “RAMON - Reference And Management Of Nomenclatures”, 2008  URL: https://ec.europa.eu/eurostat/ramon/nomenclatures/index.cfm?TargetUrl=LST_NOM_DTL&StrNom=NACE_REV2&StrLanguageCode=EN&IntPcKey=&StrLayoutCode=HIERARCHIC
  * [9] HM Hearnshaw et al.  “Are audits wasting resources by measuring the wrong things? A survey of methods used to select audit review criteria”  In _BMJ Quality & Safety_ 12.1 BMJ Publishing Group Ltd, 2003, pp. 24–28
  * [10] D Holt and T M F Smith  “Post Stratification”  In _Journal of the Royal Statistical Society, Series A_ 142.1, 1979, pp. 33–46
  * [11] Jonas Klingwort, Joep Burger and Bart Buelens  “Inferring network traffic from sensors without a sampling design”  CBS Discussion Paper, 2021
  * [12] Dennis Ougrin and Ricky Banarsee  “Clinical audit”  In _BMJ_ 333.7563 British Medical Journal Publishing Group, 2006, pp. s68–s69
  * [13] R Core Team  “R: A Language and Environment for Statistical Computing”, 2019  R Foundation for Statistical Computing  URL: https://www.R-project.org/
  * [14] Donald B Rubin  “Inference and missing data”  In _Biometrika_ 63.3 Oxford University Press, 1976, pp. 581–592
  * [15] S. Scholtus, B… Bakker and A. Delden  “Modelling Measurement Error to Estimate Bias in Administrative and Survey Variables” Discussion Paper 2015-17, Statistics Netherlands, The Hague, 2015
  * [16] Michael E Sobel and Gerhard Arminger  “Platonic and Operational True Scores in Covariance Structure Analysis: An Invited Comment on Bielby’s “Arbitrary Metrics in Multiple Indicator Models of Latent Variables””  In _Sociological Methods & Research_ 15.1-2 Sage Publications, 1986, pp. 44–58
  * [17] TheWorldBank  “Listed domestic companies”, 2019  URL: https://data.worldbank.org/indicator/CM.MKT.LDOM.NO
  * [18] N. Tollenaar et al.  “Predictieve textmining in politieregistraties: cyber- en gedigitaliseerde criminaliteit”  Rapport van het Wetenschappelijk Onderzoek- en Documentatiecentrum, Ministerie van Justitie en Veiligheid, 2018
  * [19] UN General Assembly  “Transforming our world: the 2030 Agenda for Sustainable Development”, 2015  URL: https://www.refworld.org/docid/57b6e3e44.html
  * [20] Bianca Zadrozny  “Learning and evaluating classifiers under sample selection bias”  In _Proceedings of the twenty-first international conference on Machine learning_ , 2004, pp. 114


# Generic construction of scale-invariantly coarse grained memory

Karthik H. Shankar  Center for Memory and Brain, Boston University

###### Abstract

Encoding temporal information from the recent past as spatially distributed
activations is essential in order for the entire recent past to be
simultaneously accessible. Any biological or synthetic agent that relies on
the past to predict/plan the future, would be endowed with such a spatially
distributed temporal memory. Simplistically, we would expect that resource
limitations would demand the memory system to store only the most useful
information for future prediction. For natural signals in real world which
show scale free temporal fluctuations, the predictive information encoded in
memory is maximal if the past information is scale invariantly coarse grained.
Here we examine the general mechanism to construct a scale invariantly coarse
grained memory system. Remarkably, the generic construction is equivalent to
encoding the linear combinations of Laplace transform of the past information
and their approximated inverses. This reveals a fundamental construction
constraint on memory networks that attempt to maximize predictive information
storage relevant to the natural world.

##  I Introduction

Representing the information from the recent past as transient activity
distributed over a network has been actively researched in biophysical as well
as purely computational domains Maass  _et al._ (2002); Jaeger (2001). It is
understood that recurrent connections in the network can keep the information
from distant past alive so that it can be recovered from the current state.
The memory capacity of these networks are generally measured in terms of the
accuracy of recovery of the past information Jaeger (2001); White  _et al._
(2004); Hermans and Schrauwen (2010). Although the memory capacity strongly
depends on the network’s topology and sparsity Ganguli  _et al._ (2008);
Strauss  _et al._ (2012); Legenstein and Maass (2007); Wallace  _et al._
(2013), it can be significantly increased by exploiting any prior knowledge of
the underlying structure of the encoded signal Ganguli and Sompolinsky (2012);
Charles  _et al._ (2014).

Our approach to encoding memory stems from a focus on its utility for future
prediction, rather than on the accuracy of recovering the past. In particular
we are interested in encoding time varying signals from the natural world into
memory so as to optimize future prediction. It is well known that most natural
signals exhibit scale free long range correlations Mandelbrot (1982); Voss and
Clarke (1975); West and Shlesinger (1990). By exploiting this intrinsic
structure underlying natural signals, prior work has shown that the predictive
information contained in a finite sized memory system can be maximized if the
past is encoded in a scale-invariantly coarse grained fashion Shankar and
Howard (2013). Each node in such a memory system would represent a coarse
grained average around a specific past moment, and the time window of coarse
graining linearly scales with the past timescale. Clearly the accuracy of
information recovery in such a memory system degrades more for more distant
past. In effect, the memory system sacrifices accuracy in order to represent
information from very distant past, scaling exponentially with the network
size Shankar and Howard (2013). The predictive advantage of such a memory
system comes from washing out non-predictive fluctuations from the distant
past, whose accurate representation would have served very little in
predicting the future. Arguably, in the natural world filled with scale-free
time varying signals, animals would have evolved to adopt such a memory system
conducive for future predictions. This is indeed evident from animal and human
behavioral studies that show that our memory for time involves scale invariant
errors which linearly scale with the target timescale Gibbon (1977); Rakitin
_et al._ (1998).

Our focus here is not to further emphasize the predictive advantage offered by
a scale invariantly coarse grained memory system, rather we simply assume the
utility of such a memory system and focus on the generic mechanism to
construct it. One way to mechanistically construct such a memory system is to
gradually encode information over real time as a Laplace transform of the past
and approximately invert it Shankar and Howard (2012). The central result in
this paper is that any mechanistic construction of such a memory system is
simply equivalent to encoding linear combinations of Laplace transformed past
and their approximate inverses. This result should lay strong constraints on
the connectivity structure of memory networks exhibiting the scale invariance
property.

We start with the basic requirement that different nodes in the memory system
represents coarse grained averages about different past moments. Irrespective
of the connectivity, the nodes can be linearly arranged to reflect their
monotonic relationship to the past time. Rather than considering a network
with a finite set of nodes, for analysis benefit, we consider a continuum
limit where the information from the past time is smoothly projected on a
spatial axis. The construction can later be discretized and implemented in a
network with finite nodes to represent past information from timescales that
exponentially scale with the network size.

##  II Scale Invariant Coarse Graining

Consider a real valued function ${f}(\tau)$ observed over time $\tau$. The aim
is to encode this time-varying function into a spatially distributed
representation in one dimension parametrized by $s$, such that at any moment
$\tau$ the entire past from $-\infty$ to $\tau$ is represented in a coarse
grained fashion as $\mathbf{T}(\tau,s)$

$\mathbf{T}(\tau,s)=\int_{-\infty}^{\tau}{f}(\tau^{\prime})\mathbf{W}(\tau-\tau^{\prime},s)\,\,d\tau^{\prime}.$
(1)

This is now a convolution memory model. The kernel
$\mathbf{W}{}(\tau-\tau^{\prime},s)$ is the coarse graining window function
with normalized area for all $s$,
$\int_{-\infty}^{\tau}\mathbf{W}(\tau-\tau^{\prime},s)d\tau^{\prime}=1$.
Different points on the spatial axis uniquely and monotonically represents
coarse grained averages about different instants in the past, as illustrated
in figure 1.

Figure 1: Coarse grained averages around different past instants are projected
on to different points on the spatial axis.

We require that coarse graining about any past instant linearly scales with
the past timescale. So, for any pair of points $s_{1}$ and $s_{2}$, there
exists a scaling constant $\alpha_{12}$ such that
$\mathbf{W}(\tau-\tau^{\prime},s_{1})=\alpha_{12}\mathbf{W}(\alpha_{12}(\tau-\tau^{\prime}),s_{2})$.
For the window function to satisfy this scale-invariance property, there
should exist a monotonic mapping $s(\alpha)$ from a scaling variable $\alpha$
to the spatial axis so that

$\mathbf{W}(\tau-\tau^{\prime},s(\alpha))=\alpha\mathbf{W}(\alpha(\tau-\tau^{\prime}),s(1)).$
(2)

Without loss of generality we shall pick $s(\alpha)=\alpha$ because it can be
retransformed to any other monotonic $s(\alpha)$ mapping after the analysis.
Hence with $0<s<\infty$,

$\mathbf{W}(\tau-\tau^{\prime},s)=s\mathbf{W}(s(\tau-\tau^{\prime}),1).$ (3)

##  III Space-Time Local mechanism

Equation 1 expresses the encoded memory as an integral over the entire past.
However, the encoding mechanism can only have access to the instantaneous
functional value of ${f}$ and its derivatives. The spatial pattern should self
sufficiently evolve in real time to encode eq. 1. This is a basic requirement
to mechanistically construct $\mathbf{T}(\tau,s)$ in real time using any
network architecture. Since the spatial axis is organized monotonically to
correspond to different past moments, only the local neighborhood of any point
would affect its time evolution. So we postulate that the most general
encoding mechanism that can yield eq. 1 is a space-time local mechanism given
by some differential equation for $\mathbf{T}(\tau,s)$. To analyze this, let
us first express the general space-time derivative of $\mathbf{T}(\tau,s)$ by
repeatedly differentiating eq. 1 w.r.t $\tau$ and $s$.

$\displaystyle\mathbf{T}_{(m)}^{[n]}(\tau,s)$ $\displaystyle=$
$\displaystyle\sum_{j=0}^{n-1}f^{[n-j-1]}(\tau)\mathbf{W}^{[j]}_{(m)}(0,s)$
(4) $\displaystyle+$
$\displaystyle\int_{-\infty}^{\tau}{f}(\tau^{\prime})\mathbf{W}^{[n]}_{(m)}(\tau-\tau^{\prime},s)\,\,d\tau^{\prime}.$

Here $n$ and $m$ are positive integers. For brevity, we denote the order of
time derivative within a square bracket in the superscript and the order of
space derivative within a parenthesis in the subscript.

Since ${f}(\tau)$ is an arbitrary input, $\mathbf{T}(\tau,s)$ should satisfy a
time-independent differential equation which can depend on instantaneous time
derivatives of ${f}(\tau)$. The first term in the r.h.s of eq. 4 is time-
local, while the second term involves an integral over the entire past. In
order for the second term to be time-local, it must be expressible in terms of
lower derivatives of $\mathbf{T}(\tau,s)$. Since the equation must hold for
any ${f}(\tau)$, $\mathbf{W}^{[n]}_{(m)}(\tau-\tau^{\prime},s)$ should satisfy
a linear equation.

$\sum_{n,m}C_{nm}(s)\mathbf{W}^{[n]}_{(m)}(\tau-\tau^{\prime},s)=0.$ (5)

The aim here is not to derive the time-local differential equation satisfied
by $\mathbf{T}(\tau,s)$, but just to impose its existence, which is achieved
by imposing eq. 5 for some set of functions $C_{nm}(s)$. To impose this
condition, let us first evaluate
$\mathbf{W}^{[n]}_{(m)}(\tau-\tau^{\prime},s)$ by exploiting the functional
form of the window function given by eq. 3. Defining $z\equiv
s(\tau-\tau^{\prime})$ and the function $G(z)\equiv\mathbf{W}(z,1)$, eq. 3 can
be repeatedly differentiated to obtain

$\mathbf{W}^{[n]}_{(m)}(\tau-\tau^{\prime},s)=\sum_{r=r_{o}}^{m}(n+1)!m!\frac{s^{n+1-m+r}}{r!(m-r)!^{2}}(\tau-\tau^{\prime})^{r}G^{[n+r]}(z),$

where $r_{o}=$ max[0, $m-n-1$] and the superscript on $G(z)$ represents the
order of the derivative w.r.t $z$. Now eq. 5 takes the form

$\sum_{n,m}C_{nm}(s)s^{n+1-m}\sum_{r=r_{o}}^{m}\frac{(n+1)!m!}{r!(m-r)!^{2}}z^{r}G^{[n+r]}(z)=0$
(6)

The above equation is not necessarily solvable for an arbitrary choice of
$C_{nm}(s)$. However, when it is is solvable, the separability of the
variables $s$ and $z$ implies that the above equation will be separable into a
set of linear differential equations for $G(z)$ with coefficients given by
integer powers of $z$. The general solution for $G(z)$ is then given by

$G(z)=\sum_{i,k}a_{ik}z^{k}e^{-b_{i}z},$ (7)

where $i$ and $k$ are non negative integers. The coefficients $a_{ik}$ and
$b_{i}$, and the functions $C_{nm}(s)$ cannot be independently chosen as they
are constrained through eq. 6. Once a set of $C_{nm}(s)$ is chosen
consistently with the coefficients $a_{ik}$ and $b_{i}$, the differential
equation satisfied by $\mathbf{T}(\tau,s)$ can be obtained by iteratively
substituting $\mathbf{W}^{[n]}_{(m)}(\tau-\tau^{\prime},s)$ (in the second
term of the r.h.s of eq. 4) in terms of its lower derivatives and replacing
the integral in terms of derivatives of $\mathbf{T}(\tau,s)$.

Here we shall neither focus on the choice of $C_{nm}(s)$ nor on the
differential equation for $\mathbf{T}(\tau,s)$ it yields. We shall only focus
on the set of possible window functions that can be constructed by a space-
time local mechanism. Hence it suffices to note from the above derivation that
the general form of such a window function is given by eq. 7. Since by
definition the window function at each $s$ coarse grains the input about some
past moment, we expect it to be non-oscillatory and hence restrict our focus
to real values of $b_{i}$. Further, the requirement of the window function to
have normalized area at all $s$ restricts $b_{i}$ to be positive.

##  IV Two step process

Let us consider the simplest window function, where only one of the
coefficients in the set of $a_{ik}$ and $b_{i}$ in eq. 7 are non-zero, namely
$b_{i}=b$ and $a_{ik}=b^{(k+1)}/k!$. We shall denote the corresponding window
function as $W_{\\{k,b\\}}$ to highlight its dependence on specific $k$ and
$b$. The most general window function is then simply a linear combination of
various $W_{\\{k,b\\}}$ for different values of $k$ and $b$. From eq. 7,
$W_{\\{k,b\\}}$ takes the form

$W_{\\{k,b\\}}(\tau-\tau^{\prime},s)=\frac{(bs)^{k+1}}{k!}(\tau-\tau^{\prime})^{k}e^{-bs(\tau-\tau^{\prime})}.$
(8)

It turns out that the differential equation satisfied by $\mathbf{T}(\tau,s)$
that generates this window function is simply first order in both space and
time given by

$\mathbf{T}_{(1)}^{[1]}(\tau,s)+bs\mathbf{T}_{(1)}^{[0]}(\tau,s)-\frac{(k+1)}{s}\mathbf{T}_{(0)}^{[1]}(\tau,s)=0,$
(9)

with a boundary condition $\mathbf{T}(\tau,\infty)={f}(\tau)$. This equation
can hence be evolved in real time by only relying on the instantaneous input
${f}(\tau)$ at each moment $\tau$.

For more complex window functions that are linear combinations of
$W_{\\{k,b\\}}$ for various $k$ and $b$, the order of the space and time
derivatives of $\mathbf{T}(\tau,s)$ involved in the differential equation are
not necessarily bounded when the parameters $k$ and $b$ involved in the linear
combinations of $W_{\\{k,b\\}}$ are bounded. So, it is not straight forward to
derive the mechanistic construction as a differential equation for
$\mathbf{T}(\tau,s)$. Hence the question now is, what is the mechanism to
construct a memory system with any window function?

Interestingly, there exists an alternative derivation of eq. 9 where the time
derivative and space derivative can be sequentially employed in a two step
process Shankar and Howard (2012). The first step is equivalent to encoding
the Laplace transform of the input ${f}(\tau)$ as $\mathbf{F}(\tau,s)$. The
second step is equivalent to approximately inverting the Laplace transformed
input to construct $\mathbf{T}(\tau,s)$.

${f}(\tau)\,\,\underrightarrow{\,\,\,Laplace\,\,\,}\,\,\mathbf{F}(\tau,s)\,\,\underrightarrow{\,\,\,Inverse\,Laplace\,\,\,}\,\,\mathbf{T}(\tau,s)$
$\displaystyle\mathbf{F}^{[1]}(\tau,s)$ $\displaystyle=$ $\displaystyle-
bs\mathbf{F}(\tau,s)+{f}(\tau),$ (10) $\displaystyle\mathbf{T}(\tau,s)$
$\displaystyle=$ $\displaystyle\frac{b}{k!}s^{k+1}\mathbf{F}_{(k)}(\tau,s).$
(11)

Taking ${f}(\tau)$ to be a function of bounded variation and
$\mathbf{F}(-\infty,s)=0$, eq. 10 can be integrated to see that
$\mathbf{F}(\tau,s)=\int_{-\infty}^{\tau}{f}(\tau^{\prime})e^{-bs(\tau-\tau^{\prime})}d\tau^{\prime}$.
Thus $\mathbf{F}(\tau,s)$ is the Laplace transform of the past input computed
over real time. Eq. 11 is an approximation to inverse Laplace transform
operation Post (1930). So $\mathbf{T}(\tau,s)$ essentially attempts to
reconstruct the past input, such that at any $s$,
$\mathbf{T}(\tau,s)\simeq{f}(\tau-k/bs)$. This reconstruction grows more
accurate as $k\rightarrow\infty$, and the input from each past moment is
reliably represented at specific spatial location. For finite $k$ however, the
reconstruction is fuzzy and each spatial location represents a coarse grained
average of inputs from past moments, as characterized by the window function
$W_{\\{k,b\\}}$. For further details, refer to Shankar and Howard (2012).

Since any window function is a linear combination of various $W_{\\{k,b\\}}$
for different values of $k$ and $b$, its construction is essentially
equivalent to linear combinations of the two step process given by equations
10 and 11.

| Combinations of $W_{\\{k,b\\}}$
---|---
1 | $W_{\\{100,2\\}}$
2 | $0.32W_{\\{65,1.3\\}}+0.27W_{\\{75,1.8\\}}$
| $+0.32W_{\\{75,1.26\\}}+0.09W_{\\{100,1.62\\}}$
3 | $W_{\\{8,0.16\\}}$
4 | $0.5W_{\\{35,1\\}}+0.5W_{\\{70,1\\}}$

Figure 2: For different combinations of $W_{\\{k,b\\}}$, the window functions
are plotted as a function of past time at the spatial point $s=50$.

The choice of the combinations of $W_{\\{k,b\\}}$ has strong implications on
the shape of the resulting window function. At any given $s$, $W_{\\{k,b\\}}$
is a unimodal function with a peak at $\tau-\tau^{\prime}=k/bs$ (see eq. 8).
Arbitrary combinations of $W_{\\{k,b\\}}$ could result in a spatial location
representing the coarse grained average about disjoint past moments, leading
to undesirable shapes of the window function. Hence the values of $k$ and $b$
should be carefully tuned. Figure 2 shows the window functions constructed
from four combinations of $b$ and $k$. The combinations are chosen such that
at the point $s=50$, the window function coarse grains around a past time of
$\tau^{\prime}-\tau\simeq-1$. The scale invariance property guarantees that
its shape remains identical at any other value of $s$ with a linear shift in
the coarse graining timescale. Comparing combinations 1 and 3, note that the
window function is narrower for larger $k$(=100) than for a smaller $k$(=8).
Combination 2 has been chosen to illustrate a plateau shaped window function
whose sides can be made arbitrarily vertical by fine tuning the combinations.
Combination 4 (dotted curve in fig. 2) illustrates that combining different
values of $k$ for the same $b$ will generally lead to a multimodal window
function which would be an undesirable feature.

##  V Discretized spatial axis

A memory system represented on a continuous spatial axis is not practical, so
the spatial axis should be discretized to finite points (nodes). The two step
process given by equations 10 and 11 is optimal for discretization
particularly when the nodes are picked from a geometric progression in the
values of $s$ Shankar and Howard (2013). Eq. 10 implies that the activity of
each node evolves independently of the others to construct
$\mathbf{F}(\tau,s)$ with real time input ${f}(\tau)$. This is achieved with
each node recurrently connected on to itself with an appropriate decay
constant of $bs$. Eq. 11 involves taking the spatial derivative of order $k$
which can be approximated by the discretized derivative requiring linear
combinations of activities from $k$ neighbors on either sides of any node. For
further details on implementation of the two step process on discretized
spatial axis, refer to Shankar and Howard (2013).

Figure 3: Window function $W_{\\{8,1\\}}$ at two points $s_{1}=6.72$ and
$s_{2}=2.59$ computed on a discretized spatial axis with $c=0.1$. The dotted
curves correspond to the window functions computed on the continuous spatial
axis ($c\rightarrow 0$).

By choosing the nodes along the $s$-axis from a geometric progression, the
error from the discretized spatial derivative will be uniformly spread over
all timescales, hence such a discretization is ideal to preserve scale-
invariance. Let us choose the $s$-values of successive nodes to have a ratio
$(1+c)$, where $c<1$. Figure 3 shows the window function $W_{\\{k,b\\}}$ with
$k=8$ and $b=1$ constructed from the discretized axis with $c=0.1$. The window
functions at two spatial points $s_{1}=6.72$ and $s_{2}=2.59$ are plotted to
illustrate that scale invariance is preserved after discretization. As a
comparison, the dotted curves are plotted to show the corresponding window
function constructed in the continuous $s$-axis (limit $c\rightarrow 0$). The
window function computed on the discretized axis is artificially scaled up so
that the solid and dotted curves in figure 3 are visually discernible. Note
that the discretized window function peaks farther in the past time and is
wider than the window function on the continuous spatial axis. As
$c\rightarrow 0$, the discretized window function converges on to the window
function constructed on the continuous axis, while for larger values of $c$
the discrepancy grows larger. Nevertheless, for any value of $c$, the
discretized window function always stays scale-invariant, as can be seen by
visually comparing the shapes of the window functions at $s_{1}$ and $s_{2}$
in figure 3. Now, it is straight forward to construct scale-invariant window
functions of different shapes by taking linear combinations of discretized
$W_{\\{k,b\\}}$, analogous to the construction in figure 2.

Implementing this construction on a discretized spatial axis as a neural
network has a tremendous resource conserving advantage. Since at each $s$, the
window function $W_{\\{k,b\\}}$ coarse grains the input around a past time of
$k/bs$, the maximum past timesscale represented by the memory system is
inversely related to minimum value of $s$. The geometric distribution of the
$s$ values on the discretized axis implies that if there are $N$ nodes
spanning the spatial axis for $\mathbf{T}(\tau,s)$, it can represent the
coarse grained past from timescales proportional to $(1+c)^{N}$. Hence
exponentially distant past can be represented in a coarse grained fashion with
linearly increasing resources.

##  VI Discussion and Conclusion

The formulation presented here starts from a convolution memory model (eq. 1)
and derives the form of the scale-invariant window functions (or the kernels)
that can be constructed from a space-time local mechanism. Interestingly, by
simply postulating a kernel of the form of eq. 7, Tank and Hopfield have
demonstrated the utility of such a memory system in temporal pattern
classification Tank and Hopfield (1987). In general, a convolution memory
model can adopt an arbitrary kernel, but it cannot be mechanistically
constructed from a space-time local differential equation, which means a
neural network implementation need not exist. However, the Gamma-memory model
Vries and Principe (1992) shows that linear combinations of Gamma kernels,
functionally similar to eq. 7, can indeed be mechanistically constructed from
a set of differential equations.

The construction presented here takes a complementary approach to the Gamma-
memory model by requiring scale invariance of the window function in the
forefront and then imposing a space-time local differential equation to derive
it. This sheds light on the connectivity between neighboring spatial units of
the network that is required to generate a scale invariant window function, as
described by the second part of the two step process (eq. 11). Moreover, the
linearity of the two step process and its equivalence to the Laplace and
Inverse Laplace transforms makes the memory representation analytically
tractable.

Theoretically, the utility of a scale invariantly coarse grained memory hinges
on the existence of scale free temporal fluctuations in the signals being
encoded Shankar and Howard (2013). Although detailed empirical analysis of
natural signals is needed to confirm this utility, preliminary analysis of
time series from sunspots and global temperature show that such a memory
system indeed has a higher predictive power than a conventional shift register
Shankar and Howard (2013). The predictive advantage of this memory system can
be understood as arising from its intrinsic ability to wash out non-predictive
stochastic fluctuations in the input signal from distant past and just
represent the predictively relevant information in a condensed form. Finally,
the most noteworthy feature is that a memory system with $N$ nodes can
represent information from exponentially past times proportional to
$(1+c)^{N}$. In comparison to a shift register with $N$ nodes which can
accurately represent a maximum past time scale proportional to $N$, this
memory system is exponentially resource conserving.

## Acknowldegements :

The work was partly funded by NSF BCS-1058937 and AFOSR FA9550-12-1-0369.

## References

  * Maass  _et al._ (2002) W. Maass, T. Natschläger, and H. Markram, Neural Computation 14, 2531 (2002).
  * Jaeger (2001) H. Jaeger,  _The echo state approach to analyzing and training recurrent networks_ , GMD-Report 148 (GMD - German National Research Institute for Information Technology, 2001).
  * White  _et al._ (2004) O. L. White, D. D. Lee, and H. Sompolinsky, Physical Review Letters 92, 148102 (2004).
  * Hermans and Schrauwen (2010) M. Hermans and B. Schrauwen, Neural Networks 23, 341 (2010).
  * Ganguli  _et al._ (2008) S. Ganguli, D. Huh, and H. Sompolinsky, Proceedings of the National Academy of Sciences of the United States of America 105, 18970 (2008).
  * Strauss  _et al._ (2012) T. Strauss, W. Wustlich, and R. Labahn, Neural Computation 24, 3246 (2012).
  * Legenstein and Maass (2007) R. Legenstein and W. Maass, Neural Networks 20, 323 (2007).
  * Wallace  _et al._ (2013) E. Wallace, H. R. Maei, and P. E. Latham, Neural Computation 25, 1408 (2013).
  * Ganguli and Sompolinsky (2012) S. Ganguli and H. Sompolinsky, Annual Review of Neuroscience 35, 485 (2012).
  * Charles  _et al._ (2014) A. S. Charles, H. L. Yap, and C. J. Rozell, Neural Computation 26, 1198 (2014).
  * Mandelbrot (1982) B. Mandelbrot,  _The Fractal Geometry of Nature_ (W. H. Freeman, San Fransisco, CA, 1982).
  * Voss and Clarke (1975) R. F. Voss and J. Clarke, Nature 258, 317 (1975).
  * West and Shlesinger (1990) B. J. West and M. F. Shlesinger, American Scientist 78, 40 (1990).
  * Shankar and Howard (2013) K. H. Shankar and M. W. Howard, Journal of Machine Learning Research 14, 3785 (2013).
  * Gibbon (1977) J. Gibbon, Psychological Review 84, 279 (1977).
  * Rakitin  _et al._ (1998) B. C. Rakitin, J. Gibbon, T. B. Penny, C. Malapani, S. C. Hinton, and W. H. Meck, Journal of Experimental Psychololgy: Animal Behavior Processes 24, 15 (1998).
  * Shankar and Howard (2012) K. H. Shankar and M. W. Howard, Neural Computation 24, 134 (2012).
  * Post (1930) E. Post, Transactions of the American Mathematical Society 32, 723 (1930).
  * Tank and Hopfield (1987) D. Tank and J. Hopfield, Proceedings of the National Academy of Sciences 84, 1896 (1987).
  * Vries and Principe (1992) B. d. Vries and J. C. Principe, Neural Networks 5, 565 (1992).