{'fraction_non_alphanumeric': 0.03185869698748311, 'fraction_numerical': 0.009554879130778996, 'mean_word_length': 4.6516238525032785, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Two-Stage Cox Process Model with Spatial and Nonspatial Covariates

Rich new marked point process data allow researchers to consider disparate problems such as the factors affecting the location and type of police use of force incidents, and the characteristics that impact the location and size of forest fires. We develop a two-stage log Gaussian Cox process that models these data in terms of both spatial (community-level) and nonspatial (individual or event-level) characteristics; both types of covariates are present in the examples we consider and are not easy to incorporate via existing methods. Via simulated and real data examples we find that our model is easy to interpret and flexible, accommodating multiple types of marks and multiple types of spatial covariates. In the first example we consider, our approach allows us to study the impact of community-level socioeconomic features such as unemployment as well as event-level features such as officer tenure on force used by police, illustrated through simulated examples. In our second example we consider factors that impact the locations and severity of forest fires from the Castilla-La Mancha region of Spain between 2004-2007.


Introduction
Rich marked point process data are increasingly available in a number of domains. Such data enable researchers to answer interesting new questions. However, these data often pose new statistical and computational challenges. We consider two motivating examples: police use of force where the mark is the level of force used by police, and forest fires, where the mark is the amount of burned area caused by each fire. In addition to understanding the nature and driving factors of individual use of force incidents or forest fires, it is also critical for scientists, policymakers, and the public to develop a data-driven, systematic understanding of the important factors at play in driving outcomes generally.
The first motivating dataset we describe is police use of force events in the city of Dallas, Texas. Fully understanding police use of force requires a model for a point process that has multiple types of predictors, including both spatial and nonspatial variables. Examples of spatial variables that may impact police use of force are neighborhood-level characteristics such as the unemployment rate or neighborhood diversity. Nonspatial variables include both individual-level characteristics, such as officer or citizen race and officer training, and event-level characteristics, such as presence of citizen resistance or if the officer or citizen was injured. Our modeling framework determines the potential impact of these variables on the locations of events and on the level of force used by police (the mark). Examples of simulated marks and nonspatial variables following the structure of use of force data from Dallas, Texas in plotted in Figure 1, with spatial variables from the American Communities Survey plotted at the bottom.  In the second motivating example of forest fires, potential spatial variables of interest include elevation and terrain; potential nonspatial variables of interest include the time of year that the fire occurred and the cause of the fire (intentional or not). We are interested in determining the potential effect of these variables on the location of fires and the amount of burned area at those locations (the mark). Data from this application for the Castilla-La Mancha region of Spain is shown in Figure 2. We note that our model can incorporate multiple types of marks which is illustrated by our applications (categorical for use of force, continuous for forest fires) and multiple types of spatial variables (areal for use of force, geostatistical for forest fire data).
There is considerable attention and statistical methodology for point process models involving spatial variables, such as the log Gaussian Cox process (LGCP) (cf. Møller et al., 1998;Møller and Waagepetersen, 2003). These LGCP studies have included a wide variety of applications including crime, disease cases, agriculture, and seismic activity (Shirota and Gelfand, 2017;Li et al., 2012;Brix and Møller, 2001;Siino et al., 2018). There is also some work in marked point processes incorporating spatial dependence in the marking process, where the mean and the variance of the mark is dependent on the spatial intensity function Myllymäki and Penttinen (2009). However, there has been less extensive work on incorporating both spatial and nonspatial variables into marked point process models. Furthermore, it is important to have a computationally tractable inferential procedure for this model, and we would like the resulting inferences to be easy to interpret. Here we develop a novel two-stage log Gaussian Cox process (LGCP) model for a marked point process using a Bayesian framework. Our model utilizes both event/individual-level (nonspatial) and neighborhood-level (spatial) characteristics. We use this model to determine the impact of spatial variables on the locations of events and the relationship between both kinds of variables and the mark of the point process. We note that some studies may refer to all data available only at point process locations as marks. However, we find that for clarity in our framework developed here, there is value in distinguishing between dependent and independent variables in the marking process. We refer to independent variables as nonspatial variables, and dependent variables as marks. We also allow for dependence between the location and mark determination processes. We find that inference for our model is reliable and that our approach provides an intuitive interpretation of the impact of spatial and event/individual-level on marked point processes.
We note that our model provides advantages over existing methods, especially when there is dependence between the location-determination and mark-determination stages of the process. There are multiple ways for dependence to be characterized in marked point processes, including between marks and points, points and the covariates, between marks and the covariates, and between all three (cf. Dvořák et al., 2020). Our framework allows users to capture dependence between many of these dependence structures in each stage of the model and between the two stages. We present possibilities for the location-determination stage (NHPP and LGCP) and the mark-determination stage (logistic and linear regression) but there are a vast number of possible models that could be chosen for each stage, especially when the stages are independent. Under independence of the location and mark determination processes, their analysis can also be treated independently by classical point process methods for locations and by spatial regression methods for marks. Our regression framework, which allows for simple incorporation of dependent Gaussian processes in both stages of the model, creates a method for further evaluating dependence between stages after controlling for spatial variables and other dependence structures. It is difficult to perform this extent of diagnosis and evaluation of dependence between point locations, marks, and covariates using existing methods.
Previous approaches to incorporating nonspatial variables into a point process framework have typically been in the context of public health data. Liang et al. (2008) addresses spatial misalignment between spatial and nonspatial variables with an application to disease cases. Specifically, the authors characterize spatial variation of disease occurrence to identify areas with elevated risk. Their bivariate Gaussian process approach, which we call the bivariate mark model in this paper, breaks new ground by specifying an approach for handling nonspatial covariates in point process models as well as allowing for dependence between marks. Our simulated examples suggest that their approach may pose some parameter identifiability and interpretability issues. Best et al. (2000) studies individual-level outcomes of disease cases. The authors utilize individual-level covariate information from patients, such as sex and home conditions like dampness or parental smoking, as well as a spatial random field to account for unattributed spatially varying risk. They also incorporate the level of nitrogen dioxide, which is considered to be spatially varying. The model incorporates both individual-level and spatial variables in order to compare baseline and relative risk, under given sets of conditions. As in Liang et al. (2008), the model measures both location and attribute-specific risk intensity simultaneously, which makes it difficult to describe the spatial distribution of the outcome without choosing particular individual-level attributes. Quick et al. (2015) uses both spatial and nonspatial information to study mortality rates in an LGCP framework, similar to the model proposed by Liang et al. (2008). However, instead of treating nonspatial point-level information as covariates, they treat all nonspatial variables (race, sex, education, and cause of death) as marks through a cross-tabulation of all possible combinations of categorical nonspatial variables to create 24 categorical marks, each of which is called a group. The intensity is then estimated for each of these 24 combinations/groups. The method is applied to North Carolina mortality data, where the authors analyze the relationship between variables such as education and cause of death, where both are included in this crosstabulation of marks. This is a rich model but parameter interpretation can pose some challenges, as it depends on groups that are defined by the specific combination of categorical marks, and it is difficult to incorporate continuous explanatory variables. As the number of nonspatial variables in the model increases, the number of groups would continue to grow, leading to interpretability and computational challenges.
Pinto Junior et al. (2015) studies cerebrovascular diseases (CVD) and the nonspatial/ individuallevel variables that may play a role in the intensity of cases, namely age at death, education level, gender, and marital status. These nonspatial covariates are incorporated into the intensity function in a similar way to the model proposed by Liang et al. (2008) but Pinto Junior et al. (2015 allows for additional spatial variation in the effect of the nonspatial covariates. The model includes three Gaussian processes: one for nonspatial variables, another for the interaction between spatial and nonspatial variables, and a third for an unmeasured spatial random effect. This is a very flexible model formulation. However, we are more interested in the effect of individuallevel/nonspatial variables on the mark of the point process than on how the effect of these variables change over space. The model presented by Pinto Junior et al. (2015) is not extended to marked point processes. Furthermore, their flexible model results in large computational challenges, which the authors address by discretizing the analysis into sub-regions. We summarize the main contributions of this paper below.
• We introduce a novel two-stage log Gaussian Cox process model that assigns locations based on spatial predictors in the first stage and marks based on nonspatial and/or spatial variables in the second stage. We test the inclusion of dependence structures both within and between stages of the model. The approach of Liang et al. (2008) relies on a categorical mark (or multitype) point process framework. Our two-stage marked point process model allows for both multiple types of spatial variables, including areal and geostatistical, as well as multiple types of marks, including categorical and continuous. • Current models (cf. Liang et al., 2008;Diggle et al., 2010;Pinto Junior et al., 2015) that incorporate nonspatial variables by and large assume that they are available everywhere in the spatial domain by including them in the spatial intensity function. Such an assumption imposes some challenges conceptually as well as in practice. In particular, it becomes challenging to think of a straightforward approach to simulating from such models. We describe how to simulate point process data where the intensity is a function of both spatial and nonspatial variables; this is likely to be broadly useful for practitioners and perhaps also useful in specific applications such as synthetic data creation for preserving privacy (cf. Quick et al., 2015). • We use simulation to understand the interpretability of inference based on our method, and compare it to the bivariate mark modeling approach (Liang et al., 2008). We illustrate that incorporating nonspatial variables into a spatial intensity function imposes challenges for simulation and inference. Our approach avoids using nonspatial variables in this way by separating the location and mark determination stages of the model, where nonspatial variables are only used to determine the mark of the point process. Furthermore, we study how the availability of replicates may impact parameter inference. The paper is organized as follows. In Section 2, "Marked Point Process Models with Spatial and Nonspatial Predictors", we summarize the modeling framework for marked point processes generally. We give detailed specifications of our new two-stage model, including an outline of a simulation algorithm. We also present an existing bivariate marked point process model that we will compare to the two-stage model (Liang et al., 2008). In Section 3, "Applications to Simulated Data", we outline a study of both our new two-stage model and the bivariate mark model to simulated data in order to test the ability of these models to recover simulated parameters. Next, we apply these methods, in a variety of specifications, to real data in Section 4, "Application to Castilla-La Mancha Fire Data." In the discussion section, Section 5, we summarize our results, the benefits and limitations of our model and potential for future work. We relegate to the appendix some of the details of computation and a related model, the bivariate mark model (Liang et al., 2008).

Marked Point Process Models with Spatial and Nonspatial Predictors
We describe the general setup for log Gaussian Cox Process (LGCP) models in order to motivate two models: our new two-stage marked point process model and the bivariate mark model developed by Liang et al. (2008). Point processes are defined by a set of n locations s 1 , ..., s n in a spatial region/window of interest, W . In our case, these n locations are police use of force incidents in Dallas, Texas or forest fires in the Castilla-La Mancha region of Spain. We assume the n events come from a spatial point process with intensity λ θ (s), s ∈ W , where θ consists of unknown parameters such as regression coefficients and covariance function parameters. The intensity may be determined by a vector of spatial covariates at location s, which we notate as z(s). The spatial variables for our study of police use of force in Dallas are collected at the census tract level and points are assigned the census estimate of the census tract where they are located. Spatial variables for the forest fire data are geostatistical, and are provided on a fine grid over the region. Locations of the point process are assigned values for spatial variables of the nearest grid location. The intensity may also be impacted by a vector of individual/event-level or nonspatial covariates at that event, which are denoted ν (cf. Liang et al., 2008).
For a log Gaussian Cox process (cf. Møller and Waagepetersen, 2003), the intensity function λ θ (s) may be defined as follows: log(λ θ (s)) = z(s) β + ν α + ω(s), where θ consists of β and α, which are vectors of unknown spatial and nonspatial regression coefficients respectively. The spatial intensity function is typically defined with spatial variables (z(s)) but has incorporated nonspatial variables (ν) in recent work (cf. Liang et al., 2008). The vector θ also consists of parameters of (ω(s), s ∈ W ), which is a zero-mean Gaussian process. We define the Gaussian process via a positive definite covariance function. Later we provide details about the form of the covariance function we use for each model. For estimation of the Gaussian processes in this paper, we rely on the predictive process approach, as outlined in Section 2.3.
The resulting likelihood function is given below, where the parameter vector θ typically consists of regression coefficients, such as β and α as described above, and parameters of the covariance function. Details of the estimation of the integral are given in Appendix A.2.
We denote the mark at each point process location as Y (s i ), which could be continuous or categorical. When marks are categorical, such as type of force used by police, we call these point processes multitype point processes. We also consider a continuous mark in the case of forest fire data, where we model the amount of burned area in each fire. Existing approaches are often specified specifically for categorical (multitype) point processes (cf. Waagepetersen et al., 2016;Møller et al., 1998;Brix and Møller, 2001;Liang et al., 2008). We develop a framework that is flexible and can easily accommodate multiple kinds of marks, including continuous and categorical (multitype) marks.
We use a Bayesian inferential approach for these models with computation done by Markov Chain Monte Carlo (MCMC). Both the two-stage and the bivariate mark models are fit using the NIMBLE R package, which allows for flexible Bayesian computation (de Valpine et al., 2017). We assess convergence of our MCMC-based approximations through the use of trace plots and effective sample sizes, with the latter based on rigorous Markov Chain Monte Carlo standard error estimates (Haran and Hughes, 2020;Flegal et al., 2021). Detailed prior distributions are given in the description of the two-stage and bivariate mark models. Our MCMC approach uses Metropolis-Hastings adaptive random-walk samplers with univariate normal proposal distributions for the regression coefficients and multivariate normal proposals for the Gaussian processes (de Valpine et al., 2017).

Two-Stage Marked Point Process Model
The two-stage model we develop here specifies a first stage model that determines the locations of the point process events and a second stage which determines the marks at those locations. Two-stage models, sometimes referred to as hurdle models, have been developed in the context of areal data (cf. Ver Hoef and Jansen, 2007;Michaud et al., 2014) as well as point-referenced data (cf. Recta et al., 2012).
In our model, the nonspatial covariates are not used to determine the locations of the events, but rather serve as part of the mark-determination stage. Nonspatial covariates may not be available across the spatial domain but are typically only observed where a point process event occurs, which motivates this model structure. Nonspatial variables may be associated with an event (such as arrest/hospitalization at an incident) or with an individual involved in an event (officer/citizen race/gender). The first stage of this model consists of an LGCP based on spatial variables and a Gaussian process, where the intensity takes the following form: log(λ θ (s)) = z(s) β + ω(s). For the Gaussian process, {ω(s), s ∈ W }, we assume a univariate exponential covariance function so The second stage of the model, in the case of a binary categorical mark, consists of logistic regression to determine the mark, given the locations. The second stage of the model could include any combination of nonspatial and spatial variables. In addition to logistic regression, the second stage of our two-stage model allows for multinomial models, linear regression, and other model types to accommodate marks of many forms. The logistic regression for a binary mark takes the following form, where Y (s i ) is the mark of location s i (either 0 or 1): logit(P (Y (s) = 1)) = z(s) γ + ν α.
Depending on the particulars of the application, spatial predictors can be included or excluded from the second stage of the model. This two-stage model may be more appropriate than alternatives when the relationship between covariates (especially nonspatial covariates) and marks are of interest, rather than the dependence between types of points. Dependence between types of points may be better captured through a bivariate Gaussian process approach as described in Liang et al. (2008).
Our two-stage marked point process framework can include different structures for dependence both within and between the stages of the model, which are shown in Table 1. In the first stage of the model, we investigate the cases where there is no additional spatial structure outside of spatial covariates (Model 1), spatial structure through a univariate Gaussian process (Models 2 and 3), and spatial structure that is dependent on the marking process through a bivariate Gaussian process (Model 4). In the second stage, we study the cases where there is no spatial structure in the second stage of the model (Models 1 and 2), as well as a univariate GP (Model 3) and bivariate GP (Model 4) as in the first stage.
In the most general and flexible formulation of our model, we allow for spatial dependence between the two stages of our model, or Model 4 in Table 1. We assume a Gaussian process for the error in both stages of the model, and then incorporate a dependence structure between the two Gaussian processes. Following Liang et al. (2008), the Gaussian processes {ω 1 (s), s ∈ W } and {ω 2 (s), s ∈ W } are defined such that (ω 1 (s), ω 2 (s)) ∼ M V GP (0, Γ(·, ·, θ)) where the crosscovariance matrix Γ(s i , s j , θ) is of dimension 2n × 2n, where n is the number of points. The cross-covariance matrix consists of both a univariate exponential correlation function, as defined , and the matrix Λ, as defined below. The full cross- Liang et al., 2008). If |ρ| is close to 1, this would indicate strong dependence, and close to 0 would indicate weak dependence.
This formulation leads to the following two-stage model, described in Equation 2, with the first stage determining the locations and the second stage determining the mark of the point process. The dependence introduced by the Gaussian processes allows for dependence in the spatial structure that determines both where points occur and the mark at those point locations.
In our Bayesian hierarchical framework, we assume Normal(0,100) priors for all of the regression coefficients (α, β, γ). We use an Inverse-Gamma(α = 2, β = 0.5) prior for σ, as in Liang et al. (2008). As in other studies (cf. Liang et al., 2008), we fix φ due to identifiability issues. We want the covariance function not involving σ, exp(−|s i − s j |/φ), to be small when points are far apart and large when points are close together. We find the value of φ so that the 95 th percentile of distances would have a a correlation of 0.05, and the value of φ so that the 5 th percentile of distances would have a a correlation of 0.95. We fix φ at the average of these two values.
We fit our model to data simulated from the two-stage model via Algorithm 1. In order to simulate from this model, we assume a spatial window of interest with area a (in the simulated example, it is the Dallas city boundary). For our simulated example, we also assume that W has been broken into areal units and we have spatial variables affiliated with those areal units. We use the Census data available to us from census tracts, but these variables may also be simulated. Spatial data could also be geostatistical, as with the elevation data in the forest fire dataset. We now describe in Algorithm 1 how to simulate from a two-stage marked point process with parameters α 1 , . . . , α p , the regression coefficients for the nonspatial variables, β 1 , . . . , β k , the coefficients for the spatial variables that are used in the first stage, γ 1 , . . . , γ r , the coefficients for the spatial variables that are used in the second stage, and covariance parameter σ > 0. We simulate nonspatial variables from Beta and Bernoulli distributions, but we find that the two-stage model is not sensitive to choice of distribution when other distributions are tested.
First, we calculate the maximum intensity across the space (λ max ) in order to simulate a homogeneous Poisson process based on this maximum intensity, as shown in Figure 3a and 3b. In order to calculate the maximum intensity, there needs to be at least one point per areal unit, or census tract in our case, in order to ensure that we accurately calculate the maximum intensity. Next, we thin the homogeneous points based on the intensity λ θ (s), which in our case only depends on spatial variables, to determine the final locations present in our dataset, which is shown in Figure 3c. There are no restrictions on the distributional assumptions of the nonspatial variables. For instance, unlike the bivariate mark model (Section 2.2), they need not be uniform or bounded. After the thinning procedure, both the marks and nonspatial variables are assigned to the remaining locations, as shown in Figure 3d.
Algorithm 1 Two-Stage Marked Point Process Simulation 1: Define n locations, x 1 , . . . x n ∈ W . † Simulate a Gaussian process, ω(x 1 ), . . . , ω(x n ) with parameters σ 2 , φ. 2: Evaluate the intensity λ θ (s) = z(s) β + ω(s) at x 1 , . . . x n † † and compute the maximum intensity, λ max = max(λ θ (x 1 ) . . . , λ θ (x n )). 3: Simulate a homogeneous Poisson Process with parameter λ max × a on W , where a is the area of W . 4: Affiliate each simulated point with one of x 1 , . . . , x n . ‡ This is to assign a value of ω(s) to each simulated point of the homogeneous Poisson process. 5: Calculate the probability of keeping the point at location s, λ θ (s)/ max(λ θ (s)). Thin the points by keeping each point with this probability, otherwise remove. 6: Simulate nonspatial variables according to their distribution and assign values to each point from the thinned homogeneous Poisson process. 7: Assign the location s to mark 1 with probability, p θ (s), where logit(p θ (s)) = z(s) γ + ν α. Otherwise, assign to mark 0. †: In our approach, we found it helpful to generate x 1 , . . . x n based on a grid of equally spaced points in the observation window, W . † †: When computing the maximum data using spatial variables available on the areal level, at least one point per areal unit is needed. We use a number of points per unit proportional to the area of the areal unit. ‡: Options for affiliating each simulated point with x 1 , . . . , x n include affiliation with the closest point or predictive process transformations.

Bivariate Mark Model
A bivariate marked point process model was developed by Liang et al. (2008) to study two types of cancer and their relationship to each other and with community-level and individual characteristics. The model includes mark-dependent coefficients for both spatial and nonspatial variables. Dependence between categorical marks is introduced through a dependent Gaussian process structure, where {ω 1 (s), s ∈ W } and {ω 2 (s), s ∈ W } follow the dependence structure shown in Section 2.1. Instead of the Gaussian processes being affiliated with each stage as in the two-stage model, the Gaussian processes are affiliated with each level of the mark, or level of force in our simulated example.
In general, the mark-dependent intensity function for marks k = 1, 2, notated λ θk (s, ν), takes the following form. The mark dependent coefficients β k correspond to spatial covariates at location s, z(s). The mark dependent coefficients α k correspond to nonspatial covariates, ν. Finally, {ω k (s), s ∈ W }, k = 1, 2 are the Gaussian processes affiliated with each level of the mark. The Gaussian processes are estimated in a similar fashion to the Gaussian processes described in our two-stage model, except for the cross-covariance and dependence are across levels of the mark in the bivariate mark model.
The likelihood for the bivariate mark model including both spatial and nonspatial variables and for k levels of the mark is given below, where θ k = (α k , β k , σ 1 , σ 2 , ρ):  The intensity for each mark level k, λ θk (s, v), could be specified with or without the Gaussian process ω k (s). The Gaussian process specification allows for the additional flexibility of modeling dependence between the mark levels.
The full simulation algorithm for the bivariate mark model is included in Algorithm 2 in Appendix B. The steps are depicted in Figure 4. Of particular concern, we note in the first step the need to calculate max(λ θk (s, ν)) for each of k mark levels, where λ θk (s, ν) consists of both spatial and nonspatial variables, shown in Figure 4a. This is intuitive for spatial variables but the procedure is not so clear when nonspatial variables are incorporated into the intensity function. In order to calculate λ max for each mark level, we require bounded distributions for the nonspatial variables; this is also required for model inference as shown in detail in Appendix A.1. In summary, in order to calculate the integral of the intensity over a nonspatial variable that is not categorical ( V exp(να)dv) we need to specify an upper and lower bound for the integral. This necessitates bounded distributions. Some choices for bounded distributions include uniform and Beta distributions, and other standard truncated distributions. After calculating max(λ θk (s, ν)) for each mark level k, we again simulate a homogeneous Poisson Process, but this time with a separate point process for each mark level. Next, we simulate nonspatial variables for each point. Finally, we thin the point process by calculating the intensity at each point and thinning points based on their probability in relation to the maximum intensity to obtain the final point process based on both spatial and nonspatial variables, shown in Figure 4d.

Predictive Processes
The Gaussian process, ω(s) needs to be estimated as a finite dimensional random-variable. However, estimating a Gaussian process with a n × n covariance matrix is computationally intensive. Therefore, we use the predictive process approach advocated for by Banerjee et al. (2008). We use a set of knots evenly distributed over the window of interest, as in Liang et al. (2008). These knots are shown in Figure 5 as the blue circles. If s * 1 , ..., s * n are the predictive process knots, and ω * is a realization of the Gaussian process over the knots, then we transform this realization over the knots to be realizations over both the s 1 , ..., s n observed events as well as the integration points, discussed above. The multivariate predictive process realization,w, is calculated as follows:w (s) = cov(ω(s), ω * )var −1 (ω * )ω * For the two-stage model, if m is the number of knots in the predictive process and n is the number of points in the point process, then cov(ω(s), ω * ) is the n × m covariance matrix between the two Gaussian processes. Next, var −1 (ω * ) is the m × m covariance matrix of the Gaussian process over the knots, ω * (Banerjee et al., 2008(Banerjee et al., , 2014. Then the resultingw is the n × 1 mean-zero predictive process. For the bivariate mark model, if k is the number of marks and again m is the number of knots in the predictive process and n is the number of points in the point process, then cov(ω(s), ω * ) is the nk × mk covariance matrix between the two Gaussian processes. The covariance matrix for the bivariate mark model involves both the exponential correlation function as well as the cross-covariance function. As in the two-stage model, var −1 (ω * ) is the mk × mk covariance matrix of the Gaussian process over the knots, ω * . Then the resultingw is the nk × 1 mean-zero predictive process.

Applications to Simulated Data
We studied our modeling and inferential approach via simulated examples. Our simulated data closely mirrors the police use of force data from Dallas, Texas, where we create a marked point process with both spatial, using Dallas census tracts, and nonspatial covariates of multiple types. Given several open questions about the use of force dataset, in Dallas, Texas and otherwise, in this manuscript we restrict our attention to the simulated point process data, resembling the use of force data, and the forest fire application. The areal data is not simulated, and is collected from the 2013-2017 American Communities Survey on the census tract level. For spatial variables, we use one variable that is continuous, such as median age, which we denote z 1 (s) with a corresponding coefficient β 1 . We also use another spatial variable that is bounded between 0 and 1, such as the unemployment rate or the Herfindahl index, which we denote z 2 (s) with a corresponding coefficient of β 2 . For nonspatial variables, we consider one variable that is continuous, to represent a variable such as an officer's tenure on the police force, which is denoted by ν 1 with a corresponding coefficient of α 1 . We also use a nonspatial variable that is binary, to represent an indicator variable such as male officer, which is denoted ν 2 with a corresponding coefficient of α 2 . For both the two-stage and bivariate mark models, we consider both log Gaussian Cox process and nonhomogeneous Poisson process models, or models with and without a Gaussian process respectively, in our simulation framework. For both the two-stage and bivariate mark models we analyze simulated data generated through thinning procedures, as outlined above, to generate point processes with the appropriate intensity function.

Two-Stage Model
In our first simulated example, we consider a simple model that includes spatial variables in the location determination stage of the model and only nonspatial variables in the mark determination stage of the model, and we do not include a Gaussian process in either stage of the model. Therefore, the first stage of the model is a nonhomogeneous Poisson process (NHPP), rather than an LGCP. We find that we are able to recover the logistic regression estimates for the nonspatial variables. The coefficients for the spatial variables are also recovered well when we do not include a Gaussian process.
When we add a Gaussian process to the first stage of the model to form an LGCP, we find that we are able to recover regression parameters associated with spatial and nonspatial variables again reasonably well. However, the estimate of σ from the covariance of the Gaussian process is not as accurate. It is often difficult in LGCP-type models to tease apart information about these parameters from a single realization of the point process. To study this further, and because there are situations where replications are available, we consider the situation where there are multiple realizations, or replicates, of the point process, using the same parameters. In Figure 6, the red, blue, and green lines represent the inferred values using just a single realization, that is, when we consider each of the three replicates separately. We see that they range rather widely, especially for the intercept, β 0 and the covariance parameter, σ. When we repeat this analysis with varied parameter settings, we find similar results where parameters are estimated relatively well and replicates improve our ability to accurately estimate parameters. However, when we infer parameters based on all three replicates at once, we find that we better estimate all the parameters, including σ and β 0 (orange curves in Figure 6). The likelihood function here is simply the product of the three likelihoods in Equation 3 across the individual replicates, assuming independence between replicates.

Two-stage model with cross-correlation
We also simulate from and fit the more general version of the two-stage model shown in Equation 2 to study the impact of the flexible assumption of dependence across the two stages of the model. We test all possible combinations of simulated data and models fitted to that data from Table 1, including fitting simple models to more complex data, and complex models to simple data. This allows us to test whether our modeling procedure shows an improved fit and can recover parameters from the more complicated model where there is dependence between the two stages.
We present the results in Table 2. Model comparison is conducted with the Watanabe-Akaike information criterion (WAIC) which is a Bayesian approach to balancing model fit and model complexity (Watanabe and Opper, 2010). We find that the bivariate and univariate (univariate GP in both stages) models have the lowest WAIC and therefore best model fit, regardless of the form of the simulated data. The WAIC values tend to be similar between the two models. When the second stage of the model does not include a GP (Models 1 and 2  univariate GP in both stages provides the best fit. When the second stage of the model includes a GP (Models 3 and 4), regardless of dependence in the GP, the model with a bivariate GP across stages provides the best fit. In Table 2, we include the simulated case of the bivariate GP model where ρ is 0.9. We also test the comparison of models for different values of ρ. We focus on the case where data is simulated from the two-stage model with dependent Gaussian processes across stages (model 4). We find that WAIC values indicate a similar model fit for both the model with a bivariate GP across stages and the model with a univariate GP in both stages of the model, regardless of the value of ρ that is used. However, we find that these models with spatial dependence (models 3 and 4) are distinguished effectively from models without spatial dependence in one or both stages of the model (models 1 and 2). From this analysis, we suggest that it may be difficult at times to distinguish between models 3 and 4 using WAIC, where in model 4 we add dependence between the two stages of the model, using the ρ parameter. Even if the independent model is indicated as the better fit through WAIC, there may be dependence across stages.
Nonspatial: α1, Mark1 Nonspatial: α2, Mark2 Nonspatial: α1, Mark1 Nonspatial: α2, Mark1 Spatial: β0, Mark2 Spatial: β1, Mark2 Spatial: β2, Mark2 Spatial: β0, Mark1 Spatial: β1, Mark1 Spatial: β2, Mark1  Figure 7: We simulate data from the bivariate mark model, without a Gaussian process, and compare the results for the when the nonspatial variables are simulated from uniform or nonuniform distributions. The posterior densities for all parameters are shown in black when we include nonspatial variables simulated from uniform distributions and blue for nonspatial variables simulated from nonuniform distributions. The true parameters from the simulation are shown in red.

Bivariate Mark Model
We test many different cases of the bivariate mark model with both spatial and nonspatial variables. From our results, we suspect that there are identifiability issues between the nonspatial variables and the intercept, as the intercept is not recovered well when we include nonspatial variables, even when we do not include a Gaussian process. In Step 5 of the bivariate mark model's simulation algorithm as described in Algorithm 2 in the Appendix, we simulate nonspatial variables from pre-specified distributions. We find that we are only able to recover nonspatial variables when they are simulated from uniform/equal probability distributions. As shown in the Appendix in Equation A.1, bounds on continuous distributions are required in order to calculate the integral of the intensity function. However, we find that even when we simulate from a bounded distribution that is not uniform, we cannot successfully recover the nonspatial variable coefficients. When we draw the realizations of the nonspatial variables from a uniform distribution, we are able to effectively recover the nonspatial variable coefficients, but still are unable to recover the intercept. When the nonspatial variables are binary, we also find that the only way to recover nonspatial coefficients is to simulate with equal probability from 0 and 1.
We illustrate our findings in Figure 7, where the models are estimated without a Gaussian process through nonhomogeneous Poisson processes rather than LGCPs. The black line shows the posterior distributions of the nonspatial parameters where the associated variables are simulated from a uniform distribution for the continuous variable (ν 1 ) and a Bernoulli distribution with equal probability on 0 and 1 for the binary variable (ν 2 ). We see that, with the exception of the intercept, the spatial and nonspatial variables are recovered reasonably well in most cases. We note that even though the variables are simulated from these uniform/equal probability distributions, our simulation process relies on a thinning procedure, as described in Algorithm 2 in the Appendix, and the distributions of the nonspatial variables after thinning no longer resemble uniform/equal probability distributions. The blue line in Figure 7 shows the posterior distributions when the nonspatial variables are simulated from non-uniform distributions and Bernoulli distributions with unequal probability. The nonspatial variables are not estimated well, and the spatial variables are recovered similarly well to the case above. We have tested many choices for the prior distribution of the regression coefficients, especially the intercept terms, and these results remained consistent.
We illustrate our LGCP results, with a Gaussian process added to the intensity function, in Figure 8. When we add the bivariate Gaussian process to this model, with cross-covariance Λ, we find once again that we are unable to recover the intercept parameters, β 0 for either mark. We are able to recover both spatial (β) and nonspatial (α) variable coefficients relatively well, when the ν are once again simulated from uniform/equal probability distributions. The main inferential challenge in this method is estimating the parameters of the bivariate Gaussian process. We find that we are able to recover a difference between σ 1 and σ 2 , namely that σ 2 is smaller than σ 1 , but we do not recover the true value well. We do not recover the parameter ρ well, as our point estimate for ρ is around 0.7, while the true value is 0.5.
In summary, these simulated examples have shown that our two-stage model provides accurate and intuitive estimation of parameters in both the location and the mark determination stages of the model. Replicates can be valuable, if available for an application, especially in recovering the parameter(s) of the covariance function. When incorporating dependence across the two stages of our model, it is often difficult to distinguish between the case where the Gaussian processes are independent or dependent using model choice metrics such as WAIC. Therefore, due to the computational costs of estimating the case of dependence across stages, users may choose to estimate the case where there are independent Gaussian processes in both stages. For the bivariate mark model, we have shown through the simulated examples that it is difficult to recover intercept parameters and there are constraints placed on the distributions of nonspatial  variables, given the placement of nonspatial variables in the spatial intensity function.

Application to Castilla-La Mancha Fire Data
Next, we apply the two-stage model to data on forest fires from the Castilla-La Mancha region of Spain, shown in Figure 2. The fire data used to evaluate the model is the clmfires forest fire data from the spatstat package in R (Baddeley et al., 2015). Nonspatial variables with the fire data indicate the cause of the fire and whether the fire occurred in the summer season. Spatial variables include elevation and other characteristics of the land, such as land cover and slope, which are available on a fine grid. We are also interested in characterizing which factors impact the amount of land burned in each forest fire. We consider 3,657 fires from 2004 to 2007. Data collected before 2004 were recorded with less precision than after 2004, so we exclude these cases. We also exclude fires that had "other" listed as the cause of the fire, as we are interested in comparing intentional fires to those caused by lightning or an accident. In the first stage of the model, we test the following spatial variables and their impact on determining the locations of forest fires in the Castilla-La Mancha region of Spain: indicator variable for if the location is in a forest, elevation, and the slope at that location (note that the spatial covariate data may need to be transposed after loading it from the spatstat package). In the second stage of the model, we use the following nonspatial variables: an indicator variable for if the fire was set intentionally (compared to being set by lightning or an accident) and an indicator variable for if the fire occurred in the summer season. We combine the spatial and nonspatial variables listed above to determine their potential impact on the mark of the point process, or the amount of land burned, in hectares, by a given fire. Therefore, this model is slightly modified from the two-stage model shown in Equation 2, where the second stage is now linear regression, where previously it was logistic regression in the case of simulated police use of force incidents. In the second stage, we use the log burned area as the mark and therefore response in the linear regression. This is an advantage of our two-stage model, as it allows the user to consider many different kinds of marks, including continuous and categorical, whereas the bivariate mark model proposed by Liang et al. (2008) allows only for categorical marks and therefore the bivariate mark model cannot be applied to the clmfires data.
Alba-Fernández and Ariza-López (2018) analyze the clmfires data by comparing two different fire seasons per year. The authors find it is important to consider whether fires occur in the summer season during some years, motivating the use of this nonspatial covariate in our model. Myllymäki et al. (2020) use the clmfires data to model local covariate effects. Specifically, they incorporate elevation and an indicator variable for whether the fire occurred in a forest in the spatial model, motivating our use of these variables as well. Dvořák et al. (2020) develop nonparametric methods to study the relationship between points, marks, and covariates. Through the use of test statistics, they conclude that there is dependence between the points, marks, and covariates that they consider (the mark is the area burned and the covariate is elevation). The authors consider only fires in 2007, and when studying the relationship between points and (spatial) covariates alone, they find no relationship between elevation and locations of fires, which will be tested in the first stage of our model.
We begin our analysis of the 3,657 from 2004 to 2007 in the clmfires dataset by conducting preliminary tests of correlation between different combinations of the points, marks, and covariates. We begin with tests based on random shifts, as introduced by Dvořák et al. (2020) to test the dependence between spatial covariates and the points (denoted as a Point-Covariate test, or P-C test). As in Dvořák et al. (2020), we perform 999 shifts with radius up to 150km for this dataset. The region is approximately 450km × 450km. After the variance correction which accounts for points being moved outside of the observation window, the p-values for the elevation, slope, and indicator variable for forest are 0.73, 0.052, and 0.59 respectively. When these tests are bound together by the Bonferroni test, as suggested in Dvořák et al. (2020), we find that none of the spatial covariates impact the distribution of points. This will be further investigated in the first stage of our model, where we have developed a parametric model to study the relationship between all three covariates and the distribution of points.
Next, we investigate the geostatistical marking model, which implies that the point process and marks are independent. We conduct visual tests using the functions E(r) and V (r), defined by Schlather et al. (2004) as the conditional mean and variance of the mark attached to a typical random point, given that there is a further point of the process a distance r away. Given geostatistical marking, the functions E(r) and V (r) would be constant at any distance r. For the fire data, we find that E(r) and V (r) tend to increase with distance r, providing evidence against geostatistical marking. The mark-weighted K function (cf. Guan, 2006) and corresponding simulation envelopes for random labelling (implemented through the Kmark function in R) also show evidence against random labelling, showing dependence between the points and marks. This is consistent with the findings on the subsetted fire data from Dvořák et al. (2020). This rejection of geostatistical marking suggests the need to investigate the spatial dependence between marks and points further. The Gaussian process in the second stage of our two-stage model is a useful tool to do so.
Finally, we conduct a combined test of dependence between the marks and spatial covariates through the random shift test developed by Dvořák et al. (2020). We use the Kendall's correlation coefficient test with variance correction and find the following p-values for elevation, slope and indicator variable for forest: 0.04, 0.014, and 0.002. When the tests are bound together with the Bonferroni correction, we find that with a significance level α = 0.05/3 = 0.017, only the forest indicator variable is significant when analyzing the relationship between marks and the  These tests are powerful exploratory tools to analyze the dependence in the data between points, marks, and the covariates. Our two-stage parametric model allows us to investigate these relationships in more detail. The two-stage model allows us to analyze the relationship between points and covariates in the first stage, covariates and marks in the second stage, and the relationship between all three using the dependent Gaussian processes across stages. It allows us to incorporate unobserved spatial patterns into both the location and mark determination processes through Gaussian processes. We investigate our detailed findings using the two-stage model next.
We fit four models in the two-stage framework to the fire data. The first model considers a non-homogeneous Poisson process (NHPP) in the first (location-determination) stage and standard normal error in the second (mark-determination) stage. In the second model, the first stage has a Gaussian process and the second stage has standard normal error. In the third model, there are univariate Gaussian processes in both stages of the model, with no dependence between them. In the final model, we consider dependent Gaussian processes in the first and second stages. The results for model comparison are presented in Table 3. We find that the lowest WAIC is for the model which includes a univariate GP in both stages of the model, with no dependence between the stages. This is followed by the bivariate Gaussian process model, with dependence between the stages.
Next, we analyze the effects of the variables included in both stages of the model. We show the results for Model 3, or the model with a univariate GP in both stages, in Table 4, and the results for all models in the Appendix, Table C.5. Values in bold in the table are coefficients whose credible intervals do not include 0. For Model 3, which was indicated as the best fit through WAIC, we find that the indicator variable for forest cover and elevation both have negative impacts on the spatial intensity of fires, while the slope has a positive effect on the spatial intensity. This would indicate that at locations with higher elevation, forested regions, and regions with a smaller slope, we are less likely to see fires. These trends are consistent for the first stage of the model in all model types, shown in Table C.5, except for the NHPP where the forest indicator variable is found to have no effect.
In the second stage of the model, we find for Model 3 that the forest indicator variable and slope have negative effects on the amount of area burned, and slope has no effect. For nonspatial variables, we find that the binary variable indicating that the fire was intentional has no effect, and the summer indicator variable has a negative effect. Therefore, we expect smaller fires during the summer and in regions that have forest cover and have higher elevation. The negative impact of elevation on area burned is consistent with exploratory analysis conducted by Dvořák et al. (2020). We find that the results are dependent on whether a Gaussian process is included in the second stage of the model. When we do not include a GP in the second stage (Models 1 and 2), we find that both slope and the indicator variable for whether the fire was intentionally set have impacts on the amount of area burned where these have no effect with a Gaussian Process (Models 3 and 4).
Finally, we consider the dependence parameter, ρ, in the fourth model which includes a   Table 4: We present the posterior mean as well as the 95% Credible Intervals for all covariates in our analysis of the forest fire data in the Castilla-La Mancha region of Spain. Values in bold indicate coefficients whose credible intervals do not contain 0. We present the estimates for Model 3 (univariate GP in both stages) above and the full results in the Appendix (Table C.5).
bivariate Gaussian process. The credible interval for ρ does not include 0, which indicates dependence between the two stages of the model. The ρ parameter is estimated to be negative, indicating negative spatial structure in variables unaccounted for between the location and mark determination processes. Through the WAIC, we have indicated that the model without WAIC is a better fit and our simulation studies show that WAIC may indicate the independent GP model (Model 3) even when there is dependence between the stages. Therefore, further investigation into the spatial structures in the processes that determine where forest fires happen and the severity of these fires is necessary.

Discussion
Through this study, we have demonstrated the utility of the two-stage model in estimating the relationship between spatial and nonspatial variables and marked point process outcomes. Through the two applications discussed in this paper, we have also shown how the model can be applied to multiple types of marks, such as continuous and categorical. Specifically, with the police use of force application, we are interested in nonspatial variables such as officer and citizen race/gender, officer tenure, citizen resistance, and the number of citizens and officers involved. Spatial variables include community factors, such as median income and neighborhood diversity. These spatial and nonspatial variables could be related to the mark, or the level of force that is used by the police officer. Through a simulated example created to closely resemble police use of force data, we have shown that this two-stage model provides interpretable and accurate estimates of relationships between covariates and both the intensity and mark of point processes. In the simulated example, we have shown that replicates can be helpful in recovering parameters, especially the covariance parameters of the Gaussian process. There may be opportunities in some applications to treat segments of the data as replicates, for instance annual data could be treated as replicates. In many cases, however, point processes can be expected to vary over time and space, and therefore it is difficult to consider any of the data as replicates.
For the forest fire data, we are interested in determining which factors play a role in both determining the location of forest fires and the amount of burned area at the fire locations. We show that this model provides a clear interpretation of nonspatial and spatial variables and the relationship with these marked point process outcomes by splitting the model into two stages: location determination and mark determination. We show that the factors impacting the location of fires and the factors impacting the amount of burned area at those locations are not necessarily the same, and accounting for spatial dependence when determining the impact of variables on the mark is important.
In addition to being flexible to many types of marks, our two-stage model could also be expanded to include multiple marks, with an additional stage added to the model with each additional mark of interest. Adding Gaussian processes into the location or additional mark determination stages of the model, as motivated by the application, would add computational complexity but could allow for exploration of the dependence between marks or between individual marks and the location determination process. A variety of marking models can be incorporated into the second stage, including more explicit models for the dependence between the spatial intensity and marking process, such as that of Myllymäki and Penttinen (2009).
We anticipate that this new two-stage model could be widely applicable to other types of datasets. For example, Liang et al. (2008) studies the impact of nonspatial variables on cancer. The authors parameterize their results of the bivariate mark model to show fixed effects on one mark, and differential effects for the other mark, where they hope to show if the variable has a different impact on the two cancers. In the first stage of the model, we would determine if the spatial variables contribute to an increased intensity of the outcome of interest, in Liang et al. (2008)'s case this would be cancer. In the second stage of the model, we could incorporate spatial and/or nonspatial variables to determine the perhaps differential impact of the variables on the different marks, or types of cancer. The dependence structure in the bivariate Gaussian process could be used to test for dependence between the location-determination and the marking processes.
In the future, we would like to apply this model to use of force data provided through the Police Data Initiative and Open Data city websites to conduct an in-depth analysis of the impact of both spatial and nonspatial variables on police use of force, and how these trends compare across cities and jurisdictions.

Acknowledgements
The authors would like to thank Aleksandra Slavković, Ephraim Hanks, and Corina Graif for helpful conversations. This project was supported by Award No. 2020-R2-CX-0033, awarded by the National Institute of Justice, Office of Justice Programs, U.S. Department of Justice. The opinions, findings, and conclusions or recommendations expressed in this publication/program/ exhibition are those of the author(s) and do not necessarily reflect those of the Department of Justice.

Appendix A. Computational Details
Appendix A.1. Integration and Nonspatial Variables The integral presented in Equation 4 cannot be evaluated explicitly and must therefore be approximated. If we omit the mark, the integral to be approximated is as follows: When we consider an intensity function of the form log(λ θ (s)) = z(s) β + ν α + ω(s), where there is no interaction between the spatial and nonspatial variables, then the integral can be separated by spatial and nonspatial variables, shown below.
As in Liang et al. (2008), if a nonspatial variable in ν is continuous, we must specify an upper and lower bound for the covariate, say ν u and ν l , respectively. However, if ν is categorical or discrete, it is less clear how to incorporate these nonspatial variables into the integral. Liang et al. (2008) and Banerjee et al. (2014) both indicate that integrals over discrete variables should be replaced by sums but also that when ν is categorical, than the integration is only over the spatial domain, W (Banerjee et al., 2014, p. 234). In our application, we are interested in nonspatial categorical (binary) variables such as officer gender, citizen gender, and if the officer and the citizen are the same race or not. We find that when we do not incorporate these categorical variables into the integral of the intensity function, the estimates for α are not accurate.
Therefore, we outline how to evaluate this integral while incorporating nonspatial variables. We provide an example where ν 1 is continuous over a bounded distribution from ν 1,l to ν 1,u and ν 2 is a binary indicator variable. The integral over ν 2 is the sum of the function evaluated at 0 and 1, as suggested by Banerjee et al. (2014) where the authors suggest that discrete variables are integrated with sums. This is easily extendable to additional continuous and binary variables.

. Integration Methods
We investigated many different integration methods in this study. It is important to get an accurate measure of the integral of the intensity function in order to accurately conduct inference. First, we find it is important to include at least one point in each areal unit (in our case, census tracts) as the intensity function depends on the evaluation of the function over each areal unit. The spatial variation of the intensity function depends only on the census tracts and the Gaussian process. Therefore, it is important to include more than one point in large census tracts to capture the variation in the Gaussian process. As in Liang et al. (2008), we set the number of integration points per census tracts as a minimum of one and otherwise proportional to the area. We place these points randomly in the given census tract after determining the number of points to be placed in that census tract. These integration points are shown in Figure  5 as the black crosses. We test different integration methods over the region, such as Delaunay triangulation, and we find that a quasi-Monte-Carlo estimate over the integration points produces similar estimates, and we proceed with a quasi-Monte-Carlo estimate for the integral.

Appendix B. Bivariate Mark Model Simulation Algorithm
Below, we include the detailed simulation algorithm for the bivariate marked point process (Algorithm 2 in Appendix B). The simulation algorithm for the two-stage model is shown in Section 2.1. Both algorithms are structured based on the thinning technique, where we start with a homogeneous Poisson Process and we thin to a point process with a specified intensity function. In order to do this, we need to calculate the maximum possible intensity over the region, based on our pre-specified intensity function that we are simulating from. For the two-stage model, the marks are incorporated in the last step of the simulation process. For the bivariate mark model, the intensity involves mark-dependent parameters, so the maximum intensity must be calculated for each mark. Another key difference between the bivariate mark and two-stage simulation procedures is that in the bivariate mark model, the marks and nonspatial variables are simulated before thinning, where in the two-stage models, these variables are assigned after thinning.

Algorithm 2 Bivariate Marked Point Process Model Simulation
Simulate a bivariate Gaussian process, ω 1 (x 1 ), . . . , ω 1 (x n ) and ω 2 (x 1 ), . . . , ω 2 (x n ) with parameters σ 1 , σ 2 , and ρ. 2: Decide the bounds of continuous nonspatial variables from their distributions in order to calculate the maximum contribution to the intensity function by the continuous nonspatial variables, max(ν α). 3: Calculate maximum intensity for each mark, max(λ θ1 (s)) and max(λ θ2 (s)). 4: Simulate a homogeneous Poisson Process for each mark, with parameter max(λ θk (s)) × a for mark k, on W , where a is the area of W . 5: Simulate nonspatial variables according to their distributions and assign values of the nonspatial variables to each point from the unthinned homogeneous Poisson process. 6: Affiliate each simulated point with one of x 1 , . . . , x n ‡ . This is to assign a value of ω 1 (s) and ω 2 (s) to each simulated point. 7: Thin the points by keeping the location s of mark k in the homogeneous Poisson process with probability λ θk (s)/ max(λ θk (s)), otherwise remove. †: In our approach, we found it helpful to generate x 1 , . . . x n based on a grid of equally spaced points in the observation window, W . ‡: Options for affiliating each simulated point with x 1 , . . . , x n include manual affiliation with the closest point or predictive process transformations.
Appendix C. Posterior Mean Estimates, Models 1-4, Forest Fire Data

Integration and Nonspatial Variables
The integral presented in Equation 4 in the main text cannot be evaluated explicitly and must therefore be approximated. If we omit the mark, the integral to be approximated is as follows: When we consider an intensity function of the form log(λ θ (s)) = z(s) β +ν α+ω(s), where there is no interaction between the spatial and nonspatial variables, then the integral can be separated by spatial and nonspatial variables, shown below.
As in ?), if a nonspatial variable in ν is continuous, we must specify an upper and lower bound for the covariate, say ν u and ν l , respectively. However, if ν is categorical or discrete, it is less clear how to incorporate these nonspatial variables into the integral. ?) and ?) both indicate that integrals over discrete variables should be replaced by sums but also that when ν is categorical, than the integration is only over the spatial domain, W . In our application, we are interested in nonspatial categorical (binary) variables such as officer gender, citizen gender, and if the officer and the citizen are the same race or not. We find that when we do not incorporate these categorical variables into the integral, the estimates for α are not accurate.
Therefore, we outline how to evaluate this integral while incorporating nonspatial variables. We provide an example where ν 1 is continuous over a bounded distribution from ν 1,l to ν 1,u and ν 2 is a binary indicator variable. The integral over ν 2 is the sum of the function evaluated at 0 and 1, as suggested by ?) where the authors suggest that discrete variables are integrated with sums. This is easily extendable to additional continuous and binary variables. (1)

