# Multi-resolution Spatial Regression for Aggregated Data with an Application
to Crop Yield Prediction

Harrison Zhu111Corresponding author., Adam Howes
Imperial College London
{harrison.zhu15<EMAIL_ADDRESS>Owen van Eer, Maxime Rischard
Cervest
{owen<EMAIL_ADDRESS>Dino Sejdinovic
University of Oxford
<EMAIL_ADDRESS>Seth Flaxman
Imperial College London
<EMAIL_ADDRESS>

###### Abstract

We develop a new methodology for spatial regression of aggregated outputs on
multi-resolution covariates. Such problems often occur with spatial data, for
example in crop yield prediction, where the output is spatially-aggregated
over an area and the covariates may be observed at multiple resolutions.
Building upon previous work on aggregated output regression, we propose a
regression framework to synthesise the effects of the covariates at different
resolutions on the output and provide uncertainty estimation. We show that,
for a crop yield prediction problem, our approach is more scalable, via
variational inference, than existing multi-resolution regression models. We
also show that our framework yields good predictive performance, compared to
existing multi-resolution crop yield models, whilst being able to provide
estimation of the underlying spatial effects.

##  1 Introduction

In spatial statistics, the widely-used technique of kriging is equivalent to
what is known as Gaussian process (GP) regression (Rasmussen (2003)) in the
machine learning community. This nonparametric technique has been successfully
applied to crime (Flaxman et al. (2019)), air quality monitoring (Hamelijnck
et al. (2019), remote sensing (Ton et al. (2018)) and epidemiology (Law et al.
(2018a); Bhatt et al. (2017)). Although GP regression becomes computationally
intractable for large datasets due to the $\mathcal{O}(n^{3})$ and
$\mathcal{O}(n^{2})$ complexities of inverting and storing the kernel matrix,
various approximation methods such as variational inference (Titsias (2009);
Hensman et al. (2013, 2015a); Shi et al. (2020)) and the Integrated Nested
Laplace Approximation (INLA) (Rue et al. (2016)) have been used to alleviate
this.

We focus on performing inference in the case where the outputs are available
at a lower spatial resolution than the covariates. Additionally, each of the
covariates may have a different spatial resolution. For crop yield prediction,
crop yields are only observed at the county level over an entire year due to
the scarcity of agricultural census data (Burke et al. (2021)), whereas the
covariates, such as satellite imagery, are abundantly available at daily
temporal and pixel-level (e.g. 500m by 500m squares in space) spatial
resolution. Similarly, in epidemiology, disease diagnostic data (Bhatt et al.
(2017)) is also only available at the census (or aggregated) level due to
privacy reasons. To accomplish prediction at the census level, a
straightforward approach would be to spatially average the covariates within
each census level in order to create a standard supervised learning dataset
(Bhatt et al. (2017); Mateo-Sanchis et al. (2019)). The drawback of this
approach is that it would result in a loss of within region-level variability
and the ability for disaggregation modelling (Law et al. (2018a); Lucas et al.
(2020)). For both problems of crop yield and epidemiology, high-resolution
disaggregation can be very important for policy making as it would allow for
interventions at a high precision. In this work we also present the first
attempt, to the best of our knowledge, at disaggregation for crop yield
modelling.

In the literature, this type of problem has recently been approached in many
ways under the lens of regression with aggregated outputs (Gelfand (2010); Law
et al. (2018a); Hamelijnck et al. (2019); Yousefi et al. (2019); Tanaka et al.
(2019); Zhang et al. (2020)), distribution regression (Law et al. (2018b);
Mateo-Sanchis et al. (2019); Adsuara et al. (2019); Lemercier et al. (2021))
and multiple instance learning (Maron and Lozano-Pérez (1998); Kim and De la
Torre (2010)). Distribution regression has been widely used in the context of
crop yield prediction (Thorns (2018); Mateo-Sanchis et al. (2019); Adsuara et
al. (2019); Lemercier et al. (2021)) and the reason for this is partly because
of its attractive theoretical properties (Szabó et al. (2016)) and ability to
fuse together covariates at different resolutions (Thorns (2018); Adsuara et
al. (2019)).

However, distribution regression is a computationally expensive procedure due
to its $\mathcal{O}(\sum_{i,j=1}^{n}N_{i}N_{j})$ complexity of computing the
kernel matrix, where $N_{i}$ and $N_{j}$ are the number of covariate vectors
for outputs $i$ and $j$. Typically, $N_{i}$ and $N_{j}$ may be of the order of
$10^{2}$ or $10^{3}$, depending on the resolution of the labels, and it will
become computationally intractable to perform hyperparameter search. Previous
methods that proposed to reduce this complexity include the use of random
Fourier features (Rahimi and Recht (2008)), radial basis function (RBF)
networks and Bayesian distribution regression (Law et al. (2018b)), at a
slight cost to the approximation error.

Using the framework of aggregated output regression, where Gaussian process
priors are often used for an underlying regression function, performing
inference is still an expensive procedure. Exact inference or Markov chain
Monte Carlo (MCMC) methods would still involve the
$\mathcal{O}(\sum_{i,j=1}^{n}N_{i}N_{j})$ complexity of computing the
aggregated kernel matrix and would also prohibit hyperparameter tuning. Law et
al. (2018a) and Yousefi et al. (2019) used variational inference to
approximate the Gaussian process posterior, where one compresses all the
information onto a set of $M$ inducing points, and thereby bypassing the need
for the full aggregated kernel matrix computation, reducing the computational
complexity to $\mathcal{O}(M^{3}+M^{2}\sum_{i=1}^{n}N_{i})$ and allowing for
mini-batch training.

#### Contributions:

  * •

We develop a new methodology of performing spatial regression with aggregated
outputs and multi-resolution covariates by building upon VBAgg (Law et al.
(2018a)). Our methodology could be interpreted as a generalisation of Law et
al. (2018a) in the context of using satellite imagery covariates and is an
extension of spatial kriging to the aggregated output setting due to the
addition of an extra spatial effects component.

  * •

We formulate a model, Multi-resolution Variational Bayes with Aggregated
Output Gaussian Processes (MVBAgg), that extends Law et al. (2018a) by using
multiple sets of inducing points for each resolution to perform variational
inference. We thereby also inherit the attractive properties of computational
tractability, ability for GPU-acceleration, uncertainty quantification and
disaggregation of Law et al. (2018a). We then demonstrate that it gives good
prediction performance relative to distribution regression, GPs trained on
centroid covariates and VBAgg on a crop yield dataset located in the Corn Belt
of the US (Adsuara et al. (2019); Mateo-Sanchis et al. (2019)).

##  2 Background

###  2.1 Distribution Regression:

Distribution regression is a method for nonparametric regression (Szabó et al.
(2016)), which has been recently used for crop yield prediction (Thorns
(2018); Mateo-Sanchis et al. (2019); Adsuara et al. (2019)). It assumes that
the data comes in the form of

$\textstyle\\{(\Pi_{i},y_{i})\\}_{i=1}^{n},$

where $y_{i}\in\mathbb{R}$ is a label, $\Pi_{i}$ is a distribution over
covariates $x\in\mathcal{X}\subseteq\mathbb{R}^{d}$, where
$d\in\mathbb{N}_{+}$. In practice, $\Pi_{i}$ is replaced with the empirical
measure $\hat{\Pi}_{i}=N_{i}^{-1}\sum_{j=1}^{N_{i}}\delta_{x_{i,j}}$ for some
observed covariates $\\{x_{i,j}\\}_{j=1}^{N_{i}}$, where $\delta_{x_{i,j}}$ is
a Dirac measure centred at $\\{x_{i,j}\\}$ and $N_{i}\in\mathbb{N}_{+}$.
Suppose that we want to find a predictive model $\hat{f}$ that minimises the
regularised empirical risk

$\textstyle\hat{f}:=\text{argmin}_{f\in\mathcal{H}_{\rho}}\frac{1}{n}\sum_{i=1}^{n}\left(y_{i}-f(\hat{\mu}_{i})\right)^{2}+\lambda||f||_{\mathcal{H}_{\rho}}^{2},$
(1)

where
$\hat{\mu}_{i}\vcentcolon=\int_{\mathcal{X}}k(\cdot,x)\text{d}\hat{\Pi}_{i}(x)=N_{i}^{-1}\sum_{j=1}^{N_{i}}k(\cdot,x_{i,j})$
is the empirical mean embedding of $\hat{\Pi}_{i}$ into the reproducing kernel
Hilbert space (RKHS) $\mathcal{H}_{k}$ for some kernel
$k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$. Solving (1) corresponds
to kernel ridge regression on the mean embeddings (KRRe). Predictive functions
$f$ are taken from $\mathcal{H}_{\rho}$, the RKHS of a second-level kernel
$\rho:\mathcal{H}_{k}\times\mathcal{H}_{k}\rightarrow\mathbb{R}$, and
$\lambda$ is a regularisation parameter. Given any set of covariates
$\\{x_{*,j}\\}_{j=1}^{N_{*}}$, where $N_{*}\in\mathbb{N}_{+}$, and using the
kernel trick, this yields a solution
$\hat{f}(\hat{\mu}_{*})=K_{*}\tilde{\alpha}$, for
$\hat{\mu}_{*}=N_{*}^{-1}\sum_{j=1}^{N_{*}}k(\cdot,x_{*,j})$, where

$\textstyle K_{*}$
$\textstyle\vcentcolon=(\rho(\hat{\mu}_{*},\hat{\mu}_{1}),\ldots,\rho(\hat{\mu}_{*},\hat{\mu}_{n})),$
$\textstyle\tilde{\alpha}$ $\textstyle\vcentcolon=(K+\lambda I_{n})^{-1}y,$

with $(K)_{i,j}=\rho(\hat{\mu}_{i},\hat{\mu}_{j})$. Using the kernel trick
again and with the linear kernel $\rho(f,g)=\langle
f,g\rangle_{\mathcal{H}_{k}}$, this recovers linear regression on the mean
embeddings (LRe), in which
$(K)_{i,j}=(N_{i}N_{j})^{-1}\sum_{l=1}^{N_{i}}\sum_{r=1}^{N_{j}}k(x_{i,l},x_{j,r})$
and $I_{n}$ being the $n$-dimensional identity matrix. Note that the LRe also
employs the kernel trick since mean embeddings are potentially infinite-
dimensional, i.e. no explicit coefficients of the linear model are used.
Similarly, one can also use a non-linear kernel $\rho$, which we call kernel
ridge regression on the mean embeddings (KRRe), and this has been shown to
improve its performance in certain scenarios (Muandet et al. (2012)).

A major disadvantage of distribution regression is the need to perform
$\mathcal{O}(\sum_{i=1}^{n}\sum_{j=1}^{n}N_{i}N_{j})$ number operations to
compute the matrix $K$, and a subsequent $\mathcal{O}(n^{3})$ to invert
$K+\sigma^{2}I_{n}$, therefore prohibiting the ability to efficiently tune the
regularisation parameter $\lambda$ and kernel hyperparameters of $k$ (in
contrast to the hyperparameters of $\rho$, which can be tuned more
efficiently). It is also possible to make use of RBF networks and Bayesian
distribution regression (Law et al. (2018b)) by inducing points
approximations, and RFF to approximate the kernel matrix (Rahimi and Recht
(2008)) for improving the computational scalability.

###  2.2 Gaussian Process for Aggregated Output Regression

GPs can be used to tackle aggregated output regression (Law et al. (2018a);
Hamelijnck et al. (2019); Tanaka et al. (2019); Yousefi et al. (2019)). This
is a common setting in spatial statistics, where each label is often
associated with a set of observations within a region (areal unit). Again,
suppose we have some covariate space $\mathcal{X}$ and a response space
$\mathbb{R}$. The aggregated output regression model (Tanaka et al. (2019)) is
defined as

$\textstyle
y_{i}=\int_{\mathcal{X}_{i}}f(x)\text{d}\Pi_{i}(x)+\epsilon_{i},\qquad\epsilon_{i}\sim\mathcal{N}(0,\sigma^{2}),$
(2)

where $\mathcal{X}_{i}\subseteq\mathcal{X}$ is an observation region with
distribution $\Pi_{i}$ with Lebesgue density $\pi_{i}$ and $\epsilon_{i}$ is
an independently distributed Gaussian noise with $\sigma>0$. We denote
$\mathcal{N}(a,b)$ and $N(\cdot;a,b)$ as the distribution and density function
of a Gaussian distribution with mean $a$ and covariance $b$.

In practice, we will want to approximate
$\int_{\mathcal{X}_{i}}f(x)\text{d}\Pi_{i}(x)$ using some available covariates
$\\{x_{i,j}\\}_{j=1}^{N_{i}}$ drawn from $\Pi_{i}$. In general, this integral
could be approximated by constructing quadrature weights $w_{i,j}$ and states
$\\{x_{i,j}\\}_{j=1}^{N_{i}}$ and integrating with respect to
$\hat{\Pi}_{i}=\sum_{j=1}^{N_{i}}w_{i,j}\delta_{x_{i,j}}$ such that
$\sum_{j=1}^{N_{i}}w_{i,j}=1$. Law et al. (2018a) used deterministic survey
weights and Hamelijnck et al. (2019) constructed the weights using a Gaussian
process regression network, but it is not clear how these will improve the
estimation accuracy of $\int_{\mathcal{X}_{i}}f(x)\text{d}\Pi_{i}(x)$ and so
we will assume that $w_{i,j}=N_{i}^{-1}$ for all $j=1,\ldots,N_{i}$. In the
context of classification of images, alternative aggregation methods are max-
aggregation or attention-aggregation (Kim and De la Torre (2010); Haussmann et
al. (2017); Ilse et al. (2018)).

Using this setup, we thus work with the dataset

$\textstyle\\{(\\{x_{i,j},w_{i.j}\\}_{j=1}^{N_{i}},y_{i})\\}_{i=1}^{n}.$

Using $\hat{\Pi}_{i}$, Equation 2 becomes

$\textstyle
y_{i}=\sum_{j=1}^{N_{i}}w_{i,j}f(x_{i,j})+\epsilon_{i},\qquad\epsilon_{i}\sim\mathcal{N}(0,\sigma^{2}),$
(3)

where $f$ is some underlying regression function and $\sigma>0$. We can posit
a GP prior $f\sim\mathcal{G}\mathcal{P}(m,k)$, where $m$ and $k$ are the mean
and covariance functions, in order to infer $f$. We will denote the likelihood
function of (3) as $p(y_{i}|f_{i})$, where
$f_{i}\vcentcolon=(f(x_{i,1}),\ldots,f(x_{i,N_{i}}))^{\intercal}$, and denote
the weights $w_{i}\vcentcolon=(w_{i,1},\ldots,w_{i,N_{i}})^{\intercal}$.
Similarly, denote $f_{*}\vcentcolon=(f_{*,1},\ldots,f_{*,N_{*}})^{\intercal}$
for prediction points $\\{x_{*,j},w_{*,j}\\}_{j=1}^{N_{*}}$ with true label
$y_{*}$. Suppose further that
$X_{i}\vcentcolon=(x_{i,1},\ldots,x_{i,N_{i}})^{\intercal}$,
$X_{*}\vcentcolon=(x_{*,1},\ldots,x_{*,N_{*}})^{\intercal}$ and
$k_{X_{i},X_{j}}:=(k(x_{i,l},x_{j,m}))_{l,m}$ for pairs $i,j$, and
$l=1,\ldots,N_{i}$ and $m=1,\ldots,N_{j}$,

The log-likelihood is $\log p(y_{i}|f_{i})$ is

$\textstyle\log
p(y_{i}|f_{i})=-\frac{1}{2\sigma^{2}}(y_{i}-w_{i}^{\intercal}f_{i})^{2}-\frac{1}{2}\log\left(2\pi\sigma^{2}\right).$

Since we have

$\textstyle\begin{pmatrix}y_{1}\\\ \vdots\\\ y_{n}\\\
w_{*}^{\intercal}f_{*}\end{pmatrix}=\mathcal{N}\Bigg{(}\begin{pmatrix}w_{1}^{\intercal}m_{X_{1}}\\\
\vdots\\\ w_{n}^{\intercal}m_{X_{n}}\\\
w_{*}^{\intercal}m_{X_{*}}\end{pmatrix},\begin{pmatrix}K+\sigma^{2}I_{n}&K_{*}\\\
K_{*}^{\intercal}&K_{*,*}\end{pmatrix}\Bigg{)},$

where $(K)_{ij}=(w_{i}^{\intercal}k_{X_{i},X_{j}}w_{j})_{ij}$,
$(K_{*})_{i}\vcentcolon=(w_{i}^{\intercal}k_{X_{i},X_{*}}w_{*})_{i}$ and
$K_{*,*}\vcentcolon=w_{*}^{\intercal}k_{X_{*},X_{*}}w_{*}$, by standard
identities of Gaussian distributions we have the conditional distribution

$\textstyle
w_{*}^{\intercal}f_{*}|\\{(\\{x_{i,j},w_{i.j}\\}_{j=1}^{M_{i}},y_{i})\\}_{i=1}^{n}\sim\mathcal{N}(\tilde{m}_{X_{*}},\tilde{k}_{X_{*},X_{*}}),$

with

$\textstyle\tilde{m}_{X_{*}}$
$\textstyle\vcentcolon=w_{*}^{\intercal}m_{X_{*}}+K_{*}^{\intercal}(K+\sigma^{2}I_{n})^{-1}(y-m_{X}),$
$\textstyle\tilde{k}_{X_{*},X_{*}}$
$\textstyle\vcentcolon=K_{*,*}-K_{*}^{\intercal}(K+\sigma^{2}I_{n})^{-1}K_{*}.$

With a zero mean function, the posterior mean obtained here is the exact
expression from LRe, similar to the connection between standard, non-
aggregated GP regression and kernel ridge regression (Kanagawa et al. (2018)).
With this, we are still prohibited by the large computational burden in the
kernel matrix operations.

On the other hand, if $\Pi_{i}$ is known, then we would obtain

$\textstyle\begin{pmatrix}y_{1}\\\ \vdots\\\ y_{n}\\\
\int_{\mathcal{X}_{*}}f_{*}(u)\pi_{*}(x)\text{d}x\end{pmatrix}=\mathcal{N}\left(\begin{pmatrix}\int_{\mathcal{X}_{1}}m(x)\pi_{1}(x)\text{d}x\\\
\vdots\\\ \int_{\mathcal{X}_{n}}m(x)\pi_{n}(x)\text{d}x\\\
\int_{\mathcal{X}_{*}}m(x)\pi_{*}(x)\text{d}x\end{pmatrix},\begin{pmatrix}K+\sigma^{2}I_{n}&K_{*}\\\
K_{*}^{\intercal}&K_{*,*}\end{pmatrix}\right),$

where

$\textstyle(K)_{ij}$
$\textstyle\vcentcolon=\int_{\mathcal{X}_{i}}\int_{\mathcal{X}_{j}}k(u,v)\pi_{i}(u)\pi_{j}(v)\text{d}u\text{d}v,$
$\textstyle(K_{*})_{i}$
$\textstyle\vcentcolon=\int_{\mathcal{X}_{*}}\int_{\mathcal{X}_{i}}\int_{\mathcal{X}_{i}}k(u,v)\pi_{*}(u)\pi_{i}(v)\text{d}u\text{d}v,$
$\textstyle K_{*,*}$
$\textstyle\vcentcolon=\int_{\mathcal{X}_{*}}\int_{\mathcal{X}_{*}}k(u,v)\pi_{*}(u)\pi_{*}(v)\text{d}u\text{d}v.$

The posterior distribution over the integral
$\int_{\mathcal{X}_{*}}f(x)\text{d}\Pi_{*}(x)$ would thus be

$\textstyle\int_{\mathcal{X}_{*}}f(x)\text{d}\Pi_{*}(x)\text{d}x\sim\mathcal{N}(\tilde{m},\tilde{k}),$

with

$\textstyle\tilde{m}$
$\textstyle\vcentcolon=\int_{\mathcal{X}_{*}}m(x)\pi_{*}(x)\text{d}x+K_{*}^{\intercal}(K+\sigma^{2}I_{n})^{-1}(y-m_{X}),$
$\textstyle\tilde{k}$
$\textstyle\vcentcolon=K_{*,*}-K_{*}^{\intercal}(K+\sigma^{2}I_{n})^{-1}K_{*}.$

Having learned the posterior distribution of $f$, it is then possible to make
predictions and provide statistical insights at the individual covariate level
(Law et al. (2018a)).

###  2.3 Variational Bayes with Gaussian Processes (VBAgg)

Although kernel ridge and exact Gaussian process regression provide elegant
inference procedures for learning from aggregated outputs, the computational
complexity prohibits the ability to tune hyperparameters in a data-driven way.
Following Law et al. (2018a) and Yousefi et al. (2019), which use the Sparse
Variational GP (SVGP) framework of Titsias (2009); Hensman et al. (2013) and
Hensman et al. (2015a), we will use a set of $M$ inducing points
$u=(f(z_{1}),\ldots,f(z_{M}))^{\intercal}$ at locations
$Z=(z_{1},\ldots,z_{M})^{\intercal}$, with
$z_{1},\ldots,z_{M}\in\mathbb{R}^{d}$. We pose a variational distribution
$q(f,u)\vcentcolon=p(f|u)q(u)$ and minimise $\text{KL}(q(f,u)||p(f,u|y))$ to
give the variational approximation

$\textstyle
q(u)=\text{argmin}_{q\in\mathcal{Q}}\,\text{KL}(q(u)||p(u))-\sum_{i=1}^{n}\int\log\\{p(y_{i}|f_{i})\\}q(u)p(f_{i}|u)\,\text{d}u\,\text{d}f_{i},$

that is close to the true posterior $p(f,u|y)$. We choose the variational
family
$\mathcal{Q}=\\{q(u)=N(u;\eta_{Z},\Sigma_{Z}):\eta_{Z}\in\mathbb{R}^{M},\Sigma_{Z}\in\mathbb{R}^{M\times
M}\\}$. Since $p(f_{i}|u)=N(\hat{m}_{X_{i}},\hat{k}_{X_{i},X_{i}})$ with

$\textstyle\hat{m}_{X_{i}}$
$\textstyle=m_{X_{i}}+k_{X_{i},Z}k_{Z,Z}^{-1}(u-m_{Z}),$
$\textstyle\hat{k}_{X_{i},X_{i}}$
$\textstyle=k_{X_{i},X_{i}}-k_{X_{i},Z}k_{Z,Z}^{-1}k_{Z,X_{i}},$

the approximate posterior is thus $q(f_{i})\vcentcolon=\int
p(f_{i}|u)q(u)\text{d}u=N(f_{i};\tilde{m}_{X_{i}},\tilde{k}_{X_{i},X_{i}})$
where

$\textstyle\tilde{m}_{X_{i}}$
$\textstyle=m_{X_{i}}+k_{X_{i},Z}k_{Z,Z}^{-1}(\eta_{Z}-m_{Z}),$
$\textstyle\tilde{k}_{X_{i},X_{i}}$
$\textstyle=k_{X_{i},X_{i}}-k_{X_{i},Z}(k_{Z,Z}^{-1}-k_{Z,Z}^{-1}\Sigma_{Z}k_{Z,Z}^{-1})k_{Z,X_{i}}.$

This thus allows the lower bound to be written down explicitly as

$\textstyle\mathcal{L}\vcentcolon=\text{KL}(q(u)||p(u))-\frac{1}{2}\sum_{i=1}^{n}\bigg{(}\frac{y_{i}^{2}-2y_{i}w_{i}^{\intercal}\tilde{m}_{X_{i}}+w_{i}^{\intercal}(\tilde{k}_{X_{i},X_{i}}+\tilde{m}_{X_{i}}\tilde{m}_{X_{i}}^{\intercal})w_{i}}{\sigma^{2}}-\log\left(2\pi\sigma^{2}\right)\bigg{)},$

and solved via optimisation methods in mini-batches (via a stochastic
approximation to $\mathcal{L}$. See Zhang et al. (2019) for a review). One
could also sample the optimal posterior distribution, in which
KL$(q(f,u)||p(f,u|y))=0$, directly using MCMC, as proposed in Hensman et al.
(2015b). To make a prediction for an arbitrary region with $N_{*}$ inputs
$X_{*}$, this simply requires the posterior distribution

$\textstyle
q(w_{*}^{\intercal}f_{*})=N(w_{*}^{\intercal}f_{*};w_{*}^{\intercal}\tilde{m}_{X_{*}},w_{*}^{\intercal}\tilde{k}_{X_{*},X_{*}}w_{*}).$

##  3 Aggregated Output Spatial Regression

In spatial statistics, we often make use of spatially-resolved covariates,
where we implicitly assume that they are a function of the spatial
coordinates. Suppose that $\mathcal{S}\subseteq\mathbb{R}^{2}$ represents a
spatial domain and $g:\mathcal{S}\rightarrow\mathbb{R}^{d}$, where $d>0$, is
an unknown function that represents a covariate. Then $g$, for example, could
be the solution of the heat equation if it represents temperature, although in
reality we may only observe a noisy version of $g$ collected via remote
sensing methods, such as satellite sensors and cameras. We can thus formulate
aggregated output spatial regression as

$\textstyle
y_{i}=\frac{1}{|\mathcal{S}_{i}|}\int_{\mathcal{S}_{i}}f^{1}(g(s))\text{d}s+\frac{1}{|\mathcal{S}_{i}|}\int_{\mathcal{S}_{i}}f^{0}(s)\text{d}s+\epsilon_{i},\qquad\epsilon_{i}\sim\mathcal{N}(0,\sigma^{2}),$
(4)

where $\mathcal{S}_{i}\subseteq\mathcal{S}$, and $f^{1}$ and $f^{2}$ are
functions indicating effects from remote sensing covariates and additional
spatial variability. Note that the first integral can also be written as
$\frac{1}{|\mathcal{S}_{i}|}\int_{\mathcal{T}_{i}}f^{1}(t)\text{d}(g_{\sharp}\lambda)(t)$,
where $\lambda$ is the Lebesgue measure, $g_{\sharp}\lambda$ is the push-
forward measure and $g^{-1}(\mathcal{T}_{i})=\mathcal{S}_{i}$, meaning that we
can use the aggregated output models as previously discussed. It is
interesting to note that $\lambda$ could also be other measures, such as the
spherical measure, a generalisation of the Lebesgue measure on the sphere
$\mathbb{S}^{2}\subset\mathbb{R}^{3}$ and may be a more realistic measure for
modelling Earth observations.

In the context of images, such as ones obtained via satellites, each datapoint
represents a pixel, which is a square over space, as shown in Figure 1. It
should be noted that this limits the number of points that one can use for
performing numerical integration as it would depend on the spatial resolution
of the covariates. For each pixel, suppose we have a representative latitude-
longitude pair $[s]$. For a subset with latitude and longitude
$\mathcal{S}\subseteq\mathbb{R}^{2}$ and Lebesgue measure $\lambda$, this
means that we have to discretise the domains $\mathcal{S}_{i}$ into a discrete
measure space $\widehat{\mathcal{S}}_{i}$ of the representatives $[s]$, with
measure $\hat{\lambda}_{i}(S)=\frac{1}{|\hat{\mathcal{S}}_{i}|}\sum_{s\in
S}\mathbbm{1}_{\hat{\mathcal{S}}_{i}}(s)$ for any
$S\subseteq\widehat{\mathcal{S}}_{i}$. This changes Equation 4 to

$\textstyle
y_{i}=\frac{1}{|\widehat{\mathcal{S}}_{i}|}\sum_{s\in\widehat{\mathcal{S}}_{i}}f^{1}(g(s))+\frac{1}{|\mathcal{S}_{i}|}\int_{\mathcal{S}_{i}}f^{0}(s)\text{d}s+\epsilon_{i},\qquad\epsilon_{i}\sim\mathcal{N}(0,\sigma^{2}),$
(5)

where we note that we can still approximate the second integral with
arbitrarily many points in space without discretisation as the spatial
covariates have infinite resolution. In constrast, Law et al. (2018a) only
allows for a single integral and the number of spatial points must match the
number of pixels for the spatially-resolved covariates.

Figure 1: Pixel squares of the MCD12C1 land cover mask (Friedl and Sulla-
Menashe (2015)) over Adams county, Illinois, USA. Note that some of the data
is missing, as represented by the white squares.

As it is the case that the covariate effects from $g(s)$ and spatial effects
from $s$ do not need to be treated at the same resolution, various spatially-
resolved covariates may also correspond to different resolutions. In previous
aggregated output modelling works, the covariates are usually preprocessed to
the same resolution (Mateo-Sanchis et al. (2019)). This means that when using
aggregate output methods, the integrals of the mean and covariance functions
will be integrated over a discrete measure with fewer points, giving less
accurate and higher variance approximations. On the other hand, up-sampling
may result in errors from the approximation algorithms used, introducing bias
into the integrals. This thus motivates the use of multi-resolution models,
where the integrals for different resolution are separated without up-sampling
or down-sampling by assuming an additive structure. Suppose we have $D$
covariates, each with a different spatial resolution. By treating each
resolution separately, we can modify Equation 5 and obtain

$\textstyle
y_{i}=\sum_{l=1}^{D}\frac{1}{|\widehat{\mathcal{S}}_{i,l}|}\sum_{s\in\widehat{\mathcal{S}}_{i,l}}f^{k}(g_{l}(s))+\frac{1}{|\mathcal{S}_{i}|}\int_{\mathcal{S}_{i}}f^{0}(s)\text{d}s+\epsilon_{i},\qquad\epsilon_{i}\sim\mathcal{N}(0,\sigma^{2}),$

for resolutions $l=1,\ldots,D$ and
$g_{l}:\mathcal{S}\rightarrow\mathbb{R}^{d_{l}}$ for $d_{l}>0$. Note that we
here consider an additive model across covariates corresponding to different
resolutions. This results in the additive covariance kernel structure. It is
also possible to formulate a product covariance kernel structure, but we do
not pursue it further in this work.

###  3.1 Multi-resolution Variational Bayes with Aggregated Output Gaussian
Processes (MVBAgg)

A natural extension to VBAgg is to include covariates of multiple sets and
hence observing the dataset

$\textstyle\\{(\\{\\{x_{i,j,l},w_{i,j,l}\\}_{j=1}^{N_{i,l}}\\}_{l=1}^{D},y_{i})\\}_{i=1}^{n},$

where $l=1,\ldots,D$ is the index of the covariates. Adsuara et al. (2019)
proposes incorporating sets of covariates at different resolutions, in the
context of LRe, by considering kernel regression functions in the space equal
to the direct sum of RKHSs with kernels $\rho_{1},\ldots,\rho_{D}$, where
$\rho_{l}:\mathcal{H}_{k_{l}}\times\mathcal{H}_{k_{l}}\rightarrow\mathbb{R}$.
From a GP point of view, we may learn similar functions by modifying Equation
3 to

$\textstyle
y_{i}|\\{\\{x_{i,j,l},w_{i,j,l}\\}_{j=1}^{N_{i,l}}\\}_{l=1}^{D}\sim\mathcal{N}\left(\sum_{l=1}^{D}\sum_{j=1}^{N_{i}}w_{i,j,l}f^{l}(x_{i,j,l}),\left(\sum_{l=1}^{D}\sum_{j=1}^{N_{i,l}}w_{i,j,l}^{2}\right)\sigma^{2}\right),$

where $f^{l}$ is a function induced with a GP prior
$\mathcal{G}\mathcal{P}(0,k_{l})$ with kernel
$k_{l}:\mathcal{X}_{i,l}\times\mathcal{X}_{i,l}\rightarrow\mathbb{R}$, where
$\mathcal{X}_{i,l}$ is the covariate space for the $l$th resolution in the
region indexed by $i$. We assume that the $f^{l}$’s are mutually independent.
To use a VBAgg approximation framework, we pose the variational distribution
$q(f^{1},\ldots,f^{D},u_{1},\ldots,u_{D})\vcentcolon=\prod_{l=1}^{D}p(f^{l}|u_{l})q(u_{l})$
with different inducing points
$u_{l}=(f^{l}(z_{1,l}),\ldots,f^{l}(z_{L_{l},l}))^{\intercal}$ at each sets of
locations $Z_{l}=(z_{1,l},\ldots,z_{L_{l},l})^{\intercal}$. We can thus obtain
the approximate posterior $q(f^{l})$ for each $l$ as
$q(f^{l}_{i})=N(f^{l}_{i};\tilde{m}_{X_{i,l}},\tilde{k}_{X_{i,l},X_{i,l}})$

$\textstyle\tilde{m}_{X_{i,l}}$
$\textstyle=m_{X_{i,l}}+k_{X_{i,l},Z_{l}}k_{Z_{l},Z_{l}}^{-1}(\eta_{Z_{l}}-m_{Z_{l}}),$
$\textstyle\tilde{k}_{X_{i,l},X_{i,l}}$
$\textstyle=k_{X_{i,l},X_{i,l}}-k_{X_{i,l},Z_{k}}(k_{Z_{l},Z_{l}}^{-1}-k_{Z_{l},Z_{l}}^{-1}\Sigma_{Z_{l}}k_{Z_{l},Z_{l}}^{-1})k_{Z_{l},X_{i,l}},$

where $X_{i,l}:=(x_{i,1,l},\ldots,x_{i,N_{i,l},l})^{\intercal}$. The lower
bound can then be written as

$\textstyle\mathcal{L}$
$\textstyle\vcentcolon=\sum_{l=1}^{D}\text{KL}(q(u_{l})||p(u_{l}))$
$\textstyle-\frac{1}{2}\sum_{i=1}^{n}\bigg{(}\frac{1}{\sigma^{2}}[y_{i}^{2}-2y_{i}\sum_{l=1}^{D}w_{i,l}^{\intercal}\tilde{m}_{X_{i,l}}+\sum_{l=1}^{D}w_{i,l}^{\intercal}(\tilde{k}_{X_{i,l},X_{i,l}}+\tilde{m}_{X_{i,l}}\tilde{m}_{X_{i,l}}^{\intercal})w_{i,l}]$
$\textstyle-\log\left(2\pi\sigma^{2}\right)\bigg{)},$

where $w_{i,l}:=(w_{i,1,l},\ldots,w_{i,N_{i,l},l})^{\intercal}$. Similarly,
the predictive distribution for $y_{*}$ would be

$\textstyle
q(\sum_{l=1}^{D}w_{*,l}^{\intercal}f_{*}^{l})=N(\sum_{l=1}^{D}w_{*,l}^{\intercal}f_{*}^{l};\sum_{l=1}^{D}w_{*,l}^{\intercal}\tilde{m}_{X_{*,l}},\sum_{l=1}^{D}w_{*,l}^{\intercal}\tilde{k}_{X_{*,l},X_{*,l}}w_{*,l}).$

To initialise the inducing points $Z_{k}$’s, we follow the setup in VBagg (Law
et al. (2018a)) and pick 1 point from each region $i$ using the cluster
centres with KMeans. We can also set the inducing points to be trainable or
fixed. Inheriting the nice properties of VBAgg, we can attain
$\mathcal{O}(\sum_{l=1}^{D}L_{l}^{3}+L_{l}^{2}\sum_{i=1}^{n}N_{i,l})$
complexity in time, can performing mini-batch training and make use of GPU-
acceleration. In comparison, Adsuara et al. (2019) gives
$\mathcal{O}(\sum_{l=1}^{D}\sum_{i,j=1}^{n}N_{i,l}N_{j,l})$ complexity.
Although it is possible to employ RFF to give a reduction to
$\mathcal{O}(\sum_{l=1}^{D}L_{l}^{3}+L_{l}^{2}\sum_{i=1}^{n}N_{i,l})$, where
$L_{l}$ is now the dimension of the RFF map, it still cannot perform kernel
parameter tuning without cross-validation, whereas we propose to tune the
parameters in a data-driven way via the variational lower bound. The full
inference scheme is detailed in Algorithm 1.

Algorithm 1 Training of Multi-resolution Variational Bayes with Aggregated
Output Gaussian Processes (MVBAgg)

Input
$\\{\\{x_{i,j,l},w_{i,j,l}\\}_{j=1}^{N_{i,l}}\\}_{l=1}^{D},y_{i})\\}_{i=1}^{n}$
(dataset), $\\{z_{l}\\}_{l=1}^{D}$ (landmark points),
$\\{\eta_{Z_{l}}\\}_{l=1}^{D}$, $\\{\Sigma_{Z_{l}}\\}_{l=1}^{D}$ (variational
distribution parameters), kernels $\\{k_{l}\\}_{l=1}^{D}$, noise $\sigma$,
total iterations $L$, batch size $B$

1:Initialise ELBO $\mathcal{L}\leftarrow 0$

2:for $i=1,\ldots,L$ do

3: Take a mini-batch
$\\{(\\{\\{x_{i,j,l},w_{i,j,l}\\}_{j=1}^{N_{i,l}}\\}_{l=1}^{D},y_{i})\\}_{i\in
B}$ with indexing set $B\subset\\{1,\ldots,n\\}$

4: Compute lower bound $\mathcal{L}_{B}$ with this batch

5: $\mathcal{L}\leftarrow\mathcal{L}+\frac{n}{B}\mathcal{L}_{B}$

6: if $i\mod{n/B}=0$ then

7: Backpropagate $\mathcal{L}$ and update model parameters

8: end if

9:end for

10:return
$\\{k_{l}\\}_{l=1}^{D},\sigma,\\{\eta_{Z_{l}}\\}_{l=1}^{D},\\{\Sigma_{Z_{l}}\\}_{l=1}^{D}$

Returning to the context of aggregated outputs in the spatial context,
alongside sets of remote sensing covariates, we would also include a set of
latitude and longitude coordinates that we can sample uniformly on the
continuous spatial domain. This yields the spatial GP $f^{0}$ and the
covariate GPs $f^{1},\ldots,f^{D}$. As motivated in Law et al. (2018a), the
VBAgg model is also capable of performing disaggregation by learning an
analytic form of the approximate posterior distribution on the underlying
function $f$. In a similar manner, with MVBAgg we can also study each of the
GPs at the pixel level using the posterior distribution of each $f^{l}$.

##  4 Application to Crop Yield prediction

The task of crop yield prediction has recently been attempted by scientists
using a variety of methods (Bolton and Friedl (2013); Kuwata and Shibasaki
(2015); Lobell et al. (2015); You et al. (2017); Thorns (2018); Mateo-Sanchis
et al. (2019); Adsuara et al. (2019); Yang et al. (2019); Lemercier et al.
(2021)). The wide availability of remote sensing and geophysical data from
NASA programmes such as Moderate Resolution Imaging Spectroradiometer (MODIS)
(Didan (2015)) and Soil Moisture Active Passive (SMAP) (Entekhabi et al.
(2010)) have also helped with this effort by providing modellers with pixel-
resolution information on croplands. By utilising the vast amount of data,
such as vegetation indices, that have biological links to crop health and
hence to yield (Mateo-Sanchis et al. (2019)), as covariates of the croplands,
modellers can formulate crop yield prediction as a supervised learning
problem. However, existing works on crop yield modelling using aggregated
output/distribution regression methods (You et al. (2017); Thorns (2018);
Mateo-Sanchis et al. (2019); Adsuara et al. (2019)) do not consider the model
uncertainty and interpretability, which is arguably as important as performing
crop yield prediction itself.

Our work is most closely related to kernel methods for crop yield modelling,
notably multi-source kernel ridge regression (on the mean embeddings) (Adsuara
et al. (2019)), kernel ridge regression (Mateo-Sanchis et al. (2019)) and
Bayesian distribution regression (Thorns (2018)). Another very popular
approach is to extract signals from the image histograms of the pixel values
using convolutional neural networks, as attempted in the Deep Gaussian process
for crop yield prediction work (You et al. (2017)). However, this approach is
not suitable for accounting for the underlying spatial structure and
variability due to the the use of image histograms, whereas our method keeps
the pixel-level spatial structure. Our aim is to both predict the crop yield
at the county-level and disaggregate the underlying spatial random surface at
the pixel-level

###  4.1 Data

We use a county-level crop yield dataset from Mateo-Sanchis et al. (2019) with
yield in the U.S. Corn Belt in 2015, vegetation optical depth (VOD) (Konings
et al. (2017)) derived from the SMAP products (Entekhabi et al. (2010)) and
Enhanced Vegetation Index (EVI) from MODIS (Didan (2015)). Unfortunately, the
covariates in this study are in the same spatial resolution due to the
downsampling of the EVI covariates and the data processing procedure is non-
trivial, but we can nonetheless obtain a multi-resolution dataset by adding
sets of 100 spatial covariates (latitude + longitude in EPSG:4326 projection
format) per county. Following Sanchis et al. (2019), we use only use
covariates during the period April 2015 - October 2015, which covers the
period from growth to just before harvest. For VOD and EVI respectively, this
gives 213 and 13 observations for each spatial location respectively. We
filter out regions where the geometry of the yield has not been provided,
giving 375 outputs as shown in Figure 2.

###  4.2 Experimental Setup

For our experiments, we will assess the performance of the predictive models
via 5-fold cross-validation, reporting the Root Mean Squared Error (RMSE)
$\sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_{i}-\hat{y}_{i})^{2}}$ and Mean Average
Percentage Error (MAPE)
$\frac{1}{n}\sum_{i=1}^{n}|(y_{i}-\hat{y}_{i})/y_{i}|$, averaged over 5
validation sets, with 300 and 75 training and test labels respectively per
fold. Our benchmarks are Linear Regression (LR), Gaussian process (GP)
regression, linear regression on the mean embeddings (LRe), kernel ridge
regression on the mean embeddings (KRRe) and VBAgg. We compare these to
MVBagg. For LR and GP, we use the centroid covariates, where we spatially
average the covariates to obtain a standard supervised learning dataset. For
all other models, we use the “stacked covariates” (Thorns (2018); Mateo-
Sanchis et al. (2019); Adsuara et al. (2019)), where each column represents an
observation at a certain time during the year and each row presents an
observation at a certain spatial location. For MVBAgg, VBAgg and GP, the VOD
and EVI covariates are treated in the same resolution with additive RBF
kernels, defined as:

