{'fraction_non_alphanumeric': 0.03121988822509154, 'fraction_numerical': 0.016766236269030642, 'mean_word_length': 4.686684931506849, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 43, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Data-driven chimney fire risk prediction using machine learning and point process tools

Chimney fires constitute one of the most commonly occurring fire types. Precise prediction and prompt prevention are crucial in reducing the harm they cause. In this paper, we develop a combined machine learning and statistical modeling process to predict chimney fires. Firstly, we use random forests and permutation importance techniques to identify the most informative explanatory variables. Secondly, we design a Poisson point process model and apply associated logistic regression estimation to estimate the parameters. Moreover, we validate the Poisson model assumption using second-order summary statistics and residuals. We implement the modeling process on data collected by the Twente Fire Brigade and obtain plausible predictions. Compared to similar studies, our approach has two advantages: i) with random forests, we can select explanatory variables non-parametrically considering variable dependence; ii) using logistic regression estimation, we can fit the statistical model efficiently by tuning it to focus on important regions and times of the fire data.


Introduction
During the last decade, the Dutch fire and rescue services are developing an interest in applying business intelligence to improve their strategy of fire prediction and prevention [28]. In order to prepare for risk reducing measures, such as essential public awareness campaigns and proper fire staffing and equipment arrangements, accurate predictions are required. In this paper, we focus on chimney fires, as they occur frequently, rely heavily on environmental factors and impact people's daily life. Collaborating with the Twente Fire Brigade, we conduct a complete risk prediction study for chimney fires. We define the risk prediction as an occurrence modeling problem, analyze underlying patterns and design appropriate prediction models. Our approach for chimney fire prediction is quite general and can be transferred to other similar fire types, such as kitchen fires.
The literature for fire risk prediction is mostly concerned with wildfires. Overall, the prediction methods can be divided into two categories: machine learning based approaches (e.g. [20,25,34]) and statistical approaches (e.g. [28,48,52]). Usually, machine learning based approaches do not require prior knowledge of a fire type and its corresponding dependent variables, and they can detect the dependence between the fire risk and a large number of environmental candidate variables automatically using specialised learning algorithms, such as logistic regression [32], support vector machine [36], decision trees [42], random forests [34] and neural networks [37]. Most of machine learning algorithms are applicable for discrete data, whereas continuous hazard maps of fire risk are expected. Moreover, fire occurrences are usually recorded as spatio-temporal point patterns. To solve these gaps, certain studies (e.g. [36,37,42]) discretized fire incidents into areal unit data and sometimes transformed the risk prediction from an occurrence modeling problem to a risk scale classification problem (i.e. label different amounts of fires to corresponding risk scales), however, complete and real data inference may be lost. In addition, machine learning approaches often require relatively large amounts of data in order to obtain a satisfactory performance, and due to their 'black box' behaviour, it is difficult to specify the influence of an explanatory variable on fire occurrences.
In contrast, statistical approaches are able to learn spatio-temporal point patterns directly using point process tools, and are more interpretable as they are based on a specific, parametric, mathematical model which allows for theoretical confidence intervals as well. For instance, Hering et al. [17] and Costafreda-Aumedes et al. [11] used the K-function to analyze the clustering patterns of wildfire data. Ye [54] and Boubeta et al. [6] utilized Poisson structures to generate hazard maps of fire occurrences. As for complex models, Møller and Díaz-Avalos [27] considered spatial and temporal explanatory variables in a shot-noise Cox process and fitted it with minimum contrast techniques. Pereira et al. [30] and Serra et al. [40] modeled the spatio-temporal point patterns in a log-Gaussian Cox process, which was particularly designed to simulate latent phenomena. Recently, a Bayesian framework was developed and suggested to improve fire prediction [22]. Pimont et al. [31] employed a spatiotemporal log-Gaussian Cox process model ('Firelihood') to predict wildfires and established Bayesian inference for model components using integrated nested Laplace approximation [35]. Koh et al. [21] developed a joint hierarchical model framework by combining extreme-value theory and point processes and studied summer wildfire data for the French Mediterranean basin. For chimney fires, School [39] also used a log-Gaussian Cox process to predict fire risk based on the explanatory variables selected by Pearson correlation coefficients and the random effects simulated by Gaussian random fields.
Although various point processes models have been proposed to improve the performance of fire prediction, many studies (e.g. [27,30]) determine the model structures firstly and select and fit explanatory variables afterwards. This procedure does not fully exploit information contained in variables and sometimes requires complex modeling of residuals. In addition, the selection of explanatory variables is conducted either manually [48] or using basic statistical methods [39,53]. Several studies [8,47,55] proposed regularized penalty functions in the estimation of point processes for variable selection. However, a parametric form of the intensity function including all variables needs to be specified in advance.
To conduct a data-driven fire risk prediction study, we combine machine learning and point process tools in our modeling process to leverage the advantages of both types of approaches.
In an initial study [24], we used the permutation importance techniques of random forests to select important variables, as random forests are the most accurate machine learning method for fire prediction [34]. Accordingly, we designed a generalized linear Poisson model and fitted it on areal unit data to predict chimney fire risk in Twente. In this paper, we refine this model to a spatio-temporal Poisson point process as it targets at point patterns (cf. Figure 1(b)) directly, enabling more detailed analysis. Specifically, our contributions are: i) we use machine learning algorithms to select explanatory variables of chimney fires non-parametrically while taking variable dependence into account, so that the utilization of variable information can be maximized; ii) we adhere to a statistical Poisson point process model with an interpretable parametric structure designed on explanatory variables for fire occurrence modeling, so that complete inference on the data can be carried out; iii) we fit the model parameters using the efficient logistic regression approach by tuning it to focus on important parts of point patterns and validate the model using second-order summary statistics and residuals.
The remainder of the paper is organized as follows. Section 2 introduces the fire data we use as well as some data pre-processing steps. Section 3 discusses the selection of important variables to be used in the chimney fire risk prediction. In Section 4, we motivate and fit a Poisson point process model. Section 5 tests for point interactions in our fire data and validates the Poisson model assumption of absence of such interactions. Section 6 discusses the modeling process, compares the point process model with an areal unit model [24] and explores the role of the dummy intensity in logistic regression estimation. Section 7 provides a conclusion and ideas to further refine the model.

