Variational Learning on Aggregate Outputs
with Gaussian Processes
Ho Chung Leon Law∗
Dino Sejdinovic∗
University of Oxford
University of Oxford
Seth Flaxman‡
Imperial College London

Ewan Cameron†
University of Oxford

Katherine Battle†
University Of Oxford

Tim CD Lucas†
University of Oxford

Kenji Fukumizu§
Institute of Statistical Mathematics

Abstract
While a typical supervised learning framework assumes that the inputs and theoutputs are measured at the same levels of granularity, many applications, including
global mapping of disease, only have access to outputs at a much coarser level
than that of the inputs. Aggregation of outputs makes generalization to new
inputs much more difficult. We consider an approach to this problem based on
variational learning with a model of output aggregation and Gaussian processes,
where aggregation leads to intractability of the standard evidence lower bounds. We
propose new bounds and tractable approximations, leading to improved prediction
accuracy and scalability to large datasets, while explicitly taking uncertainty into
account. We develop a framework which extends to several types of likelihoods,
including the Poisson model for aggregated count data. We apply our framework
to a challenging and important problem, the fine-scale spatial modelling of malaria
incidence, with over 1 million observations.

1

Introduction

A typical supervised learning setup assumes existence of a set of input-output examples {(x` , y` )}`from which a functional relationship or a conditional probabilistic model of outputs given inputs can belearned. A prototypical use-case is the situation where obtaining outputs y? for new, previously unseen,
inputs x? is costly, i.e., labelling is expensive and requires human intervention, but measurements
of inputs are cheap and automated. Similarly, in many applications, due to a much greater costin acquiring labels, they are only available at a much coarser resolution than the level at which
the inputs are available and at which we wish to make predictions. This is the problem of weakly
supervised learning on aggregate outputs [14, 20], which has been studied in the literature in avariety of forms, with classification and regression notably being developed separately and without
any unified treatment which can allow more flexible observation models. In this contribution, we
consider a framework of observation models of aggregated outputs given bagged inputs, which residein exponential families. While we develop a more general treatment, the main focus in the paper ison the Poisson likelihood for count data, which is motivated by the applications in spatial statistics.
In particular, we consider the important problem of fine-scale mapping of diseases. High resolution
maps of infectious disease risk can offer a powerful tool for developing National Strategic Plans,∗

Department of Statistics, Oxford, UK. <ho.law@stats.ox.ac.uk, dino.sejdinovic@stats.ox.ac.uk>
Big Data Institute, Oxford, UK. <dr.ewan.cameron@gmail.com, timcdlucas@gmail.com, katherine.battle@bdi.ox.ac.uk>‡
Department of Mathematics and Data Science Institute, London, UK. <s.flaxman@imperial.ac.uk>§
Tachikawa, Japan. <fukumizu@ism.ac.jp>†

Preprint.

allowing accurate stratification of intervention types to areas of greatest impact [5]. In low resource
settings these maps must be constructed through probabilistic models linking the limited observational
data to a suite of spatial covariates (often from remote-sensing images) describing social, economic,
and environmental factors thought to influence exposure to the relevant infectious pathways. In
this paper, we apply our method to the incidence of clinical malaria cases. Point incidence data ofmalaria is typically available at a high temporal frequency (weekly or monthly), but lacks spatial
precision, being aggregated by administrative district or by health facility catchment. The challenge
for risk modelling is to produce fine-scale predictions from these coarse incidence data, leveraging the
remote-sensing covariates and appropriate regularity conditions to ensure a well-behaved problem.
Methodologically, the Poisson distribution is a popular choice for modelling count data. In the
mapping setting, the intensity of the Poisson distribution is modelled as a function of spatial andother covariates. We use Gaussian processes (GPs) as a flexible model for the intensity. GPs are awidely used approach in spatial modelling but also one of the pillars of Bayesian machine learning,enabling predictive models which explicitly quantify their uncertainty. Recently, we have seen many
advances in variational GP posterior approximations, allowing them to couple with more complex
observation likelihoods (e.g. binary or Poisson data [21, 17]) as well as a number of effective scalable
GP approaches [24, 30, 8, 9], extending the applicability of GPs to dataset sizes previously deemedprohibitive.
Contribution Our contributions can be summarised as follows. A general framework is developed for aggregated observation models using exponential families and Gaussian processes. This
is novel, as previous work on aggregation or bag models focuses on specific types of output models
such as binary classification. Tractable and scalable variational inference methods are proposed forseveral instances of the aggregated observation models, making use of novel lower bounds on themodel evidence. In experiments, it is demonstrated that the proposed methods can scale to dataset
sizes of more than 1 million observations. We thoroughly investigate an application of the developed
methodology to disease mapping from coarse measurements, where the observation model is Poisson,
giving encouraging results. Uncertainty quantification, which is explicit in our models, is essential
for this application.
Related Work The framework of learning from aggregate data was believed to have been
first introduced in [20], which considers the two regimes of classification and regression. However,while the task of classification of individuals from aggregate data (also known as learning fromlabel proportions) has been explored widely in the literature [23, 22, 13, 18, 35, 34, 14], there
has been little literature on the analogous regression regime in the machine learning community.
Perhaps the closest literature available is [13], who considers a general framework for learningfrom aggregate data, but also only considers the classification case for experiments. In thiswork, we will appropriately adjust the framework in [13] and take this to be our baseline. A
related problem arises in the spatial statistics community under the name of ‘down-scaling’,‘fine-scale modelling’ or ‘spatial disaggregation’ [11, 10], in the analysis of disease mapping,
agricultural data, and species distribution modelling, with a variety of proposed methodologies(cf. [33] and references therein), including kriging [6]. However, to the best of our knowledge,
approaches making use of recent advances in scalable variational inference for GPs are not considered.
Another closely related topic is multiple instance learning (MIL), concerned with classification with max-aggregation over labels in a bag, i.e. a bag is positively labeled if at least one individualis positive, and it is otherwise negatively labelled. While the task in MIL is typically to predictlabels of new unobserved bags, [7] demonstrates that individual labels of a GP classifier can also be
inferred in MIL setting with variational inference. Our work parallels that approach, considering
bag observation models in exponential families and deriving new approximation bounds for somecommon generalized linear models. In deriving these bounds, we have taken an approach similar to[17], who considers the problem of Gaussian process-modulated Poisson process estimation using
variational inference. However, our problem is made more complicated by the aggregation of labels.
Other related research topics include distribution regression and set regression, as in [28, 15, 16] and[36]. In these regression problems, while the input data for learning is the same as the current setup,the goal is to learn a function at the bag level, rather than the individual level, the application of thesemethods in our setting, naively treating single individuals as “distributions”, may lead to suboptimal

2

performance. An overview of some other approaches for classification using bags of instances isgiven in [4].

2Bag observation model: aggregation in mean parameters

Suppose we have a statistical model p(y|η) for output y ∈ Y, with parameter η given by a function of
input x ∈ X , i.e., η = η(x). Although one can formulate p(y|η) in an arbitrary fashion, practitioners
often only focus on interpretable simple models, hence we restrict our attention to p(y|η) arisingfrom exponential families. We assume that η is the mean parameter of the exponential family.
Assume that we have a fixed set of points xai ∈ X such that xa = {xa1 , . . . , xaNa } is a bag of points
with Na individuals, and we wish to estimate the regression value η(xai ) for each individual. However,instead of the typical setup where we have a paired sample {(x` , y` )}` of individuals and their outputs
to use as a training set, we observe only aggregate outputs y a for each of the bags. Hence, ourtraining data is of the form
1
n Nn
n
1
({x1i }N(1)i=1 , y ), . . . ({xi }i=1 , y ),aand the goal is to estimate parameters η(xi ) corresponding to individuals. To relate the aggregate y aaand the bag xa = (xai )N
i=1 , we use the following bag observation model:
y a |xa ∼ p(y|η a ),

ηa =

Na
X

wia η(xai ),

(2)

i=1

where wia is an optional fixed non-negative weight used to adjust the scales (see Section 3 foran example). Note that the aggregation in the bag observation model is on the mean parameters
for individuals, not necessarily on the individual responses yia . This implies that each individual
contributes to the mean bag response and that the observation model for bags belongs to the same
parametric form as the one for individuals. For tractable and scalable estimation, we will use
variational methods, as the aggregated observation model leads to intractable posteriors. We consider
the Poisson, normal, and exponential distributions, but devote a special focus to the Poisson model inthis paper, and refer readers to Appendix A for other cases and experimental results for the Normal
model in Appendix H.2.
It is also worth noting that we place no restrictions on the collection of the individuals, with thebagging process possibly dependent on covariates xai or any unseen factors. The bags can also beof different sizes, with potentially the same individuals appearing in multiple bags. After we obtain
our individual model η(x), we can use it for prediction of in-bag individuals, as well as out-of-bag
individuals.

3

Poisson bag model: Modelling aggregate counts

The Poisson distribution p(y|λ) = λy e−λ /(y!) is considered for count observations, and this paper
discusses the Poisson regression with intensity λ(xai ) multiplied by a ‘population’ pai , which is aconstant assumed toP
be known for each individual (or ‘sub-bag’) in the bag. The population for a bag
a is given by pa = i pai . An observed bag count y a is assumed to follow
y a |xa ∼ Poisson(pa λa ),

λa :=

Na a
X
p
i

pa
i=1

λ(xai ).