$\textstyle
k(x,y)=\sigma^{2}\exp\left(-\frac{||x-y||_{2}^{2}}{2\ell^{2}}\right),\qquad\forall
x,y\in\mathcal{X},$

where $\sigma>0$ is the scale, $\ell>0$ is the lengthscale and $||\cdot||_{2}$
denotes the Euclidean norm. Similarly, the spatial components are specified
with a Matérn-3/2 kernel:

$\textstyle
k(x,y)=\sigma^{2}\left(1+\frac{\sqrt{3}||x-y||_{2}}{\ell}\right)\exp\left(-\frac{\sqrt{3}||x-y||_{2}}{\ell}\right).$

We fit LR and GP using maximum likelihood estimation, and LRe and KRRe using
the median heuristic for the lengthscale222As there are over 20,000 covariate
values, we take a maximum of 10 values per county, which reduces this number
to around 3,000, and then compute the subsequent median lengthscale. (Garreau
et al. (2018)) and the scale $\sigma=1$, with a regularisation parameter of
$\lambda=0.1$. In addition, LRe and KRRe use the linear and RBF second-level
kernel respectively, in which for the latter we tune the lengthscale also
using the median heuristic (of the difference of the mean embedding norms
$||\hat{\mu}_{i}-\hat{\mu}_{j}||_{\mathcal{H}_{\rho}}$ for counties $i,j$) and
set $\sigma=1$. We use the Adam optimiser with a learning rate of 0.001,
20,000 training iterations and the inducing points chosen by taking the centre
of a K-Means algorithm (Oglic and Gärtner (2017)) in each county for MVBAgg
and VBAgg. We relied on GPFlow (De G. Matthews et al. (2017)) and NumPy
(Harris et al. (2020)) to implement our models. The code is available at
https://github.com/ImperialCollegeLondon/MVBAgg. For VBAgg and MVBAgg, we used
2 Nvidia RTX 2080 GPUs, and for all the other methods, we used Intel(R)
Xeon(R) Gold 6242 CPU @ 2.80GHz CPUs, both on the NVIDIA4 GPU Compute Server
provided by the Department of Mathematics, Imperial College London.