Data
In this section, we introduce the chimney fire data, the relevant environmental variables and the data pre-processing operations. Collaborating with the Twente Fire Brigade, we collected the data of all reported chimney fire incidents occurring between January 1, 2004 and December 31, 2020, in the Twente region, in the eastern part of the Netherlands (map shown in Figure 1(a)). After a manual check that removes obvious mistakes, the dataset consists of 1759 incidents. Each incident is reported individually with its ID number, location, time and a brief description of the circumstances of the fire and the rescue process 1 . The spatial and temporal projections 2 are plotted, respectively, in Figure 1(b) and 2. It is clearly visible in Figure 1(b) that the spatial distribution of chimney fires is heterogeneous in the sense that most incidents occur in urban areas, especially in cities with a higher population, such as Almelo, Hengelo and Enschede (cf. Figure 1(a)). Apparent clustering may also arise, because neighbouring buildings tend to have identical chimney types. Moreover, locations neighbouring in space and time share similar weather conditions. In Figure 2, the temporal distribution of chimney fires is clearly periodic: chimney fires occur more frequently in winter than in summer. 0 10 20 30 2 0 0 4 J a n 2 0 0 5 J a n 2 0 0 6 J a n 2 0 0 7 J a n 2 0 0 8 J a n 2 0 0 9 J a n 2 0 1 0 J a n 2 0 1 1 J a n 2 0 1 2 J a n 2 0 1 3 J a n 2 0 1 4 J a n 2 0 1 5 J a n 2 0 1 6 J a n 2 0 1 7 J a n 2 0 1 8 J a n 2 0 1 9 J a n 2 0 2 0 J a n 2 0 2 0 D e c According to the experts of the Twente Fire Brigade, various environmental factors could influence the risk of chimney fires. Spatially, building types, population density and composition and urbanity degrees in an area may determine the baseline risk of chimney fires. Temporally, season and weather conditions may also influence chimney fire risk. Table 1 lists 27 putative explanatory variables with specific abbreviations, descriptions and sources. Among them, population and urbanity information are recorded over 6291 pre-defined 500m × 500m area boxes in Twente, whereas building information contains the precise locations, ages and functions of the houses. Temporal data derives from the daily weather information observed at the weather station at Twenthe airport. To assess the influence of small variations in weather among different parts of the Twente region, we also collected the data from two neighbouring weather stations (i.e. Heino, Hupsel) outside Twente for analysis. Their locations relative to Twente can be found on https://www.knmi.nl.

Variable Abbreviation
Description House The total number of houses IFV V σ, 2 House indu The number of houses with an industrial function IFV V σ, 3 House hotl The number of houses with a hotel function IFV V σ, 4 House resi The number of houses with a residential function IFV V σ, 5 House 20 The number of houses constructed before 1920 IFV V σ, 6 House 2045 The number of houses constructed between 1920 and 1945 IFV V σ, 7 House 4570 The number of houses constructed between 1945 and 1970 IFV V σ, 8 House 7080 The number of houses constructed between 1970 and 1980 IFV V σ, 9 House 8090 The number of houses constructed between 1980 and 1990 IFV V σ, 10 House 90 The number of houses constructed after 1990 IFV V σ, 11 House frsd The number of free standing (detached or semi-detached) houses IFV V σ, 12 Resid The total number of residents CBS V σ, 13 Resid 14 The number of residents with an age in the range of 0 till 14 CBS V σ, 14 Resid 1524 The number of residents with an age in the range of 15 till 24 CBS V σ, 15 Resid 2544 The number of residents with an age in the range of 25 till 44 CBS V σ, 16 Resid 4564 The number of residents with an age in the range of 45 till 64 CBS V σ, 17 Resid 65 The number of residents with an age of 65 or higher CBS V σ, 18 Man The number of male residents CBS V σ, 19 Woman The number of female residents CBS V σ, 20 Address The number of addresses in the neighbourhood CBS V σ, 21 Urbanity 4 The urbanity of the neighbourhood CBS V σ, 22 Town Boolean variable indicating the presence of a town CBS V τ, 1 WindSpeed Daily mean wind speed (km/h) KNMI V τ, 2 Temperature Daily mean temperature (°C) KNMI V τ, 3 WindChill Daily mean wind chill (°C) (calculated from V τ,1 , V τ,2 ) V τ, 4 Sunshine Daily sunshine duration (h) KNMI V τ, 5 Visibility 5 Daily minimum visibility KNMI

