Improved prediction accuracy for disease risk mapping using Gaussian Process stacked generalisation

Maps of infectious disease---charting spatial variations in the force of infection, degree of endemicity, and the burden on human health---provide an essential evidence base to support planning towards global health targets. Contemporary disease mapping efforts have embraced statistical modelling approaches to properly acknowledge uncertainties in both the available measurements and their spatial interpolation. The most common such approach is that of Gaussian process regression, a mathematical framework comprised of two components: a mean function harnessing the predictive power of multiple independent variables, and a covariance function yielding spatio-temporal shrinkage against residual variation from the mean. Though many techniques have been developed to improve the flexibility and fitting of the covariance function, models for the mean function have typically been restricted to simple linear terms. For infectious diseases, known to be driven by complex interactions between environmental and socio-economic factors, improved modelling of the mean function can greatly boost predictive power. Here we present an ensemble approach based on stacked generalisation that allows for multiple, non-linear algorithmic mean functions to be jointly embedded within the Gaussian process framework. We apply this method to mapping Plasmodium falciparum prevalence data in Sub-Saharan Africa and show that the generalised ensemble approach markedly out-performs any individual method.


Introduction
Infectious disease mapping with model-based geostatistics [1] can provide a powerful synthesis of the available evidence base to assist surveillance systems and support progress towards global health targets, revealing the geographical bounds of disease occurrence and the spatial patterns of transmission intensity and clinical burden. A recent review found that out of 174 infectious diseases with a strong rational for mapping only 7 (4%) have thus far been comprehensively mapped [2], with a lack of accurate, population representative, geopositioned data having been the primary factor impeding progress. In recent years this has begun to change as increasing volumes of spatially referenced data are collected from both cross-sectional household surveys and web-based data sources (e.g. Health Map [3]), bringing new opportunities for scaling up the global mapping of diseases. Alongside this surge in new data, novel statistical methods are needed that can generalise to new data accurately while remaining computationally tractable on large datasets. In this paper we will introduce one such method designed with these aims in mind.

