Gradient boosting with extreme-value theory for wildfire prediction

This paper details the approach of the team Kohrrelation in the 2021 Extreme Value Analysis data challenge, dealing with the prediction of wildfire counts and sizes over the contiguous US. Our approach uses ideas from extreme-value theory in a machine learning context with theoretically justified loss functions for gradient boosting. We devise a spatial cross-validation scheme and show that in our setting it provides a better proxy for test set performance than naive cross-validation. The predictions are benchmarked against boosting approaches with different loss functions, and perform competitively in terms of the score criterion, finally placing second in the competition ranking.


Introduction
Wildfires occur in every season of the year and are a natural phenomenon of the forest ecosystem, important for clearing out decayed vegetation and helping plants to reproduce.However, they have the potential to become conflagrations-intense, destructive fires-that may have huge environmental and ecological impacts.Apart from human casualties, these fires can lead to substantial economic losses; global insured claims due to wildfire events have increased dramatically in recent years, from below $10 billion in 2000-2009 to $45 billion in the subsequent decade1 .
Wildfires are complex dynamic processes: their occurrences and behaviour are the product of interconnected factors that include the ignition source, fuel composition, topography and the weather.For example, the wind plays a big role in the spread and ease of fire containment, but its effect is magnified in the presence of accumulated biomass on a hilly boreal forest after a prolonged dry spell.The modelling of wildfires is made even more complicated by the need to model the Wildland-to-Urban interface (Stewart et al., 2007), as 90% of fires are caused by human activity.
An important measure of wildfire impact and size is the burned area of wildfire events, commonly used by government agencies and aggregated at different spatial and temporal scales for reporting purposes, e.g., National Interagency Fire Center (2021).Although there is a positive relationship between wildfire counts and sizes, with perfect agreement of the lowest value of zero in both variables, more wildfire occurrences in a region do not necessarily translate to higher burned areas.In 2020, nearly 26,000 wildfires burned approximately 9.5 million acres (ac) in the west, compared with the over 33,000 fires that burned just under 0.7 million ac in the east.Similarly, although the numbers of wildfires have fallen since the 1990s, the average annual acreage burned since 2000 has more than doubled.
Many statistical approaches have been developed to aid in wildfire prevention and risk mitigation, with most studies modelling wildfire occurrences and sizes separately (Taylor et al., 2013;Xi et al., 2019;Pereira and Turkman, 2019;Jain et al., 2020), though models that identify latent factors affecting both have been proposed (e.g., Koh et al., 2021).Point processes are natural models for the spatiotemporal pattern of occurrences (Peng et al., 2005;Genton et al., 2006;Tonini et al., 2017;Opitz et al., 2020).Cumming (2001), Cui and Perera (2008) and Pereira and Turkman (2019) suggested modelling fire sizes with various probability distributions.As data usually show heavy-tailed behaviour, only a small fraction of wildfires account for the vast majority of the area burned.Obvious candidates to capture this stem from extreme-value theory, such as the generalized Pareto distribution (GPD) for modelling threshold exceedances (De Zea Bermudez et al., 2009;Turkman et al., 2010;Pereira and Turkman, 2019).
Both Bayesian and frequentist methods have been used for explanatory modeling, the former predominantly for hierarchical mixed effect models (Joseph et al., 2019;Pimont et al., 2021;Koh et al., 2021) and the latter within the generalized additive modelling (GAM, Wood, 2017) framework (Preisler et al., 2004;Brillinger et al., 2006;Vilar et al., 2010;Woolford et al., 2011).Covariates include weather variables such as humidity, temperature, precipitation and meteorologically-based fire danger indices such as the Canadian Fire Weather Index (van Wagner, 1977).When available, land-use or locally observed an-thropogenic variables like population density and the distance to the nearest train line are used to help to explain human-induced occurrences; spatiotemporal random effects have been incorporated as surrogates for these variables.
If accurate prediction is of primary interest, then machine learning (ML) techniques offer an attractive alternative to the statistical modelling approaches described above.Since the 1990s, the surge in the availability of data and covariates has spurred the use of these techniques to predict wildfire behaviour.Jain et al. (2020) found 127 journal papers or conference proceedings published up to the end of 2019 on ML applied to fire occurrence, susceptibility and risk; of these adversarial neural networks (ANN) were the most prominent (e.g., Dutta et al., 2013;Shidik and Mustofa, 2014;Liang et al., 2019).For wildfire occurrences, most studies focus on classification tasks instead of count modelling.Sakr et al. (2010) used metereological variables with support vector machines to predict a four-class fire risk index based on the daily number of fires in Lebanon.Dutta et al. (2013) compared ten ANN based cognitive imaging systems to determine the relationship between monthly fire incidence and climate for Australia.Mitsopoulos and Mallinis (2017) and Xie and Peng (2019) showed that ensemble learning methods like random forests and boosting trees performed well in estimating area burned or classifying wildfire sizes in Greece and Portugal, respectively.
Gradient boosting techniques (Friedman, 2001) have exploded in popularity over the last decade, in part due to the development and dissemination of open-source packages such as gbm (Greenwell et al., 2020) and xgboost (Chen and Guestrin, 2016).A key ingredient of gradient boosting is the loss function used to train these models, and choices for these functions have largely been restricted to those that emphasize good prediction of the distributional bulk instead of the tails.For example, squared loss, the default when modelling wildfire sizes, implicitly presupposes normality of the response given the covariates, which may be inappropriate if the focus is predominantly on extreme values.The Poisson loss is popular for modelling wildfire counts within a grid cell, but the zero-inflated nature and potentially heavy tails of count distributions suggest that this loss may be unsuitable.Evaluation metrics should also reflect the non-linear impact of wildfire events; e.g., in many cases, predicting a false negative occurrence is much costlier than predicting a false positive.
As ML methods are prone to overfitting, it is imperative to evaluate models with held-out datasets using robust validation schemes.A realistic approach in the forecasting context (when there are no trends) is to leave out the most recent portion of the dataset (e.g., Woolford et al., 2011;Dutta et al., 2013;Joseph et al., 2019;Koh et al., 2021).Dutta et al. (2013) explored different combinations of training-testing splits to identify the best possible paradigm to maximize the generalization capability of their ANN architecture.K-fold cross-validation is also popular (Shidik and Mustofa, 2014;De Angelis et al., 2015;Mitsopoulos and Mallinis, 2017;Xie and Peng, 2019), but it may give overly optimistic evaluations for spatially dependent data (Roberts et al., 2017).An alternative is spatial cross-validation (Pohjankukka et al., 2017), but it is still unclear how best to construct spatial folds in this context, and doing so anyway ignores relevant inter-variable or time dependencies.
Our work aims to tackle the limitations of the studies mentioned above, and does so in the context of the Extreme Value Analysis 2021 data challenge (Opitz, 2021).We develop novel gradient boosting models trained with loss functions appropriate for predicting extreme values.Our chosen model for fire counts is a discrete generalized Pareto distribution (dGPD, Shimura, 2012) relying on a covariate-dependent parameter that models a chosen high quantile of the distribution, and a shape hyperparameter selected by cross-validation.The dGPD provides added flexibility in modelling the upper tail of the count distribution, especially when compared to other more standard count distributions like the Poisson.The model for fire sizes has three components and covariate-dependent probabilities.The first component models the probability of observing no fires, and the others model the probabilities of observing medium-sized and extreme fires.We approximate the conditional distribution above a high threshold with a GPD, and the conditional distribution below the threshold with a truncated log-gamma distribution.
To improve our models, we also engineered new covariates that incorporate more spatial information into the climatic and land-use covariates provided by averaging them across neighbouring grid cells each month.With a smart imputation method for replacing missing data, we also use the wildfire counts as a covariate when predicting wildfire sizes, and vice versa.
We develop a spatiotemporal cross-validation scheme that provides a better proxy for our models' test set performance than the naive scheme.This involves fitting a space-time latent Gaussian model to pseudo-binary observations that indicate whether a grid cell was masked in a particular month, and then simulating from the fitted model to generate folds of train-test regimes.
In the remainder of the paper, we first explore the data on wildfires and their covariates, and then introduce the problem set out by the data challenge in §2.We provide general background on extreme-value theory and gradient boosting and on how to combine them in §3.Our spatiotemporal cross-validation scheme is developed in §3.4 and the specific model structure is detailed in §4.We highlight the prediction of wildfire activity components in §4.2, and compare them to related and competing approaches.We conclude with a discussion and outlook in §5.