Data Pre-processing
To unify the data to be used in the modeling process, we perform several data pre-processing operations on environmental variables. Firstly, we collect building information at the level of living units, as it may reflect the geographical information of chimneys in a more representative way. For instance, two families living in a semi-detached house actually use their individual chimneys. Secondly, some spatial variables (i.e. population and urbanity information) consist of a list of historical data, whereas others (i.e. building information) are only accessed in a single actual value. To enable similar treatment of the data, we use the averaged value over the time period of interest for the variables consisting of a list of historical data, so that for all spatial variables, a single value is accessible. Moreover, we only consider the buildings that are currently in use and assume that buildings keep their functions and types during the time period of interest. Observing that the age and type information of some buildings are missing, we assign them the label 'extra' during the counting process. Addi-tionally, we exclude the fire incidents and weather data on the leap day in leap years, while in the application of future fire risk prediction, we will use the prediction of February 28 to compensate for the missing prediction of the leap day, February 29. Finally, to obtain spatial environmental variables for every location in Twente from areal unit data (i.e. population and urbanity information) or data with precise coordinates (i.e. building information), we employ kernel smoothing. For temporal variables, we assume their values at different times of a day remain invariant and equal to the daily mean values. Based on the pre-processing steps above, we introduce the following notation for environmental variables: -V σ,i (u): the smoothed value of the i-th spatial variable at location u, where (u, t) denotes any location and time combination in the spatio-temporal domain.

Selection of Explanatory Variables
Obviously, not all putative environmental variables listed in Table 1 are required to explain the incidence of chimney fires and some of them are mutually dependent. To select the most important variables, we perform a non-parametric variable importance analysis using the permutation importance techniques of random forests [7].

Random Forests and Permutation Importance
Random forests [7] are widely used as robust classification and regression methods in many applications, such as risk prediction [51] and data mining [38]. A random forest is usually composed of hundreds or thousands of decision trees, where each tree is generated on a sampled subset of the data by repeated bagging (i.e. bootstrap sampling with replacement) and trained to fit the explanatory variables to the response variable. Afterwards, a combined result over all trees will be reported as the final output. In addition, random forests can be used to assess the importance of a variable by means of permutation importance techniques [1,7,46]. Through randomly permuting the values of a variable over the observations, the importance of a variable is defined as the mean increase of the prediction error over all trees computed on the permuted data compared to that computed on the original data.
More formally, consider a regression problem on a dataset D with n observations, and suppose that for each observation of the dataset, there are m explanatory variables. Let x j i and y i , with i = 1, ..., n and j = 1, ..., m, denote the j-th explanatory variable and the response variable respectively for the observation i; x i is the vector that collects all x j i . The construction of the random forest consists of the generation of T decision trees. To generate a tree t, a subset of D, denoted as B(t), is sampled by bagging. At each node of the tree, a number of explanatory variables are selected randomly from all variables as the candidates. Then, one of the candidates, say x j , is used to split the node into two subsets, B L (t) and , in a way that the residual sum of squares is minimized. Here,ȳ B L (t) andȳ B R (t) denote the mean of the response variables in the corresponding subsets. Such node splitting procedure is continued until the tree satisfies certain preconditions (e.g. a maximum number of levels of the tree). Similar tree generation processes are employed to construct all other trees in the forest as well. To measure the importance of an explanatory variable, say x j , in tree t, we randomly permute its values over the out-of-bag observations for tree t, denoted as oB(t) (i.e. oB(t) = D(t) \ B(t)). Writing π j for the permutation, the value of the of the j-th explanatory variable for response y i changes from x j i to x j π j (i) when i ∈ oB(t). The values of other explanatory variables are left unchanged. Then, the increase of the prediction error is whereŷ i denotes the prediction in tree t for observation i using unpermuted explanatory variables,ŷ i,π j denotes the corresponding prediction using explanatory variables with the j-th variable permuted and |oB(t)| denotes the number of out-of-bag observations for tree t. Finally, the mean increase of the prediction error over all trees, T t=1 I(x j ; t)/T , is used to illustrate the importance of the j-th variable in the forest.
In practice, the original construction algorithms of random forests tend to bias the variable selection at tree nodes to factorial variables with many categories or continuous variables with many cut points. In addition, traditional permutation importance techniques can sometimes be misled by correlated explanatory variables. To address these problems, unbiased random forests [19,44] and conditional permutation importance techniques [43] were proposed based on a conditional inference framework of recursive partitioning. Compared to earlier methods, the idea is to partition the variable space in order to obtain groups of observations with similar association patterns instead of groups of observations with merely similar values of the response variable. For instance, to measure the importance of variable x j under conditional permutation in tree t, the set of out-of-bag observations, oB(t), is partitioned into a grid where each block of it shares the same information on the remaining variables, x \ x j . Then, the permutation is applied to each block and the increase of the prediction error is computed and summed up over all blocks as the importance score of x j . With these significant increments and optimizations, the random forest approach has proved very useful to measure variable importance non-parametrically considering variable dependence [45]. Note that the notation used in this section is only valid here to explain random forests and permutation importance statistically.