Integration Methods
We investigated many different integration methods in this study. It is important to get an accurate measure of the integral of the intensity function in order to accurately conduct inference. First, we find it is important to include at least one point in each areal unit (in our case, census tracts) as the intensity function depends on the evaluation of the function over each areal unit. The spatial variation of the intensity function depends only on the census tracts and the Gaussian process. Therefore, it is important to include more than one point in large census tracts to capture the variation in the Gaussian process. As in ?), we set the number of integration points per census tracts as a minimum of one and otherwise proportional to the area. We place these points randomly in the given census tract after determining the number of points to be placed in that census tract. These integration points are shown in Figure 1 as the black crosses. We test different integration methods over the region, such as Delaunay triangulation, and we find that a quasi-Monte-Carlo estimate over the integration points produces similar estimates, and we proceed with a quasi-Monte-Carlo estimate for the integral.

Predictive Processes
The Gaussian process, ω(s) needs to be estimated as a finite dimensional randomvariable. However, estimating a Gaussian process with a n × n covariance matrix is computationally intensive. Therefore, we use the predictive process approach advocated for by ?). We use a set of knots evenly distributed over the window of interest, as in ?). These knots are shown in Figure 1 as the blue circles. If s * 1 , ..., s * n are the predictive process knots, and ω * is a realization of the Gaussian process over the knots, then we transform this realization over the knots to be realizations over both the s 1 , ..., s n observed events as well as the integration points, discussed above. The multivariate predictive process realization,w, is calculated as follows: w(s) = cov(ω(s), ω * )var −1 (ω * )ω * For the two-stage model, if m is the number of knots in the predictive process and n is the number of points in the point process, then cov(ω(s), ω * ) is the n × m covariance matrix between the two Gaussian processes. Next, var −1 (ω * ) is the m × m covariance Point Type Integration Point Knot