Data and exploratory analyses
The Extreme Value Analysis 2021 data challenge dealt with the prediction of monthly wildfire counts and burned areas at 3503 grid cells across the contiguous US over the period 1993-2015.As fuel moisture is an integral of past precipitation and evaporation mediated by soil field capacity, temporal scales longer than hourly or daily (e.g., monthly in our case) are appropriate for predicting fire risk from climatic covariates.The data comprise the monthly numbers of wildfires (CNT) and the aggregated burned area (BA) in each grid cell based on a 0.5 • × 0.5 • grid of longitude and latitude coordinates (roughly 55km ×55km) covering the study area, from March to September each year.Figure 1 shows that the grid cells with the highest averaged CNT tend to be clustered towards the west (California) and southeast (North and South Carolina) corners of the study region, while clusters in the west (near the border of Idaho and Nevada), southwest (Arizona, New Mexico and Texas), and southeast (Florida) are observed for BA.

Gulf of Mexico
Thirty-five auxiliary variables related to land cover, weather and altitude are provided at the same spatial and temporal resolution, and can be used for  modelling.Figure 2 hints at the importance of some of these variables, such as meteorological covariate 5 (clim5; potential evaporation, measured in meter water equivalent -mwe) and land cover covariate 7 (lc7; tree needleleave evergreen closed to open, measured in % of the grid cell), for predicting high BA.However it also indicates that their effect differs over space, and that the interactions between these covariates may be important for predicting large burned areas.The Rocky Mountain Area and Great Basin, two regions defined formally as 'Geographic Area Coordination Centers' by the United States Department of Agriculture, have empirical exceedance probabilities that respond differently to both covariates.For instance, the probability tends to decrease, then increase (after the 75% quantile) with potential evaporation in the Great Basin, while the opposite holds for the Rocky Mountain Area, though the associated large uncertainties suggest that there is substantial heterogeneity within each region.
The original dataset has no gaps, but missing data were artificially created to compare predictive approaches; the full dataset was split into training and testing subsets to evaluate participants' test scores.No data were masked in the odd years, but 80,000 observations of each variable were masked in the even years.Figure 3 shows that the spatial and temporal positions of test data are clustered in space, and the test grid cells for BA and CNT are correlated.This masking is reminiscent of a real-world situation in which two related and spatially dependent processes could render both CNT and BA unobservable at small spatial clusters every month (e.g., from a potentially spatially dependent measurement error), and one could only use the available covariates and responses from the surrounding non-masked regions for prediction.For example, datasets generated using satellite-based remote sensing of wildfires have known misclassification issues, very often due to cloud occlusion.
The evaluation metrics used for the competition (see Opitz, 2021) require estimates of the probability Pr(BA < u BA ) and Pr(CNT < u CNT ) for 28 thresholds u CNT and u BA .The metrics are variants of weighted ranked probability scores and put relatively strong weight on good prediction in the extremes of the distributions of counts and burned areas.As such, we expect that models that emphasize accurate modelling of the largest counts and burned areas will perform better, and we achieve this by appealing to extreme-value theory.