Variable Importance Analysis
In our study, we measure the importance of spatial and temporal variables separately on areal unit data [24]. The main reason for this is that spatial variables are accessed in a single actual value while temporal variables are accessed as historical data. Specifically, for spatial variables, we group the incidents in 6291 pre-defined 500m × 500m area boxes and merge them into a large data table consisting of 6291 rows and 23 columns. Each row indicates an area box, and the 23 columns refer to the number of incidents occurring in that box as well as the values of 22 spatial variables. In practice, we find that some boundary area boxes lie only partially in Twente. However, the population statistics for such boxes do not distinguish between Twente and the neighbouring regions. To avoid bias, we filter out the data of these boundary boxes in our variable importance analysis. For temporal variables, we group the incidents in 6205 days and merge them into a data table consisting of 6205 rows and 6 columns. Each row indicates a specific day, and the 6 columns refer to the number of incidents on that day as well as the values of 5 temporal variables. In the analysis, we construct unbiased random forests of 2000 trees and set the proportion of the number of input variables that are randomly sampled as candidates at a tree node to approximately 1/3 (i.e. spatial: 7 candidates, temporal: 2 candidates). Considering the high correlation between certain variables (e.g. the urbanity and the population), we use conditional permutation importance techniques [43] instead of the traditional ones [7] to suppress the importance scores biasing towards correlated variables. Our implementation of variable importance analysis is based on the R-package party [18,44]. We use the method proposed in [13] to compute the conditional permutation importance which provides a faster computation and shows more stable results than the original implementation [43].
The conditional permutation variable importance results for spatial and temporal variables are plotted in Figure 3 and 4 respectively. For comparison, we plot the traditional permutation results as well. In all plots, the y-axis refers to the increase of the squared prediction error when traditional or conditional permutation on a variable is applied. A large increase indicates that the variable is very important for correct predictions, while a decrease indicates that the variable has no influence on or even hampers the prediction. If we compare the traditional and conditional permutation results, we see that the bias due to correlated variables is well suppressed. For instance, the number of residents aged between 45 and 64 and the total number of residents have the largest importance scores under traditional permutation, with the number of free standing houses coming the third. However, the latter has the highest variable importance score under conditional permutation. A possible explanation could be that free standing houses are mostly occupied by families whose adult members aged between 45 and 64, which implies a strong positive correlation among all three categories. Moreover, other age groups such as the elderly or families with young children may be less inclined to use their chimney. As a result, the importance scores of the two variables concerning residents decrease a lot under conditional permutation. Similar observations hold for temporal variables as well. Both wind chill and temperature obtain high importance scores under traditional permutation although they are correlated. Conditional permutation detects the underlying correlation and suppresses temperature, as wind chill is defined in terms of both temperature and wind speed information. Since wind speed still obtains a relatively large importance score even under conditional permutation, we consider it in further modeling. Additional analysis to assess the influence of small weather variations among different parts in Twente indicates that these variations can be ignored.
Overall, according to conditional permutation variable importance, the number of   : Importance (y-axis) obtained for temporal variables using traditional (a) and conditional (b) permutation techniques.
buildings constructed between 1920 and 1945 and the number of free standing houses are the most important spatial variables, and wind chill and wind speed are the most important temporal variables. The results could have the following explanations: i) most free standing houses contain chimneys whereas other types do not; ii) chimney pipes in old buildings tend to be made of brick rather than metal, which increases the risk to catch fires; iii) strong wind can fuel a fire; iv) wind chill reflects people feeling cold thus inducing them to use their fires and chimneys.

Poisson Point Process Model
In this section, we motivate and propose a nested Poisson point process model for chimney fire prediction. We also present the model fitting and selection procedure and propose confidence intervals for both model parameters and predicted fire intensities.

Motivations
To design an appropriate prediction model structure, we investigate the relations between the selected explanatory variables and chimney fire occurrences. Firstly, we divide the houses into four house types depending on the age (whether they have been constructed between 1920 and 1945 or not) and on whether they are free standing or not, and plot the monthly intensities of chimney fires per house for different house types separately. Figure 5 shows that, generally, chimney fire occurrences in all house types are periodic, with incidents concentrated in the colder seasons. However, different house types run different risks of chimney fires: the intensities for old houses (i.e. House 2045) and free standing houses are higher than others. Moreover, note that the patterns are not perfectly periodic; the amplitudes of the peaks vary per year. To model such varying patterns, we take wind chill and wind speed into account. In addition, we plot the numerical relations between chimney fire occurrences and the number of houses of different types in Figure 6 based on areal unit data [24]. To approximate the underlying trends in scatter plots, we apply the locally estimated scatter plot smoothing method [10] and plot the estimated trends. It is interesting to find that, for each house type, chimney fire occurrences increase approximately linearly in the number of houses of that type. However, this observation only holds when the number of houses is not very large. If it is too large, the number of chimney fires tends to a saturated value. To summarize, we draw the following conclusions: i) chimneys in different house types catch fires at a type-dependent rate, and the rate is influenced by both seasonal information and temporal explanatory variables; ii) chimney fire occurrences in houses of a specific type are approximately proportional to the number of houses of that type.