Note that, by introducing unobserved counts yia ∼ Poisson(yia |pai λ(xai )), the bag observation y a
PNa a
has the same distribution as i=1
yi since the Poisson distribution is closed under convolutions.
If a bag and its individuals correspond to an area and its partition in geostatistical applications, asin the malaria example in Section 4.2, the population in the above bag model can be regarded as
the population of an area or a sub-area. With this formulation, the goal is to estimate the basic
intensity function λ(x) from the aggregated observations (1). Assuming independence given {xa }a ,the negative log-likelihood (NLL) `0 across bags is
"N
!#
Na
n
n
a
X
X
X
X
c
n
a a
a a
a
a a c
a
a
a
a
a
− log[Πa=1 p(y |x )] =p λ −y log(p λ ) =pi λ(xi ) − y log
pi λ(xi )
,
a=1

a=1

i=1

i=1

(3)
3

c

where = denotes an equality up to additive constant. During training, this term will pass informationfrom the bag level observations {y a } to the individual basic intensity λ(xai ). It is noted that once wehave trained an appropriate model for λ(xai ), we will be able to make individual level predictions, andalso bag level predictions if desired. We will consider baselines with (3) using penalized likelihoods
inspired by manifold regularization in semi-supervised learning [2] – presented in Appendix B. In thenext section, we propose a model for λ based on GPs.
3.1

VBAgg-Poisson: Gaussian processes for aggregate counts

Suppose now we model f as a Gaussian process (GP), then we have:!
NaX
y a |xa ∼ Poisson
pai λai ,λai = Ψ(f (xai )),f ∼ GP (µ, k)

(4)

i=1

where µ and k are some appropriate mean function and covariance kernel k(x, y). (For implementation, we consider a constant mean function.) Since the intensity is always non-negative, in all models,
we will need to use a transformation λ(x) = Ψ(f (x)), where Ψ is a non-negative valued function.
We will consider cases Ψ(f ) = f 2 and Ψ(f ) = ef . A discussion of various choices of this linkfunction in the context of Poisson intensities modulated by GPs is given in [17]. Modelling f witha GP allows us to propagate uncertainty on the predictions to λai , which is especially important inthis weakly supervised problem setting, where we do not directly observe any individual output yia .
Since the total number of individuals in our target application of disease mapping is typically in the
millions (see Section 4.2), we will approximate the posterior over λai := λ(xai ) using variational
inference, with details found in Appendix E.
For scalability of the GP method, as in previous literature [7, 17], we use a set of inducing points{u` }m`=1 , which are given by the function evaluations of the Gaussian process f at landmark points
W = {w1 , . . . , wm }; i.e., u` = f (w` ). The distribution p(u|W ) is thus given byu ∼ N (µW , KW W ),µW = (µ(w` ))` ,KW W = (k(ws , wt ))s,t .

(5)

The joint likelihood is given by:
Na
n Y
Y

p(y, f, u|X, W, Θ) =Poisson(y a |pa λa )p(f |u)p(u|W ), with f |u ∼ GP (µ̃u , K̃),(6)

a=1 i=1−1µ̃(z) = µz + kzW KW
W (u − µW ),−1
K̃(z, z 0 ) = k(z, z 0 ) − kzW KW
W kW z 0(7)

T

where kzW = (k(z, w1 ), . . . , k(z, w` )) , with µW , µz denoting their respective evaluations of themean function µ and Θ being parameters of the mean and kernel functions of the GP. Proceeding
similarly to [17], which discusses (non-bag) Poisson regression with GP, we obtain a lower bound of
the marginal log-likelihood log p(y|Θ):Z Z
log p(y|Θ) = log
p(y, f, u|X, W, Θ)df duZ Znp(u|W ) o≥log p(y|f, Θ)p(f |u, Θ)q(u)df du (Jensen’s inequality)
q(u)
Na
Na
X
 X
o
XZ Z n=y a logpai Ψ(f (xai ) −pai Ψ(f (xai )) p(f |u)q(u)df du
a

−

i=1

X

i=1

a

log(y !) − KL(q(u)||p(u|W )) =: L(q, Θ),(8)awhere q(u) is a variational distribution to be optimized. The general solution to the maximization overq of the evidence lower bound L(q, Θ) above is given by the posterior of the inducing points p(u|y),which is intractable. We introduce a restriction to the class of q(u) to approximate the posteriorp(u|y). Suppose that the variational distribution q(u) is Gaussian, q(u) = N (ηu , Σu ). We thenneed to maximize the lower bound L(q, Θ) over the variational parameters ηu and Σu . The resulting
q(u) gives
R an approximation to the posterior p(u|y) which also leads to a Gaussian approximationq(f ) = p(f |u)q(u)du to the posterior p(f |y), which we finally then transform through Ψ to obtain4the desired approximate posterior on each λ(xia ) (which is either log-normal or non-central χ2depending on the form of Ψ). The approximate posterior on λ will then allow us to make predictions
for individuals while, crucially, taking into account the uncertainties in f (note that even the posterior
predictive mean of λ will depend on the predictive variance in f due to the nonlinearity Ψ). We also
want to emphasis the use of inducing variables is essential for scalability in our model: we cannot
directly obtain approximations to the posterior of λ(xai ) for all individuals, since this is often large inour problem setting (Section 4.2).
As the p(u|W ) and q(u) are both Gaussian, the last term (KL-divergence) of (8) can be computed
explicitly with exact form found in Appendix E.3. To consider the first two terms, let q a (v a ) be themarginal normal distribution of v a = (f (xa1 ), . . . , f (xaNa )), where f follows the variational posteriorq(f ). The distribution of v a is then N (ma , S a ), using (7) :

−1
−1
−1
−1
a
ma = µxa + Kxa W KW
W (ηu − µW ), S = Kxa ,xa − Kxa W KW W − KW W Σu KW W KW xa(9)
In the first term of (8), each summand is of the form
Z
Z
Na
Na
X

X
y a logpai Ψ (via ) q a (v a )dv a −paiΨ (via ) q a (v a )dv a ,
i=1

(10)

i=1

in which the second term is tractable for both of Ψ(f ) = f 2 and Ψ(f ) = ef . The integral of thefirst term, however with q a Gaussian is not tractable. To solve this, we take different approaches forΨ(f ) = f 2 and Ψ(f ) = ef ; for the former, approximation by Taylor expansion is applied, while forthe latter, further lower bound is taken.
First consider the case Ψ(f ) = f 2 , and rewrite the first term of (8) as:2

y a E log kV a k, where V a ∼ N (m̃a , S̃ a ),with P a = diag pa1 , . . . , paNa , m̃a = P a 1/2 ma and S̃ a = P a 1/2 S a P a 1/2 . By a Taylor series
222
approximation for E log kV a k (similar to [29]) around E kV a k = km̃a k + trS̃ a , we obtain
Z
log

Na
X


pai (via )2 q a (v a )dv a

i=1
a>

 2m
≈ log ma> P a ma + tr(S a P a ) −P a S a P a ma + tr (S a P a )2

(ma> P a ma + tr(S a P a ))2=: ζ a . (11)with details are in Appendix E.4. An alternative approach which we use for the case Ψ(f ) = ef is totake a further lower bound, which is applicable to a general class of Ψ (we provide further details forthe analogous approach for Ψ(v) = v 2 in Appendix E.2). We use the following Lemma (proof foundin Appendix E.1):
Lemma 1. Let v = [v1 , . . . , vN ]> be a random vector with probability density q(v) with marginaldensities qi (v), and let wi ≥ 0, i = 1, . . . , N . Then, for any non-negative valued function Ψ(v),
Z
Z
N
N
X

X

log
wi Ψ(vi ) q(v)dv ≥ log
wi eξi , where ξi := log Ψ(vi )qi (vi )dvi .
i=1

i=1

Hence we obtain that
Z
log

Na
X

Na
X

a
a
pai emi ,pai evi q a (v a )dv a ≥ log

i=1

(12)

i=1

Using the above two approximation schemes, our objective (up to constant terms) can be formulated
as: 1) Ψ(v) = v 2
Ls1 (Θ, ηu , Σu , W ) :=

n
X
a=1

a a

y ζ −

Na
n X
X


a
(mai )2 + Sii/2 − KL(q(u)||p(u|W )),a=1 i=1

5(13)

Figure 1: Left: Random samples on the Swiss roll manifold. Middle, Right: Individual Average
NLL on train set for varying number of training bags n and increasing Nmean , over 5 repetitions.
Constant prediction within bag gives a NLL of 2.22. bag-pixel model gives NLL above 2.4 for thevarying number of bags experiment.

2) Ψ(v) = ev
Le1 (Θ, ηu , Σu , W ) :=

n
X

a

y log

a=1

Na
X
i=1

e

ma
i



−

Na
n X
X

a

a

emi +Sii /2 − KL(q(u)||p(u|W )).

(14)

j=1 i=1

Given these objectives, we can now optimise these lower bounds with respect to variational parameters{ηu , Σu }, parameters Θ of the mean and kernel functions, using stochastic gradient descent (SGD)on bags. Additionally, we might also learn W , locations for the landmark points. In this form, we
can also see that the bound for Ψ(v) = ev has the added computational advantage of not requiring
the full computation of the matrix S a , but only its diagonals, while for Ψ(v) = v 2 computation of ζ ainvolves full S a , which may be problematic for extremely large bag sizes.

4

Experiments

We will now demonstrate various approaches: Variational Bayes with Gaussian Process (VBAgg), a
MAP estimator of Bayesian Poisson regression with explicit feature maps (Nyström) and a neuralnetwork (NN) – the latter two employing manifold regularisation with RBF kernel (unless statedaotherwise). For additional baselines, we consider a constant within bag model (constant), i.e. λˆai = ypa
and also consider creating ‘individual’ covariates by aggregation of the covariates within a bag (bagpixel). For details of all these approaches, see Appendix B. We also denote Ψ(v) = ev and v 2 as Exp
and Sq respectively.
We implement our models in TensorFlow5 and use SGD with Adam [12] to optimise their respective
objectives, and we split the dataset into 4 parts, namely train, early-stop, validation and test set. Herethe early-stop set is used for early stopping for the Nyström, NN and bag-pixel models, while the
VBAgg approach ignores this partition as it optimises the lower bound to the marginal likelihood.
The validation set is used for parameter tuning of any regularisation scaling, as well as learning rate,
layer size and multiple initialisations. Throughout, VBAgg and Nyström have access to the same setof landmarks for fair comparison. It is also important to highlight that we perform early stopping andtuning based on bag level performance on NLL only, as this is the only information available to us.
For the VBAgg model, there are two approaches to tuning, one approach is to choose parameters
based on NLL on the validation bag sets, another approach is to select all parameters based on the
training objective L1 , the lower bound to the marginal likelihood. We denote the latter approach
VBAgg-Obj and report its toy experimental results in Appendix H.1.1 for presentation purposes.
In general, the results are relatively insensitive to this choice, especially when Ψ(v) = v 2 . To
make predictions, we use the mean of our approximated posterior (provided by a log-normal andnon-central χ2 distribution for Exp and Sq). As an additional evaluation, we report mean square error(MSE) and bag performance results in Appendix H.
5Code will be available for use.

6

4.1

Poisson Model: Swiss Roll