Extreme-value theory
The generalized Pareto distribution (GPD) arises asymptotically for excesses above a large threshold of a random variable X ∼ F , when the distribution F satisfies mild regularity conditions.Let x = sup{x : F (x) < 1}.The excess distribution above u < x can be approximated (Davison and Smith, 1990) as (1) with shape parameter ξ ∈ R and scale parameter σ = σ(u) > 0, where a + = max(a, 0).The shape parameter determines the rate of tail decay, with slow power-law decay for ξ > 0, exponential decay for ξ = 0, and polynomial decay towards a finite x , for ξ < 0. When the approximation (1) is exact asymptotically (i.e., when u → x ), we say that the random variable X lies in the maximum domain of attraction of a generalized Pareto distribution with shape parameter ξ, written X ∈ MDA ξ .

Positive discrete random variables
We say that a discrete non-negative random variable Y lies in the discrete maximum domain of attraction, Y ∈ dMDA ξ , if there exists a random variable X ∈ MDA ξ with ξ ≥ 0 such that Pr(Y ≥ k) = Pr(X ≥ k), for k = 0, 1, . . . .The random variable X is called the version of Y , and many popular discrete distributions such as the geometric, Poisson and negative binomial distributions lie in dMDA ξ .
For Y ∈ dMDA ξ and large integers u, where the last term is the probability mass function of the discrete generalized Pareto distribution (Shimura, 2012).Several studies have used this distribution to model count data; Prieto et al. ( 2014) modelled numbers of road accidents and Hitz et al. ( 2017) modelled numbers of extreme tornadoes per outbreak and multiple births.