Figure 2: Crop yield ($kg/m^{2}$) in the corn belt of the United States.

###  4.3 Results

We can see from Table 1 that the out-of-sample predictive performance for
MVBAgg is slightly better than GP, LRe and KRRe in terms of RMSE and MAPE. Due
to the larger difference in MAPE, we can interpret this as MVBAgg giving
better individual county predictions on average. In many applications, GP with
centroid covariates is the commonly-used approach (Bhatt et al. (2017); Mateo-
Sanchis et al. (2019); Martinez-Ferrer et al. (2020)) and our results suggest
that there are gains in both predictive performance and reduction in
uncertainty when using aggregated output methods.

Figure 3: Prediction versus Observed Outputs with 95$\%$ posterior credible
intervals on a test set of 20$\%$ of the points. 86.7$\%$ and 65.3$\%$ of the
observations lie in the 95$\%$ and 80$\%$ credible intervals respectively.

Method Covariates RMSE MAPE Runtime (s) LR Centroid VOD + EVI 0.321$\pm$0.0446
0.411$\pm$0.0589 0.0457$\pm$0.00521 GP Centroid VOD + EVI + space
0.0766$\pm$0.00782 0.0891$\pm$0.0113 1.26$\pm$0.0937 LRe Stacked VOD + EVI
0.0752$\pm$0.00457 0.0846$\pm$0.00732 364$\pm 1.14$ LRe Stacked VOD + EVI +
space 0.0726$\pm$0.00478 0.0798$\pm$0.00695 598$\pm 2.81$ KRRe Stacked VOD +
EVI 0.0792$\pm$0.00514 0.0888$\pm$0.00828 1140$\pm 5.46$ KRRe Stacked VOD +
EVI + space 0.0726$\pm$0.00554 0.0806$\pm$0.00826 1750$\pm 14.9$ VBAgg Stacked
VOD + EVI 0.0758$\pm$0.00411 0.0850$\pm$0.00647 288$\pm 1.14$ MVBAgg Stacked
VOD + EVI + space 0.0722$\pm$0.00454 0.0763$\pm$0.00619 478$\pm 0.485$

