text:
\section{Introduction}
\label{sec:intr}
In studies of spatial variation in disease prevalence,
 it is often necessary to
combine information from multiple prevalence surveys. This is particularly the case
in  low-resource settings, where disease registries
typically do not exist. A methodological challenge in these circumstances
is that survey designs are
severely constrained by cost constraints.
The available surveys
may therefore
be of variable quality and/or conducted at different times. In this paper, we
propose a class of generalized linear geostatistical models (GLGMs) to address two
specific issues.The first
is variation in quality, for example between randomised and non-randomised surveys, in which case our proposed methodology  assumes that
at least one of the surveys provides an
unbiased ``gold-standard''. The second is variation in the underlying prevalence
when surveys are conducted at different times. In this case, by modelling the underlying prevalence over time we
are able to use data collected at all times to estimate the underlying prevalence surface at the specific time of interest,
typically the time of the most recent survey.

Methods for the combined analysis of data from multiple surveys
have previously used meta-analysis and small area statistics
approaches; see \citet*{statmatch2001}, \citet*{cancer2005}, \citet*{multframe2006}
 and \citet*{biasmod2009}. More recently,
\citet*{combsurv2011} used
 Bayesian hierarchical models to combine smoking prevalence estimates
from multiple surveys. They noted that commercial surveys are often
ignored in constructing official estimates
because of poor information about the sampling designs used, but argued that these surveys
can nevertheless
provide useful additional information because they are
more frequently updated than official surveys.

\citet*{raghu2007} noted the potential benefits that might
accrue from spatial modelling of multiple survey
data, but to the best of our knowledge,
explicit spatial modelling of biases and/or temporal variation
in the outcome of interest has not previously been addressed, except in a few
specific applications. For example, \citet*{wanji2012}
established a logit-linear calibration relationship between
estimates of {\it Loa loa} prevalence
in part of equatorial Africa based on
two different methods,
finger-prick
blood sampling and  a short questionnaire instrument.
\citet*{crainiceanu2008}  incorporated this
calibration relationship
into
a bivariate geostatistical model for the two corresponding prevalence maps.

As discussed in \citet*{biasmod2009}, if information from multiple surveys is to be
combined, it is important to understand
the limitations of their designs in order to take account of potential
biases in the associated estimates of prevalence. As a minimal condition, the
study subjects in each survey
should be drawn from the same target population.
One potential source of bias is that
some members of the target population may be less likely than others
to be included.
Convenience samples
provide an example of this. In resource-poor settings, the relatively low cost of
convenience sampling is tempting, but its potential to produce
 biased estimates is clear. In a non-spatial context,
  \citet*{hedt2011} propose a hybrid prevalence estimator
that combines information from randomised and convenience surveys.
They demonstrate that, with suitable adjustment for the bias, their hybrid estimator
can give better prevalence estimates than would be obtained by using only the data
from the randomised surveys.

A second source of heterogeneity amongst multiple prevalence
surveys is temporal variation in prevalence. When spatially referenced
prevalence surveys are repeated over time it
is usually of interest to estimate changes in prevalence over time.
When the outcomes from consecutive surveys are correlated, there is also a potential gain in
efficiency if comparisons are made through the use of a joint model. This is especially
advantageous when the surveys do not use the same set of sampling locations, because a joint
analysis can then exploit both the temporal and spatial correlation structure of the
combined data.

In Section \ref{sec:mod_form} of the paper we propose a class of generalised linear
geostatistical models (GLGMs) for the combined analysis
of data from multiple prevalence surveys. The model allows both for biased sampling and
temporal variation in prevalence provided that one of the surveys delivers
unbiased ``gold-standard'' estimates of prevalence.   In Section \ref{sec:inference}
we describe the methods that we use to fit the model.
In Section \ref{sec:simstudy}
we report the results
of simulation experiments that illustrate how a joint model leads to
gains in efficiency of estimation and spatial prediction.
In Section \ref{sec:application} we describe an application  to malaria prevalence data from three surveys conducted
in Chikhwawa District, Southern Malawi.
Section \ref{sec:discussion} is a concluding discussion.
All computations for the paper were run on the High End Computing Cluster
at Lancaster University, using the R software environment \citep*{rsoftware}.

\section{A multivariate generalized linear geostatistical model}
\label{sec:mod_form}