Model Structure
Based on the motivations above, we design the following chimney fire risk prediction model. We suppose that each house catches a chimney fire independently at a type-dependent rate that varies in time. Under this assumption, the model structure is that of a Poisson point process defined on the Twente region and on the period of year 2004-2020, with the intensity function of the form λ(u, t) = k λ k (u, t) where λ(u, t) indicates the risk intensity at location u at time t and the subscript k indicates for houses of type k. Moreover, λ k (u, t) = h k (u)λ k (t), where h k (u) indicates the density of house type k at location u and λ k (t) indicates the risk intensity for a house of type k at time t. Such a structure reflects the conclusions reached in Section 4.1, although it relies on the condition that the density of houses at a location is not too large.
Formally, our nested Poisson point process model reads The density maps of houses of four types (i.e. h k (u)), as shown in Figure 7, are derived from V σ,6 (u) and V σ,11 (u) in Table 1 by smoothing, using a Gaussian kernel with a standard deviation of 1,000 metres. Since both wind speed, V τ,1 (t), and wind chill, V τ,3 (t), were selected as important temporal variables in Section 3.2, we propose a temporal intensity function of the form where a harmonic function 6 with order o k,1 is employed to model seasonal variations and polynomial functions with order o k,2 and o k,3 are used to model information of wind speed and wind chill. In addition, we employ another polynomial function of V τ,1 (t)V τ,3 (t) with order o k,4 to allow for interactions between wind speed and wind chill. The exponential function is applied to guarantee that the risk intensity function stays positive. For purposes of parametric representation, the temporal intensity functions λ k (t) can also be described as λ θ k (u, t), where θ k indicates the vector of coefficients in the harmonic and polynomial functions.

Model Fitting
In [24], we fitted a fire prediction model similar to (3) and (4) on areal unit data using maximum likelihood estimation. We established theoretical confidence intervals for predicted intensities based on classical theories: the Fisher information for asymptotic normality of generalized linear models [14,26] and the Delta method [50]. In principle, we could also apply this approach here to maximize the log-likelihood function of our Poisson point process where x k runs through chimney fire incidents occurring in houses of type k. However, since we have a continuous space-time domain, the integral in the log-likelihood function must be approximated numerically, e.g. using quadrature points [4]. In our case, the number of quadrature points would have to be rather large, thus, we use the logistic regression estimation [2] to fit our Poisson point process model efficiently.
The idea behind the logistic regression estimation approach is the well-known Campbell-Mecke theorem (see, e.g. [12]). Consider a point process X that is defined on a space-time domain R 2 × R and is locally finite, which means X is almost surely finite for every bounded subset of R 2 ×R. We assume that X is confined to a bounded domain W ×T ⊆ R 2 ×R and has intensity function λ(u, t), where (u, t) denotes a location and a time in the domain. For any real-valued non-negative or integrable function f (u, t) defined on W ×T , the Campbell-Mecke theorem reads where x runs through the points of X. Suppose the intensity function is parametric and depends on a parameter vector θ, thus can be described as λ θ (u, t), one might estimate both sides of the Campbell-Mecke formula and equate them to obtain an estimating equation for θ.
The logistic regression estimation is based on the function where ∇λ θ (u, t) denotes the gradient of λ θ (u, t) over θ and ρ is a positive-valued function defined on W × T . Here, one needs to assume that λ θ > 0. To estimate the integral in (6) with f (u, t) defined in (7), we use a 'dummy' point process D on W × T that is independent of X and has intensity function of ρ. Apply the Campbell-Mecke theorem (6) to D, one has Having obtained unbiased estimators for both the left and right hand side of equation (6) with f given by (7), we plug them in to obtain the estimating equation. It solves over the parameter θ. Note that the subscript x in two sums runs through the points of X and D respectively. It is interesting to observe that (9) is precisely the derivative of a logistic log-likelihood function In this way, the fitting of Poisson point processes can be transformed into a maximum likelihood estimation problem of logistic regression. It can be easily implemented based on the R-package stats [49] and treated as the fitting of generalized linear models with a logit link function. Since the log-likelihood function is concave, the existence and uniqueness of parameter θ that maximizes the log-likelihood are guaranteed under certain conditions [41].