Gradient tree boosting
The generic gradient boosting estimator (Bühlmann and Hothorn, 2007) is a sum of base procedure estimates.Regression trees (Breiman et al., 1984, CART) are popular base procedures, as they include non-linear covariate interactions by construction, and are invariant under monotone transformations of covariates, so the user need not search for good data transformations.
. ., n, be a dataset with n observations and p covariates.A binary split of the covariate space uses a splitting variable indexed by j ∈ {1, . . ., p} and a split point v ∈ R to partition the space into the pair of half-spaces {x ∈ R p : x j ≤ v} and {x ∈ R p : x j > v}, where x j is the j-th component of x.
By successive binary splits, a regression tree partitions the covariate space into a set of L disjoint regions A 1 , . . ., A L , and fits a simple model such as a constant in each region.The regions created by the splits are called nodes; a terminal node is called a leaf and an interior node is called a branch.We index leaves by l ∈ {1, . . ., L}, with leaf l representing region A l .The simplest tree is one with two leaves, known as a stump.A learning algorithm needs to decide the tree structure, i.e., the splitting variables and split points.
Suppose that L leaves with regions A 1 , . . ., A L have been chosen and we model the response as a score c l ∈ R in each region.A regression tree is a function where I is the indicator function.A gradient tree boosting model uses T such trees to model the boosting estimate where T is the space of regression trees.In this paper, θi will always represent a parameter estimate in a given model for the conditional distribution of y i given x i .The boosting estimate using stumps will be additive in the original covariates, because every base estimate is a function of a single covariate.A boosting model that has trees with at most L leaves has interactions of order at most L − 2. Thus, constraining the maximum number of nodes in the base procedure controls the complexity of the model.Gradient tree boosting learns the set of trees used in (3) by minimizing a regularized objective function in a greedy iterative fashion; at each iteration we add the tree that most improves our model according to an objective function O.More precisely, let θ(j) i be the boosting estimate for the i-th observation at the j-th iteration.We add a tree f j at each iteration to minimize where is the number of leaves in the tree f j and c ∈ R L (j) is the corresponding vector of scores.The regularization term Ω is added to penalize the complexity of each tree, and the positive parameters η and λ control the penalization.The form of Ω is simple enough to allow parallel computation (Chen and Guestrin, 2016).Using a second-order approximation for the objective (Friedman et al., 2000) gives where We then minimize (4) with respect to the tree structure and weight vector c.
Let I l = {i : x i ∈ A l } denote the instance set of leaf l.For a fixed tree structure with regions A 1 , . . ., A L (j) , the optimal weights c can easily be found and have components Plugging the weights c into (4) and removing the term that does not depend on f j gives Õ(j) = − 1 2 which can be used as a scoring function to measure the quality of the tree structure, a role similar to the impurity score in Breiman et al. (1984).
Assume that a split has been performed, and let I L and I R denote the instance sets of the left and right leaves from this split.Define I = I L ∪ I R .The loss reduction from this split is and ( 6) is used for evaluating candidate split variables and points.
As it is impossible to enumerate all possible tree structures, most existing tree boosting implementations, such as in scikit-learn (Pedregosa et al., 2011) and gbm (Greenwell et al., 2020), use greedy algorithms that start from a single leaf and iteratively add branches to the tree based on (6), until the gain for the best split is negative.Here we use the greedy algorithm implemented in xgboost (called the approximate algorithm with weighted quantile sketch Chen and Guestrin, 2016, Appendix A), that further reduces computational cost and parallelizes computations when the data do not fit into memory.This algorithm proposes candidate splitting points from the empirical quantiles of each covariate, instead of considering all possible splitting points for each variable.We also subsample a proportion s ∈ [0, 1] of the covariates at each iteration, like in the random forest algorithm (Breiman, 2001), to further prevent overfitting and to accelerate parallel computation.For the full algorithmic details, we refer the reader to Chen and Guestrin (2016).
Convexity of L is a desirable feature that would guarantee that a unique global minimum exists; if this does not hold in practice and one is concerned about the algorithm being stuck at local minimums, then a potential solution is to rerun the algorithm multiple times from different initial boosting estimates θ(0) i , i = 1, . . ., n, and select the run that gives the lowest loss.If L is convex and there are enough iterations, then the choice of the initial boosting estimate will have minimal effect on the overall prediction.A natural choice is to set θ(0 , for all i = j, where θ = arg min θ n i=1 L(y i , θ), i.e., this is the common parameter value that minimises the loss across all observations, unconditionally on the covariates.
As gradient tree boosting is an ensemble method combining predictions from many trees, the exact relationship between covariates (or their interactions) and the response is difficult to determine.Suppose we treat the covariates as random, and let X S denote the random subvector (of size l < p) of the full covariate vector X = (X 1 , . . ., X p ) T , indexed by the set S ⊂ {1, . . ., p}.Let C be the set complementary to S. For a predictive function f at a fixed point x ∈ R l , the partial dependence function (Friedman, 2001), where x 1C , . . ., x nC are realisations of X C from the training data.Metrics used to rank covariates in terms of their importance include the coverage for a chosen covariate, which is the sum of the second order gradients h i from (5), in each node which uses this covariate, standardized by dividing by the sum of the metrics for all other covariates (so the resulting metric is a proportion).The gain metric represents the fractional contribution of a chosen covariate to the model based on the total gain of all the splits involving this covariate, measured by G in (6); it is the total improvement of the model in terms of the objective function, from the branches that include the covariate.In both cases a higher proportion implies a more important predictive variable.
The loss function L in (4) strongly governs the type of models that we fit, and we discuss choices for this next.

For wildfire counts
The most common statistical model for count data is the Poisson distribution, and many existing boosting implementations use a scaled version of the simplified Poisson loss as their default loss function for counts.Given n monthly wildfire counts in a grid cell, y 1 , . . ., y n , this loss is where Stirling's approximation log(y i !) ≈ y i log(y i ) − y i is used and the boosting estimate θi models the log mean of the i-th Poisson count.The terms that do not depend on θi can be dropped when minimising this loss in an optimisation procedure.Although ( 8) is also the most common choice for modelling wildfire counts in the literature, the zero-inflation and potentially heavy tails of count distributions suggest that it may be unsuitable.
Instead, if we let α = 1/ξ > 0 and λ = σα in (2), we motivate a new loss function for counts from extreme-value theory, the discrete generalized Pareto (dGPD) loss The boosting estimate θi models the logarithm of the rescaled scale parameter and otherwise the mean does not exist.
The first and second derivatives of ( 9) with respect to the boosting estimate, i.e., are used in the second-order approximation of the objective in (4) and are essential for determining split variable and split point candidates when building trees with (6).

For wildfire sizes
The squared loss is the dominant choice in the literature on predicting burned areas in this regression context, but it implicitly presupposes normality of the response conditional on the covariates, which may be inappropriate if the distributional tail decays more slowly than exponential.Moreover, burned areas cannot be negative.Modelling log-transformed burned areas addresses the latter issue, though doing so still excludes conditional distributions with Pareto-like tails.
Another approach to modelling the full distribution of wildfire sizes is to use a mixture, by first choosing an appropriately high threshold and then fitting the burned areas below and above that threshold with different loss functions.To additionally handle the zero-inflation of wildfire sizes, we can left truncate the sizes below the high threshold at zero.This mixture approach models burned areas by splitting the distribution into three groups: zero, intermediate and extreme sizes.
To model the monthly burned area in a grid cell, y i , below a chosen threshold u > 0, the negative log-loss likelihood of a truncated distribution could be used.The probability density function of a right truncated gamma distribution is where u > 0 is the right truncation, µ > 0 is the rescaled scale parameter, k > 0 is the shape parameter and γ(k, s) = s 0 t k−1 exp(−t)dt, s > 0, is the lower incomplete gamma function.Modelling log(µ) with the boosting estimate and dropping the terms that do not depend on θi gives L trGamma (y i , θi ) =k θi + y i k/ exp( θi ) + log γ{k, ku/ exp( θi )}, + γ {k, ku/ exp( θi )}/γ{k, ku/ exp( θi )}}, where γ {k, ku/ exp( θi )} and γ {k, ku/ exp( θi )} are the first and second derivatives of γ{k, ku/ exp( θi )} with respect to θi , and have closed forms (see Supplement §6).
To model only the excesses above a threshold u, we can use the GPD in (1), and assume that ξ > 0, since burned areas tend to be heavy-tailed.If we reparameterize and model the logarithmic κ ∈ [0, 1] quantile of the excesses with the boosting estimate, i.e., θi = log[{(1 − κ) −ξ − 1}σ i /ξ], we can define the generalized Pareto (GPD) loss and obtain its derivatives where f and f are the first and second derivatives of the reparameterized probability density function f {y i , exp( θi ), ξ} given in Supplement §6.