The ingredients of a univariate
GLGM are the following. Random variables $Y_j$ and explanatory variables
$d_j$
are associated with sampling locations
$x_j$ in a region of interest $A \subseteq \mathbb{R}^2$.  Each $d_j$ is a vector of
length $p \geq 1$.
Conditional on the realisation of a zero-mean
latent Gaussian process $S(x)$ and a set of mutually independent zero-mean
latent Gaussian variables $Z_j$, the $Y_j$ follow a classical
generalized linear model \citep*{GLM1989}, hence:
\begin{itemize}
\item[(i)] the $Y_j$ are mutually independent conditional on the $S(x_j)$
and $Z_j$, with conditional
expectations $\mu_j = m_j g^{-1}(\eta_j)$, where $m_j$ is a known
scalar and $g(\cdot)$ a known {\it link function};
\item[(ii)] $\eta_j = d_j^\top \beta + S(x_j) + Z_j$;
\item[(iii)] the conditional distribution of the $Y_j$ falls within the exponential family.
\end{itemize}
In the remainder of the paper, we assume that the conditional distributions in (iii)
are binomial, with the $y_j$ representing the number of positives amongst $m_j$ individuals
sampled at location $x_j$. We also adopt the  standard logistic link function,
$g(\mu) = \log \{\mu/(1-\mu)\}$, but other link functions could also be used.
We specify the Gaussian process $S(x)$ to have
covariance function $\ensuremath{\mathrm{Cov}}\{S(x),S(x')\}=\sigma^2\rho(x,x';\phi)$,
and the mutually independent $Z_j$ to have variance $\tau^2$.
The $Z_j$
have a dual interpretation as either non-spatial extra-binomial
variation or spatial variation at scales smaller  than the
smallest distance between sampling locations; the two interpretations
can only be disentangled unambiguously
if repeated measurements are taken at coincident locations. Finally, we write
$d_j = d(x_j)$  to
emphasise its spatial context.
\begin{figure}[ht]
\vspace{-0.2in}
\begin{center}
\includegraphics[scale=0.6]{dag.pdf}
\caption{Representation of the multivariate generalized linear geostatistical model \eqref{eq:GLGM} as a directed acyclic graph; $S_1$ and $S_2$ represent prevalences at
times $t_1$ and $t_2>t_1$; $B_2$ represents bias; $Y_1$, $Y_2$ and $Y_{2}^*$
are observed prevalences from unbiased surveys at times $t_1$ and $t_2$, and
a biased survey at time $t_2$, respectively. The target for predictive inference is
$S_2$. \label{fig:graph}}
\end{center}
\end{figure}

To extend the  model to accommodate multiple surveys taken at possibly different times, some of which may be biased,
let $i=1,\dots,r$ denote the index of the survey and $x_{ij}: j = 1,\ldots,n_{i}$
the corresponding set of sample locations. We replace the single process $S(x)$ by
a set of $r$ processes $S_i(x)$ which relate to the true prevalence at different
times. We assume that at least one of the surveys is known to be
unbiased, define ${\cal B}$ to be the
index set of the potentially biased surveys
and introduce an additional set of  latent
Gaussian processes $B_i(x): i \in {\cal B}$
to represent the spatially
varying  biases. Finally, we assume that data from different surveys
are generated by conditionally
independent univariate GLGMs, with link functions
\begin{eqnarray}
\label{eq:GLGM}
g_i(\mu_{ij}/m_{ij}) &=&  \eta_{ij} = d(x_{ij})^\top \beta_1 + S_{i}(x_{ij}) +
        I(i \in {\cal B})[B_i(x_{ij}) + d(x_{ij})^\top \beta_i]+  Z_{ij},\nonumber \\
            && j=1,\ldots,n_{i}; i=1,\ldots,r.
\end{eqnarray}
On the right-hand-side of (\ref{eq:GLGM}), we assume that the marginal properties of each
$S_i(x)$  are the same
as previously specified for $S(x)$, and add a set of cross-covariance
functions, $\ensuremath{\mathrm{Cov}}\{S_{i}(x),S_{i'}(x')\}=\sigma^2 \alpha_{ii'} \rho(x,x';\phi)$,
where $-1 < \alpha_{ii'} < 1$. The parameters $\alpha_{ii'}$
capture the temporal correlation between
the true prevalence surfaces at different times, hence if surveys $i$ and $i^\prime$
are taken at the same time, $S_i(x) = S_{i^\prime}(x)$ for all $x$
and $\alpha_{ii'}=1$. Note that if $r>2$, some combinations
of $\alpha_{ii'}$ result in a non-positive-definite variance matrix. If $r$ is small, this can be handled
by setting the likelihood to zero for all such combinations. When $r$ is large the issue can be
avoided
by imposing a spatio-temporally continuous parametric structure. This has the incidental benefit of  making the model more
parsimonious. One such example would be an exponentially decaying cross-covariance structure
 with $\alpha_{ii'} = \exp\{-|t_i - t_{i'}|/\psi\}$, where $t_i$ is the time at which the
$i$th survey is taken.

The processes $B_i(x)$ in (\ref{eq:GLGM})
are assumed to be independent, with
zero mean and covariance functions
$\ensuremath{\mathrm{Cov}}\{B_i(x),B_i(x')\}=\nu_i^2\rho(x,x';\delta_i)$.

Finally, the random variables $Z_{ij}$ are again assumed to be mutually independent and Normally distributed
with common mean 0 and variances $\tau_i^2$.

As already noted, when all surveys are taken at the same time, $S_i(x)=S_1(x)$ for all $i$, which formally
corresponds to $\alpha_{ii^\prime}=1$ for all $(i,i^\prime)$. When
all surveys are unbiased but are taken at different times, ${\cal B}$ is the empty set and
the terms $[B_i(x_{ij}) + d(x_{ij})^\top \beta_i]$ in (\ref{eq:GLGM}) are omitted; formally, this
 corresponds to $\nu_i^2=0:i=2,\ldots,r$. If it is appropriate to use different
explanatory variables to model the true prevalence and the bias, this is
accommodated by setting some elements of the $\beta_i$ to zero.
The dependence structure of the model is illustrated by the directed acyclic graph in Figure \ref{fig:graph} for the special case of two gold-standard surveys conducted at two different times and a biased survey at the second time period. This scenario corresponds to case study analysed in Section \ref{sec:application}, where the aim is predictive inference for $S_2(x)$. In this case, the potential gains in efficiency
by jointly modelling the  data from all  three surveys stem from the direct
links between  $S_2$ and both $Y_2$ and $Y_2^*$ and the indirect link
between $S_2$ and $Y_1$ via $S_1$.

\section{Inference}
\label{sec:inference}

In this section, we focus on the case $r=2$.
The generalization to more than two surveys
is straightforward. We set $B_1(x)=0$,  write $B(x)$ in place of $B_2(x)$ and
write the parameters of this bivariate version of (\ref{eq:GLGM})
as $\beta=(\beta_1,\beta_2)$ and
 $\theta = (\sigma^2, \nu^2,\tau_1^2,\tau_2^2,\phi,\delta,\alpha)$.

\subsection{Likelihood}
\label{subsec:lik}

Let $y_{i} = (y_{i1},\ldots,y_{in_{i}})^\top$ denote the outcome data from surveys
$i=1,2$ and let $D_i$ be the $n_i$ by $p$ matrix whose $j$th row contains the
values $d(x_{ij}) = (d_1(x_{ij}),\ldots,d_p(x_{ij}))^\top$. Similarly, let
$T_i$ denote the vector of the $n_i$ values of
the linear predictor for survey $i$,
hence $T_i = D_i \{\beta_1 + I(i=2)\beta_2\} + W_i$,
where $W_i = (W_{i1},\ldots,W_{in_i})^\top$
and
\begin{equation}
W_{ij} = S_{i}(x_{ij}) + I(i=2)B(x_{ij}) + Z_{ij}.
\label{eq:random_part}
\end{equation}
Now, let $T$ denote
the $(n_1+n_2)$-element vector $T = (T_1^\top,T_2^\top)^\top$ and
$D$ the $(n_1+n_2)$ by $2p$ matrix,
\begin{equation}
D = \left[
\begin{array}{ll}
D_1 & 0 \\
D_2 & D_2
\end{array}
\right].
\label{eq:mean}
\end{equation}
Also, write $R_{ii'}(\phi)$ for the $n_i$ by $n_{i'}$ matrix with
$(h,k)$th element $\rho(x_{ih},x_{i^\prime k};\phi)$ and
$R_b(\delta)$ for the $n_2$ by $n_2$ matrix with
$(h,k)$th element $\rho(x_{2h},x_{2k};\delta)$.
Then,
\begin{equation}
T \sim {\rm MVN}(D \beta, V(\theta))
\label{eq:MVN}
\end{equation}
where
\begin{equation}
V(\theta) =
\left[
\begin{array}{ll}
\sigma^2 R_{11}(\phi) + \tau_1^2 I & \sigma^2 \alpha R_{12}(\phi) \\
\sigma^2 \alpha R_{21}(\phi)       & \sigma^2 R_{22}(\phi) + \nu^2 R_b(\delta) + \tau_2^2 I
\end{array}
\right].
\label{eq:covmat}
\end{equation}
The conditional distribution of $Y$ given $T=t$ is a product of independent binomial
probability mass functions. We write this as
\begin{equation}
f(y|t) = \prod_{i=1}^{2}\prod_{j=1}^{n_{i}}f(y_{ij}|t_{ij}).
\label{eq:conditional}
\end{equation}
Combining (\ref{eq:mean}), (\ref{eq:MVN}), (\ref{eq:covmat}) and (\ref{eq:conditional}) then
gives the likelihood function as the  high-dimensional integral
\begin{equation}
\label{eq:likelihood}
L(\beta,\theta) = \int  h(t; D\beta, V(\theta))f(y|t) \: dt,
\end{equation}
where $h(\cdot|\mu,V)$ is the density function of a multivariate Normal
distribution with mean $\mu$ and covariance matrix $V$.

\subsection{Conditional simulation}
\label{subsec:simre}

We propose to use Monte Carlo methods to evaluate the high-dimensional
integral in (\ref{eq:likelihood}). These methods require us to  simulate
from the conditional distribution of the spatial random effect
$T$ given the data $Y=y$. Using Bayes' formula,
this conditional density is
\begin{equation}
\label{eq:conditional_dens}
\pi(t |y) \propto h(t|D \beta, V(\theta))f(y|t).
\end{equation}
To simulate from (\ref{eq:conditional_dens}),
\citet*{mcmc2006} propose a Langevin-Hastings (LH) Markov chain
Monte Carlo (MCMC) algorithm. This operates by updating a
linear transformation of $T$, chosen to make the components of $T|y$ approximately
independent.
\citet*{mcmc2006}
use a Gaussian approximation to the distribution of $T|y$, with
mean $D \beta$ and covariance matrix
\begin{equation}
\label{eq:sigmatilde}
\tilde{V} = \{V(\theta) +\Lambda(\hat{t})\}^{-1}.
\end{equation}
In (\ref{eq:sigmatilde}), $\Lambda(t)$ is a diagonal matrix with entries
$-\partial^2/\partial t_{i}^2 \log f(y|t)$ and
$\hat{t}$ is a typical value of $T$ such as
the mode of $f(y|t)$.  For
the binomial model with logistic link, this gives
 $\Lambda(\hat{t}) = \ensuremath{\mathrm{diag}}\{y_{i}(1-y_{i}/m_{i})\}$.
\citet*{mcmc2006} demonstrate that updating the centred random variable
$
\tilde{T} = \tilde{V}^{-1/2}(T-D \beta)
$
gives better mixing and convergence properties than the analogous
MCMC algorithms based on either $T$ or on $\bar{T} = V^{-1/2}(T-D \beta)$, as suggested by \citet*{waag2002}.

\subsection{Monte Carlo Maximum Likelihood: estimation and spatial prediction}
\label{subsec:mmle}

The Monte Carlo Maximum Likelihood (MCML) method \citep*{geyer1992, geyer1994, geyer1996, geyer1999}
uses conditional simulations of $T$ given $Y$ to obtain a computationally efficient
approximation to the intractable
likelihood function.
From \eqref{eq:likelihood}, the likelihood function can be written as
\begin{eqnarray}
\label{eq:likelihood2}
L({\beta},{\theta}) &=& \int  h({t}|D{\beta}, V({\theta}))f({y}| {t}) \: d{t}
                                             = \int  \frac{h({t}|D{\beta}, V({\theta}))f({y}| {t})}
                                                                {\tilde{f}({y},{t})}\tilde{f}({y},{t})\: d{t} \nonumber \\
                             & \propto &  \int \frac{h({t}|D{\beta}, V({\theta}))f({y}| {t})}
                                                          {\tilde{h}({t})f({y}| {t})}
                                                          \tilde{f}({t} | {y})\: d{t}
                                                          =
                                                          E_{\tilde{f}}\left[\frac{h({t}|
                                                          D{\beta}, V({\theta}))}{\tilde{h}({t})}\right].
\end{eqnarray}
In (\ref{eq:likelihood2}), $\tilde{f}({t}, {y}) = f({y} | {t}) \tilde{h}({t})$,
where $\tilde{h}({t})$ is any density function with support in $\mathbb{R}^{n}$,
and $E_{\tilde{f}}$ denotes
 expectation with respect to $\tilde{f}(\cdot | {y})$. MCML estimates are
then obtained by maximizing
\begin{equation}
\label{eq:mml}
L_{m}({\beta},{\theta}) = \frac{1}{m}\sum_{h=1}^{m}\frac{h({t}_{h}|
                                                          D{\beta}, V({\theta}))}{\tilde{h}({t}_{h})},
\end{equation}
where ${t}_{1},\ldots,{t}_{m}$  are samples from $\tilde{f}(\cdot | {y})$.

The accuracy of the approximation for a given value of $m$ depends critically on
the choice of $\tilde{h}(\cdot)$. A suitable choice is
$h({t}_{h}|D{\beta}_{0}, V({\theta}_{0}))$, where ${\beta}_{0}$ and ${\theta}_{0}$
are as close as possible to the maximum likelihood estimates,
$\hat{{\beta}}$ and $\hat{{\theta}}$. In practice, we embed the maximisation of
$L_m({\beta},{\theta})$ within the following
iterative procedure  as suggested in \citet*{geyer1992} and \citet*{geyer1994}:
 let $(\hat{\beta}_1,\hat{\theta}_1)$ denote the values that
maximise $L_m({\beta},{\theta})$ using an initial guess
at suitable values  $({\beta}_{0},{\theta}_{0})$; repeat the maximisation
with $(\hat{\beta}_1,\hat{\theta}_1)$ replacing $({\beta}_{0},{\theta}_{0})$; continue
until convergence.

For the numerical maximization of \eqref{eq:mml} we use a similar procedure to
the one presented in \citet*{mmle2004}. Write
 $V({\theta}) = \sigma^2 V({\psi})$ where
$\psi = ( \nu^2/\sigma^2,\tau_1^2/\sigma^2,\tau_2^2/\sigma^2,\phi,\delta,\alpha)^\top$
For a given value of ${\psi}$,
the first and second derivatives
of \eqref{eq:mml} with respect to ${\beta}$ and $\sigma^2$
are analytically tractable
and we use an iterative Newton-Raphson algorithm. We then plug into
\eqref{eq:mml}
the values $\hat{\beta}(\psi)$ and $\hat{\sigma}(\psi)^2$  and maximize
with respect to ${\psi}$ using direct numerical optimization with a
further re-parameterisation to remove any restrictions on the permissible ranges of
the parameters; we use a log-transformation
for all elements of $\psi$ except $\alpha$, for which we use
$\log\{(1+\alpha)/(1-\alpha)\}$ to correspond to the
range $-1 < \alpha < 1$. We also consider a variety of starting values
 to guard against false convergence to either a local maximum or an arbitrary point
on a
plateau of the likelihood surface.

We now consider the prediction of
${T}^* = (T(x_{n+1}),\ldots,T(x_{n+q}))^\top$ at $q$ additional
prediction locations that are not included in any of the prevalence
surveys. This requires all relevant explanatory variables to
be available at the prediction locations.  We include the mutually
independent random variables $Z_{ij}$ in (\ref{eq:random_part})
as part of our target for prediction.  Note that in a linear
Gaussian geostatistical model, the $Z_{ij}$ would be conflated with Normally
distributed  measurement errors, whereas in a GLGM for
prevalence survey data the
analogue of measurement
error is binomial sampling variation and is formally distinguishable
from the extra-binomial variation induced by the $Z_{ij}$.

\citet*{em2002} gives approximate expressions for
 the minimum mean square predictor ${\rm E}[{T}^*|{y}]$ and its variance
 using samples from the conditional distribution
of $T|y$ generated by conditional simulation. For prediction of non-linear functionals of
the prevalence surface,
we first use our MCMC algorithm to generate samples $t_h: h=1,...,m$  from the conditional distribution of $T|y$,
then simulate samples $t_h^*:h=1,...,m$ directly from the multivariate Normal
conditional distribution of ${T}^* |T= {t}_{h}$.  This has expectation
\begin{equation}
\label{eq:condpredmean}
D^*{\beta} + C^\top V^{-1}({t}_{h}-D{\beta}),
\end{equation}
where $D^*$ is the matrix of covariates at the prediction locations,
and covariance matrix
\begin{equation}
\label{eq:conpredvar}
V^* - C^\top V^{-1} C,
\end{equation}
where
 $V^*$ is the covariance matrix of ${T}^*$ and $C$ is the cross-covariance matrix between ${T}$ and ${T}^*$.
Finally,  we transform the sampled values $t_h^*$ to predicted prevalences,
$$p_h^* = g^{-1}({t}_{h}^*) = (g^{-1}(t_{n+1,h}^*),\ldots,g^{-1}(t_{n+q,h}^*))^\top,$$
where $g^{-1}(\cdot)$ is the inverse link function. Typically, the prediction locations will form a fine grid to cover the area of interest, $A$,
so as to approximate
a set of predicted surfaces, ${\cal P}^* = \{p_h^*(x): x \in A\}$ which can then be summarised according to the
needs of each application. For example, we might want to map pointwise means, or selected quantiles,
or predictive probabilities of the exceedance of policy-relevant thresholds.

\section{Simulation study}
\label{sec:simstudy}

We have conducted a simulation study of our proposed methodology
with three aims:  to show that the parameters in (\ref{eq:GLGM}) are identifable;
 to illustrate the
finite sample
properties of the MCML estimators; and to demonstrate the potential gains in predictive performance
that can be obtained by combining data from unbiased and biased surveys.

\subsection{Identifiability and finite sample properties}
\label{subsec:ident}

For this part of the simulation study we simulated data from two surveys, the first of which was unbiased, the
second biased. We specified the covariance structure of the model to correspond to the MCML estimates that were obtained in the analysis of malaria prevalence data
to be reported in Section \ref{sec:application}. We also
used the same sample sizes as in the malaria application, hence
 $n_{1} = 425$ (to correspond to the second of the two randomised surveys)
and $n_{2} = 249$ (to correspond to the convenience survey), and the same binomial denominators $m_{ij}$.
We did not use covariates but constant means $\beta_{1}$ for the first survey and $\beta_{1}+\beta_{2}$ for the second survey.
\begin{table}[hbtp]
\centering
\caption{Estimated means and relative biases (RB) of the MCML estimators for the covariance parameters,
and ordered eigenvalues (EV) of their correlation matrix under three scenarios.  \label{tab:ident}}
\vspace{0.6cm}
\begin{tabular}{lr|rrr|rrr|rrr}
   & & \multicolumn{3}{c|}{(1)} & \multicolumn{3}{c|}{(2)} & \multicolumn{3}{c}{(3)} \\
& True value & Mean & RB &  EV & Mean & RB & EV & Mean & RB & EV\\
  \hline
$\beta_1$ & 1.000 & 0.997 & -0.003 & 1.677 & 1.011 & 0.011 & 1.677 & 0.998 & -0.002 & 1.811 \\
 $\beta_2$ & -1.000 &  -1.011 & 0.011 & 1.287 & -1.013 & 0.013 & 1.425 & -0.980 & -0.020 & 1.472 \\
 $\sigma^2$ & 2.186 & 2.132 & -0.025 & 1.173 & 2.093 & -0.042 & 1.298 & 2.005 & -0.083 & 1.141 \\
 $\tau^2$ & 0.558 & 0.465 & -0.166 & 0.903 & 0.476 & -0.148 & 0.840 & 0.486 & -0.130 & 0.835 \\
  $\nu^2$ & 0.672 & 0.900 & 0.339 & 0.772 & 1.011 & 0.504 & 0.715 & 1.193 & 0.776 & 0.806 \\
 $\phi$ & 0.017 & 0.016 & -0.045 & 0.695 & 0.016 & -0.033 & 0.577 & 0.016 & -0.085 & 0.503 \\
  $\delta$ & 0.004 & 0.005 & 0.249 & 0.492 & 0.006 & 0.496 & 0.468 & 0.008 & 1.037 & 0.433 \\
   \hline
\end{tabular}

\end{table}
We generated the sampling locations for the unbiased survey as an independent random sample from
the uniform distribution in the rectangle $[34.700, 34.900]\times[-16.170, -15.880]$. The usefulness of
the data from the biased survey may depend on the degree of overlap between the two sets of sampling locations. For this reason we
generated the sampling locations for the biased survey from each of three inhomogeneous Poisson processes,
 with intensity $\lambda(x) = \exp\{-\|x-x_{0}\|/0.02\}$ and $x_{0}$ set to
each of the three following locations:  $(34.800, -16.025)$,  the center of the the rectangle;
 $(34.700, -16.170)$, the lower left corner of the rectangle; $(34.600, -16.315)$, a point outside the rectangle. Figure \ref{fig:ident} shows an example of simulated locations for the biased survey under each of these three scenarios.
\begin{figure}
\begin{center}
\vspace*{-0.5in}
\includegraphics[scale=0.6]{ident.pdf}
\caption{example of simulated locations from a biased survey under three different scenarios as defined in Section \ref{subsec:ident}; the dashed lines encompass the region within which locations of an unbiased survey are uniformly generated.\label{fig:ident}}
\end{center}
\end{figure}
For each simulation we computed the mean and relative bias of the MCML estimates of the covariance parameters
and the eigenvalues  of their correlation matrix, based on 1000 replications of each of the three scenarios.
The results are shown in Table \ref{tab:ident}.  The estimates
of $\beta_{1}$, $\beta_{2}$, $\sigma^2$, $\tau^2$ and $\phi$ are approximately unbiased
under all three scenarios whereas the estimates of $\nu^2$ and $\delta$,
which relate to the process $B(x)$, become increasingly biased as
the overlap between the two sampled  areas decreases.
Under all  three scenarios, the
smallest eigenvalue of the correlation matrix
corresponds to about $6\%$ of its total variation as measured by the sum of
the eigenvalues. Also, the off-diagonal elements of the
 correlation matrix, whose entries are never greater than
 $0.47$ in absolute value, which represents the correlation
between the estimates of $\tau^2$ and $\phi$  in the third scenario.

The overall conclusion from this part of the simulation study is that all
of the model parameters are identifiable, and that the parameter estimates
are approximately unbiased provided that there is a substantial
overlap in the spatial  coverage of the unbiased and biased surveys. This is
as expected,  because without such overlap the two surveys can only
estimate the properties of the sum, $S(x)+B(x)$, in the area covered
by the biased survey.

\subsection{Quality variation and temporal variation}
\label{subsec:tqv}

In this part of the simulation study we focus on predictive performance. Our main objective is to indicate to what extent the inclusion of data from
a biased survey can improve predictive inference, under circumstances similar to those that hold in
our malaria application. A secondary objective, as suggested by a reviewer,
is to demonstrate the unreliability of a naive analysis that ignores bias and temporal variation. We therefore
conducted three analyses of each simulated data-set as follows.
\begin{itemize}
\item Joint (J). The combined data are analysed using the bivariate GLGM as specified in Section \ref{sec:mod_form}.
\item First-survey-only (FSO). Only the data from the first, unbiased survey are used.
\item ``Naive'' (N). The data from the two surveys are analysed using a GLGM that does not account for bias or temporal variation.
\end{itemize}

We consider a {\it quality variation} (QV) scenario, in which one survey is unbiased and
the other biased, and a {\it temporal variation} (TV) scenario, in which both surveys
are unbiased but at different times, with predictions required for the first time period.

The following features are common to both scenarios. The processes
$S_i(x):i=1,2$ have mean $\beta_1=1$,
variance $\sigma^2=1$ and correlation function
$\rho(u) = \exp(-u/\phi)$ with $\phi=0.15$.
 Locations of unbiased surveys are uniformly generated in the unit square centred on $x_{0} = (1/2, 1/2)$. Both surveys have the same number of
sampling locations, $n_{1}=n_{2}=300$.
The binomial denominators
at each sampling location are all set equal to 1.
Our primary focus is on prediction of prevalence at
$x_0$ but we also consider estimation of the parameters
$\beta_{1}$, $\log \sigma^2$
and  $\log \phi$ that define the model for the underlying prevalence process $S_1(x)$.

In the QV scenario, $S_1(x)=S_2(x)$ for all $x$ and the process $B_2(x)$ has mean $\beta_2=-1$ and correlation function
$\rho(u) = \exp(-u/\delta)$ with $\delta=0.15$.
Locations from the biased survey are generated from a Poisson process with intensity $\lambda(x) = \exp\{-\|x-x_{0}\|/0.15\}$ so that points closer to $x_{0}$ are more likely to be sampled, as might occur when using a convenience sampling strategy and
$x_0$ is the location of a health-care facility. Finally,
we consider four values, $\nu^2=0.5, 1, 2, 4 $, for the variance of the
process $B_2(x)$,  corresponding to increasingly severe
 spatial variation in the bias.

In the TV scenario, the cross-correlation function between
$S_1(x)$ and $S_2(x)$ is $\alpha \exp(-u/\phi)$. We consider
three values, $\alpha=0.2,0.5,0.8$, to correspond to weak, moderate
and strong correlation between the two prevalence surfaces.

 The results are summarised in Tables \ref{tab:qv} and \ref{tab:tv}.  These show
 estimates of the root-mean-square-error (RMSE)
 and coverage of nominal $95\%$ confidence intervals
(CIC) for MCML estimates of the parameters $\beta_{1}$, $\log \sigma^2$
and  $\log \phi$,
and for the minimum mean square error predictors of
$S_{1}(x_{0})$ and $\beta_{1}+S_{1}(x_{0})$. Each
entry is calculated from 1000 independent replicates of the simulation
model.

\begin{table}[htp]
\caption{Estimated RMSE, bias, SD and $95\%$ CIC for the MCML estimates of $\beta_{1}$, $\log\sigma^2$, $\log\phi$, for the minimum mean square error predictor of $S_{1}(x_{0})$ at location $x_{0}$ and $\beta_{1}+S(x_{0})$, under QV scenarios. \label{tab:qv}}
\centering
\small
\vspace{0.6cm}
\begin{tabular}{rcrrrr|rrrr}
Model & Parameter & \multicolumn{4}{c}{RMSE} &  \multicolumn{4}{c}{CIC}\\
   & & \multicolumn{4}{c|}{$\nu^2$} &\multicolumn{4}{c}{$\nu^2$} \\
& & 0.5 & 1 & 2 & 4 & 0.5 & 1 & 2 & 4 \\
  \hline
J & & 0.37 & 0.36 & 0.36 & 0.36 & 0.95 & 0.95 & 0.95 & 0.94 \\
  FSO & $\beta_{1}$& 0.37 & 0.36 & 0.35 & 0.35 & 0.95 & 0.95 & 0.95 & 0.95 \\
  N & & 0.50 & 0.52 & 0.51 & 0.48 & 0.82 & 0.79 & 0.83 & 0.77 \\
  \hline
  J & & 0.94 & 0.60 & 0.63 & 1.33 & 0.99 & 0.94 & 0.95 & 0.99 \\
  FSO & $\log\sigma^2$& 1.08 & 1.14 & 1.04 & 0.99 & 0.97 & 0.97 & 0.97 & 0.97 \\
N & & 0.78 & 0.49 & 0.52 & 0.54 & 0.99 & 0.96 & 0.95 & 0.91 \\
  \hline
  J & & 0.84 & 0.81 & 0.98 & 0.95 & 0.92 & 0.90 & 0.92 & 0.92 \\
 FSO & $\log\phi$& 1.45 & 1.42 & 1.32 & 1.44 & 0.94 & 0.94 & 0.92 & 0.94 \\
  N & & 0.69 & 0.68 & 0.66 & 0.70 & 0.92 & 0.90 & 0.88 & 0.83 \\
  \hline
 J & & 0.80 & 0.79 & 0.88 & 0.86 & 0.95 & 0.95 & 0.95 & 0.95 \\
  FSO & $S_{1}(x_{0})$& 0.91 & 0.86 & 0.92 & 0.85 & 0.95 & 0.95 & 0.95 & 0.94 \\
  N & & 0.93 & 0.96 & 1.13 & 1.21 & 0.80 & 0.80 & 0.77 & 0.76 \\
  \hline
    J & & 0.72 & 0.73 & 0.82 & 0.83 & 0.95 & 0.95 & 0.95 & 0.95 \\
  FSO & $\beta_{1}+S_{1}(x_{0})$ & 0.83 & 0.80 & 0.84 & 0.81 & 0.95 & 0.94 & 0.95 & 0.94 \\
  N & & 1.10 & 1.15 & 1.33 & 1.35 & 0.78 & 0.78 & 0.75 & 0.74 \\
   \hline
\end{tabular}
\end{table}
\begin{table}[htp]
\caption{Estimated RMSE, bias, SD and $95\%$ CIC for the MCML estimates of $\beta_{1}$, $\log\sigma^2$, $\log\phi$, for the minimum mean square error predictor of $S_{1}(x_{0})$ at location $x_{0}$ and $\beta_{1}+S(x_{0})$, under TV scenarios. \label{tab:tv}}
\centering
\small
\vspace{0.6cm}
\begin{tabular}{rcrrr|rrr}
Model & Parameter & \multicolumn{3}{c}{RMSE} & \multicolumn{3}{c}{CIC} \\
   & & \multicolumn{3}{c|}{$\alpha$} & \multicolumn{3}{c}{$\alpha$} \\
   & & 0.2 & 0.5 & 0.8 & 0.2 & 0.5 & 0.8 \\
  \hline
J & & 0.35 & 0.35 & 0.35 & 0.95 & 0.93 & 0.94 \\
 FSO & $\beta_{1}$& 0.36 & 0.36 & 0.35 & 0.94 & 0.93 & 0.93 \\
  N & & 0.63 & 0.63 & 0.64 & 0.29 & 0.38 & 0.42 \\
  \hline
  J & & 0.60 & 0.68 & 0.69 & 0.94 & 0.94 & 0.93 \\
  FSO & $\log\sigma^2$& 1.04 & 1.46 & 1.47 & 0.97 & 0.96 & 0.96 \\
  N & & 1.18 & 0.95 & 0.87 & 0.91 & 0.92 & 0.95 \\
  \hline
  J & & 0.93 & 0.92 & 0.91 & 0.91 & 0.93 & 0.93 \\
  FSO & $\log\phi$ & 1.47 & 1.56 & 1.75 & 0.92 & 0.93 & 0.94 \\
  N & & 1.32 & 1.09 & 1.02 & 0.92 & 0.94 & 0.94 \\
  \hline
  J & & 1.37 & 1.30 & 1.28 & 0.95 & 0.95 & 0.94 \\
  FSO & $S_{1}(x_{0})$& 1.37 & 1.30 & 1.31 & 0.94 & 0.93 & 0.92 \\
  N & & 1.39 & 1.32 & 1.28 & 0.87 & 0.91 & 0.92 \\
  \hline
    J & & 1.34 & 1.28 & 1.26 & 0.96 & 0.96 & 0.94 \\
  FSO & $\beta_{1}+S_{1}(x_{0})$ & 1.35 & 1.27 & 1.27  & 0.95 & 0.95 & 0.94 \\
  N & & 1.50 & 1.39 & 1.36 & 0.86 & 0.90 & 0.91 \\
   \hline
\end{tabular}
\end{table}

Overall, J outperforms FSO, which in turn outperforms N. Under the QV scenario, the main benefits of J are in the prediction of $S_{1}(x_{0})$ for values of $\nu^2$ smaller than $4$. The N approach yields much higher values of RMSE for the estimates of $\beta_{1}$, $S_{1}(x_{0})$ and $\beta_{1}+S_{1}(x_{0})$ and
very poor CIC.
 Under the TV scenario, the biggest gains achieved by J over FSO are in estimating
 the parameters $\log \sigma^2$ and $\log \phi$. Both J and FSO
perform similarly with respect to prediction of $S_{1}(x_{0})$ and $\beta_{1}+S_{1}(x_{0})$. The N approach, which in this scenario
consists of combining the data
under the assumption that  $S_{1}(x)=S_{2}(x)$ for all $x$,
i.e. $\alpha=1$, has the worst performance.

\section{Application: malaria prevalence mapping}
\label{sec:application}

 In this Section, we
use our proposed methodology to construct malaria prevalence maps
for an area of Malawi by combining information from three surveys.
All three surveys
were directed at the same target population, covering a 400 square km area
within Chikhwawa District, Southern Malawi.
 Two of the surveys were ``rolling'' Malaria
Indicator Surveys (MIS) \citep*{roca2012},
that used two different  practical strategies to obtain random, and
therefore unbiased, samples
from the population at risk.
The third was a facility-based  survey  that used
 a convenience sampling strategy, in which
 recruitment took place at a central child-vaccination
clinic at the main hospital in the centre of the study area.
We refer to this as the Easy Access Group (EAG)
study.

\subsection{Data}
\label{subsec:data}

\begin{figure}
\begin{center}
\includegraphics[scale=0.52]{coords_rMIS.pdf}
\includegraphics[scale=0.52]{coords_eMIS.pdf}
\includegraphics[scale=0.52]{coords_EAG.pdf}
\caption{Sampled locations for (a) rMIS, (b) eMIS and (c) EAG.
Coordinates are of individual houses  for rMIS and eMIS,
and of the villages for EAG; in (c), the radius of each circle is proportional to the number of the sampled households from the respective village and the black solid point corresponds to CDH village. The light blue lines represent waterways, with the thicker line corresponding to Shire river. \label{fig:coords}}
\end{center}
\end{figure}
Two population-level continuous malaria indicator surveys were conducted over the
period May
2010 to April 2012. Both surveys recruited
children aged less than five years in a sample of 50 village communities in order to
monitor the malaria intervention coverage and childhood burden of malaria in
a designated area
containing the sampled villages, which was chosen to represent the catchment area of
the Chikhwawa District Hospital (CDH). The two surveys differed in
the  sampling strategy used, as described below. We
refer to these two surveys  as the rMIS, covering the period May 2010 to April
2011 and the eMIS, covering the period May 2011 to April 2012. Throughout the two-year period
seven or eight villages were randomly selected
per month so as to sample all 50 villages twice yearly, once during the high-transmission season and once
during the low-transmission season. Within sampled villages, selection of households
was as follows. In the
rMIS, households were randomly selected within each village from a list of households, with
sampling probability proportional to village population size, based on a population enumeration
exercise. In the eMIS, a more economical
 ``spin-the-bottle'' method was used to identify a random set of households within villages. A
bottle was placed in the center of a village and used
to select random directions. A virtual line was drawn in each chosen direction to the border of the
village, the households that intersected this line were counted, and from
these a random household number was chosen as the starting point. The number of houses selected within each village was proportional to the estimated village population size.
Figures \ref{fig:coords} (a)-(b) show the sampled locations for the rMIS and the eMIS.

The third survey is a continuous facility-based MIS in children attending the immunization clinic at the CDH, conducted from May 2011 to April 2012. The objective of this study was to determine if estimates of uptake of control interventions and the burden of malaria from convenience sampling
were comparable to those from a randomised MIS conducted
 within the same catchment area of CDH. Children from 3 months of age who attended the vaccination clinic, and any accompanying sibling below 5 years, were recruited.
Between  30 and 50 children were recruited per month. Village of origin was extracted by direct questioning. If the village was not one of the 50 eMIS/rMIS villages for which the location was already known, its coordinates were determined retrospectively. The results for villages within 15km of CDH were extracted to make the
catchment area of the EAG comparable to that of the rMIS and eMIS.
Malaria control efforts by the national control program during the first period included a
district-wide household indoor residual spraying campaign between February and April 2011.
Practical difficulties resulted in this campaign
being conducted at the end of the rainy season rather than, as would have been
ideal, before the
start of the rainy season. This will have reduced its potential impact. Insecticide-treated net
control efforts were stable over the three months of the campaign, with distribution to women attending antenatal clinics
and mother and child clinics.

\subsection{Results}
\label{subsec:results}

The response from each child was a binary indicator of the outcome of
a rapid diagnostic test (RDT) used to test for the presence of malaria from a finger-prick blood sample.
Six
explanatory variables were considered, as
defined in Table \ref{tab:explanatory_variables}. Socio-Economic-Status (SES), an indicator of household wealth taking discrete values from 1 (poor) to 5 (wealthy),
was derived by an application of principal component analysis as discussed in \citet*{vyas2006}.

\begin{table}[htp]
\caption{Explanatory variables used in the analysis of the Chikhwawa malaria prevalence surveys}
\begin{center}
\begin{tabular}{ll} \hline
1 & intercept \\
2 & at least one treated bed-net in the household (yes/no)\\
3 & indoor residual spraying in the past two months (yes/no) \\
4 & high-transmission season (January-June/July-December) \\
5 & distance from the closest waterway (km) \\
6 & Socio-Economic-Status (SES, 1 to 5)\\ \hline
\end{tabular}
\label{tab:explanatory_variables}
\end{center}
\end{table}

It was thought  that health facility utilization might be  associated with SES
as previously observed in \cite{gahutu2011}, where children with relatively
high SES were  more likely to attend
a CDH.  Table \ref{tab:SES} shows the average SES observed in each
of our  three surveys. Enrolled children in
the EAG study show a higher average SES then those in the two other surveys. Additionally, Table \ref{tab:SESandRDT} shows that the relationship
between SES and the distribution of the number of  RDT positive results
per household differs
 between the two gold-standard surveys and the convenience survey.
We therefore allowed SES to have a direct effect on the spatially
structured bias  of the EAG survey in addition to its possible association with prevalence.
\begin{table}[ht]
\centering
\tabcolsep = 0.1cm
\caption{Mean and standard deviation (SD) of SES in the three surveys.\label{tab:SES}}
\vspace{0.5cm}
\begin{tabular}{rrrr}
  \hline
  & \multicolumn{3}{c}{SES}\\
 & rMIS & eMIS & EAG\\
  \hline
Mean & 2.76 & 2.50 & 3.45\\
  SD & 1.45 & 1.37 & 1.39 \\
   \hline
\end{tabular}
\vspace{1cm}
\caption{Distribution (percentage) of the number of positive RDTs
per household for each value of  SES, in the
convenience survey (EAG, left-columns) and in
the gold-standard surveys (rMIS and eMIS, right columns)
 \label{tab:SESandRDT}}
\vspace{0.5cm}
\begin{tabular}{crrrrrr|rrrrr}
 & & \multicolumn{5}{c}{SES (EAG)} & \multicolumn{5}{c}{SES (rMIS and eMIS)} \\
 &  & 1 & 2 & 3 & 4 & 5 & 1 & 2 & 3 & 4 & 5 \\
  \hline
  & 0 &75.76 & 80.56 & 77.50 & 79.10 & 89.04 & 54.58 & 63.09 & 71.72 & 73.29 & 83.74\\
 RDT & 1 & 21.21 & 19.44 & 20.00 & 20.90 & 10.96 & 40.49 & 33.56 & 25.25 & 22.60 & 15.45 \\
  positives& 2 & 3.03 & 0.00 & 2.50 & 0.00 & 0.00 & 4.58 & 3.35 & 3.03 & 3.42 & 0.81\\
  & 3 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.35 & 0.00 & 0.00 & 0.69 & 0.00  \\
  \hline
\end{tabular}
\end{table}

\begin{table}[ht]
\caption{Monte Carlo maximum likelihood estimates and $95\%$ confidence intervals. \label{tab:estim}}
\begin{center}
\begin{tabular}{lrr}
  \hline
 Term & Estimate & 95\% confidence interval \\
  \hline
$\beta_{1}$ & -0.272 & (-1.382,  0.862) \\
$\beta_{2}$ & -0.439 & (-0.623, -0.277) \\
$\beta_{3}$ & -0.399 & (-0.621, -0.189)\\
$\beta_{4}$ & 0.415 & (0.206, 0.598)\\
$\beta_{5}$ & -0.373 & (-0.970,  0.116) \\
$\beta_{6}$ & -0.151 & (-0.233, -0.072)\\
$\beta_{7}$ & -0.096 & (-0.222,  0.021) \\
$\sigma^2$ & 2.186 & (0.955,  3.155)\\
$\tau^2$ & 0.558 & (0.089,  1.231) \\
$\nu^2$ & 0.672 & (0.525,  0.802)\\
$\alpha$ & 0.859 & (0.483,  0.924)\\
$\phi$ & 0.017 & (0.006,  0.032) \\
$\delta$ & 0.004 & (0.001,  0.025)\\
   \hline
\end{tabular}
\end{center}
\end{table}

The resulting model for the combined data therefore included
seven regression parameters, $\beta_1, \beta_2, ... , \beta_7$. Let
 $\beta^\top = (\beta_{1}, \dots, \beta_{6})$ and denote
by  $d(x_{ij})$  the vector of covariates associated with location $x_{ij}$.
Use $i=1,2,3$ to denote
rMIS, eMIS and EAG, respectively. Then, the linear predictor is
\begin{eqnarray*}
\eta_{ij} &=& d(x_{ij})^\top\beta + S_{i}(x_{ij}) + I(i = 3)[B(x_{ij})+\beta_{7}\text{SES}_{ij}]+Z_{ij}, \\\
&& i=1,\ldots, 3; j = 1,\ldots, n_{i},
\end{eqnarray*}
where $n_{1} = 475$, $n_{2}=425$ and $n_{3}=249$.
Note that in the joint model for $S_1(x)$, $S_2(x)$ and
$S_3(x)$, $\alpha_{23} = 1$ because the EAG study took place
over the same period as the eMIS.
We therefore use $\alpha$ to denote $\alpha_{12}$ and set $S_{3}(x) = S_{2}(x)$.
We also assume equal variances $\tau^2$
 for the nugget term $Z$ across all three
surveys. Finally, we define $\ensuremath{\mathrm{Cov}}\{S_{1}(x), S_{2}(x')\} = \sigma^2\alpha\exp\{-\|x-x'\|/\phi)$ where $\sigma^2 > 0$, $\phi > 0$ and $-1 < \alpha < 1$. \par
Table \ref{tab:estim} shows the Monte Carlo maximum likelihood estimates
of the model parameters
together with $95\%$ confidence intervals. Each evaluation of the log-likelihood used $5000$ simulated values, obtained by conditional simulation of $110000$ values
and sampling every $20$th realization after
discarding a burn-in of $10000$ values. Figure \ref{fig:diag} shows two diagnostic plots for the average random effect: convergence of the MCMC algorithm appears to be
satisfactory.

The confidence intervals in Table \ref{tab:estim} were
calculated using the following parametric bootstrap procedure.
\begin{figure}[h]
\begin{center}
\includegraphics[scale=0.48]{diag.pdf}
\caption{Diagnostic plots for the convergence of the posterior average of the random effect. Left panel: correlogram of the 5000 simulated values. Right panel: empirical cumulative density function for the first (black line) and second (red line) 2500 simulated values. \label{fig:diag}}
\end{center}
\end{figure}
Using the parameter estimates in Table \ref{tab:estim} we simulated $1000$ data-sets
 from the model,
applied to each simulated data set the Monte Carlo maximum likelihood method
with $5000$ conditional simulations, and
 computed the empirical quantiles of the $1000$
resulting estimates of each parameter.
Although this procedure introduces additional Monte Carlo error, it allows us to compute confidence intervals without relying on
questionable Normal approximations for the distribution of
the Monte Carlo maximum likelihood estimates. \par
From Table \ref{tab:estim}, we see
 that the ownership of at least one treated bed net, the presence of residual indoor spraying and an increase in SES are all associated with a reduction in the prevalence of a positive RDT. The distance from the closest waterway is not significant, although the sign of the regression coefficient suggests
 that prevalence  decreases with increasing distance. The period January to June,
which is known to be a period of high malaria transmission,
 is associated with a significant increase in prevalence, by an estimated factor of
$\exp(0.415) \approx 1.51$.

\begin{figure}[h]
\begin{center}
\includegraphics[scale=0.91]{bias.pdf}
\end{center}
\caption{Predicted bias surface $B(x)$ (a) without the interaction term of SES and (b) including the effect of SES on the spatial bias. \label{fig:bias}}
\end{figure}

The regression coefficient $\beta_7$, which represents the
additional effect of SES on the bias of the EAG data,
is not significant, but its inclusion nevertheless
makes a noticeable difference to the
 predicted bias surface.  Figures \ref{fig:bias}(a) and \ref{fig:bias}(b) show
the minimum mean square error predictions of the bias with and without
including the regression on SES.

\begin{figure}[h]
\centerline{\includegraphics[scale=0.7]{pred.png}}
\caption{Predictions of (a) $d(x_{2j})^\top\beta$ and (b) $d(x_{2j})^\top\beta + S_{2}(x_{2j})$ at observed locations; (c) predicted surface of the unexplained spatial variation $S_{2}(x)$. The same colour scale has been used for the point predictions.\label{fig:pred}}
\end{figure}

The estimate $\hat{\alpha}=0.859$, albeit with a wide confidence interval,
 indicates a
strong correlation between prevalences in the two sampling
periods, 2010-2011 and 2011-2012.

 Figures \ref{fig:pred}a and b show
 the contributions of the linear regression and of the unexplained spatial variation
to the predicted  log-odds of prevalence at each of the observed locations. Figure \ref{fig:pred}c shows the unexplained component, $\hat{S}(x)$,
of the predicted prevalence as a spatially continuous surface.  The clear and
substantial difference between adjacent areas to the east and west of the river Shire strongly suggests the existence of one, or
more, social or environmental risk-factors that are not captured
by the available explanatory variables.

Figure \ref{fig:std_errors} shows pairwise scatter plots to compare the
prediction standard deviations for $S_{2}(x)$ at the sampling locations.
Figure \ref{fig:std_errors} (a) shows that analysing rMIS and eMIS data in the joint model
for temporal variation results in substantially better precision than using only the eMIS;
Figures \ref{fig:std_errors} (c) and (d)
show the further, but more modest, gains resulting from addition of the
data from the EAG; in contrast, Figure \ref{fig:std_errors} (b) suggests little or no benefit
from adding the EAG data to the eMIS data, with predictive standard deviations decreasing at
some locations but increasing at others.

\begin{figure}[h]
\begin{center}
\includegraphics[scale=0.8]{std_errors.pdf}
\end{center}
\caption{Scatter plots of the prediction standard errors for $S_{2}(x)$
at sampled locations $x$, using models fitted to the data from: (a) eMIS against eMIS and rMIS; (b) eMIS against eMIS and EAG; (c) eMIS, rMIS and EAG against eMIS; (d) eMIS, rMIS and EAG against eMIS and rMIS.\label{fig:std_errors}. The solid line represents the identity line.}
\end{figure}

\section{Discussion}
\label{sec:discussion}

We have developed a class of multivariate GLGMs for the combined data from
multiple spatially referenced surveys, and
associated Monte Carlo methods for
 maximum likelihood estimation and spatial prediction within the
proposed class of models.

The model as defined by (\ref{eq:GLGM}) is
the minimally parameterised model that captures the essential features of our motivating
application: variation in
data-quality arising from non-randomised sampling;
variation in prevalence over time; binomial and extra-binomial
sampling variation.  We have shown that
all of the
model parameters are identifiable from surveys of comparable size to the ones available
to us for the application. If substantially
larger data-sets were available, it would be of interest to extend the model in various
ways, for example by relaxing the assumption of common parameters for
 the prevalence surfaces $S_i(x)$ at different times or by allowing cross-correlation
between the $S_i(x)$ and their paired bias surfaces $B_i(x)$.
Additionally, if a large number of surveys were conducted at irregularly spaced
time-points within partly overlapping time periods, the use of a structured spatio-temporally
continuous process $S(x,t)$, as mentioned in Section \ref{sec:mod_form}, would
be more appealing than a discrete set of processes $S_i(x)$ at specific times $t_i$.

The Monte Carlo maximum likelihood
estimation procedure is computationally intensive, primarily because of the need to
use parametric bootstrapping to compute standard errors reliably. For this reason, we are currently developing a much faster
Monte Carlo method for approximate evaluation of the likelihood function.

In our application to malaria prevalence surveys, we
combined data from three surveys, two of which
were unbiased and conducted in two consecutive years,whilst the third
was a potentially biased convenience survey conducted over the same time-period as the
second unbiased survey. We obtained substantial gains in the precision of
spatial predictions by combining the data from the two unbiased surveys
and further, but smaller, gains from combining the data from all three surveys.

One of the limitations of our approach is that it assumes that
at least one of the available surveys represents an unbiased gold-standard. This is a
reasonable assumption when, as in our application, at least one of the surveys uses a
properly randomised sampling scheme. When we cannot assume that one of the surveys
is unbiased by design, it is difficult to see how any method could deliver reliable
predictions without additional assumptions that would be difficult or impossible
to validate
empirically.

The problem that we have addressed in this paper is related to, but distinct from,
the problem of preferential sampling as formulated in \cite{prefer2010}. In both
settings, the goal is to predict the realisation of a latent spatial process $S(x)$
using data obtained by a potentially biased sampling scheme.  In preferential
sampling, the bias arises from a direct relationship between the value of $S(x)$ and
the probability that the location $x$ will be sampled. In the present paper, the
bias is a function of the location $x$ itself, rather than of the value of $S(x)$. In
the context of disease prevalence mapping, a further distinction is between
properties of a location and properties of a person who happens to live at that location.  Thus,
in our application a relationship between a child's location and the
likelihood that they would present at the CDH would not, in itself, result in bias. Rather,
the bias surface $B(x)$ allows for the possibility that the sub-population
of children who present at the CDH differs from the general population
with respect to their exposure to unmeasured risk-factors
for malaria.

Our approach is of potentially wide application to disease monitoring and control
in low-resource settings, where registry data are typically not available.
 The ability to combine data from surveys that vary in their level of bias and
timing can inform more accurate, local-area burden maps,
allowing for improved risk stratification of high burden areas and identification of
transmission hot-spots. For example, although substantial progress has been made over the past decade with malaria control by homogeneous scaling up of interventions at national level, it is increasingly recognized by funders and policy makers that a more targeted approach focused on high-burden areas or hot-spots may be more
cost-effective. Furthermore, apart from its potential to optimize the use of available data,
our approach can also inform improved prospective data collection,
by using the fitted model
in simulation studies to identify
efficient prospective hybrid sampling approaches
that combine convenience and random sampling strategies
in ways that acknowledge and exploit spatial and/or temporal heterogeneity as
revealed by analyses of the kind described in Section \ref{sec:application}.

In conclusion, our proposed approach provides a way of making use of
mixed source prevalence data to improve estimates of spatial predictions.
These are urgently needed to
support control programmes and develop more accurate local spatio-temporal risk stratification maps that can
inform more targeted control efforts. Malaria is one of a number of diseases
that bring a high public health burden in
low-resource settings, whilst exhibiting
 highly heterogeneous distributions
 across space and time. Control of such diseases needs methods of
 continuous monitoring of prevalence
and evaluation of control measures that make  the best possible use of
limited resources, and
will  therefore benefit greatly from
the ability to combine national household surveys with more local convenience sampling strategies
without compromising the validity of the resulting prevalence estimates.

\section*{Acknowledgements}
 We thank the participants and staff involved in the data collection of the presented household prevalence surveys, especially Dr Roca-Feltrer who oversaw the survey field teams between 2010 and 2011

\section*{Funding}
Emanuele Giorgi holds an ESRC-NWDTC funded Ph.D. studentship. Dr Sanie Sesay holds an MCDC funded Ph.D. studentship. Dr Dianne Terlouw acknowledges support from the ACT consortium
for the presented household prevalence surveys from Malawi. Prof Peter Diggle is supported
by the UK Medical Research Council (grant number G0902153).

\bibliographystyle{biometrika.bst}


text:
\section{Introduction}\label{ref3}

For every $Q=1,2,3,\ldots$, let
$\mathcal{F}(Q)=\set{0=s_0<s_1<s_2<\cdots<s_{n(Q)-1}}$ be a finite subset of the half-open real unit interval $[0,1)$. Assume that $\mathcal{F}(1)\subset\mathcal{F}(2)\subset\cdots$, with union dense in $[0,1)$. Setting $s_{n(Q)}=1$, the \newword{normalized gap} at $s_i\in\mathcal{F}(Q)$ is
\[
\nngg_Q(s_i)=n(Q)(s_{i+1}-s_i).
\]
If, for every $z\in\mathbb{R}_{\ge0}$, the limit
\begin{equation}\label{eq1}
H(z)=\lim_{Q\to\infty}\frac{\sharp\set{0\le i<n(Q):\nngg_Q(s_i)\le z}}{n(Q)}
\end{equation}
exists, then we say that the sequence of the $\mathcal{F}(Q)$'s has \newword{cumulative gap distribution}~$H$.

Two extreme cases of this setting are the Heaviside distribution (the distribution of a random variable which is $1$ almost surely), which is easily realizable via ``evenly spaced'' $\mathcal{F}(Q)$'s, and the exponential distribution $1-\exp(-z)$, which is almost surely induced whenever the points of $\mathcal{F}(Q)$ are given by i.i.d.~random variables uniformly distributed on~$[0,1)$.

Throughout this paper rational numbers $s=p/q$ are always written in reduced form (i.e., $q>0$ and $p,q$ relatively prime). The \newword{Farey sequence of order $Q$} is the set $\mathcal{F}_{\rm 1\mskip-4mu l}(Q)$ of all rational numbers in $[0,1)$ whose denominator is $\le Q$ (we'll explain the subscript ${\rm 1\mskip-4mu l}$ shortly). All intervals $[p_i/q_i,p_{i+1}/q_{i+1}]$ between consecutive points of $\mathcal{F}_{\rm 1\mskip-4mu l}(Q)\cup\set{1}$ are \newword{unimodular} (i.e., $\det\smatrix{p_{i+1}}{p_i}{q_{i+1}}{q_i}=1$), and hence have length $(q_{i+1}q_i)^{-1}$, which is bounded from below by $Q^{-2}$. It is a classical fact that the number of intervals in $\mathcal{F}_{\rm 1\mskip-4mu l}(Q)$ is asymptotic, for $Q\to\infty$, to $\bigl(2\,\zeta(2)\bigr)^{-1}Q^2$.
This immediately implies that the normalized gaps are bounded from below by $\bigl(2\,\zeta(2)\bigr)^{-1}=3/\pi^2=0.30396\ldots$, so that $H(z)$ exists and has value $0$ for $z\le 3/\pi^2$. This remark is just a fraction of Hall's classical result~\cite{hall70}, according to which the limit~\eqref{eq1} exists for the Farey sequence and agrees with the cumulative distribution $H_{\rm 1\mskip-4mu l}(z)$ of the random variable
\begin{equation*}
Z_1(x,y)=\frac{1}{2\,\zeta(2)\,xy}.
\end{equation*}
The latter is defined on the space $(\Omega_1,P_1)$, where
$\Omega_1$ is the triangle $\set{(x,y)\in\mathbb{R}_{>0}^2:x,y\le 1<x+y}$ and $P_1$ is the Lebesgue measure, normalized by $P_1(\Omega_1)=1$.

Explicit computation gives
\begin{equation*}
H_{\rm 1\mskip-4mu l}(z)=\begin{cases}
0, & \text{if $z\le 3/\pi^2$;}\\
2-6\bigl(1+\log(\pi^2z/3)\bigr)/(\pi^2z),
& \text{if $3/\pi^2<z\le 12/\pi^2$;} \\
2-6/(\pi^2z)-\sqrt{1-12/(\pi^2z)}+ \\
\quad 12\log(1/2+\sqrt{1/4-3/(\pi^2z)})/(\pi^2z),
& \text{if $12/\pi^2 < z$.}
\end{cases}
\end{equation*}
Differentiating, we obtain the density distribution function
\begin{equation*}
h_{\rm 1\mskip-4mu l}(z)=\begin{cases}
0, & \text{if $z\le 3/\pi^2$;}\\
6\log(\pi^2z/3)/(\pi^2z^2),
& \text{if $3/\pi^2<z\le 12/\pi^2$;} \\
-12\log\bigl(1/2+\sqrt{1/4-3/(\pi^2z)}\bigr)/(\pi^2z^2),
& \text{if $12/\pi^2 < z$.}
\end{cases}
\end{equation*}

We plot $h_{\rm 1\mskip-4mu l}(z)$ in Figure~\ref{ref1}; the two points of nondifferentiability correspond to the hyperbola $\set{Z_1(x,y)=z}$ hitting~$\Omega_1$ in the upper right corner (at $z=3/\pi^2$) and in the midpoint of the hypothenuse (at $z=12/\pi^2$).
\begin{figure}[h!]
\includegraphics[width=8cm,height=5cm]{h}
\caption{}
\label{ref1}
\end{figure}

In this paper we are interested in the statistics of gap distribution when the denominators of rational points are weighted by a fixed function $u:[0,1]\to\mathbb{R}_{>0}$ such that:
\begin{itemize}
\item[(i)] $u$ is continuous;
\item[(ii)] $[0,1]$ can be partitioned into finitely many subintervals $[s,s']$ (with $s,s'$ not necessarily rational), overlapping only at endpoints, and such that $u$ is $C^1$ with bounded first derivative on each $(s,s')$.
\end{itemize}
We call such a function a \newword{unit} (the name stems from algebraic logic~\cite{pantiravotti}). The \newword{$u$-denominator} of $s=p/q\in\mathbb{Q}\cap [0,1]$ is then $\den_u(s)=u(s)\,q$; the constant function~${\rm 1\mskip-4mu l}$ is a unit, and $\den_{\rm 1\mskip-4mu l}(s)$ is the ordinary denominator of $s$. We can then form the \newword{$u$-weighted Farey sequence} $\mathcal{F}_u(Q)$, whose elements are all the rational points in $[0,1)$ of $u$-denominator $\le Q$.

We can look at things projectively, by defining the \newword{cone} over $[0,1)$ by $\Cone[0,1)=\mathbb{R}_{>0}\cdot\bigl\{\cvect{s}{1}:s\in [0,1)\bigr\}\subset\mathbb{R}^2$. Every function $f$ on $[0,1)$ gives rise to its \newword{homogeneous correspondent} $\bb{f}:\Cone[0,1)\to\mathbb{R}$ by $\bb{f}(\bb{s})=
y\, f(x/y)$, where $\bb{s}=\cvect{x}{y}$.
By defining the homogeneous correspondent of the rational point $s=p/q$ to be $\bb{s}=\cvect{p}{q}$ (we are trying to use consistently boldface type for projective objects, and lightface for affine ones), we immediately see that $\den_u(s)=\bb{u}(\bb{s})$, and the points in $\mathcal{F}_u(Q)$ correspond bijectively to the primitive integer points in $\Cone[0,1)\cap\set{\bb{u}\le Q}$.

Summing up, we are looking at the limiting gap distribution ---call it $H_u$--- of the projections of the above primitive integer points on $[0,1)\times\set{1}$.
An explicit expression for $H_u$ as a mean over the Hall distribution $H_{\rm 1\mskip-4mu l}$ is obtained by Boca, Cobeli and Zaharescu in~\cite[Theorem~0.2]{boca_et_al00}; up to changes in parametrization and notation, it is formula~\eqref{eq2} below. Their proof uses incomplete Kloosterman sums and analytic number theory, spanning several pages of delicate computation. Building upon work of Athreya and Cheung~\cite{athreyache14}, we exploit here the properties of the horocycle flow on the space of rank-2 unimodular lattices to provide a short and reasonably self-contained proof of~\eqref{eq2}. Our key technical tool, and the main novelty of this paper, is the use of a ``sliding Poincar\'e section'' for the flow (see \S\ref{ref21}).

We remark that the applicability of dynamical equidistribution to the statistics of primitive lattice points inside an arbitrary star-shaped domain ---even in a higher dimensional setting--- was already pointed out by Marklof in~\cite{marklof13} (paragraph starting at the bottom of p.~50). Although in dimension greater than $1$ the resulting limit distributions are not as visualizable as those in the classical case, these techniques grant the transfer of much information; see~\cite{marklofstr10}, \cite{marklofstr11}, \cite{strombergsson11}. The crux of the matter lies ultimately in a very general result on the equidistribution of Farey sequences on large closed horospheres~\cite[Theorem~6]{marklof10}, in which the test function has an explicit dependence on the sequence points.

Let us state our results; throughout this paper $u$ is a fixed unit. We set
\begin{align*}
v(s)&=u(s)^{-1},\\
C&=\int_0^1 v(s)^2\,\mathrm{d} s \in\mathbb{R}_{>0},\\
m(x)&=C^{-1} v(s)^2.
\end{align*}
The explicit expression for $H_u$, to be proved in Theorem~\ref{ref20}, is then
\begin{equation}\label{eq2}
H_u(z)=\int_0^1 H_{\rm 1\mskip-4mu l}\bigl(m(s)\, z\bigr)\, m(s)\,\mathrm{d} s.
\end{equation}
Differentiating under the integral sign is safe here, and we obtain
\begin{equation}\label{eq3}
h_u(z)=\int_0^1 h_{\rm 1\mskip-4mu l}\bigl(m(s)\, z\bigr)\, m(s)^2\,\mathrm{d} s.
\end{equation}

\begin{example}
Let\label{ref5} $u$ be the unit
\[
u(x)=\begin{cases}
(5x+1)/2, &\text{if $0\le x\le 1/5$};\\
(x-2/5)^2+24/25, &\text{otherwise}.
\end{cases}
\]
In Figure~\ref{ref2} we plot the graph of $u$, as well as the histogram of the gap distribution of $\mathcal{F}_u(400)$ against the expression for $h_u(z)$ in~\eqref{eq3}.
\begin{figure}[h!]
\includegraphics[height=4.5cm]{figura1}
\includegraphics[height=4.5cm]{hu_and_hu_histogram}
\caption{}
\label{ref2}
\end{figure}
\end{example}

In Theorem~\ref{ref20} we will express $H_u$ as the cumulative distribution function of the random variable $Z(x,y)=C/(2\,\zeta(2)\,xy)$, with the Hall triangle $\Omega_1$ replaced by a pentagon endowed with an appropriate probability measure; see Figure~\ref{ref29} for the case of the unit of Example~\ref{ref5}. As a consequence, we will show in Theorem~\ref{ref30} that $h_u$ is piecewise-smooth, with finitely many points of nondifferentiability which can be explicitly determined. In the case of the above unit, $C=1.14002\ldots$ and the points of nondifferentiability are
\[
\frac{3C}{\pi^2}\biggl\{\frac{1}{4},\frac{576}{625},1,\frac{1089}{625},
\frac{2304}{625},4,\frac{4356}{625}\biggr\}.
\]
The first, second, and fifth of these points (at $0.08663\ldots$, $0.31935\ldots$ and $1.27743\ldots$, respectively)
are clearly visible in Figure~\ref{ref2}, while the others are quite hidden.

It is a pleasure to thank Jayadev Athreya for introducing me to the study of gap distribution via homogeneous dynamics, and for many clarifying and stimulating discussions on these topics.

\section{Basics}\label{ref24}

We first prove a fact which is interesting in its own right.

\begin{theorem}
There\label{ref36} exists $Q'$ such that, for every $Q\ge Q'$, all intervals $[p/q,p'/q']$ between successive elements of $\mathcal{F}_u(Q)\cup\set{1}$ are unimodular.
\end{theorem}
\begin{proof}
Recall that the \newword{Ford circle} $C_{p/q}$ at the rational number $p/q$ is the circle of radius $1/(2q^2)$, lying in the upper-half plane and tangent to the real axis at $p/q$~\cite{ford38}.
The circles at the points of $\mathcal{F}_u(Q)\cup\set1$ are then precisely those touching the real axis at points $s\in[0,1]\cap\mathbb{Q}$ and having center on or above the graph of $u(s)^2/(2Q^2)$; see Figure~\ref{ref37} for the case of the unit of Example~\ref{ref5} and $Q=5$.
\begin{figure}[h!]
\includegraphics[height=4.5cm]{figura5}
\caption{}
\label{ref37}
\end{figure}
By the basic properties of the construction, two distinct circles $C_{p/q}$ and $C_{p'/q'}$ are either wholly external to one another, or tangent, and are tangent iff $[p/q,p'/q']$ is unimodular.

Suppose now that the statement of the theorem is false. Then there exists an infinite sequence $Q_0<Q_1<Q_2<\cdots$ such that, for every $i=0,1,2,\ldots$, some interval $[p(i)/q(i),p'(i)/q'(i)]$ between successive elements of $\mathcal{F}_u(Q_i)\cup\set1$ is not unimodular. Without loss of generality  $q(i)\ge q'(i)$ and $q(i)>1$. By~\cite[Theorem~4]{ford38} there exists precisely one fraction $p''/q''$ such that $p(i)/q(i)<p''/q''$, $q(i)>q''$, and $[p(i)/q(i),p''/q'']$ is unimodular. The point $p''/q''$ must necessarily lie between $p(i)/q(i)$ and $p'(i)/q'(i)$, since otherwise $C_{p'(i)/q'(i)}$ and $C_{p''/q''}$ would intersect. Since $p(i)/q(i)$ and $p'(i)/q'(i)$ are consecutive in $\mathcal{F}_u(Q_i)\cup\set1$, the rational $p''/q''$ does not belong to $\mathcal{F}_u(Q_i)\cup\set1$, and therefore
\[
\frac{u(p''/q'')^2}{2Q_i^2}>\frac{1}{2(q'')^2}>
\frac{1}{2(q(i))^2}\ge\frac{u(p(i)/q(i))^2}{2Q_i^2}.
\]
Hence
\begin{multline*}
\frac{u(p''/q'')^2-u(p(i)/q(i))^2}{Q_i^2(p''/q''-p(i)/q(i))}>
\frac{(q'')^{-2}-(q(i))^{-2}}{p''/q''-p(i)/q(i)}=
\frac{q(i)^2-(q'')^2}{q''q(i)}\ge\\
\frac{q(i)^2-(q(i)-1)^2}{(q(i)-1)q(i)}=
\frac{2q(i)-1}{q(i)^2-q(i)}>\frac{1}{q(i)}.
\end{multline*}
As $p(i)/q(i)\in\mathcal{F}_u(Q_i)$, we have $Q_i\ge u(p(i)/q(i))q(i)\ge
\min(u)\,q(i)$, and hence
\[
\frac{u(p''/q'')^2-u(p(i)/q(i))^2}{p''/q''-p(i)/q(i)}>
\frac{Q_i^2\min(u)}{Q_i}=Q_i\min(u).
\]
Now the last term tends to infinity, but this is impossible, since $u$ being a unit immediately implies that the set
\[
\biggl\{\biggl\lvert\frac{u(s)^2-u(s')^2}{s-s'}\biggr\rvert:
s\not= s'\in[0,1]\biggr\}
\]
is bounded.
\end{proof}

For the rest of this paper, and without loss of generality, we assume that $Q$ is so large to satisfy the statement of Theorem~\ref{ref36}.

We define
\[
\nabla(1)=\biggl\{
\begin{pmatrix}
x\\ y
\end{pmatrix}
\in\Cone[0,1):
\bb{u}
\begin{pmatrix}
x\\ y
\end{pmatrix}
\le 1\biggr\};
\]
it constitutes a star-shaped sector, bounded by the lines $\set{x=0}$, $\set{x=y}$, and by the curve $\set{v(s)\cvect{s}{1}:0\le s<\le 1}$.

\begin{lemma}
We\label{ref25} have:
\begin{itemize}
\item[(i)] the area of $\nabla(1)$ is $C/2$;
\item[(ii)] $n(Q)$ is asymptotic to $CQ^2/\bigl(2\,\zeta(2)\bigr)$;
\item[(iii)] as $Q\to\infty$, the probability on $[0,1)$
\[
\frac{1}{n(Q)}\sum\set{\delta_s:s\in\mathcal{F}_u(Q)}
\]
(where $\delta_s$ is the Dirac mass at $s$) converges weakly${}^*$ to $m(s)\,\mathrm{d} s$.
\end{itemize}
\end{lemma}
\begin{proof}
(i) By elementary calculus, the sector swept by the line segment $\set{w\,v(s)\cvect{s}{1}:0<w\le 1}$ in time $\,\mathrm{d} s$ has area $\,\mathrm{d} A=2^{-1}\,v(s)^2\,\mathrm{d} s$; hence
$$
\area\bigl(\nabla(1)\bigr)=\int_0^1 2^{-1}\,v(s)^2\,\mathrm{d} s = C/2.
$$

(ii)-(iii) Fix a subinterval $[a,b)$ of $[0,1)$. The cardinality of
$\mathcal{F}_u(Q)\cap [a,b)$
is equal to the number of primitive points of the lattice $Q^{-1}\cvect{\mathbb{Z}}{\mathbb{Z}}$ inside the sector $\set{w\,v(s)\cvect{s}{1}:0<w\le1\text{ and }a\le s< b}$. For $Q\to\infty$, the number of such points is asymptotic to $Q^2$ times the area of the sector divided by $\zeta(2)$~\cite[Theorem~459]{hardywri85}, and it follows that the ratio between the number of primitive points inside the sector and the total number of primitive points in $\nabla(1)$ is asymptotic to the ratio of the relative areas. By the first part of the proof we then have
\[
\lim_{Q\to\infty}\frac{\sharp\bigl(\mathcal{F}_u(Q)\cap [a,b)\bigr)}{n(Q)}=
2/C\int_a^b 2^{-1} v(s)^2\,\mathrm{d} s=
\int_a^b m(s)\,\mathrm{d} s.
\]
\end{proof}

We recall a few basic facts about the horocycle flow; see~\cite[Chapter~IV]{bekkamayer} or~\cite[Chapter~11]{einsiedlerward} for a full treatment. The group $\PSL_2\mathbb{R}$ acts on the upper halfplane $\mathcal{H}\subset\mathbb{C}$ on the left: if
$A=\smatrix abcd$ and $\alpha\in\mathcal{H}$, then
\[
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
\ast\alpha
=\frac{a\alpha+b}{c\alpha+d}.
\]
We identify the unit tangent space at $\alpha$ with $\set{\tau\in\mathbb{C}:\abs{\tau}=\im\alpha}$; the above action extends then to a left action of $\PSL_2\mathbb{R}$ on the unit tangent bundle $T^1\mathcal{H}$ via $A\ast(\alpha,\tau)=(A\ast\alpha,(c\alpha+d)^{-2}\tau)$. This latter action is transitive with trivial stabilizers, so we get a bijection
\[
\PSL_2\mathbb{R}\ni A \mapsto A\ast (i,i)\in T^1\mathcal{H},
\]
whose inverse is given by the Iwasawa decomposition
\[
\begin{pmatrix}
1 & \re\alpha \\
0 & 1
\end{pmatrix}
\begin{pmatrix}
(\im\alpha)^{1/2} & 0 \\
0 & (\im\alpha)^{-1/2}
\end{pmatrix}
\begin{pmatrix}
\cos\theta & -\sin\theta \\
\sin\theta & \cos\theta
\end{pmatrix}
\leftarrow\!\mapstochar
(\alpha,\tau),
\]
where $\theta=-\arg(\tau/i)/2$.

The above bijection is preserved under the quotient by the left action of the modular group $\Gamma=\PSL_2\mathbb{Z}$, so we have an induced bijection
\begin{equation*}
\Phi:\Gamma\backslash\PSL_2\mathbb{R} \to T^1M,
\end{equation*}
where $M=\Gamma\backslash\mathcal{H}$ is the modular surface. We keep on writing $\xi=(\alpha,\tau)$ for the points of $T^1M$, without explicit mention of the action of~$\Gamma$.

We can extend $\Phi$ to the space $X_2$ of all unimodular lattices in $\mathbb{R}^2$. To this purpose, it is expedient to identify $\mathbb{R}^2$ with $\mathbb{C}$ as real vector spaces, with the standard basis $\rvect{\bb{e}_1}{\bb{e}_2}$ corresponding to $\rvect{i}{1}$ (we always drop reference to $\rvect{\bb{e}_1}{\bb{e}_2}$, by writing $\cvect{x}{y}$ for $\rvect{\bb{e}_1}{\bb{e}_2}\cvect{x}{y}$).
We also set $S=\smatrix{0}{-1}{1}{0}$ and remark that, for $A\in\PSL_2\mathbb{R}$, the matrix $SA^{-1} S^{-1}$ is the transpose $A^T$ of $A$.

\begin{lemma}
Let\label{ref16} $\Xi:\Gamma\backslash\PSL_2\mathbb{R}\to X_2$ and $\Psi:T^1M\to X_2$ be defined by
\begin{align*}
\Xi:\Gamma A & \mapsto
SA^{-1} S^{-1}
\begin{pmatrix}
\mathbb{Z} \\
\mathbb{Z}
\end{pmatrix}
=A^T
\begin{pmatrix}
\mathbb{Z} \\
\mathbb{Z}
\end{pmatrix},\\
\Psi:(\alpha,\tau) & \mapsto
(i/\tau)^{1/2}
\begin{pmatrix}
\alpha & 1
\end{pmatrix}
\begin{pmatrix}
\mathbb{Z} \\
\mathbb{Z}
\end{pmatrix}.
\end{align*}
Then $\Xi,\Psi$ are well-defined bijections and the square
\[
\begin{xy}
\xymatrix{
\Gamma\backslash\PSL_2\mathbb{R} \ar[r]^\Phi \ar[d]_\Xi &
T^1M \ar[d]^\Psi \\
X_2 \ar@{=}[r] & X_2
}
\end{xy}
\]
whose bottom row results from the identification of $\mathbb{R}^2$ with $\mathbb{C}$ as above, commutes.
\end{lemma}
\begin{proof}
The definition of $\Xi$ is clearly unambiguous, and so is that of $\Psi$, since $(i/\tau)^{1/2}$ is determined up to sign. One checks easily that $\Xi$ and $\Psi$ are bijections and that $\Xi^{-1}\Psi\Phi=\id$.
\end{proof}

The group $\PSL_2\mathbb{R}$ acts on the right on all spaces in Lemma~\ref{ref16}. Namely, for $A,B,R\in\PSL_2\mathbb{R}$ we have
\begin{align*}
\Gamma A\ast R &= \Gamma AR,\\
\bigl(A\ast(i,i)\bigr)\ast R &= AR\ast(i,i),\\
B\begin{pmatrix}
\mathbb{Z}\\
\mathbb{Z}
\end{pmatrix}\ast R &= R^T B
\begin{pmatrix}
\mathbb{Z}\\
\mathbb{Z}
\end{pmatrix},\\
\begin{pmatrix}
i & 1
\end{pmatrix}
B\begin{pmatrix}
\mathbb{Z}\\
\mathbb{Z}
\end{pmatrix}\ast R &=
\begin{pmatrix}
i & 1
\end{pmatrix}
R^T B
\begin{pmatrix}
\mathbb{Z}\\
\mathbb{Z}
\end{pmatrix}.
\end{align*}
By construction, the square in Lemma~\ref{ref16} is equivariant under this action. In particular, the (unstable, time-reversing) horocycle flow $h_t$ is defined on $T^1M$ and~$X_2$ by
\begin{align*}
h_t:A\ast(i,i)&\mapsto A
\begin{pmatrix}
1 & 0\\
-t & 1
\end{pmatrix}
\ast(i,i),\\
h_t:B\begin{pmatrix}
\mathbb{Z}\\
\mathbb{Z}
\end{pmatrix}
&\mapsto
\begin{pmatrix}
1 & -t\\
0 & 1
\end{pmatrix}
B\begin{pmatrix}
\mathbb{Z}\\
\mathbb{Z}
\end{pmatrix}.
\end{align*}

\section{A sliding section for the horocycle flow}\label{ref21}

We say that the lattice $\Lambda\in X_2$ \newword{contains a vertical vector} if it contains a vector of the form $\cvect{0}{d}$ for some $d>0$.
The set
\[
\mathcal{S}=\set{\Lambda\in X_2:\Lambda\text{ contains a vertical vector}}
\]
is a 2-dimensional immersed submanifold of $X_2$, dense in $X_2$.
We can easily parametrize $\mathcal{S}$: indeed, every $\Lambda\in\mathcal{S}$ is of the form
\[
\Lambda=\begin{pmatrix}
d^{-1} & 0\\
c & d
\end{pmatrix}
\begin{pmatrix}
\mathbb{Z} \\ \mathbb{Z}
\end{pmatrix},
\]
with $d>0$ uniquely determined, and $c$ determined up to translation by integer multiples of $d$.
This can be described as follows: let $\mathbb{Z}$ act on the real upper halfplane $\mathbb{R}\times\mathbb{R}_{>0}$ by $k\ast(c,d)=(c+kd,d)$. Then the map
\begin{align*}
\tilde\phi:\mathcal{S}&\to \mathbb{Z}\backslash(\mathbb{R}\times\mathbb{R}_{>0}),\\
\Lambda&\mapsto
\mathbb{Z}\ast(c,d),
\end{align*}
is a homeomorphism, so that $\mathcal{S}$ is an immersed cylinder.

For every $D>0$, the set $\mathcal{S}\restriction D$ of all elements of $\mathcal{S}$ that contain a vertical vector of length $\le D$ is a Poincar\'e section for the horocycle flow. Indeed, all elements of $X_2$ ---except the codimension-$1$ set of lattices of the form $\smatrix{a}{b}{0}{a^{-1}}\cvect{\mathbb{Z}}{\mathbb{Z}}$ for some $a<D^{-1}$--- enter $\mathcal{S}\restriction D$ countably many times under the action of the flow, both in the past and in the future.
Up to the bijection $\Phi\circ\Xi^{-1}$ the cylinder
$\mathcal{S}$ is identifiable with $\Gamma\backslash\set{(\alpha,(\im\alpha)\,i):\alpha\in\mathcal{H}}$. The section
$\mathcal{S}\restriction D$ is then the set of all elements of $T^1M$ that have a lift to $T^1\mathcal{H}$ of the form $\bigl(\alpha,(\im\alpha)\,i\bigr)$ for some $\alpha$ having imaginary part $\ge D^{-2}$.

We now want ---this being the key idea in this paper--- to let $D$ vary with time.

\begin{definition}
Given\label{ref14} $Q$, for $t$ in $[0,Q^2)$ we define $\mathcal{S}_t=\mathcal{S}\restriction v(Q^{-2}t)$.
Let $\Lambda_0(Q)=\smatrix{Q}{0}{0}{Q^{-1}}\cvect{\mathbb{Z}}{\mathbb{Z}}$; we safely assume $Q^{-1}\le v(0)$, so that $\Lambda_0(Q)\in\mathcal{S}_0$. Writing $\Lambda_t(Q)$ for $h_t\bigl(\Lambda_0(Q)\bigr)$, we say that $\Lambda_0(Q)$ \newword{hits} the sliding section $\mathcal{S}_t$ at time~$t$ if $\Lambda_t(Q)\in\mathcal{S}_t$.
\end{definition}

For ease of notation, whenever $Q$ is understood we write $\Lambda_t$ for $\Lambda_t(Q)$.

\begin{lemma}
The\label{ref26} $h_t$-orbit of $\Lambda_0$ is periodic of period $Q^2$.
The hitting times $0=t_0<t_1<t_2<\cdots$ are precisely the multiples
$t_i=Q^2s_i$ of the elements $s_i\in\mathcal{F}_u(Q)$.
\end{lemma}
\begin{proof}
The first statement is clear.
As remarked in the proof of Lemma~\ref{ref25}, the elements of $\mathcal{F}_u(Q)$ are in 1-1 correspondence
with the primitive points of the lattice $Q^{-1}\cvect{\mathbb{Z}}{\mathbb{Z}}$ which are inside~$\nabla(1)$.
Let
\[
\nabla(Q)=\biggl\{
\begin{pmatrix}
Q^2x\\
y
\end{pmatrix}:
\begin{pmatrix}
x\\
y
\end{pmatrix}
\in\nabla(1)\biggr\};
\]
then
\[
\frac{p}{q}\leftrightarrow
\begin{pmatrix}
Qp\\ Q^{-1} q
\end{pmatrix}
\]
is a 1-1 correspondence between the points of
$\mathcal{F}_u(Q)$ and the primitive points of $\Lambda_0$
that are in $\nabla(Q)$.

Now, to say that $\Lambda_t\in\mathcal{S}_t$ amounts to saying that the ray $\mathbb{R}_{>0}\cvect t1$ passes through one of these primitive points. Therefore $t$ is a hitting time iff $Qp=tQ^{-1} q$ (i.e., $t=Q^2p/q$) for some $p/q\in\mathcal{F}_u(Q)$.
\end{proof}

Let $t_i=Q^2s_i=Q^2p_i/q_i$ be a hitting time. We choose a lift $\phi(\Lambda_{t_i})$ of $\tilde\phi(\Lambda_{t_i})=\mathbb{Z}\ast(c,d_i)$ to $\mathbb{R}_{>0}^2$ as follows: of all possible choices for $c$ we pick the largest one --- call it $c_i$---
such that
\[
0< c \le v\bigl(Q^{-2}(t_i+(cd_i)^{-1})\bigr),
\]
and set $\phi(\Lambda_{t_i})=(c_i,d_i)$.

\begin{lemma}
We\label{ref22} have
\begin{align*}
(c_i, d_i)&=Q^{-1}(q_{i+1},q_i),\\
t_{i+1}-t_i&=\frac{1}{c_id_i}.
\end{align*}
\end{lemma}
\begin{proof}
By construction, and according to Lemma~\ref{ref26},
\[
h_{t_i}^{-1}(\Lambda_{t_i})=\Lambda_0=
\begin{pmatrix}
Qp_{i+1} & Qp_i\\
Q^{-1} q_{i+1} & Q^{-1} q_i
\end{pmatrix}
\begin{pmatrix}
\mathbb{Z}\\ \mathbb{Z}
\end{pmatrix};
\]
this is justified by our standing assumption that $Q$ is so large that all intervals determined by $\mathcal{F}_u(Q)$ are unimodular, as guaranteed by Theorem~\ref{ref36}. We thus get
\[
\Lambda_{t_i}=\begin{pmatrix}
1 & -Q^2p_i/q_i\\
0 & 1
\end{pmatrix}
\begin{pmatrix}
Qp_{i+1} & Qp_i\\
Q^{-1} q_{i+1} & Q^{-1} q_i
\end{pmatrix}
\begin{pmatrix}
\mathbb{Z}\\ \mathbb{Z}
\end{pmatrix}=
\begin{pmatrix}
Qq_i^{-1} & 0\\
Q^{-1} q_{i+1} & Q^{-1} q_i
\end{pmatrix}
\begin{pmatrix}
\mathbb{Z}\\ \mathbb{Z}
\end{pmatrix},
\]
so that $d_i=Q^{-1} q_i$.
Let
\[
s'=\frac{p'}{q'}=\frac{p_{i+1}+p_i}{q_{i+1}+q_i}
\]
be the Farey mediant of $s_i$ and $s_{i+1}$. By definition of $\mathcal{F}_u(Q)$, we have $u(s_{i+1})q_{i+1}\le Q$ and $u(s')q'>Q$.
Therefore
\[
0<Q^{-1} q_{i+1}\le
v(s_{i+1})
=v\bigl(Q^{-2}t_i+(q_{i+1}q_i)^{-1}\bigr)=
v\bigl(Q^{-2}(t_i+(Q^{-1} q_{i+1}d_i)^{-1})\bigr),
\]
while
\begin{multline*}
Q^{-1} q_{i+1}+d_i=Q^{-1} q'>
v(s')
=v\bigl(s_i+(q'q_i)^{-1}\bigr)\\
=v\bigl(Q^{-2}(t_i+(Q^{-1} q' Q^{-1} q_i)^{-1})\bigr)=
v\bigr(Q^{-2}(t_i+((Q^{-1} q_{i+1}+d_i)d_i)^{-1})\bigr).
\end{multline*}
Therefore $c_i=Q^{-1} q_{i+1}$ as claimed, and $t_{i+1}-t_i=Q^2(s_{i+1}-s_i)=Q^2(q_{i+1}q_i)^{-1}=(c_id_i)^{-1}$.
\end{proof}

\begin{definition}
We\label{ref27} denote the minimum and maximum of $v$ on $[0,1]$
by $l$ and $L$, respectively.
For $w\in[l,L]$ we set:
\begin{align*}
\Omega_w &=\text{the triangle $\set{(x,y)\in\mathbb{R}_{>0}^2:x,y\le w<x+y}$},\\
P_w&=\text{the Lebesgue measure on $\Omega_w$, normalized by $P_w(\Omega_w)=1$},\\
\Omega&=\textstyle{\bigcup_w}\Omega_w=\text{the pentagon }\set{0<x,y\le L}\cap\set{l<x+y},\\
P&=\text{the probability measure on $\Omega$ defined by }\int_0^1 P_{v(s)}\,m(s)\,\mathrm{d} s.
\end{align*}
Given $Q$, we write $\phi(\Lambda_{t_i})=(c_i,d_i)$ for any hitting time $t_i$ of $\Lambda_0$ to the sliding section $\mathcal{S}_t$.
Then
\[
P(Q)=\frac{1}{n(Q)}\sum\delta_{(c_i,d_i)}
\]
is a point-process probability measure
on $\mathbb{R}_{>0}^2$.
\end{definition}

The following is our main theorem.

\begin{theorem}
As\label{ref19} $Q$ goes to infinity, $P(Q)$ converges weakly${}^*$ to $P$.
\end{theorem}

We will prove Theorem~\ref{ref19} in \S\ref{ref23}. In Figure~\ref{ref29} we plot the support of $P(400)$ for the unit of Example~\ref{ref5}.
\begin{figure}[h!]
\includegraphics[width=8cm]{figura3}
\caption{}
\label{ref29}
\end{figure}

\begin{theorem}
The\label{ref20} gap distribution $H_u(z)$ of the weighted Farey sequence $\mathcal{F}_u(Q)$
is the cumulative density distribution of the random variable
\begin{align*}
Z:(\Omega,P) &\to\mathbb{R}_{>0},\\
(x,y) &\mapsto \frac{C}{2\,\zeta(2)\,xy}.
\end{align*}
Also,
\[
P(Z\le z)=\int_0^1 H_{\rm 1\mskip-4mu l}\bigl(m(s)\, z\bigr)\, m(s)\,\mathrm{d} s,
\]
so formula~\eqref{eq2} holds.
\end{theorem}
\begin{proof}
Fix $z>0$. Then
\begin{align*}
&\lim_{Q\to\infty}\frac{1}{n(Q)}
\sharp\set{0\le i<n(Q):\nngg_Q(s_i)\le z}\\
=&
\lim_{Q\to\infty}\frac{1}{n(Q)}
\sharp\set{i:\bigl(2\,\zeta(2)\bigr)^{-1} CQ^2(s_{i+1}-s_i)\le z}&&\text{by Lemma~\ref{ref25}(ii)}\\
=&\lim_{Q\to\infty}\frac{1}{n(Q)}
\sharp\set{i:\bigl(2\,\zeta(2)\bigr)^{-1} C(c_id_i)^{-1}\le z} &&\text{by Lemma~\ref{ref22}}\\
=&\lim_{Q\to\infty}\bigl(P(Q)\bigr)(Z\le z)\\
=&\;P(Z\le z)&&\text{by Theorem~\ref{ref19}}.
\end{align*}

For the second claim we observe that $Z(x,y)\le z$ iff $Z(w^{-1} x,w^{-1} y)\le w^2z$, so that $P_w(Z\le z)=P_1(Z\le w^2z)$ for every $w\in[l,L]$. We then compute
\begin{align*}
P(Z\le z)&=\int_0^1 P_{v(s)}(Z\le z)\,m(s) \,\mathrm{d} s\\
&=\int_0^1 P_1(Z\le v(s)^2\,z)\,m(s) \,\mathrm{d} s\\
&=\int_0^1 P_1(Z_1\le C^{-1} v(s)^2\,z)\,m(s) \,\mathrm{d} s\\
&=\int_0^1 H_{\rm 1\mskip-4mu l}(m(s)\, z)\, m(s) \,\mathrm{d} s.
\end{align*}
\end{proof}

From here on $\lambda$ denotes the $1$-dimensional Lebesgue measure.

\begin{theorem}
(i)
The\label{ref30} measure $P$ is absolutely continuous w.r.t.~$\,\mathrm{d} x\wedge\,\mathrm{d} y$. More precisely,
let
\[
p(x,y)=2C^{-1}(v_*\lambda)\bigl([\max(x,y),x+y)\cap[l,L]\bigr),
\]
where $v_*\lambda$ is the pushforward of $\lambda$ by $v$.
Then
\[
\,\mathrm{d} P=p(x,y)\,\mathrm{d} x\wedge\,\mathrm{d} y.
\]

(ii) The probability density function $h_u$ is piecewise-smooth, with finitely many nondifferentiability points. These points are those in the set
\[
\frac{C}{2\,\zeta(2)}\bigl(\set{u(s)^2:s\in E}\cup\set{4u(s)^2:s\in E}\bigr),
\]
where $E\subset[0,1]$ contains $0$, $1$, and all points at which $v$ (equivalently, $u$) is nondifferentiable or has a local maximum or minimum.
\end{theorem}
\begin{proof}
(i) By definition, $\,\mathrm{d} P_w=2w^{-2}{\rm 1\mskip-4mu l}_{\Omega_w}(x,y)\,\mathrm{d} x\wedge\,\mathrm{d} y$. Thus, by Fubini,
\begin{align*}
\int f\,\mathrm{d} P &= \int_0^1\biggl[\int f\,\mathrm{d} P_{v(s)}\biggr]\,m(s)\,\mathrm{d} s\\
&=\int_0^1\biggl[\int f\,2\,v(s)^{-2}{\rm 1\mskip-4mu l}_{\Omega_{v(s)}}(x,y)
\,\mathrm{d} x\wedge\,\mathrm{d} y\biggr]\,C^{-1} v(s)^2\,\mathrm{d} s\\
&=2C^{-1} \int f\biggl[\int_0^1 {\rm 1\mskip-4mu l}_{\Omega_{v(s)}}(x,y)\,\mathrm{d} s\biggr]\,\mathrm{d} x\wedge\,\mathrm{d} y\\
&=2C^{-1} \int f\,\lambda\set{s:x,y\le v(s)<x+y}\,\mathrm{d} x\wedge\,\mathrm{d} y\\
&=2C^{-1} \int f\,\lambda(v^{-1}[\max(x,y),x+y))\,\mathrm{d} x\wedge\,\mathrm{d} y\\
&=\int f\,p(x,y)\,\mathrm{d} x\wedge\,\mathrm{d} y.
\end{align*}

(ii) This is best conveyed in geometrical language. As the sliding triangle $\Omega_{v(s)}$ moves through $\Omega$, it deposits mass. By our assumptions about $u$, this process is smooth except at the points $s\in E$, where the triangle starts or stops moving, reverses direction, or changes speed abruptly. We thus get a finite set $\set{\Omega_{v(s)}:s\in E}$ of triangles along whose borders the density $p(x,y)$ is singular. Now, the probability distribution function $h_u$ is singular at $z$ precisely when the hyperbola $\set{Z(x,y)=z}$ ---in its downward movement as $z$ goes from $0$ to infinity--- touches one of these triangles, say $\Omega_{v(s)}$, either in the upper right corner or in the midpoint of the hypothenuse. In the first case we have $Z(v(s),v(s))=z$ (i.e., $z=\bigl(2\,\zeta(2)\bigr)^{-1} Cu(s)^2$), while in the second we have $Z(v(s)/2,v(s)/2)=z$ (i.e., $z=2\,\zeta(2)^{-1} Cu(s)^2$).
\end{proof}

\section{Proof of Theorem~\ref{ref19}}\label{ref23}

Fix a sequence of real numbers $\varepsilon_0,\varepsilon_1,\varepsilon_2,\ldots$, strictly decreasing and converging to $0$. For every $k$, fix $Q_k$ so large that every $v$-image of an interval in $\mathcal{F}_u(Q_k)$ has length $\le\varepsilon_k$. Given $k$ and $0\le h<n(Q_k)$, we let $I(k,h)$ be the $h$-th interval
(closed to the left and open to the right) of the partition of $[0,1)$ determined by $\mathcal{F}_u(Q_k)$. We also let $r(k,h)$ be a point in the topological interior of $I(k,h)$ such that
\begin{equation}\label{eq5}
\lim_{Q\to\infty}
\frac{\sharp\bigl(\mathcal{F}_u(Q)\cap I(k,h)\bigr)}{\sharp\mathcal{F}_u(Q)}=
\int_{I(k,h)} m(s)\,\mathrm{d} s=
m(r(k,h))\,\lambda(I(k,h));
\end{equation}
such a point exists by Lemma~\ref{ref25}(iii) and the intermediate value theorem for integrals of continuous functions.

For short, we write
\begin{align*}
\mathcal{F}_u(Q,k,h)&=\mathcal{F}_u(Q)\cap I(k,h),\\
\mathcal{F}_{u(r)}(Q,k,h)&=\mathcal{F}_{u(r(k,h))}(Q)\cap I(k,h).
\end{align*}
Note that $u(r(k,h))$ and $v(r(k,h))$ may denote either a number or ---as above--- the constant function whose value is that number; the context always makes the meaning clear.

\begin{lemma}
For\label{ref15} $Q\to\infty$, the cardinalities of $\mathcal{F}_u(Q,k,h)$ and of $\mathcal{F}_{u(r)}(Q,k,h)$ are asymptotically equal.
\end{lemma}
\begin{proof}
As in Lemma~\ref{ref25}, the ratio of the two cardinalities is asymptotic to the ratio of the areas of the two sectors
$\set{w\,v(s)\cvect{s}{1}:0<w\le 1\text{ and }s\in\ I(k,h)}$
and $\set{w\,v(r(k,h))\cvect{s}{1}:0<w\le 1\text{ and }s\in\ I(k,h)}$.
The first sector has area
$$
\int_{I(k,h)} 2^{-1} v(s)^2\,\mathrm{d} s,
$$
while the second has area
$$
2^{-1}  v(r(k,h))\bigl(\lambda(I(k,h))\,v(r(k,h))\bigr).
$$
Taking into account the definition of $r(k,h)$ and the fact that $v(s)^2=C\,m(s)$, one checks immediately that the two areas agree.
\end{proof}

Let $l(k,h)$, $L(k,h)$ be the infimum and the supremum of $v$ on $I(k,h)$, respectively. Then the triangle
$$
\Omega_{I(k,h)}^o=\set{x,y\le l(k,h)}\cap\set{L(k,h)< x+y}
$$
lies inside $\Omega_{v(r(k,h))}$. In turn, the latter lies inside the pentagon
$$
\Omega_{I(k,h)}=\set{0<x,y\le L(k,h)}\cap\set{l(k,h) < x+y};
$$
see Figure~\ref{ref33}.
\begin{figure}[h!]
\includegraphics[width=8cm]{figura4}
\caption{}
\label{ref33}
\end{figure}

For $Q\ge Q_k$ and $t$ varying in $Q^2I(k,h)$, we are aiming at freezing the sliding section $\mathcal{S}_t$ to the fixed section $\mathcal{S}_{Q^2r(k,h)}$. Lemma~\ref{ref15} guarantees that for large $Q$'s this has no impact on the number of hits.
The following Lemma~\ref{ref34} copes with the change of the coordinate function. Indeed, replacing $u$ with the constant function $u(r(k,h))$ we are forced to assign the lattice
\begin{equation}
\label{ref18}
\Lambda_t(Q)=
\begin{pmatrix}
d^{-1} & 0\\
c & d
\end{pmatrix}
\begin{pmatrix}
\mathbb{Z} \\ \mathbb{Z}
\end{pmatrix}
\end{equation}
(with $t\in Q^2I(k,h)$ and $d\le v(r(k,h))$) the coordinates $\phi'(\Lambda_t(Q))=(c',d)$, where $c'$ is the largest element in $c+d\mathbb{Z}$ such that $c'\le v(r(k,h))$; assuming that $t$ is a hitting time for the sliding section as well, the two lifts, $\phi(\Lambda_t)$ and $\phi'(\Lambda_t)$, of $\tilde\phi(\Lambda_t)=\mathbb{Z}\ast(c,d)$ may differ in their first component.

\begin{lemma}
Let\label{ref34} $Q\ge Q_k$, $t\in Q^2I(k,h)$, and let $\Lambda_t=\Lambda_t(Q)$ be as in~\eqref{ref18}. Then:
\begin{itemize}
\item[(i)] if $t$ is a hitting time for the sliding section, then $\phi(\Lambda_t)\in\Omega_{I(k,h)}$;
\item[(ii)] if $t$ is a hitting time for the fixed section, then $\phi'(\Lambda_t)\in\Omega_{v(r(k,h))}$;
\item[(iii)] if $\Sigma\ast(c,d)=\Sigma\ast(c',d)$ with $(c',d)\in\Omega_{I(k,h)}^o$, then $t$ is a hitting time for both sections and $\phi(\Lambda_t)=\phi'(\Lambda_t)=(c',d)$.
\end{itemize}
\end{lemma}
\begin{proof}
(i) Say $Q^{-2}t=s_i\in I(k,h)$, with $s_{i+1}$ the element following $s_i$ in $\mathcal{F}_u(Q)$. By Lemma~\ref{ref22},
$\phi(\Lambda_t)=(Q^{-1} q_{i+1},Q^{-1} q_i)$. Since $Q^{-1} q_i\le v(s_i)$ and $Q^{-1} q_{i+1}\le v(s_{i+1})$, both $Q^{-1} q_i$ and $Q^{-1} q_{i+1}$ are $\le L(k,h)$. On the other hand, as the Farey mediant $s'$ of $s_i$ and $s_{i+1}$ does not belong to $\mathcal{F}_u(Q)$, we have $Q^{-1} q_{i+1}+Q^{-1} q_i=Q^{-1}(q_{i+1}+q_i)>v(s')\ge l(k,h)$.
(ii) follows from (i), applied to the constant unit $u(r(k,h))$.
(iii) Assume $(c',d)\in\Omega^o_{I(k,h)}$. Then $d\le l(k,h)\le v(Q^{-2}t), v(r(k,h))$, so $t$ is a hitting time for both sections. Since $c'+d>L(k,h)$ and $(c'-d)+d=c'\le l(k,h)$, the only lift of $\tilde\phi(\Lambda_t(Q))$ which is contained in $\Omega_{I(k,h)}$ is $(c',d)$. By~(i) and~(ii) both lifts $\phi(\Lambda_t)$ and $\phi'(\Lambda_t)$ must necessarily be equal to $(c',d)$.
\end{proof}

Recall from~\S\ref{ref21} that $\mathcal{S}\restriction v(r(k,h))$
is a transverse section for the horocycle flow. The Liouville measure on $X_2$ ---namely, the pushforward via $\Xi$ of the unique $\PSL_2\mathbb{R}$ right invariant measure on $\Gamma\backslash\PSL_2\mathbb{R}$, see the diagram in Lemma~\ref{ref16}--- decomposes locally as a product of a transverse measure for the horocycle foliation (the foliation whose leaves are the $h_t$-orbits) and the linear measure on the leaves.
Let $F$ be the map $\Phi\circ\Xi^{-1}\circ\phi^{-1}$ from $\Omega_{v(r(k,h))}$ to $T^1M$; then
$$
F(x,y)=(xy^{-1}+y^{-2} i,y^{-2} i)=(\xi +\nu i,\nu i).
$$
Write $U$ for the infinitesimal generator of the horocycle flow, and $\omega=\nu^{-2}\,\mathrm{d}\xi\wedge\,\mathrm{d}\nu\wedge\,\mathrm{d}\theta$ for the riemannian volume form.
Then the transverse measure corresponds on
$\Omega_{v(r(k,h))}$ to the pullback via $F$ of the contraction $\iota_U\omega$ of $\omega$ w.r.t.~the vector field~$U$. A straightforward computation shows now that $F^*(\iota_U\omega)=2\,\mathrm{d} x\wedge\,\mathrm{d} y$.

A key result by Sarnak~\cite{sarnak81} states that closed horocycles on $X_2$ become equidistributed as their length goes to infinity. This has been refined by Hejhal~\cite{hejhal96}, who shows that equidistribution still holds if the orbits are restricted to a constant fraction of their full length. In our case, this implies that for every pair $(k,h)$ the probability
$$
\frac{1}{\sharp\mathcal{F}_{u(r)}(Q,h,k)}
\sum\set{\delta_{\phi'\Lambda_{Q^2s}}:s\in\mathcal{F}_{u(r)}(Q,k,h)}
$$
converges weakly${}^*$ to the normalization of $2\,\mathrm{d} x\wedge\,\mathrm{d} y$ to a probability on $\Omega_{v(r(k,h))}$, namely to $\,\mathrm{d} P_{v(r(k,h)}$. Summing up,
\begin{equation}\label{eq7}
\lim_{Q\to\infty}\frac{\sum\set{f(\phi'(\Lambda_{Q^2s}(Q))):s\in\mathcal{F}_{u(r)}(Q,k,h)}}{\sharp\mathcal{F}_{u(r)}(Q,k,h)}
=\int f\,\mathrm{d} P_{v(r(k,h))},
\end{equation}
for every continuous function $f:\mathbb{R}_{>0}^2\to\mathbb{R}$ of compact support.

Given such a function $f$, let us write for short
\begin{align*}
\phi\Lambda_{Q^2s}&=\phi\bigl(\Lambda_{Q^2s}(Q)\bigr),\text{ and analogously for $\phi'$},\\
\alpha(Q,k,h)&=
\frac{\Sigma\set{f(\phi'\Lambda_{Q^2s}):s\in\mathcal{F}_{u(r)}(Q,k,h)}-\set{f(\phi\Lambda_{Q^2s}):s\in\mathcal{F}_u(Q,k,h)}}{\sharp\mathcal{F}_{u(r)}(Q,k,h)},\\
a(Q,k)&=\sum_{h=0}^{n(Q_k)-1}\Biggl(\alpha(Q,k,h)\,\frac{\sharp\mathcal{F}_u(Q,k,h)}{\sharp\mathcal{F}_u(Q)}\Biggr).
\end{align*}

\begin{lemma}
Given\label{ref31} $k$ there exists $N(k)$ such that, for every $Q$ sufficiently large, both $\abs{\alpha(Q,k,h)}$ and $\abs{a(Q,k)}$ are $\le N(k)$. Also, $\lim_{k\to\infty}N(k)=0$.
\end{lemma}
\begin{proof}
Take first $R=\max(f)-\min(f)$. By Lemma~\ref{ref34}(iii) $\abs{\alpha(Q,k,h)}$ is bounded from above by
\begin{multline*}
R\Biggl(
\frac{\sharp\set{s\in\mathcal{F}_{u(r)}(Q,k,h):\phi'\Lambda_{Q^2s}\notin\Omega_{I(k,h)}^o}}
{\sharp\mathcal{F}_{u(r)}(Q,k,h)}\\
+
\frac{\sharp\set{s\in\mathcal{F}_u(Q,k,h):\phi\Lambda_{Q^2s}\notin\Omega_{I(k,h)}^o}}
{\sharp\mathcal{F}_{u(r)}(Q,k,h)}
\Biggr).
\end{multline*}
Possibly increasing $R$ to take care of the asymptotic $\sharp\mathcal{F}_{u(r)}(Q,k,h)\sim\sharp\mathcal{F}_u(Q,k,h)$ given by Lemma~\ref{ref15}, the above sum is bounded by
\[
2R
\frac{\sharp\set{s\in\mathcal{F}_{u(r)}(Q,k,h):\phi'\Lambda_{Q^2s}\notin\Omega_{I(k,h)}^o}}
{\sharp\mathcal{F}_{u(r)}(Q,k,h)}.
\]
By~\eqref{eq7} this last expression tends, for $Q\to\infty$, to
\[
2R\frac{\area\bigl(\Omega_{v(r(k,h))}\setminus\Omega_{I(k,h)}^o\bigr)}{\area\bigl(\Omega_{v(r(k,h))}\bigr)}.
\]
By elementary geometric considerations, the area ratio is bounded by $4\varepsilon_k$. Taking $N(k)=8R\varepsilon_k$ gives the bound on $\abs{\alpha(Q,k,h)}$, and the bound on $\abs{a(Q,k)}$ follows easily.
\end{proof}

We can now conclude the proof of Theorem~\ref{ref19}. Since $P$ is defined via a Riemann integral, we have
\begin{align*}
\int f \,\mathrm{d} P &= \int_0^1\biggl[\int f\,\mathrm{d} P_{v(s)}\biggr] m(s)\,\mathrm{d} s\\
&=\lim_{k\to\infty}\sum_{h=0}^{n(Q_k)-1}
\Biggl(\biggl[\int f\,\mathrm{d} P_{v(r(k,h))}\biggr]\,\lambda(I(k,h))\,m(r(k,h))\Biggr)\\
&=\lim_{k\to\infty}\sum_{h=0}^{n(Q_k)-1}\lim_{Q\to\infty}
\Biggl(\frac{\sum\set{f(\phi'\Lambda_{Q^2s}):s\in\mathcal{F}_{u(r)}(Q,k,h)}}{\sharp\mathcal{F}_{u(r)}(Q,k,h)}
\cdot
\frac{\sharp\mathcal{F}_u(Q,k,h)}{\sharp\mathcal{F}_u(Q)}\Biggr).
\end{align*}
Summing and subtracting
$$
\frac{\sum\set{f(\phi\Lambda_{Q^2s}):s\in\mathcal{F}_u(Q,k,h)}}{\sharp\mathcal{F}_{u(r)}(Q,k,h)}
$$
to the left factor of the parenthesized product above, we get
\begin{multline*}
\int f \,\mathrm{d} P =\\
\lim_{k\to\infty}\sum_{h=0}^{n(Q_k)-1}\lim_{Q\to\infty}
\Biggl(
\frac{\sum\set{f(\phi\Lambda_{Q^2s}):s\in\mathcal{F}_u(Q,k,h)}}{\sharp\mathcal{F}_u(Q)}
+\alpha(Q,k,h)\frac{\sharp\mathcal{F}_u(Q,k,h)}{\sharp\mathcal{F}_u(Q)}
\Biggr)\\
=\lim_{k\to\infty}\lim_{Q\to\infty}
\Biggl(\int f\,\mathrm{d} P(Q)+a(Q,k)\Biggr)
\end{multline*}

By Lemma~\ref{ref31},
$$
\int f\,\mathrm{d} P(Q) - N(k) \le \int f\,\mathrm{d} P(Q) + a(Q,k),
$$
and hence
\begin{multline*}
\biggl(\limsup_{Q\to\infty}\int f\,\mathrm{d} P(Q)\biggr)-N(k)=
\limsup_{Q\to\infty}\biggl(\int f\,\mathrm{d} P(Q)-N(k)\biggr)\\
\le
\limsup_{Q\to\infty}\biggl(\int f\,\mathrm{d} P(Q)+a(Q,k)\biggr)=
\lim_{Q\to\infty}\biggl(\int f\,\mathrm{d} P(Q)+a(Q,k)\biggr).
\end{multline*}
Taking the limit for $k\to\infty$, we get
$$
\limsup_{Q\to\infty}\int f\,\mathrm{d} P(Q) \le \int f\,\mathrm{d} P.
$$
By a dual argument we obtain
$$
\int f\,\mathrm{d} P \le \liminf_{Q\to\infty} \int \,\mathrm{d} P(Q),
$$
thus concluding the proof of Theorem~\ref{ref19}.

redpajama_set_name: RedPajamaArXiv