Modeling Chimney Fire Data
We use logistic regression estimation to fit the intensity functions of our four house types separately. Consider (4) for fixed k and write θ k for the vector of coefficients in the harmonic and polynomial functions. Also consider the chimney fires occurring in houses of type k as a point process X k with intensity function λ θ k (u, t). Denote the corresponding dummy point process by D k with intensity function ρ k (u, t). Assuming we have m parameters in θ k , the estimating equations read where C p (x) denotes the p-th temporal covariate (i.e. the harmonic components, wind chill terms with different orders and the interaction term of wind chill and wind speed in (4)) at point x. Furthermore, we perform a three-step implementation to estimate θ k . Firstly, we specify the intensity function ρ k (u, t) of the dummy point process D k . Considering that, spatially, most chimney fires occur in urban areas, and temporally, there are more chimney fires in winter than in summer, we tune ρ k (u, t) to concentrate on important regions and seasons so as to fit the model efficiently. Specifically, we use the density map of houses of the given type k (cf. Figure 7) to distinguish urban areas from rural areas and design a sine function to assign higher ρ k (u, t) to winter seasons so as to comply with the observations in Figure 2. Thus, the intensity function ρ k (u, t) of D k is given by where r k is a multiplication factor 7 used to guarantee that ρ k (u, t) is at least four times λ θ k (u, t) as suggested in [2] and h k (u) is the density of houses of type k at location u. Secondly, we generate a realisation from the dummy point process D k using the R-package spatstat [3]. Thirdly, based on the observation of the target point process X k and the realisation of the dummy point process D k , we use the R-package stats [49] to estimate the model parameters θ k for house type k. In addition, to determine the optimal function orders (i.e. o k,1 , o k,2 , o k,3 , o k,4 in (4)) simultaneously, following [24], we apply a grid search to select the combination that yields the smallest Akaike information criterion [9] given a set of ranges 8 for them. Considering that wind chill already contains some information of wind speed, we also perform likelihood ratio tests over the best models with and without wind speed and the interaction term. Similar to [24], we separate the fire data into two sets: the data on the period 2004-2019 and the data on the year 2020, for fitting and predicting respectively. Performing the model selection discussed above on the fitting data, we obtain corresponding optimal prediction models for all four house types and find that wind chill alone is mostly sufficient to obtain accurate risk predictions, except for house type 4, where the interaction term provides improved performance. Finally, the optimal temporal risk intensity functions obtained for the four house types read as follows: + θ 1,10 V τ,3 (t) + θ 1,11 V 2 τ,3 (t)), λ 2 (t) = exp(θ 2,1 + θ 2,2 cos( 2π 365 t) + θ 2,3 sin( 2π 365 t) + θ 2,4 cos( 4π 365 t) + θ 2,5 sin( 4π 365 t) 7 Specifically, we set r k to 60, 20, 20 and 8 for the four house types respectively. where t refers to a specific day during 2004-2020, thus t ∈ [0, 6205]. We provide the estimates of all model parameters (i.e. θ) in Table 2. Additionally, we also plot the spatial and temporal predictions and the actual realisations of year 2020 in Figure 8. Overall, both the spatial and temporal predictions capture the correct trends. Specifically, in the spatial domain, our model learns the urbanity features of the Twente region which are highly relevant for the fire data. Since the spatial prediction depends on the distribution of houses of different types which can only be accessed as actual data, our model displays similar spatial patterns for different years. In the temporal domain, our model captures the periodic pattern on an annual basis and adds weather variable dependent information, which explains the noisier aspect of Figure 8(b) compared to Figure 8(a). A comparison with the areal unit model and the ρ-tuning experiments will be specifically discussed in Section 6. In summary, given appropriate building information, our prediction model can detect areas with higher risks of chimney fires in Twente and, additionally, estimate the risks for specific days based on the weather forecast.

Confidence Interval
For risk prediction problems, it is important to quantify the uncertainty of the estimates. Baddeley et al. [2] demonstrated asymptotic normality for their logistic regression estimators for stationary point processes in space when the domain increases. Recently, Choiruddin et al. [8] proved asymptotic normality in an increasing-domain setting for spatial point process models that involve explanatory variables. In a companion paper, we study consistency and asymptotic normality for the logistic regression estimators for spatio-temporal Poisson point processes in an infill asymptotics setting. Specifically, under appropriate conditions, the maximiserθ of (9), under true value θ o , is approximately normally distributed with mean θ o and covariance matrix G given by the inverse of the Godambe matrix [16]. A plug-in estimator for G isĜ assuming we have m parameters in θ. Approximate confidence intervals for the model parameters are readily obtained. For chimney fire modeling, the approximate 95% confidence intervals of the model parameters are listed in Table 2. The Delta method [50] can then be used to caclulate approximate confidence intervals for predicted fire intensities. We visualize the temporal risk predictions for the year 2020 in Figure 8(b).

Model Validation
In this section, we validate the independence assumption underlying Poisson point processes by analyzing the second-order properties. We also validate the first-order structure of our model by visualizing the residuals.

Second-order Property Analysis
To assess the validity of the Poisson assumption, we perform a second-order analysis on the chimney fire data using two summary statistics: the pair correlation function and the K-function (see, e.g. [23]). The pair correlation function can be used to detect inter-point interaction for specific distances, whereas the K-function is a cumulative statistic that can be used to test the existence of point interactions in data.

The Pair Correlation Function
Consider a point process X defined on a space-time domain W × T ⊂ R 2 × R that is locally finite. The second-order product density λ (2) (x, y) of the pair (x, y), where x, y ∈ W × T , is the infinitesimal probability that X places points at x and y. The pair correlation function for the pair (x, y) is then defined as where λ(x) and λ(y) denote the intensities at x and y, and we use the convention that a/0 = 0 for all a ≥ 0. For purposes of visualisation, we consider the pair correlation functions for the projections in space and time domain separately. Under appropriate stationarity and isotropy assumptions, in d dimensions (spatially, d = 2; temporally, d = 1), (15) is a function of r = ||x − y|| only 10 and can be estimated bŷ where s d denotes the surface area of the unit sphere in R d , = denotes the summation over all pairs of distinct points, k b (·) is an one-dimensional smoothing kernel with bandwidth b, 1/|W ∩ W x−y | is the spatial edge correction factor [29] (temporally, 1/|T ∩ T x−y |), and |A| is the volume of A ⊂ R d . If there is no interaction in the point pattern, the pair correlation function is 1 at any distance. A value ofĝ(r) larger than 1 suggests clustering in the sense that points are more likely to occur from one another at distance r, whereas a value smaller than 1 suggests inhibition that points tend to mutually avoid each other at this distance. Moreover, if the intensity function λ(x) is unknown, a plug-in estimator can be used. In practice, we employ (16) with a Gaussian smoothing kernel as k b and plug in the fitted intensity functions based on data of all 17 years. We set the possible interaction intervals for the spatial and temporal domain to 10, 000 metres and 100 days, respectively, and the smoothing bandwidths to 500 metres and 10 days. The estimated pair correlation functions are plotted in Figure 9. Overall, both the spatial and temporal results show a curve that converges to 1. At smaller range, there is some attraction between the points spatially. Temporally, possibly because of the precision of the observations (daily basis), there is some repulsion and a weak peak at an interval of about twenty days. The latter might be due to an underestimate of the intensity function on a specific day.   Figure 1(b) and 2) of chimney fire incidents during 2004-2020.