We first demonstrate our method on the swiss roll dataset6 , illustrated in Figure 1 (left). To make thisan aggregate learning problem, we first construct n bags with sizes drawn from a negative binomial
distribution Na ∼ N B(Nmean , Nstd ), where Nmean
Pnand Nstd represents the respective mean andstandard deviation of Na . We then randomly select a=1 Na points from the swiss roll manifold tobe the locations, giving us a set of colored locations in R3 . Ordering these random locations by theirz-axis coordinate, we group them, filling up each bag in turn as we move along the z-axis. The aimof this is to simulate that in real life the partitioning of locations into bags is often not independentof covariates. Taking the colour of each location as the underlying rate λai at that location, we
PNa asimulate yia ∼ P oisson(λai ), and take our observed outputs to be y a = i=1
yi ∼ P oisson(λa ),
P
Na
a
a
a
where λ =i=1 λi . Our goal is then to predict the underlying individual rate parameter λi ,agiven only bag-level observations y . To make this problem even more challenging, we embedthe data manifold into R18 by rotating it with a random orthogonal matrix. For the choice of k for
VBAgg and Nyström, we use the RBF kernel, with the bandwidth parameter learnt. For landmark locations, we use the K-means++ algorithm, so that landmark points lie evenly across the data manifold.
Varying number of Bags: n To see the effect of increasing number of bags available for training,
we fix Nmean = 150 and Nstd = 50, and vary the number of bags n for the training set from100 to 350 with the same number of bags for early stopping and validation. Each experiment isrepeated for 5 runs, and results are shown in Figure 1 for individual NLL on the train set. Again weemphasise that the individual labels are not used in training. We see that all versions of VBAgg
outperform all other models, in terms of MSE and NLL, with statistical significance confirmed by asigned rank permutation test (see Appendix H.1.1). We also observe that the bag-pixel model haspoor performance, as a result of losing individual level covariate information in training by simply
aggregating them.
Varying number of individuals per bag: Nmean
To study the effect of increasing bag
sizes (with larger bag sizes, we expect "disaggregation" to be more difficult), we fix the numberof training bags to be 600 with early stopping and validation set to be 150 bags, while varying thenumber of individuals per bag through Nmean and Nstd in the negative binomial distribution. To
keep the relative scales between Nmean and Nstd the same, we take Nstd = Nmean /2. The results
are shown in Figure 1, focusing on the best performing methods in the previous experiment. Here, we
observe that VBAgg models again perform better than the Nyström and NN models with statisticalsignificance as reported in Appendix H.1.1, with performance stable as Nmean increases.
Discussion
To gain more insight into the VBAgg model, we look at the calibration of
our two different Bayesian models: VBAgg-Exp and VBAgg-Square. We compute their respective
posterior quantiles and observe the ratio of times the true λai lie in these quantiles. We present
these in Appendix H.1.1. The calibration plots reveal an interesting nature about using the twodifferent approximations for using ev versus v 2 for Ψ(v). While experiments showed that the twomodel perform similarly in terms of NLL, the calibration of the models is very different. While the
VBAgg-Square is well calibrated in general, the VBAgg-Exp suffers from poor calibration. This is
not surprising, as VBAgg-Exp uses an additional lower bound on model evidence. Thus, uncertainty
estimates given by VBAgg-Exp should be treated with care.
4.2

Malaria Incidence Prediction

We now demonstrate the proposed methodology on an important real life malaria prediction problem
for an endemic country from the Malaria Atlas Project database7 . In this problem, we would liketo predict the underlying malaria incidence rate in each 1km by 1km region (referred to as a pixel),
while having only observed aggregated incidences of malaria y a at much larger regional levels, which
are treated as bags of pixels. These bags are non-overlapping administrative units, with Na pixels per
bag ranging from 13 to 6,667, with a total of 1,044,683 pixels. In total, data is available for 957 bags8 .
6The swiss roll manifold function (for sampling) can be found on the Python scikit-learn package.
Due to confidentiality reasons, we do not report country or plot the full map of our results.
8
We consider 576 bags for train, 95 bags each for validation and early-stop, with 191 bags for testing, withdifferent splits across different trials, selecting them to ensure distributions of labels are similar across sets.
77

Figure 2: Triangle denotes approximate start and end of river location, crosses denotes non-train set
bags. Malaria incidence rate λai is per 1000 people. Left, Middle: log(λ̂ai ), with constant model(Left), and VBAgg-Obj-Sq (tuned on Ls1 ) (Middle). Right: Standard deviation of the posterior v in(9) with VBAgg-Obj-Sq.
Along with these pixels, we also have population estimates pai (per 1000 people) for pixel i in bag a,spatial coordinates given by sai , as well as covariates xai ∈ R18 , collected by remote sensing. Some
examples of covariates includes accessibility, distance to water, mean of land surface temperatureand stable night lights. It is clear that rather than expecting malaria incidence rate to be constantthroughout the entire bag (as in Figure 2), we expect pixel incidence rate to vary, depending on social,economic and environmental factors [32]. Our goal is therefore to build models that can predict
malaria incidence rates at a pixel level.
P
We assume a Poisson model on each individual pixel, i.e. y a ∼ P oisson( i pai λai ), where λai is
the underlying pixel incidence rate of malaria per 1000 people that we are interested in predicting.
We consider the VBAgg, Nyström and NN as prediction models and use a kernel given as a sum ofan ARD (automatic relevance determination) kernel on covariates and a Matérn kernel on spatial
locations for the VBAgg and Nyström methods, learning all kernel parameters (the kernel expression
is provided in Appendix G). We use the same kernel for manifold regularisation in the NN model.
This kernel choice incorporates spatial information, while allowing feature selection amongst othercovariates. For choice of landmarks, we ensure landmarks are placed evenly throughout space byusing one landmark point per training bag (selected by k-means++). This is so that the uncertainty
estimates we obtain are not too sensitive to the choice of landmarks. In this problem, no individuallevel labels are available, so we report Bag NLL and MSE (on observed incidences) on the test bagsin Appendix G over 10 different re-splits of the data. Although we can see that Nyström is the best
performing method, the improvement over VBAgg models is not statistically significant. On theother hand, both VBAgg and Nyström models statistically significantly outperform NN, which also
has some instability in its predictions, as discussed in Appendix G.1. However, a caution should be
exercised when using the measure of performance at the bag level as a surrogate for the measure of
performance at the individual level: in order to perform well at the bag level, one can simply utilise
spatial coordinates and ignore other covariates, as malaria intensity appears to smoothly vary between
the bags (Left of Figure 2). However, we do not expect this to be true at the individual level.
To further investigate this, we consider a particular region, and look at the predicted individual malaria
incidence rate, with results found in Figure 2 and in Appendix G.1 across 3 different data splits,
where the behaviours of each of these models can be observed. While Nyström and VBAgg methods
both provide good bag-level performance, Nyström and VBAgg-Exp can sometimes provide overlysmooth spatial patterns, which does not seem to be the case for the VBAgg-Sq method (recall that
VBAgg-Sq performed best in both prediction and calibration for the toy experiments). In particular,
VBAgg-Sq consistently predicts higher intensity along rivers (a known factor [31]; indicated bytriangles in Figure 2) using only coarse aggregated intensities, demonstrating that prediction of(unobserved) pixel-level intensities is possible using fine-scale environmental covariates, especiallyones known to be relevant such as covariates indicated by the Topographic Wetness Index, a measureof wetness, see Appendix G.2 for more details.
In summary, by optimising the lower bound to the marginal likelihood, the proposed variational
methods are able to learn useful relations between the covariates and pixel level intensities, while
avoiding the issue of overfitting to spatial coordinates. Furthermore, they also give uncertainty
estimates (Figure 2, right), which are essential for problems like these, where validation of predictions
is difficult, but they may guide policy and planning.
85Conclusion

Motivated by the vitally important problem of malaria, which is the direct cause of around 187
million clinical cases [3] and 631,000 deaths [5] each year in sub-Saharan Africa, we have proposed a
general framework of aggregated observation models using Gaussian processes, along with scalable
variational methods for inference in those models, making them applicable to large datasets. The
proposed method allows learning in situations where outputs of interest are available at a much coarser
level than that of the inputs, while explicitly quantifying uncertainty of predictions. The recent uptakeof digital health information systems offers a wealth of new data which is abstracted to the aggregate
or regional levels to preserve patient anonymity. The volume of this data, as well as the availability ofmuch more granular covariates provided by remote sensing and other geospatially tagged data sources,
allows to probabilistically disaggregate outputs of interest for finer risk stratification, e.g. assisting
public health agencies to plan the delivery of disease interventions. This task demands new highperformance machine learning methods and we see those that we have developed here as an important
step in this direction.

Acknowledgement
We thank Kaspar Martens for useful discussions, and Dougal Sutherland for providing the codebase in which this work was based on. HCLL is supported by the EPSRC and MRC through the
OxWaSP CDT programme (EP/L016710/1). HCLL and KF are supported by JSPS KAKENHI
26280009. EC and KB are supported by OPP1152978, TL by OPP1132415 and the MAP databaseby OPP1106023. DS is supported in part by the ERC (FP7/617071) and by The Alan Turing Institute(EP/N510129/1). The data were provided by the Malaria Atlas Project supported by the Bill andMelinda Gates Foundation.

9

References[1] LU Ancarani and G Gasaneo. Derivatives of any order of the confluent hypergeometric function f
1 1 (a, b, z) with respect to the parameter a or b. Journal of Mathematical Physics, 49(6):063508,
2008.
[2] Mikhail Belkin, Partha Niyogi, and Vikas Sindhwani. Manifold regularization: A geometricframework for learning from labeled and unlabeled examples. Journal of machine learning
research, 7(Nov):2399–2434, 2006.
[3] Samir Bhatt, DJ Weiss, E Cameron, D Bisanzio, B Mappin, U Dalrymple, KE Battle, CL Moyes,
A Henry, PA Eckhoff, et al. The effect of malaria control on plasmodium falciparum in africa
between 2000 and 2015. Nature, 526(7572):207, 2015.
[4] Veronika Cheplygina, David M.J. Tax, and Marco Loog. On classification with bags, groups
and sets. Pattern Recognition Letters, 59:11 – 17, 2015.
[5] Peter W Gething, Daniel C Casey, Daniel J Weiss, Donal Bisanzio, Samir Bhatt, Ewan Cameron,
Katherine E Battle, Ursula Dalrymple, Jennifer Rozier, Puja C Rao, et al. Mapping plasmodium
falciparum mortality in africa between 1990 and 2015. New England Journal of Medicine,375(25):2435–2445, 2016.
[6] Pierre Goovaerts. Combining areal and point data in geostatistical interpolation: Applications
to soil science and medical geography. Mathematical Geosciences, 42(5):535–554, Jul 2010.
[7] Manuel Haußmann, Fred A Hamprecht, and Melih Kandemir. Variational bayesian multiple
instance learning with gaussian processes. In Proceedings of the IEEE Conference on Computer
Vision and Pattern Recognition, pages 6570–6579, 2017.
[8] James Hensman, Nicolo Fusi, and Neil D Lawrence. Gaussian processes for big data. 2013.
[9] James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable Variational GaussianProcess Classification. In Guy Lebanon and S. V. N. Vishwanathan, editors, Proceedings ofthe Eighteenth International Conference on Artificial Intelligence and Statistics, volume 38
of Proceedings of Machine Learning Research, pages 351–360, San Diego, California, USA,
09–12 May 2015. PMLR.
[10] Richard Howitt and Arnaud Reynaud. Spatial disaggregation of agricultural production datausing maximum entropy. European Review of Agricultural Economics, 30(3):359–387, 2003.
[11] Petr Keil, Jonathan Belmaker, Adam M Wilson, Philip Unitt, and Walter Jetz. Downscalingof species distribution models: a hierarchical approach. Methods in Ecology and Evolution,
4(1):82–94, 2013.
[12] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint
arXiv:1412.6980, 2014.
[13] Dimitrios Kotzias, Misha Denil, Nando De Freitas, and Padhraic Smyth. From group to
individual labels using deep features. In Proceedings of the 21th ACM SIGKDD International
Conference on Knowledge Discovery and Data Mining, pages 597–606. ACM, 2015.
[14] H. Kueck and N. de Freitas. Learning about individuals from group statistics. In UAI, pages
332–339, 2005.
[15] H. C. L. Law, C. Yau, and D. Sejdinovic. Testing and learning on distributions with symmetricnoise invariance. In NIPS, 2017.
[16] Ho Chung Leon Law, Dougal Sutherland, Dino Sejdinovic, and Seth Flaxman. Bayesian
approaches to distribution regression. In International Conference on Artificial Intelligence andStatistics, pages 1167–1176, 2018.
[17] Chris Lloyd, Tom Gunter, Michael Osborne, and Stephen Roberts. Variational inference
for gaussian process modulated poisson processes. In International Conference on Machine
Learning, pages 1814–1822, 2015.
[18] Vitalik Melnikov and Eyke Hüllermeier. Learning to aggregate using uninorms. In Joint
European Conference on Machine Learning and Knowledge Discovery in Databases, pages
756–771. Springer, 2016.
[19] Krikamol Muandet, Kenji Fukumizu, Bharath Sriperumbudur, and Bernhard Schölkopf. Kernelmean embedding of distributions: A review and beyonds. arXiv preprint arXiv:1605.09522,
2016.
10