Table 1: Cross-validated predictive and average run-time performance of the
models on the crop yield data. We also report the standard error over 5-fold
cross-validation of the RMSE and MAPE, as indicated as se-RMSE and se-MAPE
respectively, and well as the standard error over the run-time.
Figure 4: Absolute prediction error ($kg/m^{2}$) of MVBAgg on 75 counties that
were held out as test set. Figure 5: The spatial effects component of MVBAgg
at the individual level after using the whole dataset for training.

LRe, the existing benchmark, requires hyperparameter tuning via cross-
validation over a grid of configurations without the median heuristic. We also
find that with our computational hardware, a full inference with VBAgg and
MVBAgg is faster than (See Table 1) fitting LRe or KRRe on a single set of
hyperparameters, which confirms that our method is indeed more scalable.

As shown in Figure 4 of the prediction on 75 counties held out as test set,
our model is unable predict well on Lincoln County, Minnesota which is an
outlier as can be seen in Figure 2. However, in Figure 3, we can see that most
of the 95$\%$ posterior credible intervals for the test points contain the
true crop yield. Finally, by looking at the posterior mean of the spatial
component at the pixel-level in Figure 5, it is also interesting to see that
we are able to uncover the additional spatial heterogeneity: for example,
regions in Nebraska exhibit within-county variability as seen via the colour
gradient within each county. This type of disaggregated information could be
useful for agricultural policy makers to study where to impose interventions
to improve productivity.