For wildfire size classification
Adopting the mixture modelling approach to wildfire sizes requires us to model the probability that a fire belongs to each size component 1, . . ., C, where C is the number of components.Let y i = (y i,1 , . . ., y i,C ) denote the vector of wildfire size component indicators, where y i,c = 1 if the i-th fire size is in component c, and otherwise y i,c = 0. We can model the probability of each class with the boosting estimate using the softmax function σ : The generalization of the logistic (Cox, 1958) loss to multiple classes is the cross-entropy loss, which can be reweighted to give where the vector of component probabilities is modelled by applying (11) to the boosting estimate θi = ( θi,1 , . . ., θi,C ) T , and the weights w 1 , . . ., w n could be chosen to improve predictions in unbalanced classification tasks.

A spatiotemporal cross-validation scheme
The use of k-fold cross-validation generally presupposes independent replicates, so it would produce optimistic predictive performance estimates in our setting because data points that are geographically closer will have stronger dependencies.To address this, we first study the spatiotemporal process leading to grid cells being masked, which we call the masking process.Figure 3 hints at either a common or two inter-correlated spatially dependent latent processes governing the observed masking processes for CNT and BA, which we model with a common latent Gaussian process (Rasmussen and Williams, 2005).We then fit a Bernoulli model to observations arising from the masking process, and simulate observations from the model to generate cross-validation folds.
We consider only the months m with masked observations, and let M denote the number of those months and D denote the number of grid cells.Let R CNT d,m and R BA d,m denote the binary 0-1 observations indicating whether the grid cell d ∈ {1, . . ., D} (with centroid s d ) at month m ∈ {1, . . ., M } was masked for the CNT and BA responses, respectively.Our hierarchical model for the masking processes is where expit(x) = {1 + exp(−x)} −1 is the inverse logit function.
We fit this model using the integrated nested Laplace approximation (INLA, Rue et al., 2009), which is an approximate Bayesian inference technique wellsuited for latent Gaussian models.The parameter β governs the degree of latent sharing between the two masking processes and we use a flat and independent zero-centered Gaussian hyperprior for it.Similar frameworks were used by Koh et al. (2021) for the joint modelling of different wildfire risk components and by Diggle et al. (2010) and Pati et al. (2011) to model preferential sampling.The spatial process g t is independently replicated in time and each replicate has a Gaussian process prior GP with a Matérn covariance structure governed by the parameter vector ζ.We represent these Gaussian processes via a numerically convenient Gauss-Markov random field approximation, constructed by solving a stochastic partial differential equation (Lindgren et al., 2011).Supplement §6 details the full procedure.
We generate samples from this Bayesian model by first sampling parameters from the posterior distribution, and then generating observations from the Bernoulli model with the sampled parameters.We do this for all months, including in those where observations were masked; if a location was already part of the test set, i.e., if it was already masked, then we removed it from the validation set.Five samples were generated to obtain five folds for our crossvalidation scheme.Figure 4 shows two samples from this model for March 1993.The degree of spatial and inter-variable dependencies resemble those of the masking processes in Figure 3, and the numbers of grid cells masked and chosen for validation in each month are also similar.The triplet (2.5% posterior quantile, posterior mean, 97.5% posterior quantile) for the scaling parameter β is (0.28, 0.42, 0.58).