[20] David R Musicant, Janara M Christensen, and Jamie F Olson. Supervised learning by training on
aggregate outputs. In Data Mining, 2007. ICDM 2007. Seventh IEEE International Conferenceon, pages 252–261. IEEE, 2007.
[21] H. Nickisch and CE. Rasmussen. Approximations for binary gaussian process classification.
Journal of Machine Learning Research, 9:2035–2078, October 2008.
[22] Giorgio Patrini, Richard Nock, Tiberio Caetano, and Paul Rivera. (Almost) no label no cry. In
NIPS. 2014.
[23] Novi Quadrianto, Alex J Smola, Tiberio S Caetano, and Quoc V Le. Estimating labels fromlabel proportions. JMLR, 10:2349–2374, 2009.
[24] Joaquin Quiñonero Candela and Carl Edward Rasmussen. A unifying view of sparse approximate gaussian process regression. J. Mach. Learn. Res., 6:1939–1959, December 2005.
[25] Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In NIPS,
pages 1177–1184, 2007.
[26] Carl Edward Rasmussen and Christopher KI Williams. Gaussian processes for machine learning,
2006.
[27] Alex J Smola and Peter L Bartlett. Sparse greedy gaussian process regression. In Advances inneural information processing systems, pages 619–625, 2001.
[28] Zoltán Szabó, Bharath K Sriperumbudur, Barnabás Póczos, and Arthur Gretton. Learning theory
for distribution regression. The Journal of Machine Learning Research, 17(1):5272–5311, 2016.
[29] Yee W Teh, David Newman, and Max Welling. A collapsed variational bayesian inference
algorithm for latent dirichlet allocation. In Advances in neural information processing systems,
pages 1353–1360, 2007.
[30] Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In
David van Dyk and Max Welling, editors, Proceedings of the Twelth International Conferenceon Artificial Intelligence and Statistics, volume 5 of Proceedings of Machine Learning Research,
pages 567–574, Hilton Clearwater Beach Resort, Clearwater Beach, Florida USA, 16–18 Apr
2009. PMLR.
[31] DA Warrel, T Cox, J Firth, and Jr E Benz. Oxford textbook of medicine, 2017.
[32] Daniel J Weiss, Bonnie Mappin, Ursula Dalrymple, Samir Bhatt, Ewan Cameron, Simon I
Hay, and Peter W Gething. Re-examining environmental correlates of plasmodium falciparum
malaria endemicity: a data-intensive variable selection approach. Malaria journal, 14(1):68,
2015.
[33] António Xavier, Maria de Belém Costa Freitas, Maria do Socorro Rosário, and Rui Fragoso.
Disaggregating statistical data at the field level: An entropy approach. Spatial Statistics, 23:91 –
108, 2018.
[34] Felix X Yu, Krzysztof Choromanski, Sanjiv Kumar, Tony Jebara, and Shih-Fu Chang. On
learning from label proportions. arXiv preprint arXiv:1402.5902, 2014.
[35] Felix X Yu, Dong Liu, Sanjiv Kumar, Tony Jebara, and Shih-Fu Chang. propto svm for learningwith label proportions. arXiv preprint arXiv:1306.0886, 2013.
[36] Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Ruslan Salakhutdinov,and Alexander Smola. Deep sets. In NIPS, 2017.

11

A

Aggregated Exponential Family Models

Consider an observation model of the formp(y|θ) = expyθ − c(θ)τ

h(y, τ ),(15)

where response y is one-dimensional, θ is a natural parameter corresponding to the statistic y, τ
is a dispersion parameter, and h(y, τ ) is base measure. For simplicity, we will assume that naturalparameters corresponding to the other parts of the sufficient statistic are fixed and folded into the basemeasure. Let η be the corresponding mean parameter, i.e.
Z
η = Eθ y = yp(y|θ)dy
and θ = F (η) be the link function mapping from mean to the natural parameters and G(θ) its inverse.
We wish to model the mean parameter η = η(x) using a Gaussian process on a domain X together
with a function Ψ which transforms the GP value to the natural parameter space, i.e.
f ∼ GP(µ, k).

η(x) = Ψ(f (x)),(16)

For example, the mean parameter for some models is restricted to the positive part of the real line,
while the GP values cover the whole real line. We will consider the following examples:• Normal (with fixed variance). F = G = idenity and Ψ can be identity, too, as there are no
restrictions on the mean parameter space.
• Poisson. F (η) = log η, G(θ) = eθ . Ψ should take a positive value, so we consider
Ψ(v) = ev or Ψ(v) = v 2 .
• Exponential. p(y|η) = exp(−y/η)/η and θ = −η, F (η) = −1/η, G(θ) = −1/θ. Ψshould take a positive value, so we consider Ψ(v) = ev or Ψ(v) = v 2
Note that the link function F is concave for all the examples above.
A.1

Bag model

We will consider the aggregation in the mean parameter space. Namely, let y 1 , . . . , y n be n independent aggregate responses for each of the n bags of covariates xa = {xa1 , . . . , xaNa }, a = 1, . . . , n.
We assume the following aggregation model:y a ∼ p(y|ηa ),

ηa =

Na
X
i=1

wia ηia =NaX

wia Ψ(f (xai )),a = 1, . . . , n.

(17)

i=1

where wia are fixed weights to adjust the scales among the individuals and the bag (e.g., adjusting for
population size).
We also can model individual (unobserved) variables yia (i = 1, . . . , Na ), which follow:
yia ∼ p(y|ηia ),ηia = Ψ(f (xai )),

i = 1, . . . , Na , a = 1, . . . , n.

(18)Note that we consider aggregation in mean parameters of responses, not in the responses themselves.
If we consider a case where underlying individual responses yia aggregate to y a as a weighted
sum, the form of the bag likelihood and individual likelihood would be different unless we restrict
attention to distribution families which are closed under both scaling and convolution. However,when aggregation occurs in the mean parameter space, the form of the bag likelihood and individual
likelihood is always the same. This corresponds to the following measurement process:• Each individual has a mean parameter ηia - if it were possible to sample a response for thatparticular individual, we would obtain a sample yia ∼ p(·|ηia )• However, we cannot sample the individual and we can only observe a bag response. But inthat case, only a single bag response is taken and depends on all individuals simultaneously.
Each individual contributes in terms of an increase in a mean bag response, but this measurement process is different from the two-stage procedure by which we aggregate individual
responses.
12

A.2

Marginal likelihood and ELBO

Let Y = (y 1 , . . . , y n ) (bag observations). With the inducing points u = f (W ), the marginal
likelihood is
Z Z Ynp(Y ) =p(y a |η a )p(f |u)p(u)dudf.
(19)
a=1

The evidence lower bound can be derived as
Z Z nYnp(u) op(f |u)q(u)dudf
log p(Y ) = log
p(y a |η a )
q(u)
a=1
Z Z
n
nY
p(u) o≥log
p(y a |η a )p(f |u)q(u)dudf
q(u)
a=1
Z X
Z  X
n


X
ya
=
F
wia Ψ(f (xai )) q(f )df − c F
wia Ψ(f (xai )) q(f )df
τ
a=1
i
i
Z
q(u)
du,
(20)
− q(u) log
p(u)R
where q(f ) = p(f |u)q(u)du.
By setting the variational distribution q(u) as Gaussian, the third term is tractable. The first andsecond terms are however tractable only in limited cases. The cases we develop are the Poisson bag
model, described in the main text, as well as the normal bag model and the exponential bag model,
described below.
A.3

Normal bag model

F is identity and c(θ) = θ2 /2, which makes both the first and the second terms tractable with the
choice of Ψ(v) = v. Moreover, the viewpoints of aggregating in the mean parameters and in theindividual responses are equivalent for this model and we can also allow different variance parameters
for different bags (and individuals).
aaaConsider a bag a of items {xai }N
i=1 . Each item xi is assumed to have a weight wi . At the individuallevel, we model the (unobserved) responses yia as2yia |xai ∼ N wia µai , (wia ) τia(21)

where µai = µ(xai ), thus µai is a mean parameter per unit weight corresponding to the item xai and
it is assumed to be a function of both xai . Similarly, τia is a variance parameter per unit weight. At
the bag level, we consider the following model for the observed aggregate response y a , assuming
conditional independence of individual responses given covariates xa = {xa1 , . . . , xaNa }:
ya =

Na
X

yia , i.e. y a |xa ∼ N (wa µa , (wa )2 τ a ),

µa =

i=1

Na
X
wa
i=1

i a a
µ ,τ =
wa i

PNa

a 2 a
i=1 (wi ) τi
a
(w )2

(22)

where µa and τ a are the mean and variance parameters per unit weight of the whole bag a andPNa awa = i=1
wi is the total weight of bag a. Although we can take τia to also be a function of
the covariates, here for simplicity, we take τia = τ a to be constant per bag (note the abuse ofnotation). We can now compute the negative log-likelihood (NLL) across bags (assuming conditional
independence given the xa ):!  a PNa a a 2 