##  5 Conclusion

We developed a framework for spatial regression of aggregated outputs on
multi-resolution covariates. Using Gaussian processes, we derived a multi-
resolution version of VBAgg (Law et al. (2018a)), giving at least as good
performance on a real world crop yield dataset as existing multi-resolution
crop yield prediction models. In addition, our model is capable of performing
inference on its model hyperparameters at scale and in a data-driven manner
due to the variational formulation and the ability for GPU-acceleration,
whereas kernel ridge regression requires tuning via cross-validation or the
median-heuristic. We were also able to reconstruct the underlying fine-grained
spatial effects by looking at the spatial component of MVBAgg and prediction
uncertainty and we believe that this information would be particularly useful
for organisations making policy interventions. We hope that this model could
be useful for agriculture domain-experts who have access to other interesting
crop yield datasets.

For future work, it would be interesting to consider our model for other
problems such as disease mapping, as attempted previously by Law et al.
(2018a) with VBAgg and Lucas et al. (2020). For improved performance with
fewer observations of covariates in each county, it may be interesting to use
control variates to give estimators with lower variance for the required
integrals (Hammersley and Handscomb (1964)). Our model also makes the
assumption that data for different resolutions are not dependent on each
other, which may not be the case, and it would be useful to explore methods
that can account for dependence without data up-sampling or down-sampling. It
also may be of interest to explore alternative stochastic process priors for
the underlying function other than GPs. Lastly, while our modelling approach
is flexible for aggregation spatially, the dimensionality of the covariate
space at each spatial location may be high and so it may be useful to consider
learning low-dimensional projections within the model (Guhaniyogi and Dunson
(2016); Delbridge et al. (2020)).