Fitting procedure
We fit our gradient tree boosting models with the approach outlined in §3.2, minimizing the loss functions described in §3.3.
We use the Poisson and dGPD loss functions to fit models on the full CNT distribution.We experimented with different high thresholds u but achieved the best prediction when we modelled the full distribution with the dGPD, i.e. u = 0.For all observations with a masked count response but an unmasked zero valued burnt area, we artificially set their corresponding counts to zero, and use these observations to fit the CNT models.
For wildfire sizes, we fit models with the squared loss on log(1+BA i ), i = 1, . . ., n, and call the best fitted model from this category 'log-Normal'.We also consider mixture models which require split modelling of the distribution.For these, we first choose a sufficiently high (95% empirical quantile) threshold at 200ac.We then use the fire sizes above the threshold to fit a model with the GPD loss, and the log-transformed positive sizes below the threshold with the truncated gamma loss.Lastly, we fit a multi-class classifier to the wildfire size component indicators y i = (y i,1 , . . ., y i,3 ), defined in §3.3.3, using the crossentropy loss from (12); here, y i,1 = 1 if we observe no fire, y i,2 = 1 if BA i is a medium fire (between 0 to 200ac), and For the classifier, we also used all the observations with a masked burned area but an unmasked zero valued count, i.e., we artificially set the corresponding burned area to be zero for these observations.Given the covariates at the i-th observation x i , we combine the three model components to get the prediction of the cumulative conditional probability for each observation We also engineered new covariates to improve our predictions.To incorporate more spatial information from our covariates (other than the longitude and latitude coordinates), we took the average value of the covariate across neighbouring grid cells for each month; this smooths the climatic variables across space.We also allow land-use covariates at neighbouring grid cells to help predictions at each grid cell.
The relationship between CNT and BA is positive, though a high CNT does not imply a high BA.Nevertheless, it is still natural to consider the other response as a covariate when modelling a given response, or at least to use this information whenever possible.As the test grid cells for BA and CNT are correlated, there are instances where the BA response was masked but the CNT wasn't, and vice versa; for 39% of masked BA observations, their corresponding CNT observations were unmasked.Using the CNT/BA covariate to predict BA/CNT thus raises the question of how best to impute its value if it was masked for a given observation.The default way to handle a missing covariate is to impute its mean across all observations, such as in algorithms from xgboost or gbm, though this will be sub-optimal for predictions on a spatially heterogenous dataset.Instead, we use an imputation method which fits a model for the covariate and then imputes the best estimate from this model.
When modelling BA, we first fit a gradient boosting model with the dGPD loss on the CNT response, without using BA as a covariate so as to prevent data leakage, and then use the estimated parameters to find the estimated mean CNT for each observation using (10); we then impute this estimate whenever CNT was masked.When modelling CNT, we first fit a gradient boosting model, without using CNT as a covariate, with the cross-entropy loss on the wildfire size component indicators y i , i = 1, . . ., n, and then impute the estimates of the probabilities Pr(y i,1 = 1), Pr(y i,2 = 1) and Pr(y i,3 = 1) from the fitted model, whenever BA was masked; when BA wasn't masked, we use the observed indicators y i as a covariate.
To assess models using these covariates, a cross-validation scheme should also reflect this imputation procedure; thus, it becomes even more important to appropriately model the inter-variable dependence between CNT and BA of the masking processes, i.e., the parameter β in our spatiotemporal crossvalidation scheme described in §3.4.
Our models have hyperparameters that must be tuned by cross-validation.They include the regularization parameters λ and η in (5), the proportion s of covariates subsampled at each iteration, the maximum number of leaves for each tree L, and the number of trees T (see §3.2).Other hyperparameters from the loss functions include k, ξ, and α, which govern the shape and tails of the fitted conditional distribution, and weights w i (i = 1, . . ., n) which govern the importance of each observation in the cross entropy loss.Some of our models assume common shape parameters ξ and α governing the tails of wildfire sizes and counts across the whole sample space; our preliminary analysis suggests that this assumption is reasonable in space, as the frequentist estimates of the shape parameters from pooled data are relatively homogeneous across the wildfire coordination regions.
We use the Bayesian optimization procedure outlined in Snoek et al. (2012) to choose the parameters (excluding the number of trees M ) minimizing the average evaluation metric, calculated on the five cross-validation folds generated in §3.4.This procedure treats the objective function as random and first places a Gaussian process prior on it.After gathering function evaluations, the prior is updated to form the posterior distribution over the objective function, which is then used to construct an infill sampling criterion.For the mixture model, we implement separate Bayesian optimization procedures for each of the three model components.Given the other parameters, we choose M with the one-standard-error rule (Hastie et al., 2009); i.e., we select the largest M within one standard error of the parameter that achieves the minimum in terms of the evaluation metric.Figure 5 shows the evaluation metric on the validation folds as a function of M for the wildfire size classifier in a mixture model.
We fitted our models by combining our own R routines with the xgboost package.We implemented the described Bayesian optimisation procedure with the rBayesianOptimization package.