N
n
a

y − i=1 wi µi X
1X
`0 = − log [Πna=1 p(y a |xa )] =log 2πτ a(wia )2 +(23)
PNa a 2 a2(wi ) τ


a=1

i=1

i=1

where µai = f (xai ) is the function we are interested in, and τ a are the variance parameters to belearnt.
13

We can now consider the lower bound to the marginal likelihood as below (assuming wia = 1 here tosimplify notation, while the analogous expression with non-uniform weights is straightforward):
Z Z
log p(y|Θ) = log
p(y, f, u|X, W, Θ)df du
!!
PNa
Z Z Y
n
f (xai ))2(y a − i=11p(u|W )
√
= log
exp −
p(f |u)q(u)df du
a
a
2N
τ
q(u)
2πN
τ
a
a
a=1
!
)
( n
PNa
Z Z
Y
f (xai ))2 p(u|W )(y a − i=1
1√p(f |u)q(u)df du
≥
log
exp −
2Na τ aq(u)2πNa τ a
a=1

2 
P
PNa
Na


a
a
a 2
a
Z
Z


(y ) − 2yi=1 f (xi )i=1 f (xi ) +1X=−
p(f |u)q(u)df du
2 aNa τ a


Z
1X
q(u)
−
log(2πNa τ a ) − q(u) log
du.
(24)
2 a
p(u|W )
R
Using again a Gaussian distribution for q(u), we have q(f ) = p(f |u)q(u)du, which is a normal
distribution and let q a (f a ) be its marginal normal distribution of f a = (f (xa1 ), . . . , f (xaNa )) withmean and covariance given by ma and S a as before in (9).
Then all expectations with respect to q(f ) are tractable and the ELBO is simplyn1X
L(q, θ) = −
2 a=1

(

 )
(y a )2 − 2y a 1> ma + 1> S a + ma (ma )> 1
1X−
log(2πNa τ a )
a
Na τ
2 a

− KL(q(u)||p(u|W )).
A.4

(25)

Exponential bag model

In this case, P
we have F (η) = −1/η. We can apply the similar argument as in Lemma 1. For any
αi > 0 with i αi = 1, by the concavity of F ,
!
!
Z
Z
X
X
F
wi Ψ(vi ) q(vi )dvi = F
αi wi /αi Ψ(vi ) q(vi )dvi
i

i

≥

Z X

αi F (wi /αi Ψ(vi )) q(vi )dvi

i

=

X

Z
αi

F (wi /αi Ψ(vi )) q(vi )dvi .

i

For F (η) = −1/η, the last line is equal toX α2 Ziiwi

1q(vi )dvi .
Ψ(vi )WhenRusing a normal q, this is tractable for several choices of Ψ including ev and v 2 . If we let1
ξi := Ψ(vq(vi )dvi , and maximize
i)
X
ξi
αi2
wi
i
P
under the constraint i αi = 1, we obtain(wi /ξi )αi = P.
` (wi /ξi )14

Finally, we have a lower bound!
P
ZX(wi /ξi )iFw Ψ(vi ) q(vi )dvi ≥ − P i
2i (wi /ξi )
i
where

(26)

Z

1
q(vi )dvi .
Ψ(vi )
which is tractable for a Gaussian variational family. Also with an explicit form of Ψ, it is easy to takethe derivatives of the resulting lower bound with respect to the variational parameters in q(v).
ξi =BAlternative approaches

Constant For the Poisson model, we can take λai = λac , a constant rate across the bag, then:ya
λ̂ac = aphence the individual level predictive distribution is the form yia ∼ P oisson(λ̂ac ), and for unseen bag
Pn
1
a
r
r bag
r, λ̂bagc = Pn pa
a=1 y , with predictive distribution given by y ∼ P oisson(p λ̂c ).
a=1