##  6 Acknowledgements

The authors would like to thank Hitesh Kumar and Jean-François Ton for useful
feedback on a previous draft of the paper. HZ, AH were supported by the EPSRC
Centre for Doctoral Training in Modern Statistics and Statistical Machine
Learning (EP/S023151/1) and the Department of Mathematics of Imperial College
London. HZ was supported by Cervest Limited. AH was supported by the Bill and
Melinda Gates Foundation. DS was supported in part by the Alan Turing
Institute. SF was supported by the EPSRC (EP/V002910/1). We would also like to
thank Andy Thomas for his endless support with using the NVIDIA4 GPU Compute
Server. Many thanks to Gustau Camps-Valls and Anna Mateo-Sanchis for sharing
the data used for their work.

## References

  * (1)
  * Adsuara et al. (2019) Jose E Adsuara, Adrián Pérez-Suay, Jordi Muñoz-Marí, Anna Mateo-Sanchis, Maria Piles, and Gustau Camps-Valls. 2019.  Nonlinear distribution regression for remote sensing applications.  _IEEE Transactions on Geoscience and Remote Sensing_ 57, 12 (2019), 10025–10035.
  * Bhatt et al. (2017) Samir Bhatt, Ewan Cameron, Seth R Flaxman, Daniel J Weiss, David L Smith, and Peter W Gething. 2017\.  Improved prediction accuracy for disease risk mapping using Gaussian process stacked generalization.  _Journal of The Royal Society Interface_ 14, 134 (2017), 20170520\.
  * Bolton and Friedl (2013) Douglas K. Bolton and Mark A. Friedl. 2013.  Forecasting crop yield using remotely sensed vegetation indices and crop phenology metrics.  _Agricultural and Forest Meteorology_ 173 (May 2013), 74–84.  https://doi.org/10.1016/j.agrformet.2013.01.007
  * Burke et al. (2021) Marshall Burke, Anne Driscoll, David B. Lobell, and Stefano Ermon. 2021.  Using satellite imagery to understand and promote sustainable development.  _Science_ 371, 6535 (March 2021), eabe8628.  https://doi.org/10.1126/science.abe8628
  * De G. Matthews et al. (2017) Alexander G De G. Matthews, Mark Van Der Wilk, Tom Nickson, Keisuke Fujii, Alexis Boukouvalas, Pablo León-Villagrá, Zoubin Ghahramani, and James Hensman. 2017.  GPflow: A Gaussian process library using TensorFlow.  _The Journal of Machine Learning Research_ 18, 1 (2017), 1299–1304.  Publisher: JMLR. org.
  * Delbridge et al. (2020) Ian Delbridge, David Bindel, and Andrew Gordon Wilson. 2020\.  Randomly projected additive Gaussian processes for regression. In _International Conference on Machine Learning_. PMLR, 2453–2463.
  * Didan (2015) Kamel Didan. 2015\.  _MOD13Q1 MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid V006_.  NASA EOSDIS Land Processes DAAC.  https://doi.org/10.5067/MODIS/MOD13Q1.006
  * Entekhabi et al. (2010) Dara Entekhabi, Eni G Njoku, Peggy E O’Neill, Kent H Kellogg, Wade T Crow, Wendy N Edelstein, Jared K Entin, Shawn D Goodman, Thomas J Jackson, Joel Johnson, and others. 2010.  The soil moisture active passive (SMAP) mission.  _Proc. IEEE_ 98, 5 (2010), 704–716.  Publisher: IEEE.
  * Flaxman et al. (2019) Seth Flaxman, Michael Chirico, Pau Pereira, and Charles Loeffler. 2019.  Scalable high-resolution forecasting of sparse spatiotemporal events with kernel methods: A winning solution to the NIJ “Real-Time Crime Forecasting Challenge”.  _The Annals of Applied Statistics_ 13, 4 (2019), 2564–2585.  https://doi.org/10.1214/19-AOAS1284
  * Friedl and Sulla-Menashe (2015) M Friedl and D Sulla-Menashe. 2015.  MCD12C1 MODIS/Terra+ Aqua Land Cover Type Yearly L3 Global 0.05 Deg CMG V006 [Data set].  _NASA EOSDIS Land Processes DAAC_ (2015).
  * Garreau et al. (2018) Damien Garreau, Wittawat Jitkrittum, and Motonobu Kanagawa. 2018\.  Large sample analysis of the median heuristic.  _arXiv:1707.07269 [math, stat]_ (Oct. 2018).  http://arxiv.org/abs/1707.07269 arXiv: 1707.07269.
  * Gelfand (2010) Alan E. Gelfand (Ed.). 2010\.  _Handbook of spatial statistics_.  CRC Press, Boca Raton.  OCLC: ocn176924834.
  * Guhaniyogi and Dunson (2016) Rajarshi Guhaniyogi and David B Dunson. 2016.  Compressed Gaussian process for manifold regression.  _The Journal of Machine Learning Research_ 17, 1 (2016), 2472–2497.  Publisher: JMLR. org.
  * Hamelijnck et al. (2019) Oliver Hamelijnck, Theodoros Damoulas, Kangrui Wang, and Mark Girolami. 2019.  Multi-resolution Multi-task Gaussian Processes. In _Advances in Neural Information Processing Systems_ , H. Wallach, H. Larochelle, A. Beygelzimer, F. d\textquotesingle Alché-Buc, E. Fox, and R. Garnett (Eds.), Vol. 32. Curran Associates, Inc.  https://proceedings.neurips.cc/paper/2019/file/0118a063b4aae95277f0bc1752c75abf-Paper.pdf
  * Hammersley and Handscomb (1964) John Michael Hammersley and David Christopher Handscomb. 1964\.  _Monte Carlo Methods_.  Methuen.  Google-Books-ID: Kk4OAAAAQAAJ.
  * Harris et al. (2020) Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant. 2020.  Array programming with NumPy.  _Nature_ 585, 7825 (Sept. 2020), 357–362.  https://doi.org/10.1038/s41586-020-2649-2 Publisher: Springer Science and Business Media LLC.
  * Haussmann et al. (2017) Manuel Haussmann, Fred A. Hamprecht, and Melih Kandemir. 2017\.  Variational Bayesian Multiple Instance Learning with Gaussian Processes. In _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_. IEEE, Honolulu, HI, 810–819.  https://doi.org/10.1109/CVPR.2017.93
  * Hensman et al. (2013) James Hensman, Nicolò Fusi, and Neil D. Lawrence. 2013\.  Gaussian Processes for Big Data. In _Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence_ _(UAI’13)_. AUAI Press, Arlington, Virginia, USA, 282–290.  event-place: Bellevue, WA.
  * Hensman et al. (2015a) James Hensman, Alexander Matthews, and Zoubin Ghahramani. 2015a.  Scalable variational Gaussian process classification. In _Artificial Intelligence and Statistics_. PMLR, 351–360.
  * Hensman et al. (2015b) James Hensman, Alexander G de G Matthews, Maurizio Filippone, and Zoubin Ghahramani. 2015b.  MCMC for variationally sparse Gaussian processes.  _arXiv preprint arXiv:1506.04000_ (2015).
  * Ilse et al. (2018) Maximilian Ilse, Jakub Tomczak, and Max Welling. 2018\.  Attention-based deep multiple instance learning. In _International conference on machine learning_. PMLR, 2127–2136.
  * Kanagawa et al. (2018) Motonobu Kanagawa, Philipp Hennig, Dino Sejdinovic, and Bharath K. Sriperumbudur. 2018.  Gaussian Processes and Kernel Methods: A Review on Connections and Equivalences.  _arXiv:1807.02582 [cs, stat]_ (July 2018).  http://arxiv.org/abs/1807.02582 arXiv: 1807.02582.
  * Kim and De la Torre (2010) Minyoung Kim and Fernando De la Torre. 2010.  Gaussian processes multiple instance learning. In _ICML_.
  * Konings et al. (2017) Alexandra G. Konings, Maria Piles, Narendra Das, and Dara Entekhabi. 2017.  L-band vegetation optical depth and effective scattering albedo estimation from SMAP.  _Remote Sensing of Environment_ 198 (Sept. 2017), 460–470.  https://doi.org/10.1016/j.rse.2017.06.037
  * Kuwata and Shibasaki (2015) Kentaro Kuwata and Ryosuke Shibasaki. 2015.  Estimating crop yields with deep learning and remotely sensed data. In _2015 IEEE International Geoscience and Remote Sensing Symposium (IGARSS)_. IEEE, 858–861.
  * Law et al. (2018a) Ho Chung Law, Dino Sejdinovic, Ewan Cameron, Tim Lucas, Seth Flaxman, Katherine Battle, and Kenji Fukumizu. 2018a.  Variational learning on aggregate outputs with Gaussian processes. In _Advances in Neural Information Processing Systems_. 6081–6091.
  * Law et al. (2018b) Ho Chung Leon Law, Dougal Sutherland, Dino Sejdinovic, and Seth Flaxman. 2018b.  Bayesian approaches to distribution regression. In _International Conference on Artificial Intelligence and Statistics_. PMLR, 1167–1176.
  * Lemercier et al. (2021) Maud Lemercier, Cristopher Salvi, Theodoros Damoulas, Edwin Bonilla, and Terry Lyons. 2021\.  Distribution Regression for Sequential Data. In _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ _(Proceedings of Machine Learning Research, Vol. 130)_ , Arindam Banerjee and Kenji Fukumizu (Eds.). PMLR, 3754–3762.  http://proceedings.mlr.press/v130/lemercier21a.html
  * Lobell et al. (2015) David B. Lobell, David Thau, Christopher Seifert, Eric Engle, and Bertis Little. 2015.  A scalable satellite-based crop yield mapper.  _Remote Sensing of Environment_ 164 (July 2015), 324–333.  https://doi.org/10.1016/j.rse.2015.04.021
  * Lucas et al. (2020) Tim C.D. Lucas, Anita K. Nandi, Suzanne H. Keddie, Elisabeth G. Chestnutt, Rosalind E. Howes, Susan F. Rumisha, Rohan Arambepola, Amelia Bertozzi-Villa, Andre Python, Tasmin L. Symons, Justin J. Millar, Punam Amratia, Penelope Hancock, Katherine E. Battle, Ewan Cameron, Peter W. Gething, and Daniel J. Weiss. 2020.  Improving disaggregation models of malaria incidence by ensembling non-linear models of prevalence.  _Spatial and Spatio-temporal Epidemiology_ (July 2020), 100357.  https://doi.org/10.1016/j.sste.2020.100357
  * Maron and Lozano-Pérez (1998) Oded Maron and Tomás Lozano-Pérez. 1998.  A framework for multiple-instance learning. In _Advances in neural information processing systems_. 570–576.
  * Martinez-Ferrer et al. (2020) Laura Martinez-Ferrer, Maria Piles, and Gustau Camps-Valls. 2020\.  Crop Yield Estimation and Interpretability With Gaussian Processes.  _IEEE Geoscience and Remote Sensing Letters_ (2020), 1–5.  https://doi.org/10.1109/LGRS.2020.3016140
  * Mateo-Sanchis et al. (2019) Anna Mateo-Sanchis, Maria Piles, Jordi Muñoz-Marí, Jose E. Adsuara, Adrián Pérez-Suay, and Gustau Camps-Valls. 2019.  Synergistic integration of optical and microwave satellite data for crop yield estimation.  _Remote Sensing of Environment_ 234 (Dec. 2019), 111460\.  https://doi.org/10.1016/j.rse.2019.111460
  * Muandet et al. (2012) Krikamol Muandet, Kenji Fukumizu, Francesco Dinuzzo, and Bernhard Schölkopf. 2012.  Learning from Distributions via Support Measure Machines. In _NIPS_. 10–18.  http://papers.nips.cc/paper/4825-learning-from-distributions-via-support-measure-machines
  * Oglic and Gärtner (2017) Dino Oglic and Thomas Gärtner. 2017.  Nyström Method with Kernel K-means++ Samples as Landmarks. In _Proceedings of the 34th International Conference on Machine Learning_ _(Proceedings of Machine Learning Research, Vol. 70)_ , Doina Precup and Yee Whye Teh (Eds.). PMLR, 2652–2660.  http://proceedings.mlr.press/v70/oglic17a.html
  * Rahimi and Recht (2008) Ali Rahimi and Benjamin Recht. 2008.  Random features for large-scale kernel machines. In _Advances in neural information processing systems_. 1177–1184.
  * Rasmussen (2003) Carl Edward Rasmussen. 2003\.  Gaussian processes in machine learning. In _Summer School on Machine Learning_. Springer, 63–71.
  * Rue et al. (2016) H. Rue, A. Riebler, S. H. Sørbye, J. B. Illian, D. P. Simpson, and F. K. Lindgren. 2016\.  Bayesian computing with INLA: A review.  _Annual Reviews of Statistics and Its Applications_ 4 (2016), 395–421.
  * Sanchis et al. (2019) Anna Mateo Sanchis, Jose Adsuara, Maria Piles, Adrián Perez-Suay, Jordi Muñoz-Marí, and Gustau Camps-Valls. 2019.  Multisensor distribution regression for crop yield estimation.. In _Geophysical Research Abstracts_ , Vol. 21.
  * Shi et al. (2020) Jiaxin Shi, Michalis Titsias, and Andriy Mnih. 2020\.  Sparse orthogonal variational inference for gaussian processes. In _International Conference on Artificial Intelligence and Statistics_. PMLR, 1932–1942.
  * Szabó et al. (2016) Zoltán Szabó, Bharath K Sriperumbudur, Barnabás Póczos, and Arthur Gretton. 2016.  Learning theory for distribution regression.  _The Journal of Machine Learning Research_ 17, 1 (2016), 5272–5311.  Publisher: JMLR. org.
  * Tanaka et al. (2019) Yusuke Tanaka, Toshiyuki Tanaka, Tomoharu Iwata, Takeshi Kurashima, Maya Okawa, Yasunori Akagi, and Hiroyuki Toda. 2019.  Spatially aggregated Gaussian processes with multivariate areal outputs. In _Advances in Neural Information Processing Systems_. 3000–3010.
  * Thorns (2018) Daniel Thorns. 2018\.  _Distribution Regression for Crop Yield Prediction_.  Master’s thesis. Department of Statistics, Oxford University, United Kingdom.
  * Titsias (2009) Michalis Titsias. 2009\.  Variational learning of inducing variables in sparse Gaussian processes. In _Artificial Intelligence and Statistics_. 567–574.
  * Ton et al. (2018) Jean-Francois Ton, Seth Flaxman, Dino Sejdinovic, and Samir Bhatt. 2018.  Spatial mapping with Gaussian processes and nonstationary Fourier features.  _Spatial Statistics_ 28 (Dec. 2018), 59–78.  https://doi.org/10.1016/j.spasta.2018.02.002
  * Yang et al. (2019) Qi Yang, Liangsheng Shi, and Lin Lin. 2019.  Plot-scale rice grain yield estimation using UAV-based remotely sensed images via CNN with time-invariant deep features decomposition. In _IGARSS 2019 - 2019 IEEE International Geoscience and Remote Sensing Symposium_. IEEE, Yokohama, Japan, 7180–7183.  https://doi.org/10.1109/IGARSS.2019.8898061
  * You et al. (2017) Jiaxuan You, Xiaocheng Li, Melvin Low, David Lobell, and Stefano Ermon. 2017.  Deep gaussian process for crop yield prediction based on remote sensing data. In _Thirty-First AAAI Conference on Artificial Intelligence_.
  * Yousefi et al. (2019) Fariba Yousefi, Michael T Smith, and Mauricio Álvarez. 2019\.  Multi-task Learning for Aggregated Data using Gaussian Processes. In _Advances in Neural Information Processing Systems_ , H. Wallach, H. Larochelle, A. Beygelzimer, F. d\textquotesingle Alché-Buc, E. Fox, and R. Garnett (Eds.), Vol. 32. Curran Associates, Inc.  https://proceedings.neurips.cc/paper/2019/file/64517d8435994992e682b3e4aa0a0661-Paper.pdf
  * Zhang et al. (2019) Cheng Zhang, Judith Butepage, Hedvig Kjellstrom, and Stephan Mandt. 2019.  Advances in Variational Inference.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 41, 8 (Aug. 2019), 2008–2026.  https://doi.org/10.1109/TPAMI.2018.2889774
  * Zhang et al. (2020) Yivan Zhang, Nontawat Charoenphakdee, Zhenguo Wu, and Masashi Sugiyama. 2020.  Learning from Aggregate Observations. In _Advances in Neural Information Processing Systems_ , H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin (Eds.), Vol. 33. Curran Associates, Inc., 7993–8005.  https://proceedings.neurips.cc/paper/2020/file/5b0fa0e4c041548bb6289e15d865a696-Paper.pdf