Results
The benchmark models are described in Opitz (2021).That for CNT corresponds to a generalized linear model (GLM, ?) with a Poisson response distribution and log link linear predictor using all the original covariates.The benchmark for BA first fits a generalized linear model with Gaussian response and log-link using all of the original covariates.The probability predictions, Pr(BA i ≤ u BA ), are obtained by combining the log-Gaussian BA model with the estimated probability that CNT i =0, obtained from the benchmark Poisson model for CNT.
We relied on our cross-validation scheme devised in §3.4,using the evaluation metrics used from the competition outlined in §2, to choose which model predictions to use for the competition.After the competition, we had access to the truth and could calculate how the predictions of every model would have performed on the test set.
Table 1 shows that incorporating the engineered CNT and BA covariates improves the scores of all models by up to 7%.According to our cross-validation Table 1: The averaged rescaled evaluation score for all models, according to the data-driven cross-validation (CV) scheme outlined in §3.4, the 5-fold CV scheme with random partitioning, and the true score.The bold figures highlight the best model chosen by our cross-validation approach.The asterisks indicate models not using the CNT/BA covariate.
scheme, the best model for wildfire sizes is the mixture model, and for the counts it is the dGPD model.The best mixture model and dGPD models from the Bayesian optimization procedure have α = 52 and ξ = 0.8.This implies a fat tail for the size distribution, but a thinner (Gumbel-like) tail for the wildfire count distribution, though the parameter α in the dGPD loss provides additional flexibility to the model that gives slightly better predictions than the Poisson model.All our gradient boosting models outperform the benchmark by around 10-50%.
Our cross-validation scheme tends to perform better than the 5-fold crossvalidation scheme as a proxy for the true test set performance; the scores from the 5-fold cross-validation scheme are generally too optimistic compared to the true test error, especially for the wildfire size models.This optimism is especially pronounced when evaluating models that use the engineered CNT and BA covariates.Our scheme is better able to capture the inter-variable dependence between the CNT and BA masking processes, giving a better reflection of how the models that incorporate the engineered covariates would perform when predicting responses on the test set.
Figure 6 shows that the gain and coverage metrics introduced in §3.2 give similar orderings for the importance of covariates when predicting the probability of being in a given wildfire size component with the best mixture model.As hypothesized in §2, clim5, lc7 and the spatial covariates of longitude (lon) and latitude (lat) are among the five most important variables for both metrics.The other variables (e.g., clim4, clim7, clim9, year, lc16, etc.) are relevant, but each is less than half as important as clim5, the most important covariate.
To evaluate the marginal effect of clim5 on the CNT response in the best dGPD model, we transformed the partial dependence estimate ( 7) by ( 10) to get the predicted mean count mi , and evaluated and plotted it with clim5 in the set of interest S in (7).As it is computationally infeasible to evaluate all data points x iC in our setting with over 500, 000 observations, we subsampled 10, 000 observations to obtain our Monte Carlo estimates.
Figure 7 shows that the marginal effect of clim5 on mi tends to be negative, especially above −0.03mwe.Figure 8 displays the joint marginal effects of clim5 and land cover covariate 12 (lc12; grassland, in %) on the predicted mean count, i.e., with clim5 and lc12 in the set of interest S in ( 7).The figure hints at interaction between the two covariates; increasing lc12 tends to decrease the response CNT slightly if clim5 is low, but not when clim5 is high.Although the partial dependence plot is useful for showing the overall marginal trend of a covariate on the response across all observations considered, it is important to be honest about the uncertainty associated with the Monte Carlo estimate in (7).The estimates in Figures 7 and 8 are associated with high uncertainty throughout (not shown for the latter); our dataset is very heterogeneous and it is not possible to quantify the marginal effects of covariates with less uncertainty.
Our chosen models perform competitively when compared to the other teams' submission in the data challenge (Opitz, 2021), placing second out of 28 teams in the final ranking; the other top three teams used other popular prediction techniques such as random forests, hierarchical Bayesian modelling and ANN models with adapted loss functions (Opitz, 2022).

Discussion
We have implemented novel gradient boosting models for wildfire activity that are trained with loss functions motivated by extreme-value theory.Compared to models trained on the Poisson loss, our chosen model for wildfire counts has an additional parameter α that governs the tail of the count distribution, which, after tuning by cross-validation, enables the model to give better predictions.Our chosen model for burned areas has specific components for extreme fire sizes.According to the given score criteria which put more weight on large fire sizes, this model improves on the models that do not discriminate between extreme and non-extreme fire sizes.
As the use of other data sources, e.g., covariates not included in the provided dataset, was strictly prohibited for the competition, we were not able to leverage other spatial information, such as the Geographic Area Coordination Center of a grid cell.Each center follows its own governing jurisdictions that could affect its fire mitigation and suppression strategies.Including this information, either by incorporating additional covariates, or by having a sep-arate gradient boosting model for each coordination center, would improve predictions.
Our mixture model has the threshold fixed at u = 200ac.However, u could have also been allowed to be an additional parameter chosen by crossvalidation, or as a spatio-temporally varying threshold u d,m , d ∈ {1, . . ., D}, m ∈ {1, . . ., M }, that could be either independent (Opitz et al., 2018) or dependent (Velthoen et al., 2021) on covariates.However, implementing these approaches would significantly raise the computational cost here, especially if the model for the threshold also has hyperparameters that need to be tuned, as the optimization of all of our model components would need to be done jointly because they all rely on the threshold.Future work could investigate whether these potential approaches that increase complexity do indeed improve wildfire prediction.
Our best chosen dGPD model from the Bayesian optimization procedure has α = 52, which implies a thin tail for the CNT distribution that is not too different from the Poisson distribution.This explains the slight improvement of the model using the dGPD compared to the Poisson loss.Had the tail of the CNT distribution been fatter, as in the case of another data application (e.g., insurance claim counts), we would have also noticed a larger improvement in predictions by the dGPD model.
We have implemented a spatial cross-validation scheme for our context which partly fixes the optimism when using traditional k-fold cross-validation to evaluate complex models with engineered covariates over a spatially heterogeneous but dependent dataset.One should always have the real-world prediction scenario in mind when choosing a cross-validation scheme, and we appealed to tools from spatial statistics to aid the validation of model predictions on our test data.
Apart from cross-validation approaches, model comparison using other predictive scores, e.g., Continuous Ranked Probability Scores (CRPS, Matheson and Winkler, 1976) or tail-weighted CRPS (Gneiting and Ranjan, 2011), could be used to compare simulated predictive distributions of burned areas or counts.These scores, along with graphical summaries for validating models, are an important future topic of research, and have already been explored in the extreme wildfire prediction context by Joseph et al. (2019), Pimont et al. (2021) and Koh et al. (2021).Due to the time constraints of the competition however, this is out of the scope of this paper.
Our gradient boosting models have hyperparameters from the loss functions that govern the tails of the predictive distributions: ξ, and α, which, once chosen by cross-validation, are fixed in the model.To allow more flexible modelling of the distributional tails, one could incorporate them as an additional boosting estimate in §3.2 which would allow these parameters to depend on the covariates.The boosting estimate θi , gradient g i and hessian h i in (4) would then be two-dimensional vectors.The described approach is similar to the recent work by Velthoen et al. (2021).Another avenue for future work is to assess how the the tail indices affect the cross-validation score when the other hyperparameters are fixed.
Apart from better predictions, our models improve decision support in wildfire management.The partial dependence plots in Figures 7 and 8 allow marginal and interaction effects of covariates to be assessed, though one should be aware of the large uncertainty associated with these estimates.Importance metrics like the gain and coverage in Figure 6 could be used for covariate selection and could prompt national wildfire predictive services to rethink designs of fire danger warning systems (e.g., indices) across the contiguous US. + γ {k, ku/ exp( θi )}.