bag-pixel: Bag as Individual Another baseline is to train a model from the weighted average of
PNa pai a
PNa wia aathe covariates, given by xa = i=1
i=1 wa xi in the normal
pa xi in the Poisson case, and x =case. The purpose of this baseline is to demonstrate that modelling at the individual level is importantduring training. Since we now have labels and covariates at the bag level, we can consider thefollowing model:
y a |xa ∼ P oisson(pa λ(xa ))aawith λ(x ) = Ψ(f (x )) for the Poisson model. For the normal model, we have:y a |xa ∼ N ormal(wa µ(xa ), (wa )2 τ )aawhere µ(x ) = f (x ) and τ is a parameter to be learnt (assuming constant across bags). Now weobserve that these models are identical to the individual model, except for a difference in indexing.
Hence, after learning the function f at the bag level, we can transfer the model to the individuallevel. Essentially here we have created fake individual level instances by aggregation of individual
covariates inside a bag.
Nyström: Bayesian MAP for Poisson regression on explicit feature maps Instead of the posterior based on the model (6), we can also consider an explicit feature map in order to directly
construct a MAP estimator. While this method does not provide posterior uncertainty over λai , it does
provide an interesting connection to the settings we have considered and also manifold-regularizedneural networks, as discussed below. Let Kzz be the covariance function defined on covariates−1{z1 , . . . zn }, and consider its low rank approximation Kzz ≈ kzW KW
W kWz with landmark pointsmT
W = {w` }`=1 and kzW = (k(z, w1 ), . . . , k(z, w` )) . By using landmark points W , we haveavoided computation of the full kernel matrix, reducing computational complexity. Under this setup,− 12
we have that Kzz ≈ Φz Φ>
z , with Φz = kzW KW W being the explicit (Nyström) feature map. Using
this explicit feature map Φ, we have the following model:
fia = φai β,β ∼ N (0, γ 2 I)
!
Na
X
a a
a
a
y |x ∼ Poissonpi λ(xi ) ,λ(xai ) = Ψ(fia ),i=1

where γ is a prior parameter and φai is the corresponding ith row of Φxa . We can then consider a
MAP estimator of the model coefficients β:β̂ = argmaxβ log[Πna=1 p(y a |β, xa )] + log p(β).
(27)
This essentially recovers the same model as in (3) with the standard l2 loss regularising the complexityof the function. This model can be thought of in several different ways, for example as a weight spaceview of the GP ([26] for an overview), or as a MAP of the Subset of Regressors (SoR) approximation[27] of the GP when σ = 1. Additional we may include manifold regulariser as part of the prior, seediscussion below about neural network.
15

NN: Manifold-regularized neural networks The next approach we consider is a parametric model
for f as in [13], and search the best parameter to minimize negative log-likelihood `0 across bags. This
paper considers a neural network with parameters θ for the model f , and uses the back-propagation tolearn θ and hence individual level model f . However, since we only have aggregated observations at
the bag level, but lots of individual covariate information, it is useful to incorporate this informationalso, by enforcing smoothness on the data manifold given by the unlabelled data. To do this, following[13] and [22], we pursue a semisupervised view of the problem and include an additional manifold
2
regularisation term [2] (rescaling with Ntotal
during implementation):`1 =Ntotal N
total
X
X

(fu − fw )2 kL (xu , xw ) = f > L f(28)

w=1 u=1

where we have suppressed the bag index, Ntotal represents the total number of individuals, kL (·, ·) issome user-specified kernel9 , f = [f1 , . . . , fNtotal ]> , L is the Laplacian defined as L = diag(KL 1> ) −KL , where 1 is just [1, . . . , 1] and KL is a kernel matrix. Although this term involves calculation of akernel matrix across individuals, in practice we consider stochastic gradient descent (SGD) and also
random Fourier features [25] or Nyström approximation (see Appendix C), with scale parameter λ1
to control the strength of the regularisation. Similarly, one can also consider manifold regularisation
at the bag level, if bag-level covariates/embeddings are available, for further details, see Appendix D.
In fact, the same regularisation can be applied to the MAP estimation with the explicit feature maps.
This is equivalent to having a prior β ∼ N (0, σ 2 I + (λ1 Φ> LΦ)−1 ) that is data dependent andincorporates the structure of the manifold 10 .
For implementation, we consider a one hidden layer neural network, with also an output layer, for afair comparison to the Nyström approach. For activation function, we consider the Rectified Linear
Unit (ReLU).
MAP estimation of GP We introduce p(f, u) = p(f |u)p(u|W ) and consider the posterior givenby p(u|f, y, w, θ), where here the conditional distribution f |u is given by:
f |u ∼ GP (µ̃u , K̃),

(29)

−1
−1
µ̃(z) = µz + kzW KW
K̃(z, z 0 ) = k(z, z 0 ) − kzW KW
W (u − µW ),W kWz0T
where kzW = (k(z, W1 ), . . . , k(z, W` )) . Using Bayes rule, we obtain:

log[p(u|f, y, w)]= log[p(y|f, u)p(f, u|X, W )]= log[p(y|f )p(f |u, X)p(u|W )]
n
n
n
X
X
X
=
y a log(pa λa ) +pa λa −
log(y a !) + log(p(f |u, X)) + log(p(u|W ))
a=1

a=1

a=1

where p(f |u, X) ∼ N (µ̃u , K̃) given by above, and p(u|W ) ∼ N (µW , ΣW W ), i.e.
1
log p(f |u, X)+log p(u|W ) = − (log(|K̃||ΣW W |)+(f −µ̃u )> K̃ −1 (f −µ̃u )+(u−µW )> Σ−1
W W (u−µW )2(30)Here, we can not perform SGD, as the latter terms does not decompose into a sum over the data.
More importantly, here we require the computation of K̃, which contains the kernel matrix K, evenafter the use of landmarks. This direct approach is not feasible for large number of individuals, which
is true in our target application, and hence we do not pursue this method, and consider Nyström andNN as baselines.

C

Random Fourier Features on Laplacian

Here we discuss using random Fourier features [25] to reduce computational cost in calculation of
the Laplacian defined as L = diag(K1> ) − K, where 1 is just [1, . . . , 1] and K. Suppose the kernel
9

In practice, this does not have to be a positive semi-definite kernel, it can be derived from any notion ofsimilarity between observations, including k-nearest neighbours.
10
In order to guarantee positive definiteness of Laplacian, one can add I, where  > 0.

16

is stationary i.e. kw (x − y) = k(x, y) (some examples include the gaussian and matern kernel),then using random Fourier features, we obtain K ≈ ΦΦ> , where Φ ∈ RbN ×m , bN denotes the total
number of individuals in the batch and m denotes the number of frequencies. Now we have:f > Lf ≈ f > diag(ΦΦ> 1> )f − f > ΦΦ> f = f > diag(ΦΦ> 1> )f − ||Φ> f ||22(31)
In both terms, we can avoid computing the kernel matrix, by carefully selecting the order of computation. Note another option is to consider Nyström approximation with landmark points {z1 , . . . zm },−1then K ≈ Knm Kmm
Kmn , where Kmm denotes the kernel matrix on landmark points, while Knm−12
is the kernel matrix between landmark and data. Then Φ = Knm Kmm.

D

Bag Manifold regularisation

Suppose we have bag covariates sa (note these are for the entire bag), and also some summaryPNah(xai ), with some userstatistics of a bag, e.g. mean embeddings [19] given by H a = N1a i=1
defined h. Then similarly to individual level manifold regularisation, we can consider manifold
regularisation at the bag level (assuming a seperable kernel for simplicity), i.e.
n XnX(F l − F m )2 ks (sl , sm )kh (H l , H m ) = F> Lbag F
(32)
`2 =
l=1 m=1
a

= N1l

PNa

a
a
a
where Fi=1 fi , ks is a kernel on bag covariates s , kµ is a kernel on H , Lbag is the bag
level Laplacian with the corresponding kernel, and F = [F 1 , . . . , F n ]> . Combining all these terms,
we have the following loss function to minimise:1
λ1
λ2` = `0 + 2 `1 + 2 `2
(33)
b
bN
bN
where b is the mini-batch size in SGD, BN is the total number of individuals in each mini-batch, λ1
and λ2 are parameters controlling the strength of the respective regularisation.

E

Additional details for Poisson variational derivation

E.1

Log-sum lemma

Lemma 2. Let v = [v1 , . . . , vN ]> be a random vector with probability density q(v), and let wi ≥ 0,
i = 1, . . . , N . Then, for any non-negative valued function Ψ(v),
Z
N
N
X

X

log
wi Ψ(vi ) q(v)dv ≥ log
wi eξi ,
i=1

i=1

where

Z
ξi :=

log Ψ(vi )qi (vi )dvi .

PN
Proof. Let α1 , . . . , αN be non-negative numbers with i=1 αi = 1. It follows from Jensen’s
inequality that
Z
N
X

log
wi Ψ(vi ) q(v)dv =
i=1

Z
log

N
X


i
αi w
Ψ(v
)
q(v)dv ≥
i
αi

i=1
N
X
i=1

Z
αi


wi
log Ψ(vi ) q(vi )dvi + log
=
αi




N
X

αi ξi +

i=1

N
X
i=1

17

αi log

wi
.
αi

(34)

By Lagrange multiplier method, maximizing the last line with respect to α gives
w i eξ i
αi = PN.
ξj
j=1 wj e
Plugging this to (34) completes the proof.
E.2

A lower bound of marginal likelihood for Ψ(f ) = ef and Ψ(f ) = f 2

Using Lemma 2, we obtain that
Z
log

N
X

pai Ψ(via ) q(v a )dv a ≥ log


i=1

N
X


pai Ψ(ξia ) ,

(35)

i=1

where
ξia =

Z

log Ψ(via )qia (via )dvia .

The above lower bound is tractable for the popular functions Ψ(v) = v 2 and Ψ(v) = ev under thenormal variational distributions q a (v a ) ∼ N (ma , S a ) . In particular,ZvaΨ(v) = e : ξi =via qia (via )dvia = mai ,


 a
Z
ma
Sii
Ψ(v) = v 2 : ξia =
log(via )2 qia (via )dvia = −G − ia + log
− γ,
2Sii
2
where γ is the Euler constant andG(t) = 2t

∞
X

j!
tj
(2)
(3/2)
j
j
j=0

is the partial derivative of the confluent hypergeometric function [17, 1]. However, in this work wefocus on the Taylor series approximation for Ψ(v) = v 2 , as implementation of the above bound uses
a large look-up table and involves linear interpolation. Furthermore, it is suggested in experimentsthat the secondary lower bound proposed above in Lemma 2 can lead to poor calibration, for moredetails, refer to Section 4.
E.3

KL Term

Since q(u) and p(u|W ) are both normal distribution, the KL divergence is tractable:
KL(q(u)||p(u|W )) =

E.4

o
1n
|KW W |−1
−1
T r[KW−m+(µW −ηu )T KW(36)W Σu ]+log
W (µW −ηu )2|Σu |Taylor series approximation in the variational method

We consider the integral
Z
log

N
X


pai (via )2 q a (v a )dv ai=1
2

where q a is N (ma , S a ). We note that this can be written as E log kV a k , where V a ∼ N (m̃a , S̃ a ),
2with P a = diag pa1 , . . . , paNa , m̃a = P a 1/2 ma and S̃ a = P a 1/2 S a P a 1/2 . Note that kV a kfollows a non-central chi-squared distribution. We now resort to a Taylor series approximation for
18

2

2

2

E log kV a k (similar to [29]) around E kV a k = km̃a k + trS̃ a , resulting in




2
2
E log kV a k= log E kV a k



2


a 2
a 2
2
2


kV k − E kV k
3
 kV a k − E kV a k

2
2
−
+E 
+O
kV a k − E kV a k


2
2
a
2
E kV ka
2 E kV k  
2a> a a
2
m̃
S̃
m̃
+
tr
S̃ a


2
.
≈ log km̃a k + trS̃ a −2
2
km̃a k + trS̃ a
2

As commented in [29], approximation is very accurate when E kV a k is large, but the caveat is that22
the Taylor series converges only for kV k ∈ (0, 2E kV k ) so this approach effectively ignores thetail of the non-central chi-squared.

F

Code

All of our models were implemented in TensorFlow, and code will be published and available for use.

G

Additional Malaria Experimental Results

Here we provide additional experimental results for the malaria dataset. In table 1, we provide
results for bag level performance for NLL and MSE with 10 different test sets (after retrial of theexperiments, splitting the data across train, early-stop, validation and testing). Statistical significancewas not establish for the best performing Nyström method versus the VBAgg methods, this is shownin Table 2. We further provide additional prediction/uncertainty patches for 3 different splits tohighlight the general behaviour of the trained models, with further explanation and details below.
It is also noted in all cases λai is the incidence rate per 1000 people. For VBAgg and Nyström, we
use an additive kernel, between an ARD kernel and a Matern kernel:
!
!
!
√
√
18
1X 1
3||sx − sy ||23||sx − sy ||22(xk − yk ) +γ2 1 +exp −k((x, sx ), (y, sy )) = γ1 exp −
2
`k
ρ
ρ
k=1
(37)
where x, y are covariates, and sx , sy are their respective spatial location. Here, we learn any scaleparameters and weights during training. For the NN, we also use this kernel as part of manifold
regularisation, however we use an RBF kernel instead of an ARD kernel, due to parameter tuning
reasons (we can no longer learn these scales).
Pn
bag
1
a
For constant model, bag rate predictions are computed by, pa λ̂bag
c ,where λ̂c = Pn pa
a=1 y .
a=1
This essentially takes into account of population.
Table 1: Results for the Poisson Model on the malaria dataset with 10 different re-splits of train,
early-stopping, validation and test. Approximately, 191 bags are used for test set. Bag performance is
PNa a ameasured on a test set, with MSE computed between log(y a ) and log( i=1pi λ̂i ). Brackets include
standard deviation.
Bag NLL
Bag MSE (Log)Constant
173.1 (31.2)
4.08 (0.13)
Nyström-Exp
88.1 (25.1)
1.31 (0.15)
VBAgg-Sq-Obj
94.1 (34.0)
1.21 (0.05)
VBAgg-Exp-Obj 97.2 (39.6)
1.04 (0.11)
VBAgg-Sq
97.6 (39.0)
1.38 (0.18)
VBAgg-Exp
99.2 (39.8)
1.21 (0.19)
NN-Exp164.4 (127.8) 1.82 (0.29)
19

Table 2: p-values from a Wilcoxon signed-rank test for Nyström-Exp versus the methods below for
Bag NLL and MSE for the malaria dataset. The null hypothesis is Nyström-Exp performs equal orworse than the considered method on the test bag performance.
NLL
MSE
Constant
0.0009766 0.0009766
NN-Exp0.002930.0009766
VBAgg-Sq-Obj
0.1162
0.958
VBAgg-Sq
0.13770.1611VBAgg-Exp-Obj 0.08008
1.0
VBAgg-Exp
0.09668
0.958

Table 3: p-values from a Wilcoxon signed-rank test for VBAgg-Sq versus the methods below for Bag
NLL and MSE for the malaria dataset. The null hypothesis is VBAgg-Sq performs equal or worse
than the considered method on the test bag performance.
NLL
MSE
Constant
0.0009766 0.0009766
NN-Exp0.01855
0.001953
VBAgg-Sq-Obj
0.62340.9861
Nyström-Exp
0.88380.8623
VBAgg-Exp-Obj 0.6875
1.0
VBAgg-Exp
0.34770.9346

G.1

Predicted log malaria incidence rate for various models

Constant: Bag level observed incidences This is the baseline with λ̂ai being constant throughout
the bag, as shown in Figure 3. For training, we only use 60% of the data.

Figure 3: Predicted λ̂ai on log scale using constant model, for 3 different re-splits of the data.× denote
non-train set bags.

VBAgg-Sq-Obj This is the VBAgg model with Ψ(v) = v 2 and tuning of hyperparameters is
performed based on training objective, the lower bound to the marginal likelihood, we ignore earlystop and validation set here. The uncertainty of the model seems reasonable, and we also observe thatin general the areas that are not in the training set have higher uncertainties. Furthermore, in all cases,
malaria incidence was predicted to be higher near the river, as discussed in Section 4.2.
20

Figure 4: Top: Predicted λ̂ai on log scale for VBAgg-Sq-Obj. Bottom: Standard deviation of theposterior v in (9) with VBAgg-Sq-Obj.
VBAgg-Sq This is the VBAgg model with Ψ(v) = v 2 and tuning of hyperparameters is performed
based on NLL at the bag level. Predicted incidence are similar to the VBAgg-Sq-Obj model.
The uncertainty of the model is less reasonable here, this is expected behaviour, as we are tuning
hyperparameters based on NLL here. In the first patch, the same parameters was chosen as VBAggSq-Obj.

Figure 5: Top: Predicted λ̂ai on log scale for VBAgg-Sq. Bottom: Standard deviation of the posterior
v in (9) with VBAgg-Sq.
VBAgg-Exp-Obj This is the VBAgg model with Ψ(v) = ev and tuning of hyperparameters
is performed based on training objective, the lower bound to the marginal likelihood, we ignore
early-stop and validation set here. Predicted incidence seem to be stable in general, though somesmoothness is observed. The uncertainty of the model is also not very reasonably here, but thisbehaviour was observed in the Toy experiments, and likely due to an additional lower bound.
21

Figure 6: Top: Predicted λ̂ai on log scale for VBAgg-Exp-Obj.Bottom: Standard deviation of theposterior v in (9) with VBAgg-Exp-Obj.

VBAgg-Exp This is the VBAgg model with Ψ(v) = ev and tuning of hyperparameters is performed
based on NLL. For details, see discussion above for the VBAgg-Exp-Obj model.

Figure 7: Top: Predicted λ̂ai on log scale for VBAgg-Exp. Bottom: Standard deviation of theposterior v in (9) with VBAgg-Exp.

Nyström-Exp This is the Nyström-Exp model, it is clear that while it performs best in terms of bag
NLL, sometimes prediction are too smooth in the pixel space, this is because it optimises directly bag
NLL. This pattern might be seen to be unrealistic, and may cause useful covariates to be neglected.
22

Figure 8: Predicted λ̂ai on log scale for Nyström-Exp.
NN-Exp We can see that the model is not very stable, this can be potentially due to the model
does not have an inbuilt spatial smoothness function unlike other methods. It only uses manifold
regularisation for training. Also, the maximum predicted pixel level intensity rate λ̂ai is over 1000 insome cases, this is clearly physically impossible given λai is rate per 1000 people.

Figure 9: Predicted λ̂ai on log scale for NN-Exp.
G.2

Remote Sensing covariates that provide the existence of a river

Here, we provide figures for some covariates that give information that there is a river as indicated by
the triangles in Figure 2.

Figure 10: Topographic wetness index, measures the wetness of an area, rivers are wetter than others,
as clearly highlighted.
23

Figure 11: Land Surface Temperature at night, river is hotter at night, due to river being able to retain
heat better.

H

Additional Toy Experimental Results

In this section, we provide additional experimental results for the Normal and Poisson model. In
particular, we provide results on test bag level performance, and provide also prediction, calibration
and uncertainty plots.
For the VBAgg model, during the tuning process, it is possible to choose tuning parameters (e.g.
learning rate, multiple-initialisations, landmark choices) based on NLL with an additional validation
set or on the objective L1 on the training set. To compare the difference, we denote the model tuned
on NLL as VBAgg and the model tuned on L1 as VBAgg-Obj. Intuitively, as VBAgg-Obj attempts toobtain as tight a bound to the marginal likelihood, we would expect better performance in calibration,i.e. more accurate uncertainties.
For calibration plots, we compute the α quantiles of the approximated posterior distribution andconsider the ratio of times the underlying rate parameter λai (or µai for the normal model) appear
inside the quantiles of the posterior distribution. If the model provides good uncertainties/calibration,
we should expect to see the quantiles to match with the observed ratio.
In the case of Ψ(v) = v 2 , the approximated posterior distribution is simply a non-central χ2
distribution, while for Ψ(v) = ev , this is a log-normal distribution. For the Normal Model, it issimply a normal distribution, as we do not have any transformations. Calibration plots can be foundin Figure 20 and Figure 21 for the Normal Model, with Figure 14 and Figure 15 for the Poisson
Model.
For uncertainty plots, we plot the standard deviation of the posterior of v ∼ N (ma , S a ) (i.e. before
transformation through Ψ), as this provides better interpretability. Uncertainty plots can be found in
Figure 17 and 23
To demonstrate statistical significance of our result, we aggregate the repetitions in each experiment
for each method and consider a one sided rank permutation test (Wilcoxon signed-rank test) to see
whether VBAgg is statistically significant better than other approaches for individual NLL and MSE.
H.1
H.1.1

Poisson Model
Swiss Roll Dataset

We provide additional results here for the experimental settings that we consider.
24

2.20

Bag Avg NLL

Individual Avg NLL

2.22

2.18
2.16
2.14
150

200

250

Number of Bags

300

350

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1

100

150

100

150

200

250

300

350

200

250

300

350

Number of Bags

250
200

Bag MSE

Individual MSE

100

4.84
4.82
4.80
4.78
4.76
4.74
4.72
4.70

150
100
50

100

150

200

250

Number of Bags

VBAgg-Exp-Obj
VBAgg-Sq-Obj

300

VBAgg-Exp
VBAgg-Sq

350
NN-Exp
NN-Sq

Number of Bags

Nystrom-Exp
Nystrom-Sq

Figure 12: Varying number of bags over 5 repetitions.Left Column: Individual average NLL and
MSE on train set. Right Column: Bag average NLL and MSE on test set (of size 500). Constantprediction NLL and MSE is 2.23 and 0.85 respectively. bag-pixel model prediction NLL is above 2.4
and MSE is above 3.0, hence not shown on graph.

The varying number of bags experimental results is found in Figure 12, with the corresponding tableof p-values in Table 4, 5 demonstrating statistical significance of the VBAgg-Exp and VBAgg-Sq
method. Similarly, the varying number of individuals per bag through Nmean experimental result
can be found in Figure 13, with the corresponding table of p-values in Table 6, 7. The comparison
between VBAgg-Exp and VBAgg-Sq was found to be non-significant.

Table 4: p-values from a Wilcoxon signed-rank test for VBAgg-Sq versus the methods below for
the varying number of bags experiment for the Poisson model. The null hypothesis is VBAgg-Sq
performs equal or worse than NN or Nyström in terms of individual NLL or MSE on the train set.
NLL
MSE
NN-Exp6.98e−06 0.00025
Nyström-Exp 0.00048
0.00015

Table 5: p-values from a Wilcoxon signed-rank test for VBAgg-Exp versus the methods below for
the varying number of bags experiment for the Poisson model. The null hypothesis is VBAgg-Exp
performs equal or worse than NN or Nyström in terms of individual NLL or MSE on the train set.
NLL
MSE
NN-Exp2.48e−06 2.48e−05
Nyström-Exp 0.0005
0.00025
25

Bag Average NLL

Individual Average NLL

5.4
2.17
2.16
2.15
2.14
300

400

Nmean

500

5.0
4.8

600

0.40

600

0.35

500

0.30

400

Bag MSE

Individual MSE

200

5.2

0.25
0.20

200

0.10

100
200
VB-Exp
VB-Sq

300

400

Nmean

VBAgg-Exp-Obj
VBAgg-Sq-Obj

500

0

600

300

Nmean

400

500

600

200

300

Nmean

400

500

600

300

0.15
0.05

200

Nystrom-Exp
NN-ExpFigure 13: Varying number of individuals per bag Nmean over 5 repetitions.Left Column: Individual
average NLL and MSE on train set. Right Column: Bag average NLL and MSE on test set (of size
500). Constant prediction NLL and MSE is 2.23 and 0.85 respectively.

Table 6: p-values from a Wilcoxon signed-rank test for VBAgg-Sq versus the methods below for
the varying number of individuals per bag experiment for the Poisson model. The null hypothesis is
VBAgg-Sq performs equal or worse than NN or Nyström in terms of individual NLL or MSE on the
train set.
NLL
MSE
NN-Exp1.81e−05 9.53e−06
Nyström-Exp 0.062
0.041

Table 7: p-values from a Wilcoxon signed-rank test for VBAgg-Exp versus the methods below for
the varying number of individuals per bag experiment for the Poisson model. The null hypothesis is
VBAgg-Exp performs worse than NN or Nyström in terms of individual NLL or MSE on the train set.
NLL
MSE
NN-Exp6.68e−05 0.00016
Nyström-Exp 0.049
0.062

Calibration Plots for the Swiss Roll Dataset In Figure 14 and 15, we provide calibration results
for both experiments that we have considered. See top of Appendix H for a further details. It is clearthat while VBAgg-Sq-Obj and VBAgg-Sq provides good calibration in general, this is not the casefor VBAgg-Exp-Obj and VBAgg-Exp. This is not surprising as the VBAgg-Exp methods uses anadditional lower bound.
26

Absolute Error in Coverage

0.95

Absolute Error in Coverage

0.95

Absolute Error in Coverage

Absolute Error in Coverage
Absolute Error in Coverage
Absolute Error in Coverage

0.95

Number of Bags 100
0.75
0.50
0.25
0.00

0.70

0.75

0.80

0.85

Claimed Coverage

0.90

Number of Bags 200
0.75
0.50
0.25
0.00

0.70

0.75

0.80

0.85

Claimed Coverage

0.90

Number of Bags 300
0.75
0.50
0.25
0.00

0.70

0.75

0.80

0.85

Claimed Coverage

VBAgg-Exp
VBAgg-Sq

0.90

Number of Bags 150
0.75
0.50
0.25
0.00

0.70

0.75

0.80

0.85

Claimed Coverage

0.90

0.95

0.90

0.95

0.90

0.95

Number of Bags 250
0.75
0.50
0.25
0.00

0.70

0.75

0.80

0.85

Claimed Coverage

Number of Bags 350
0.75
0.50
0.25
0.00

0.70

0.75

VBAgg-Exp-Obj
VBAgg-Sq-Obj

0.80

0.85

Claimed Coverage

Figure 14: Absolute Error in coverage from 70% to 95% for the increasing number of bags experiment
for the Poisson Model. Shaded regions highlight the standard deviation. Perfect coverage would
provide a straight line at 0 error.

Nmean = 300

0.8Absolute Error in Coverage

Absolute Error in Coverage

Nmean = 150

0.6
0.4
0.2
0.0

0.70

0.75

0.80

0.85

Claimed Coverage

0.90

0.8
0.6
0.4
0.2
0.0

0.95

0.70

0.75

0.8
0.6
0.4
0.2
0.0

0.70

0.75
VBAgg-Exp
VBAgg-Sq

0.80

0.85

Claimed Coverage

0.85

0.90

0.95

0.90

0.95

Nmean = 600
Absolute Error in Coverage

Absolute Error in Coverage

Nmean = 450

0.80

Claimed Coverage

0.90

0.95

0.8
0.6
0.4
0.2
0.0

0.70

0.75

0.80

0.85

Claimed Coverage

VBAgg-Exp-Obj
VBAgg-Sq-Obj

Figure 15: Absolute Error in coverage from 70% to 95% for the increasing number of individuals
per bag Nmean and Nstd for the Poisson Model. Shaded regions highlight the standard deviation.
Perfect coverage would provide a straight line at 0 error.
27

Prediction and uncertainty plots In Figure 16 and 17, we provide some prediction plots fordifferent models, and uncertainties for VBAgg models.

NN-ExpNystrom-Exp
15

15

7

10

10
6

6

5

5
5

0

0

4

5

5

10

5

0

5

10

0

5

4

10
3

10

7

3

20
15
10
5

10

5

0

5

10

0

20
15
10
5

Figure 16: Individual predictions on the train set for the swiss roll dataset with 150 bags for NN and
Nyström model. Here Nmean = 150, with Nstd = 50.

VBAgg-Obj-Exp

10

5

0

5

10

VBAgg-Obj-Exp uncertainty
15
10
5
0
5
10
20
15
10
05

7

15
10
5
0
5
10
20
15
10
05

7

6
5
4
3
10

VBAgg-Obj-Sq

10

5

0

5

10

5

0

5

10

15
10
5
0
5
10
20
15
10
05

0.06
0.05
0.04
0.03
0.02

VBAgg-Obj-Sq uncertainty
6
5
4
3
10

5

0

5

10

15
10
5
0
5
10
20
15
10
05

0.250
0.225
0.200
0.175
0.150
0.125
0.100
0.075

Figure 17: Predictions and uncertainty on the swiss roll dataset with 150 bags for the VBAgg-Obj
models. Here Nmean = 150, with Nstd = 50. For uncertainty, we plot the standard deviation of theposterior of v, coming from v a ∼ N (ma , S a ) in (9).

28

H.2

Normal Model

H.2.1

Swiss Roll Dataset

In this section, we provide some experimental results for the Normal model, where throughout weassume τia = τ , same for all individuals.
We consider the same swiss roll dataset as in the Poisson model, here the colour of each point to be
the underlying mean µai . We then consider yia ∼ N (µa , τ ) with τ = 0.1, hence bag observations are
PNa a
PNa a
given by y a = i=1
yi ∼ N (µa , Na τ ) with µa = i=1µi . Here, the goal is to predict µai and τ ,given bag observations y a only. The results for the experiments are shown below in Figure 18 and
Figure 19, which shows the VBAgg outperforming the NN and Nyström model. To show statisticalsignificance, we also report the corresponding table of p-values in Table 8 and Table 9. Furthermore,
we would also like to point out that the VBAgg is well calibrated as shown in Figure 20.
VBAgg
VBAgg-Obj
Nystrom
NN

0.55
0.60
0.65
0.70

VBAgg
VBAgg-Obj
Nystrom
NN

2.00
1.95

Bag Avg NLL

Individual Avg NLL

0.50

0.75

1.90
1.85
1.80
1.75

0.80

1.70
150

200

250

Number of Bags

0.0045350

VBAgg
VBAgg-Obj
Nystrom
NN

0.0040Individual MSE

300

0.00350.00300.00250.0020

100

150

200

250

Number of Bags

1.75

300

350

VBAgg
VBAgg-Obj
Nystrom
NN

1.50

Bag MSE

100

1.25
1.00
0.75
0.50

0.00150.25

0.0010
100

150

200

250

Number of Bags

300

350

100

150

200

250

Number of Bags

300

350

Figure 18: Varying number of bags over 5 repetitions for the Normal model.Left Column: Individual
average NLL and MSE on train set. Right Column: Bag average NLL and MSE on test set (of size
500). Constant model individual MSE is 0.04.

Table 8: p-values from a Wilcoxon signed-rank test for VBAgg versus the methods below for thevarying number of bags experiment for the Normal model. The null hypothesis is VBAgg performs
equal or worse than NN or Nyström in terms of individual NLL or MSE on the train set.
NLL
MSE
NN
5.96e−07 4.79e−09
Nyström 4.01e−08 6.52e−09

Table 9: p-values from a Wilcoxon signed-rank test for VBAgg versus the methods below for thevarying number of individuals per bag Nmean experiment for the Normal nodel. The null hypothesis
is VBAgg performs worse than NN or Nyström in terms of individual NLL or MSE on the train set.
NLL
MSE
NN
4.77e−06 4.77e−06
Nyström 4.77e−06 4.77e−06
Calibration Plots for the Swiss Roll Dataset In Figure 20 and 21, we provide calibration results
for both experiments that we have considered. See top of Appendix H for further details. It is clear
29

2.6

Bag Average NLL

Individual Average NLL

0.65
0.70
0.75
0.80
200

300

400

Nmean

500

600

2.4
2.2
2.0
1.8
1.6

200

300

Nmean

400

500

600

200

300

Nmean

400

500

600

7
6
5

Bag MSE

Individual MSE

0.004
0.003
0.002

4
3
2
1

0.001
200
VBAgg

300

400

Nmean

VBAgg-Obj

500
Nystrom

600

0NN

Figure 19: Varying number of individuals per bag Nmean over 5 repetitions.Left Column: Individual
average NLL and MSE on train set. Right Column: Bag average NLL and MSE on test set (of size
500). Constant model individual MSE is 0.039.
that VBAgg-Obj has better calibration in general, this is not surprising as it is tuned based on thecorrect objective, rather than NLL.

30

0.2
0.1
0.0

0.70

0.75

0.80

0.85

Claimed Coverage

0.90

0.95

Number Of Bags 200

0.4

VBAgg
VBAgg-Obj

0.3
0.2
0.1
0.0

0.70

0.75

0.80

0.85

Claimed Coverage

0.90

0.95

Number Of Bags 300

0.4

VBAgg
VBAgg-Obj

0.3
0.2
0.1
0.0

0.70

0.75

0.80

0.85

Claimed Coverage

0.90

0.95

Absolute Error in Coverage

0.3

Absolute Error in Coverage

VBAgg
VBAgg-Obj

Absolute Error in Coverage

Absolute Error in Coverage
Absolute Error in Coverage
Absolute Error in Coverage

Number Of Bags 100

0.4

Number Of Bags 150

0.4

VBAgg
VBAgg-Obj

0.3
0.2
0.1
0.0

0.70

0.75

0.80

0.85

Claimed Coverage

0.90

0.95

Number Of Bags 250

0.4

VBAgg
VBAgg-Obj

0.3
0.2
0.1
0.0

0.70

0.75

0.80

0.85

Claimed Coverage

0.90

0.95

Number Of Bags 350

0.4

VBAgg
VBAgg-Obj

0.3
0.2
0.1
0.0

0.70

0.75

0.80

0.85

Claimed Coverage

0.90

0.95

Figure 20: Absolute Error in coverage from 70% to 95% for the increasing number of bags experiment
for the Normal Model. Shaded regions highlight the standard deviation. Perfect coverage would
provide a straight line at 0 error.

Nmean = 150
0.5
0.4
0.3
0.2
0.1
0.0

0.70

0.75

0.80

0.85

Claimed Coverage

0.90

VBAgg
VBAgg-Obj

0.6

Absolute Error in Coverage

Absolute Error in Coverage

Nmean = 300
VBAgg
VBAgg-Obj

0.6

0.5
0.4
0.3
0.2
0.1
0.0

0.95

0.70

0.75

Nmean = 450

0.4
0.3
0.2
0.1
0.70

0.75

0.80

0.85

Claimed Coverage

0.90

0.95

0.90

0.95

VBAgg
VBAgg-Obj

0.6

Absolute Error in Coverage

Absolute Error in Coverage

0.5

0.0

0.85

Nmean = 600
VBAgg
VBAgg-Obj

0.6

0.80

Claimed Coverage

0.5
0.4
0.3
0.2
0.1
0.0

0.70

0.75

0.80

0.85

Claimed Coverage

0.90

0.95

Figure 21: Absolute Error in coverage from 70% to 95% for the increasing number of individuals per
bag Nmean and Nstd for the Normal Model. Shaded regions highlight the standard deviation. Perfectcoverage would provide a straight line at 0 error.

31

Prediction and uncertainty plots Here, we provide some prediction plots for different models.

NN

Nystrom
15
10

1.0

15
10

0.8

5

0.6

0

10

10

5

0

5

10

0

0.8

5
0.6

5

1.0

0
0.4

5

0.2

20
15
10
5

10

10

5

0

5

10

0

0.4

0.2

20
15
10
5

Figure 22: Individual predictions on the train set for the swiss roll dataset with 150 bags for NN and
Nyström model. Here Nmean = 150, with Nstd = 50.

1.0

VBAgg-Obj uncertainty

VBAgg-Obj

0.25

15

15
0.8

10
5

10

0
5

0.20

5

0.6

0
0.4

0.15

5
0.10

10

10
0.2

10

5

0

5

10

0

20
15
10
5

10

5

0

5

10

0

20
15
10
5

0.05

Figure 23: Predictions and uncertainty on the swiss roll dataset with 150 bags for the VBAgg-Obj
model. Here Nmean = 150, with Nstd = 50. For uncertainty, we plot the standard deviation of theposterior of v, coming from v a ∼ N (ma , S a ) in (9).

H.2.2

Elevators Dataset

For a real dataset experiment, we consider the elevators dataset11 , which is a large scale regression
dataset12 containing 16599 instances, with each instance ∈ R17 . This dataset is obtained from thetask of controlling F16 aircraft, with the label y being a particular action taken on the elevators of theaircraft ∈ R. For the model formulation we assume each label follows a normal distribution, i.e.
yl ∼ N (µl , τ ), where τ is a fixed quantity to be learnt. In practice, we can imagine the action taken
may differ according to the operator.
In order formulate this dataset in an aggregate data setting, we sample bag sizes from a negative
binomial distribution as before, with Nmean = 30 and Nstd = 15, and also take wia = 1. To placeobservations into bags, similar to the swiss roll dataset, we consider a particular covariate, and placeinstances into bags based on the ordering of the covariate. We now have the bag-level model given byy a ∼ N (µa , Na τ ), with individual model yia ∼ N (µai , τ ) and it is our goal to predict µai (and alsoinfer τ ), given only y a . After the bagging process, we obtain approximately 225 bags for training,and 33 bags each for early stopping, validation and testing (for bag level performance). Further,in order to neglect variables that do not provide signal, we use an ARD kernel for the VBAgg and1112This dataset is publicly available at http://sci2s.ugr.es/keel/dataset.php?cod=94
We have removed one column that is almost completely sparse.

32

Table 10: Results for the Normal Model on the elevators dataset with 50 repetitions. Indiv represents
individuals on train set here, while bag performance is measured on a test set. Numbers in brackets
denotes p-values from a Wilcoxon signed-rank test for VBAgg versus the method. The null hypothesis
is VBAgg performs equal or worse than NN or Nyström in terms of individual NLL or MSE on the
train set. It is also noted MSE is computed on the observed yia or y a , rather than the unknown µai orµa .
Indiv NLL
Bag NLL Indiv MSE
Bag MSE
Constant
N/A
N/A
0.010
0.366
VBAgg
-1.69
0.003
0.00180.052
VBAgg-Obj -1.71-0.02
0.00180.052
Nyström
−1.57(1.5e−13)
0.003
0.0024 (8.9e−16) 0.041
NN-1.64 (0.0001258) 0.082
0.0021 (8.8e−10) 0.041
Nyström model, as below:d1X 1(xk − yk )2
kard (x, y) = γscale exp −
2
`k

!
(38)

k=1

and learn kernel parameters γscale and {`k }dk=1 . We repeat this process and splitting of the dataset
50 times and report individual NLL results, and also MSE results in Table 10. From the results,
we observe that the VBAgg model performs better the Nyström and NN model, with statisticalsignificance.

33