2/33
Owing to both a long history of published research in the field and a widespread appreciation amongst endemic countries for the value of cross-sectional household surveys as guides to intervention planning, malaria is an example of a disease that has been comprehensively mapped.
Over the past decade, volumes of publicly-available malaria prevalence data-defined as the proportion of parasite positive individuals in a sample-have reached sufficiency to allow for detailed spatio-temporal mapping [4]. From a statistical perspective, the methodological mainstay of these malaria prevalence mapping efforts has been Gaussian process regression [5][6][7][8].
Gaussian processes are a flexible semi-parametric regression technique defined entirely through a mean function, µ(·), and a covariance function, k(·, ·). The mean function models an underlying trend, such as the effect of environmental/socio-economic factors, while the covariance function applies a Bayesian shrinkage to residual variation from the mean such that points close to each other in space and time tend towards similar values. The resulting ability of Gaussian processes to strike a parsimonious balance in the weighting of explained and unexplained spatio-temporal variation has led to their near exclusive use in contemporary studies of the geography of malaria prevalence [1,4,7,[9][10][11].
Outside of disease mapping, Gaussian processes have been used for numerous applications in machine learning, including regression [1,5,6], classification [5], and optimisation [12]; their popularity leading to the development of efficient computational techniques and statistical parametrisations. A key challenge for the implementation of Gaussian process models arises in statistical learning (or inference) of any underlying parameters controlling the chosen mean and covariance functions. Learning is typically performed using Markov Chain Monte Carlo (MCMC) or by maximizing the marginal likelihood [5], both of which are made computationally demanding by the need to compute large matrix inverses returned by the covariance function.
The complexity of this inverse operation is O(n 3 ) in computation and O(n 2 ) in storage in the naive case [5], which imposes practical limits on data sizes [5]. MCMC techniques may be further confounded by mixing problems in the Markov chains. These challenges have necessitated the use of highly efficient MCMC methods, such as Hamiltonian MCMC [13], or posterior approximation approaches, such as the integrated nested Laplace approximation [14], expectation propagation [5,15,16], and variational inference [17,18]. Many of these methods adopt finite dimensional representations of the covariance function yielding sparse precision PLOS 3/33 matrices, either by specifying a fully independent training conditional (FITC) structure [19] or by identifying a Gaussian Markov Random Field (GMRF) in approximation to the continuous process [20].
Alongside these improved methods for inference, recent research has focussed on model development to increase the flexibility and diversity of parametrisations for the covariance function, with new techniques utilising solutions to stochastic partial differential equations (allowing for easy extensions to non-stationary and anisotropic forms [20]), the combination of kernels additively and multiplicatively [21], and various spectral representations [22].
One aspect of Gaussian processes that has remained largely neglected is the mean function which is often-and indeed with justification in some settings-simply set to zero and ignored.
However, in the context of disease mapping, where the biological phenomena are driven by a complex interplay of environmental and socioeconomic factors, the mean plays a central role in improving the predictive performance of Gaussian process models. Furthermore, it has also been shown that using a well-defined mean function can allow for simpler covariance functions (and hence simpler, scalable inference techniques) [23].
The steady growth of remotely-sensed data with incredible spatio-temporal richness [24] combined with well-developed biological models [25] has meant that there is a rich suite of environmental and socio-economic covariates currently available. In previous malaria mapping efforts these covariates have been modelled as simple linear predictors [7,9,10] that fail to capture complex non-linearities and interactions, leading to a reduced overall predictive performance. Extensive covariate engineering can be performed by introducing large sets of non-linear and interacting transforms of the covariates, but this brute force combinatorial problem quickly becomes computationally inefficient [4,24].
In the field of machine learning and data science there has been great success with algorithmic approaches that neglect the covariance and focus on learning from the covariates alone [26,27].
These include tree based algorithms such as boosting [28] and random forests [29], generalized additive spline models [30,31], multivariate adaptive regression splines [32], and regularized regression models [33]. The success of these methods is grounded in their ability to manipulate the bias-variance tradeoff [34], capture interacting non-linear effects, and perform automatic PLOS 4/33 covariate selection. The technical challenges of hierarchically embedding these algorithmic methods within the Gaussian process framework are forbidding, and many of the approximation methods that make Gaussian process models computationally tractable would struggle with their inclusion. Furthermore, it is unclear which of these approaches would best characterize the mean function when applied across different diseases and settings. In this paper we propose a simplified embedding method based on stacked generalisation [35,36] that focuses on improving the mean function of a Gaussian process, thereby allowing for substantial improvements in the predictive accuracy beyond what has been achieved in the past.

Results
The results of our analysis are summarised in figure 1 where pairwise comparisons of MSE versus MAE versus correlation are shown. We found a consistent ranking pattern in the generalisation performance of the methods across the Eastern, Southern and Western African regions (Figures 1(a),1(b) and 1(d)), with the stacked Gaussian process approach presented in this paper outperforming all other methods. The constrained weighted mean stacked approach was the next best method followed by the standard Gaussian process (with a linear mean) and Gradient boosted trees. Random forests, multivariate adaptive regression splines and generalised additive splines all had similar performance and the worst performing method was the elastic net regularised regression. For the North Eastern region (Figure 1(c)), again the stacked Gaussian process approach was the best performing method, but the standard Gaussian process performed better than the constrained weighted mean stacked approach, though only in terms of MAE and MSE.

Discussion
All the level zero generalisation methods used in this paper have been previously applied to a diverse set of machine learning problems and have track records of good generalisability [26].
For example, in closely related ecological applications, these level zero methods have been shown to far surpass classical learning approaches [37]. However, as introduced by Wolpert [35], rather than picking one level zero method, an ensemble via a second generaliser has the ability to improve prediction beyond that achievable by the constituent parts [38]. Indeed, in all previous applications [35,36,39,40] ensembling by stacking has consistently produced the best predictive models across a wide range of regression and classification techniques. The most popular level one generaliser is the constrained weighted mean with convex combinations. The key attraction of this level one generaliser is the ease of implementation and theoretical properties [38,39].
In this paper we show that, for disease mapping, stacking using Gaussian processes is more predictive and generalises better than both any single method in isolation, and the more common stacking approach using a constrained weighted mean.
The key benefit of stacking is summarised in equation 6 where the total error of an ensemble model can be reduced by using multiple, very different, but highly predictive models. However, stacking using a constrained weighted mean only ensures that the predictive power of the covariates are fully utilised, and does not exploit the predictive power that could be gained from characterising any residual covariance structure. The standard Gaussian process suffers from the inverse situation where the covariates are underexploited, and predictive power is instead gained from leveraging residual spatio-temporal covariance. In a standard Gaussian process the mean function is usually paramaterised through simple linear basis functions [41] that are often unable to model the complex non linear interactions needed to correctly capture the true underlying mean. This inadequacy is best highlighted by the poor generalisation performance of the elastic net regularised regression method across all regions. The trade off between the variance explained by the covariates versus that explained by the covariance function will undoubtedly vary from setting to setting. For example in the Eastern, Southern and Western African regions, the constrained weighted mean stacking approach performs better than the standard Gaussian process, and the level 0 gradient boosted trees generaliser performs almost as well as the standard Gaussian process. For these regions, this shows a strong influence of the covariates on the underlying process. In contrast, for the North Eastern African region, the standard Gaussian process does better than both the constrained weighted mean approach (in terms of error not correlation) and all of the Level 0 generalisers, suggesting a weak influence of the covariates. However, for all zones, the stacked Gaussian process approach is consistently the best approach across all predictive metrics. By combining both the power of Gaussian processes to characterise a complex covariance structure, and multiple algorithmic approaches to fully exploit the covariates, the stacked Gaussian process approach combines the best of both worlds, and predicts well in all settings.
This paper introduces one way of stacking that is tailored for spatio-temporal data. However the same principles are applicable to purely spatial or purely temporal data, settings in which Gaussian process models excel. Additionally, there is no constraint on the types of level 0 generalisers than can be used; dynamical models of disease transmission e.g. [ implemented (e.g. [40]).
Gaussian processes are increasingly being used for expensive optimisation problems [44] and Bayesian quadrature [45]. In current implementations both of these applications are limited to low dimensional problems typically with n < 10. Future work will explore the potential for stacking to extend these approaches to high dimensional settings. The intuition is that the level 0 generalisers can accurately and automatically learn much of the latent structure in the data, including complex features like non-stationarity, which are a challenge for Gaussian processes. Learning this underlying structure through the mean can leave a much simpler residual structure [23] to be modelled by the level one Gaussian process.
In this paper we have focused primarily on prediction, that is neglecting any causal inference and only searching for models with the lowest generalisation error. Determining causality from the complex relationships fitted through the stacked algorithmic approaches is difficult, but empirical methods such as partial dependence [28] or individual conditional expectation [46] plots can be used to approximate the marginal relationships from the various covariates. Similar statistical techniques can also be used to determine covariate importance.
Increasing volumes of data and computational capacity afford unprecedented opportunities to scale up infectious disease mapping for public health uses [47]. Maps of diseases and socioeconomic indicators are increasingly being used to inform policy [4,48], creating demand for methods to produce accurate estimates at high spatial resolutions. Many of these maps can subsequently be used in other models but, in the first instance, creating these maps requires continuous covariates, the bulk of which, come from remotely sensed sources. For many indicators, such as HIV or Tuberculosis, these remotely sensed covariates serve as proxies for complex phenomenon, and as such, the simple mean functions in standard Gaussian processes are insufficient to predict with accuracy and low generalisation error. The stacked Gaussian process approach introduced in this paper provides an intuitive, easy to implement method that predicts accurately through exploiting information in both the covariates and covariance structure.

Gaussian process regression
We define our response, y s,t = {y (s,t) [1] , ..., y (s,t)[n] }, as a vector of n empirical logit transformed malaria prevalence surveys at location-time pairs, (s, t) denoting a corresponding n×m design matrix of m covariates (see Section ). The likelihood of the observed response is P(y s,t |f s,t , X s,t , θ), which we will write simply as P(y|f (s, t), θ), suppressing the spatio-temporal indices for ease of notation. Naturally Following Bayes theorem the posterior distribution resulting from this hierarchy becomes: where the denominator in Equation 2 is the marginal likelihood, P(y). , for a given θ is also Gaussian with form: where Σ y|(s,t),θ = Σ θ + 1σ 2 e . For specific details on the parametrisation of Sigma see Appendix When examining the conditional expectation in equation 5 and splitting the summation into terms µ (s ,t )|θ and Σ (s ,t ),(s,t)|θ Σ −1 y|(s,t),θ y − µ (s,t)|θ , it is clear that the first specifies a global underlying mean while the second augments the residuals from that mean by the covariance function. Clearly if the mean function fits the data perfectly the covariance in the second term of the expectation would drop out and conversely if the mean function is zero, then only the covariance function would model the data. This expectation therefore represents a balance between the underlying trend and the residual correlated noise.
In most applications of Gaussian process regression a linear mean function (µ θ = X s,t β) is used, where β is a vector of m coefficients. However, when a rich suite of covariates is available this linear mean may be sub-optimal, limiting the generalisation accuracy of the model. To improve on the linear mean, covariate basis terms can be expanded to include parametric nonlinear transforms and interactions, but finding the optimal set of bases is computationally demanding and often leaves the researcher open to data snooping [49]. In this paper we propose using an alternative two stage statistical procedure to first obtain a set of candidate non-linear mean functions using multiple different algorithmic methods fit without reference to the assumed spatial covariance structure and then include those means in the Gaussian process via stacked generalisation.

Stacked generalisation
Stacked generalisation [35], also called stacked regression [36], is a general ensemble approach to combining different models. In brief stacked generalisers combine different models together to produce a meta-model with equal or better predictive performance than the constituent PLOS 11/33 parts [39]. In the context of malaria mapping our goal is to fuse multiple algorithmic methods with Gaussian process regression to both fully exploit the information contained in the covariates and model spatio-temporal correlations.
To present stacked generalisation we begin by introducing standard ensemble methods and show that stacked generalisation is simply a special case of this powerful technique. To simplify notation we suppress the spatio-temporal index and dependence on θ. Consider L models, with outputsỹ i (x), i = 1, ..., L. The choice of these models is described in appendix 1.1. We denote the true target function as f (x) and can therefore write the regression equation as The average sum-of-squares error for model i is defined as Our goal is to estimate an ensemble model across all L models, denoted as M (ỹ 1 , ..,ỹ L ). The simplest choice for C is an average across all models M (ỹ 1 , ..,ỹ L ) =ỹ avg (x) = 1 L L i=1ỹ i (x). However this average assumes that the error of all models are the same, and that all models perform equally well. The assumption of equal performance may hold when using variants of a single model (i.e. bagging) but is unsuitable when very different models are used. Therefore a simple extension would be to use a weighted mean across models M (ỹ 1 , ..,ỹ L ) =ỹ wavg (x) = L i=1 β iỹi (x) subject to constraints β > 0 ∀i, L i=1 β = 1 (convex combinations). These constraints prevent extreme predictions in well predicting models and impose the sensible inequalityỹ min (x) ≤ỹ wavg (x) ≤ỹ max [36]. The optimal βs can be found by quadratic programming or by Bayesian linear regression with a Dirichlet/categorical prior on the coefficients. One particularly interesting result of combining models using this constrained weighted mean is the resulting decomposition of error into two terms [38] Equation 6 is a reformulating of the standard bias-variance decomposition [34] where first term describes the average error of all models, and the second (termed the ambiguity) is the spread of each member of the ensemble around the weighted mean, and measures the disagreement among models. Equation 6 shows that combining multiple models with low error but with large disagreements produces a lower overall error. It should be noted that equation 6 makes the assumption that y(x) = f (x).

12/33
Combination of models in an ensemble as described above can potentially lead to reductions in errors. However the ensemble models introduced so far are based only on training data and therefore neglect the issue of model complexity and tell us nothing about the ability to generalise to new data. To state this differently, the constrained weighted mean model will always allocate the highest weight to the model that most over fits the data. The standard method of addressing this issue is to use cross validation as a measure of the generalisation error and select the best performing of the L models. Stacked generalisation provides a technique to combine the power of ensembles described above, but also produces models that can generalise well to new data. The principle idea behind stacked generalisation is to train L models (termed level 0 generalisers) and generalise their combined behaviour via a second model (termed the level 1 generaliser). Practically this is done by specifying a K-fold cross validation set, training all L level 0 models on these sets and using the cross validation predictions to train a level 1 generaliser. This calibrates the level 1 model based on the generalisation ability of the level 0 models. After this level 1 calibration, all level 0 models are refitted using the full data set and these predictions are used in the level 1 model without refitting (This procedure is more fully described in algorithm 1 and the schematic design shown in Appendix figure 3). The combination of ensemble modelling with the ability to generalise well has made stacking one of the best methods to achieve state-of-the art predictive accuracy [36,39,50].
Defining the most appropriate level 1 generaliser based on a rigorous optimality criteria is still an open problem, with most applications using the constrained weighted mean specified above [36,39]. Using the weighted average approach can be seen as a general case of cross validation, where standard cross validation would select a single model by specifying a single β i as 1 and all other β i s as zero. Additionally, it has been shown that using the constrained weighted mean method will perform asymptotically as well as the best possible choice among the family of weight combinations [39].
Here we suggest using Gaussian process regression as the level 1 generaliser. Revisiting equation 4 we can replace µ (s ,t )|θ with a linear stacked function µ (s ,t )|θ = L i=1 β iỹi (s , t ) across L level 0 generalisers, where the subscript denotes predictions from the ith Level 0 generaliser (see algorithm 1. We also impose inequality constraints on β i such that β > 0 ∀i, L i=1 β = 1. This constraint allows the βs to approximately sum to one and helps computational tractability.

13/33
It should be noted that empirical analysis suggests that simply imposing β > 0 is normally sufficient [36].
The intuition in this extended approach is that the stacked mean of the Gaussian process uses multiple different methods to exploit as much predictive capacity from the covariates and then leaves the spatio-temporal residuals to be captured through the Gaussian process covariance function. In Appendix 1.3 we show that this approach yields all the benefits of using constrained weighted mean (equation 6) but allows for a further reduction in overall error from the covariance function of the Gaussian process.

Data, Covariates and Experimental Design
The hierarchical structure most commonly used in infectious disease mapping is that shown in Equation 1. In malaria studies our response data are often discrete random variables representing the number of individuals testing positive for the Plasmodium falciparum (Pf ) malaria parasite, N + , out of the total number tested, N , at a given location. If the response is aggregated (e.g.) from the individual household level to a cluster or enumeration area level usually the centroid of the component sites is used as the spatial position datum. The ratio of N + to N is defined as the parasite rate or prevalence. To aid computation, the response data transformed to the empirical logit of observed prevalence [1,4], which for N 20 can be well approximated by a Gaussian likelihood. Pre-modelling standardisation of the available prevalence data for age and diagnostic type has also been performed on the sample used here, as described in depth in [4,7].
Our analysis is performed over Sub-Saharan Africa with the study area and dataset partitioned into 4 epidemiologically-distinct regions [7]-Eastern Africa, Western Africa, North Eastern Africa, and Southern Africa-each of which was modelled separately (see Figure 2. The covariates (i.e., independent variables) used in this research consist of raster layers spanning the entire continent at a 2.5 arc-minute (5 km x 5 km) spatial resolution. The majority of these raster covariates were derived from high temporal resolution satellite images that were first gap filled [51] to eliminate missing data (resulting primarily from persistent cloud cover over equatorial forests) and then aggregated to create a dynamic (i.e. temporally varying) dataset for every month throughout the study period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). The list of covariates is presented in  Table 1 and detailed information on individual covariates can be found here [24,25,51]. The set of monthly dynamic covariates was further expanded to include lagged versions of the covariate at 2 month, 4 month and 6 month lags. The main objective of this study was to judge the predictive performance of the various generalisation methods and therefore no variable selection or thinning of the covariate set was performed. It should be noted however that many of the level 0 generalisers performed variable selection automatically (e.g. elastic net regression).
The level 0 generalisers used were gradient boosted trees [28,52], random forests [53], elastic net regularised regression [33], generalised additive splines [26,31] and multivariate adaptive regression splines [32]. The level 1 generalisers used were stacking using a constrained weighted mean and stacking using Gaussian process regression. We also fitted a standard Gaussian process for benchmark comparisons with the level 0 and 1 generalisers. Stacked fitting was performed following algorithm 1. Full analysis and K-fold Cross validation was performed 5 times and then averaged to reduce any bias from the choices of cross validation set. The averaged cross validation results were used to estimate the generalisation error by calculating the mean squared error (MSE), mean absolute error (MAE) and the correlation.

16/33
Tables and Algorithms Input X as a n × m design matrix 3: Input y as a n vector of responses 4: Input v cross validation folds 5: choose l, L(y, X) models level zero generalisers 6: define n × l matrix P matrix of predictions 7: for i ← 1, l do 8: fit L i (y, X) 9: predict P ·,i = L i (y, X) 10: ., X gv } and {y g 1 , .., y gv } training set 11: add remaining samples to {X /g 1 , .., X /gv } and {y /g 1 , .., y /gv } testing set 12: define n × l matrix H matrix cross validation of predictions 13: for i ← 1, l do 14: for j ← 1, v do 15: fit L i (y g j , X g j ) 16: predict H /g j ,i = L i (y /g j , X /g j ) 17: choose L * (y, H) model level one generaliser 18: fit L * (y, H) 19: return L * (y, P ) final prediction output

S1 -Generalisation methods
We choose 5 level 0 generalisers to feed into the level 1 generaliser. We chose these 5 due to (a) ease of implementation through existing software packages, (b) the differences in their approaches and (c) a proven track record in predictive accuracy.

Gradient boosted trees and Random forests
Both Gradient boosted trees and random forests produce ensembles of regression trees [54].
Regression trees partition the space of all joint covariate variable values into disjoin regions R j (j = {1, 2, ..., J}) which are represented as terminal nodes in a decision tree. A constant γ j is assigned to each region such that the predictive rule is x ∈ R j → f (x) = γ j [26]. A tree is therefore formally expressed as T (x, θ) = J j=1 γ j I(x ∈ R j ). Gradient boosted trees model the target function by a sum of such trees induced by forward stagewise regression Solving equation 7 is done by functional gradient decent with regularisation performed via shrinkage and cross validation [28]. In our implementations we also stochastically sampled covariates and samples in each stagewise step [52].
Random forests combine trees by bagging [26] where bootstrap samples (B) of covariates and data are used to create an average ensemble of trees: The optimal number of bootstraped trees B is found by cross validation [29].
coarse grid search.

Elastic net regularised regression
Elastic net regularised regression [33] is a penalised linear regression where the coefficients of the regression are found by Where the subscripts on the penalty terms represent the 1 (i.e. lasso) and 2 (i.e. ridge) norms. These norms induce shrinkage in the coefficients of the linear model allowing for better generalisation.
Equation 10 is fitted by numerically maximising the likelihood and optimal parameters for λ 1 , λ 2 are computed by cross validation over the full regularisation path. Fitting was done using the H2O package in R.

Generalised additive splines
Generalised additive splines extend standard linear models by allowing the linear terms to be modelled as nonlinear flexible splines [55].
Where f denotes the second derivative and penalises non smooth functions (that can potentially overfit). Fitting was done using generalised cross validation with smoothness selection done via restricted maximum likelihood. Fitting was done using the mgcv package in R.

27/33
Multivariate adaptive regression splines Multivariate adaptive regression splines [32] build a model of the form Where β i s are coefficients and V i s are basis functions that can either be constant, be a hinge function of the form max(0, x − const), max(0, const − x), or the product of multiple hinge functions. Fitting is done using a two stage approach with a forward pass adding functions and a backward pruning pass via generalised cross validation. Fitting was done using the earth package in R.

S1 -Details of the Gaussian Process
For the Gaussian process spatial component the covariance function is chosen to be Matérn of smoothness ν = 1: with κ = √ 2/ρ an inverse range parameter (for range, ρ), τ a precision (i.e., inverse variance) parameter, and K 1 the modified Bessel function of the second kind and order 1. Typically, θ is defined with elements {log κ, log τ } to ensure positivity via the exponential transformation.
For computational efficiency and scalability we follow the stochastic partial differential equation (SPDE) approach [56] to approximate the continuous Gaussian process in Equation 1 with a discrete Gauss-Markov random field (GRMF) of sparse precision matrix, Q θ [= Σ −1 θ ], allowing for fast computational matrix inversion [14]. To find the appropriate GMRF structure, the Matérn [5] spatial component is parametrised as the finite element solution to the SPDE, is the Laplacian (for Cartesian coordinates s (1) , s (2) , s (3) ) and W(s) is a spatial white noise process [56].
To extend this spatial process to a spatio-temporal process, temporal innovations are modelled by first order autoregressive dynamics:

28/33
where |φ| < 1 is the autoregressive coefficient and ω(s i , t) is iid Normal. Practically the spatiotemporal process is achieved through a Kronecker product of the (sparse) spatial SPDE precision matrix and the (sparse) autoregressive temporal precision matrix.
In this paper for readability and consistency with standard notation equations 4 and 5 are parameterised through the covariance matrix. As specified above in our GMRF implementation we instead specified the precision (inverse covariance) matrix. Using the precision matrix the conditional predictive distribution takes the form Where Equation 18 A is introduced as a sparse observation matrix that maps the finite dimensional GMRF at locations (s • , t • ) to functional evaluations at any spatio-temporal locations e.g prediction locations (s , t ) or data locations (s, t), provided these locations are within a local domain.
1.3 S1 -Bias variance derivation for a Gaussian process stacked generaliser Theorem 1. Consider a function f from R N to R for which a sample D = {x i , y i } exists, where y i = f (x i ) and i = {1, ..., n}. Fit L level 0 generalisers, M 1 (x), .., M L (x), trained on data D. Next define two ensembles of the L models, the first using a weighted mean, , and the second as the mean of a Gaussian process,

29/33
Consider two ensembles,M , of the M 1 (x), .., M L (x) models With convex combination constraints β i ≥ 0∀i and L i=1 β i = 1. Define the squared error between the target function and each individual level zero generaliser Define the squared error between the target function and the ensemble as e(x) = (f (x) −M (x)) 2 . Following from [38] define the ambiguity of a given model As derived in [38] the ensemble ambiguity from using a constrained weighted mean ensemble (model 1 above) subject to convex combinations is defined as the weighted sum of the individual ambiguities.ā where¯ (x) = L i=1 β i i (x). Therefore the ensemble squared error when using a constrained weighted mean is e cwm (x) =¯ (x) −ā(x).
Using the Gaussian process (model 2 above) the ensemble ambiguity is defined as.

30/33
In equations 24 and 25 above the covariance matrices operate on the ambiguities, and as such are suffixed with asterixes to distinguish those from equation 21. Additionally β in equations 24 and 25 are also subject to convex combination constraints. Substituting 23 into equation 25 a(x) =¯ (x) − e(x) + Σ 2 * Σ −1 1 * e(x) The right hand side of equation 28 is identical to that derived using a constrained weighted mean in equation 23, but the left hand side error term, e gp (x) is augmented by (I − Σ 2 * Σ −1 1 * ). Clearlyā(x) ≤¯ (x) ∀x, withā(x) =¯ (x) only when the ensemble equals the true target function, M = f (x), and e(x) = 0. It follows that the left hand side of equation 28 I−Σ 2 * Σ 1 * e(x) ≥ 0 ∀x.
Therefore from the contribution of the precision or covariance stacking using a Gaussian process approach always has a lower error than stacking via a constrained weighted mean, with the error terms being equal when the contribution of the covariance is zero. That is (I − Σ 2 * Σ 1 * )e gp (x) ≤ e cwm (x) ∀x

Alternative stacking designs
In section we presented the rational for stacking and introduced a basic design (design 1 in figure 3) where multiple level 0 generalisers are stacked through a single level 1 generaliser.
For this design we proposed a Gaussian process or constrained weighted mean as the level 1 generaliser. In figure 3 we suggest two alternative stacking designs. In design 2, multiple level 0 generalisers are fitted and then passed through individual level 1 generalisers before being finally combined in a level 2 generaliser. An example of this scheme would be to fit multiple level 0 generalisers using different algorithmic methods (see Appendix 1.1) and then feed each of these into a Gaussian process regression. This design allows for the Gaussian processes to learn individual covariance structures for each level 0 algorithmic method (as opposed to a joint structure as in design 1). These level 1 Gaussian processes can then be combined through a constrained weighted mean level 2 generaliser.

31/33
In design 3, a single level 0 generaliser is used and fed into multiple level 1 generalisers before being combined via a level 2 generaliser. An example of this scheme would be to fit a single level 0 method, such as a linear mean or random forest, and then feed this single generaliser into multiple level 1 Gaussian processes. These multiple Gaussian processes can learn different aspects of the covariance structure such as long range, short range or seasonal interactions.
These level 1 generalisers can then be combined, as in design 2, through a constrained weighted mean level 2 generaliser.