Priors and SPDE triangulation
This section details the prior and SPDE triangulation specifications of the spatiotemporal Bernoulli model for the masking processes.The fixed effect coefficient β CNT 0 and β BA 0 in our model was assigned a flat Gaussian prior with zero mean and precision 0.001.The prior for the scaling parameter β is a zero-centered Gaussian distribution with precision ω = 1/20.Lastly, we assigned a log-gamma hyperprior with mean unity and precision 0.0005 to the variance hyperparameter φ.
The parameter vector ζ consists of r and σ parameters which were assigned penalized complexity priors (Fuglstad et al., 2018).Such priors penalize the distance of the prior of a model component towards a simpler baseline at a constant rate (Simpson et al., 2017).
The discretization points when triangulating the spatial domain are chosen as the nodes of a finite element representation which enables efficient inference for random effects representing spatial variation.Our spatial triangulation mesh in Figure 9 has 508 nodes.It is sparser in the extended zone around the study area to ensure that the SPDE boundary conditions have negligible influence on the study area.
For computational reasons, we only use observations of the masking processes from the first ten months, and subsample 30% of the data when fitting our model.

Fig. 1 :
Fig. 1: Maps of CNT (top) and BA (bottom) averaged across all months and years for each grid cell.

Fig. 2 :
Fig. 2: Empirical frequency Pr(BA > 200ac) with 95% error bars, as a function of the potential evaporation (in m, left) and tree needleleave evergreen closed to open (in %, right) covariates, grouped by observations within four empirical quantile ranges, for grid cells within the Rocky Mountain Area (blue) and Great Basin (red) coordination centers.The two regions are highlighted in blue and red in the map inserted in the right panel.

Fig. 3 :
Fig. 3: The test set grid cells (in red) for CNT (left) and BA (right) in March 1994.The number at the top right indicates the number of masked grid cells.

Fig. 4 :
Fig. 4: The first (top) and second (bottom) validation folds (in red) for CNT (left) and BA (right) for March, 1993, from our spatiotemporal cross-validation scheme, generated from the Bayesian spatial model.The number at the top right indicates the sum of grid cells chosen.

Fig. 5 :
Fig. 5: The average rescaled evaluation score across all five folds of our spatiotemporal cross-validation scheme for a mixture model, as a function of the number of trees M .The shaded region shows the pointwise one-standard-error bound.The red point shows the minimum average validation score and the blue point shows the M chosen by the one-standard-error rule.

Fig. 6 :
Fig.6: The coverage (left) and gain (right) metrics for the top 20 covariates in the wildfire size classifier submodel of the best mixture model for BA (without using the engineered covariates).More information about the covariates is given inOpitz (2021).

Fig. 7 :Fig. 8 :
Fig.7: The partial dependence plot with clim5 in the set of interest S, transformed so the y-axis shows the predicted mean CNT from (10), for the observations in the Rocky Mountain Area (left) and Great Basin (right) regions.The empirical 95% and 5% pointwise quantiles from the Monte Carlo estimates are shown by the dotted lines.

Fig. 9 :
Fig. 9: Triangulation mesh of the spatial study region (blue contours) for the SPDE approach.Neumann boundary conditions are set on the exterior (black) boundary to obtain a unique solution.The blue points represent the centroids of the grid cells.The finite element solution defines a Gauss-Markov random vector with one variable in each node.