The K-function
Since the pair correlation function shows some interactions at small distances, we perform a joint space-time K-function test on the chimney fire data. The K-function, also known as Ripley's reduced second moment function, of a stationary point process X is defined as the expected number of other points within a given distance r of an arbitrary point divided by the intensity. In an inhomogeneous and space-time version, it takes the form of where B denotes any Borel set, x and y are distinctive points with x(u), y(u) and x(t), y(t) representing their spatial and temporal positions, and 1(·) is the indicator function. Under appropriate weak stationarity assumptions [15], a point-wise unbiased estimator of the Kfunction isK with notation as in the previous section. If no interactions exist, K inhom (r, v) = 2πr 2 v. Moreover, when the intensity function λ(x) is unknown, a plug-in estimator can be used. We compute the inhomogeneous K-function over data from all 17 years. To obtain an elegant view of it, we set the testing sequences of spatial and temporal distances to 10, 000 metres and 100 days and divide both of them into 100 step size pairs of 100 metres and 1 day. We calculate the K-function using the fitted intensity functions to estimate λ(·) in (18) and plot the result in Figure 10(a). Same with the pair correlation function, we find some evidence of clustering. To test whether the interaction is statistically significant, we implement a Monte Carlo test. We generate 99 realisations of a Poisson point process with the fitted intensity functions, and compute their K-functions and that of the actual chimney fire data using kernel estimators to estimate λ(·) in (18). We plot the local envelopes in Figure 10(b). As the empirical K-function of the actual chimney fire data lies completely within the envelope, there is no reason to look beyond a Poisson model.

Residual Analysis
To further verify the validity of our point process model, we perform a residual analysis to check whether typical patterns still exist after the modeling of Poisson point process with specific form defined in (3), (4) and (13). We fit the model using the data of all 17 years and compute spatial and temporal residuals separately. As suggested in [5], the spatial residuals are defined as the difference between the smoothed fitted fire risks and the smoothed actual realisations using a Gaussian kernel with a standard deviation of 1, 000 metres. Temporally, we define the residuals as the difference between predicted and actual monthly counts. The results of the residual analysis are plotted in Figure 11. In the spatial domain, no specific pattern is shown, except that our point process model overestimates for big city centres (e.g. Enschede and Hengelo) whereas underestimates a bit for small city centres. It is reasonable, as our model merely depends the spatial distribution of chimney fires on the density of corresponding house types and we observe a saturation pattern of chimney fire occurrences when the density of a house type increases (cf. Figure 6). A possible solution for that is to apply a piece-wise function instead of the linear function to model the relation between chimney fire occurrences and the density of a house type. In the temporal domain, no specific pattern is evident. 2 0 0 4 J a n 2 0 0 5 J a n 2 0 0 6 J a n 2 0 0 7 J a n 2 0 0 8 J a n 2 0 0 9 J a n 2 0 1 0 J a n 2 0 1 1 J a n 2 0 1 2 J a n 2 0 1 3 J a n 2 0 1 4 J a n 2 0 1 5 J a n 2 0 1 6 J a n 2 0 1 7 J a n 2 0 1 8 J a n 2 0 1 9 J a n 2 0 2 0 J a n 2 0 2 0 D e c (b) temporal residuals (monthly counts) Figure 11: Spatial (a) and temporal (b) residuals based on the chimney fire data during 2004-2020. The unit in the spatial plot is metre −2 ; the unit in the temporal plot is month −1 .

Discussion
In this section, we discuss the data-driven modeling process designed for chimney fire prediction, which can be generalized to similar fire types. We also compare the areal unit model in [24] with our point process model. Additionally, we present the influence of tuning ρ to important regions and times in the logistic regression estimation.