Fig. 1. Integration Points and Predictive Process Knots
matrix of the Gaussian process over the knots, ω * (??). Then the resultingw is the n × 1 mean-zero predictive process.
For the bivariate model, if k is the number of marks and again m is the number of knots in the predictive process and n is the number of points in the point process, then cov(ω(s), ω * ) is the nk × mk covariance matrix between the two Gaussian processes. The covariance matrix for the bivariate model involves both the exponential correlation function as well as the cross-covariance function. As in the two-stage model, var −1 (ω * ) is the mk × mk covariance matrix of the Gaussian process over the knots, ω * . Then the resultingw is the nk × 1 mean-zero predictive process.

Bivariate GP Simulation Algorithm
Below, we include the detailed simulation algorithm for the bivariate marked point process (Algorithm 1). The simulation algorithm for the two-stage model is shown in Section 2.1 of the main paper. Both algorithms are structured based on the thinning technique, where we start with a homogeneous Poisson Process and we thin to a point process with a specified intensity function. In order to do this, we need to calculate the maximum possible intensity over the region, based on our pre-specified intensity function that we are simulating from. For the two-stage model, the marks are incorporated in the last step of the simulation process. For the bivariate model, the intensity involves mark-dependent parameters, so the maximum intensity must be calculated for each mark. Another key difference between the bivariate and two-stage simulation procedures is that in the bivariate model, the marks and nonspatial variables are simulated before thinning, where in the two-stage models, these variables are assigned after thinning.
The bivariate marked point process has parameters α 1 , . . . , α p for each mark k (α k ), the regression coefficients for the nonspatial variables, β 1 , . . . , β k for each mark k (β k ), the coefficients for the spatial variables, and covariance parameters σ 1 , σ 2 , > 0 and ρ ∈ (−1, 1). As shown in Section 3.2 of the main paper, we need to simulate nonspatial variables from uniform distributions/equal probability distributions for the bivariate model.
Algorithm 1 Bivariate Marked Point Process Simulation 1: Define n locations, x 1 , . . . , x n ∈ W † . Simulate a bivariate Gaussian process, ω 1 (x 1 ), . . . , ω 1 (x n ) and ω 2 (x 1 ), . . . , ω 2 (x n ) with parameters σ 1 , σ 2 , and ρ. 2: Decide the bounds of continuous nonspatial variables from their distributions in order to calculate the maximum contribution to the intensity function by the continuous nonspatial variables, max(ν α). 3: Calculate maximum intensity for each mark, max(λ θ1 (s)) and max(λ θ2 (s)). 4: Simulate a homogeneous Poisson Process for each mark, with parameter max(λ θk (s)) × a for mark k, on W , where a is the area of W . 5: Simulate nonspatial variables according to their distributions and assign values of the nonspatial variables to each point from the unthinned homogeneous Poisson process. 6: Affiliate each simulated point with one of x 1 , . . . , x n ‡ . This is to assign a value of ω 1 (s) and ω 2 (s) to each simulated point. 7: Thin the points by keeping the location s of mark k in the homogeneous Poisson process with probability λ θk (s)/ max(λ θk (s)), otherwise remove. †: In our approach, we found it helpful to generate x 1 , . . . x n based on a grid of equally spaced points in the observation window, W . ‡: Options for affiliating each simulated point with x 1 , . . . , x n include manual affiliation with the closest point or predictive process transformations.

Two-stage Model with cross-correlation, Dependence Tests
As discussed in Section 3.1 of the main paper, we simulate data with many different true values of ρ, shown in Table 1. We find the WAIC indicates similar model fit between the model with a bivariate GP and the model with a univariate GP in both stages. The ρ parameter is also consistently overestimated when fitting the bivariate model, regardless of the true value of ρ. Table 2 includes results for models 1-4, as applied to the forest fire data from the Castilla-La Mancha region of Spain.