Modelling Process
In this paper, we design a two-step data-driven modeling process for chimney fire prediction.
In the first step, we combine machine learning and statistics, where the former is utilized to select the most important environmental variables and the latter appropriately approximates and explains the relations between explanatory variables and fire incidence. We employ parametric functions to model such relations and fit them using the logistic regression approach, so that the influence of a variable on chimney fires can be analyzed by its coefficient in the function model. This combination offers significant advantages. From the machine learning side: i) compared to naive statistics used for the selection of explanatory variables, random forests can detect both linear and non-linear correlation between a candidate variable and chimney fire occurrences; ii) random forests also consider the dependence among candidate variables non-parametrically; iii) by applying the conditional permutation importance technique instead of the traditional one, the bias towards correlated variables is suppressed as well. From the statistics side: i) classic statistical methods for regression modeling with specific mathematical expressions are more interpretable and based on natural assumptions; ii) statistical methods also allow for significance testing and uncertainty quantification.
In the second step, we use summary statistics, the pair correlation function and the Kfunction, to detect point interactions. We also perform a residual analysis to validate the model structure. Although, in our study, a Poisson point process model is already proved sufficient for chimney fire prediction, a hierarchical modeling where random effects are incor-porated to deal with latent effects, might be more appropriate for different fire types.

Model Comparison
To compare the performance of the areal unit model in [24] and our point process model, we plot spatial and temporal predictions for the year 2020 for the former model in Figure 12. Since the spatial predictions from this areal unit model are computed for 500m× 500m area boxes, to comply with the scaling level of point process model, we average the predictions by the volume of these area boxes. Generally, comparing Figure 12 to Figure 8, we can find that both the areal unit model and the point process model capture similar spatial and temporal patterns from the fire data. However, the spatial predictions from the point process model behave smoother than those from the areal unit model, as we apply the kernel smoothing in the former to obtain the density of a house type at a location, which helps establish spatially continuous predictions. The temporal predictions from the two models are almost identical for every day, because our prediction model (cf. (3), (4) and (13)) mainly learns temporal information and merely distributes temporal predictions to spatial locations according to the densities of certain house types.
More specifically, in terms of the total chimney fire risks predicted for the year 2020, the areal unit model predicts 92.86 and the point process model predicts 92.95. The observed count of 81 lies in the 95% confidence intervals of both models. However, the areal unit model overestimates many more fire risks for certain area boxes of big city centres (e.g. Enschede). Such phenomena not only result from the different smoothing strategies the two models employ, but also derive from the fact that we fit the fire data over the whole Twente region, which means that detailed information between the number of houses of a type and chimney fire occurrences can be lost. An additional explanation is offered by the saturation effect suggested by Figure 6. From a computational perspective, the logistic regression approach employed for our point process model enables to capture spatial information efficiently and obtains more reasonable spatial predictions.  Figure 12: Spatial (a) and temporal (b) predictions based on the areal unit model in [24] for the year 2020. The unit in the spatial plot is metre −2 ; the unit in the temporal plot is day −1 .

ρ Tuning in the Logistic Regression Estimation
The role of the dummy point process D is to estimate the integral in (6). In [2], a rule of thumb was proposed that ρ(u, t) should be at least four times λ(u, t); additionally, a data-driven selection of ρ was suggested to obtain relatively precise estimate under a limited number of dummy points.
Based on such considerations, in (12), we tune ρ to concentrate on urban areas and winter. Here, we perform two experiments to assess the influence of ρ tuning in the spatial and temporal domains separately by fitting our prediction model with and without ρ tuning on the data of all 17 years and plotting the difference of the fitted values. We also set the expectations of the number of dummy points as equal both with and without ρ tuning. In addition, since the realisations of the dummy point process generated in different runs contain randomness, we perform each experiment 60 times and plot the averaged spatial and temporal differences. Figure 13(a) shows that, spatially, a concentration on urban areas reduces the bias caused by the saturation effect shown in Figure 6. Figure 13(b) shows that, temporally, with a concentration on winter seasons, the fitted risk intensities increase in March and April, which are known as the weather tipping seasons and are found more likely to catch chimney fires in [39]. To summarize, tuning ρ to important regions and times is helpful to estimate the fire risk intensities more efficiently and to capture important underlying patterns.  Figure 13: Spatial (a) and temporal (b) difference of the fitted intensities using tuned ρ compared to uniform ρ with an equal expected number of dummy points in the logistic regression estimation (averaged on 60 random runs). The unit in the spatial plot is metre −2 ; the unit in the temporal plot is day −1 .

Conclusion
In this paper, we proposed a data-driven modeling process for chimney fire risk prediction, which can be generalized to studies on other fire types. Firstly, we applied random forests and permutation importance techniques to identify the important explanatory variables nonparametrically from a large number of environmental variables. From a practical perspective, our results indicate that pre-war detached or semi-detached houses run a higher risk of chimney fires and therefore public awareness campaigns should preferentially target the owners of such houses. Secondly, based on the observed relations between the selected variables and chimney fire occurrences, we designed a Poisson point process model to predict fire risk to learn the observed spatio-temporal point pattern directly. We applied the logistic regression estimation for model fitting and provided asymptotic confidence intervals. Additionally, to validate our model assumption, we performed a second-order property analysis and a residual study to test point interactions in the fire data. Last but not least, we reviewed our modeling process, compared the areal unit model and the point process model and discussed the ρ tuning required for the logistic regression estimation.
For future work, the first interesting direction would be to extend the data-driven modeling process of chimney fire risk prediction to similar fire types. Additionally, we are investigating the asymtotics of the logistic regression estimation under the framework of infill asymptotics using limit results for U-statistics [33].