Evaluating the INLA-SPDE approach for Bayesian modeling of earthquake damages from geolocated cluster data

Modeled damage estimates are an important source of information in the hours to weeks following major earthquake disasters, but often lack sufficient spatial resolution for highlighting specific areas of need. Using damage assessment data from the 2015 Gorkha, Nepal Earthquake, this paper evaluates a Bayesian spatial model (INLA-SPDE) for interpolating geolocated damage survey data onto 1 km grid cells. The proposed approach uses a combination of geospatial covariate data and Gaussian spatial process random effects modeling to estimate the percentage of structures attaining complete damage states from sparse survey clusters. Model performance is evaluated across fifty iterations of 100, 250, and 1000 simulated survey clusters and compared to observed damage assessments and model predictions using more traditional fragility-based methods. Results show strong model fit to observed values, with mean absolute errors of .17, .13, and .11 and correlation coefficients of .75, .82, and .85 for increasing numbers of survey clusters. These results show improvements over traditional damage estimation methods with a small percentage of the damage surveys that were available within several weeks after the Gorkha event. Thus, with sufficient rapid damage assessment mobilization, the proposed model is able to provide improved damage estimates in the time frame required to deliver a Post Disaster Needs Assessment—even in cases where no additional damage data is available.


Introduction
In the hours to weeks following major earthquake disasters, detailed information on the spatial variability, extent, and severity of damages is often sparse. This lack of consistent and verifiable post-disaster information poses major challenges for rapid emergency response efforts (Comfort, Ko, and Zagorecki, 2005;Goodchild and Glennon, 2010;Lallemant, Soden, et al., 2017). Until field surveyors can be mobilized at scale, disaster response decisions are informed by coarse modeled damage estimates, scattered eye-witness reports, and any remotely-sensed damage assessments that might be available (Goodchild and Glennon, 2010;Xie et al., 2016;Lallemant, Soden, et al., 2017;Monfort, Negulescu, and Belvaux, 2019). Understanding how to best leverage and synthesizing these diverse sets of impact data is an active area of research, highlighted by work like Loos et al. (2018) which proposes a spatial integration framework for combining modeled damage estimates, field surveys, remotesensing proxies, and auxiliary shaking estimates into a single estimate of damage. Without a clear standard for collecting, organizing, and utilizing disaster impact data, continuing to develop approaches for modeling earthquake damages with various sources of input data is an important pursuit.
Common model-based approaches for estimating earthquake damages rely on building stock records and functions that describe the probability of damage at various shaking levels for different housing typologies (Whitman et al., 1997;Yeh, Loh, and Tsai, 2006;Kircher, Whitman, and Holmes, 2006;Robinson et al., 2018). These relationships can be derived empirically, semi-empirically, or analytically depending on the availability of relevant local seismological studies (Jaiswal et al., 2011;Porter, 2014). While modeled areal damage estimates are successfully applied for rapid order-of-magnitude impact estimates in systems like PAGER, their accuracy at high spatial resolutions is strongly dependent on the quality of input data (Erdik et al., 2011;Jaiswal et al., 2011;Lallemant, Soden, et al., 2017). On the other hand, field-based engineering assessments provide detailed ground-truthed damage data, but are relatively sparse in the weeks following a major earthquake. Over 60,000 rapid visual assessments were collected in several weeks following the 2015 Gorkha, Nepal Earthquake to inform the post-disaster needs assessment (PDNA), yet still represented less than 10% of affected structures. As a result, the damage statistics that inform requests for disaster aid are often based on modeled estimates. The timeline for delivery of a PDNA to funding stakeholders (approximately a month) is simply too short to fully survey the damages (Lallemant, Soden, et al., 2017).
With significant focus in the Sendai Framework for Disaster Risk Reduction placed on leaving no one behind and accounting for vulnerable populations, improving the capacity of damage models to capture spatial heterogeneities is an important goal. Accordingly, this paper evaluates a spatial interpolation model that uses the Integrated Nested Laplace Approximation (INLA) and Stochastic Partial Differential Equation (SPDE) approaches to estimate the probability of complete structural damage on a high resolution grid from geolocated clusters of surveyed damages and geospatial covariates. The INLA methodology efficiently implement approximate Bayesian inference for a subset of models that can be defined with latent Guassian Markov random fields (Rue, Martino, and Chopin, 2009). The SPDE approach provides a computationally convenient way to implement models with spatial effects using INLA methods (Lindgren, Rue, and Lindström, 2011). The INLA-SPDE approach has been applied to a wide variety of development and demographic contexts (Tatem et al., 2014;Bhatt et al., 2015;Bosco et al., 2017;C. E. Utazi et al., 2018), leveraging the correlations of particular response variables with geographic, environmental, or socio-demographic variables for which higher resolution data is available.
This work is closely related to the spatial integration modeling framework developed in Loos et al. (2018) that combines limited field surveys with a modeled damage estimate, InSAR-based coherence differences, and auxiliary shaking estimates to predict mean damage grades and associated uncertainties at high resolution. Using the 2015 Nepal Earthquake as an example, the authors find improved damage estimates over more traditional modeling approaches-even in cases with low numbers of field surveys. In contrast to this approach, the model proposed herein focuses only relies only on field surveyed damage levels and geospatial covariates, leaving out other types of modeled damage estimates. The information encoded in a traditional fragility damage model (i.e. the distribution of building types and their respective collapse rates at various ground motions) is instead accounted for at the covariate level by selecting data layers that reflect the spatial distribution of different building types. This approach is potentially beneficial in scenarios where fragilities are poorly calibrated but housing types are predominately spatially homogeneous. To evaluate the INLA-SPDE damage model, this study also uses data from the Gorkha Earthquake due to the availability of data on post-event earthquake damage states.
The rest of the paper is organized as follows. Section 2 reviews the fragility-based damage modeling approaches used in many applications. This section provides relevant background for understanding the damages statistics in the Gorkha Earthquake PDNA and motivates the proposed model. Section 3 covers the input data sources used in the model and their relevance to damage estimation. The details of the INLA-SPDE approach and model evaluation are described in section 4, with model results for different numbers of simulated survey clusters presented in section 5. The final section discusses model performance, limitations, and some recommendations for future work.

Motivation
Earthquake-related building damages are a function of both shaking intensity at a given location and the seismic resistance of exposed structures. Fragility curves capture this relationship, specifying the probability of of a structure exceeding a certain damage state at a given shaking level (Porter, 2014). A standard approach for estimating damages combines fragility curves, ground motions, and data on housing type distributions to estimate the percentages of structures attaining specific damage states for a given area (Whitman et al., 1997;Yeh, Loh, and Tsai, 2006;Kircher, Whitman, and Holmes, 2006;Robinson et al., 2018). The damage estimates included in Nepal's PDNA used this type of methodology, drawing average Modified Mercalli Intensity values from the most recent USGS Shakemap, housing counts for four different building typologies from Nepal's 2011 Housing Census, and fragility curves from Guragain (2015) (Government of Nepal, 2015).
A slightly modified version of the PDNA estimate is implemented here to illustrate the performance for the 2015 Gorkha event. The methodological details provided in the PDNA are not sufficient to reproduce the exact estimates. Consequently, instead of estimating both fully damaged and partially damaged structures, only the most severe damage state is estimated (also termed 'complete'). This is assumed to be similar to the 'fully damaged' classification in the PDNA and is the quantity of interest for the spatial interpolation model. Nepal-specific fragility curves from Guragain (2015) are used for stone and brick buildings with mud or concrete mortar, the predominate building types in rural areas. For concrete and wood buildings, fragility curves are drawn from HAZUS-MH (FEMA, 2013). Ground motion and housing typology data are drawn from the same sources as the PDNA, a USGS Shakemap and the 2011 Housing Census, respectively. The PDNA estimates were calculated at the district level, but are reproduced here for village development committees (VDCs)-the finest spatial unit for which housing typology data is available. It is not clear why the PDNA used district-level estimates when VDC level data was available.  The mean absolute error (MAE) and root mean squared error (RMSE) across all VDCs are 42% and 52%, respectively, suggesting that traditional stone and mud structures performed better than estimated in many VDCs. This margin of error might be acceptable for an event-level 'order of magnitude' estimate produced within minutes to hours of an earthquake, but it is not clear whether these modeled estimates should be considered accurate enough to inform aid allocation decisions. Model approaches that are not linked to ground-truth damage assessments are reliant on the accuracy of uncertain ground motions and fragility curves.

The INLA-SPDE Approach
This study proposes a damage model that integrates randomly sampled clusters of post-earthquake damage assessments with gridded geospatial covariates to predict earthquake damages on a uniform grid. Consider a study region A ∈ R 2 discretized into a uniformly spaced grid with n p grid points s 1 , ..., s np . For any given location i in the study region, let Y i represent the number of structures at a specified damage state given N i total structures. In this application, quantities Y i and N i are assumed to be observed for a given set of survey clusters n c with known locations and otherwise unknown at the grid point level. To estimate probabilities p i at the grid point level, the model is defined by: where F i is the vector of covariate values, θ is a vector of hyperparameters, and x is the spatial latent Gaussian field that defines spatial random effects x i . In this application, spatial random effects are specified with a Matérn covariance function defined by a smoothness parameter v, scaling parameter k, and marginal variance σ 2 η (Lindgren, Rue, and Lindström, 2011). A single latent Gaussian field is used for both the observed point-referenced data and the estimations on the prediction grid. Covariate values are defined across the survey area such that values exist for survey and prediction points.
Approximate Bayesian inference via the INLA-SPDE approach is used to fit the model in Equation 1.1 Rue, Martino, and Chopin, 2009. Principally, the INLA approach is an computationally efficient alternative to Markov chain Monte Carlo methods that produces numerical approximations of posterior parameters. For the vector of model hyperparameters θ = (F i , k, σ 2 η ), the joint posterior is given by: where p(θ) is the joint prior distribution of model parameters. The SPDE methodology (Lindgren, Rue, and Lindström, 2011) is used for estimating the Gaussian field (η). This approaches defines a triangular mesh across the study region A using basis functions that provide a sparse representation of a Gaussian field with Matérn covariance at each triangulation node. A projector matrix is then used to linearly map values from triangulation nodes to points of interest inside the mesh. For full details on the SPDE approach, readers are referred to (Lindgren, Rue, and Lindström, 2011).

Geospatial Covariates
A set of gridded geospatial covariates are used to guide interpolation between damage levels at observed survey clusters. Instead of modeling a damage trend from ex-ante secondary damage estimates along the lines of Loos et al. (2018), covariates are chosen to capture ground motion and the spatial variability of housing types with different fragilities. Ground motions are included via a U.S. Geological Survey (USGS) Shakemap for Modified Mercalli Intensity. USGS Shakemaps are produced globally in near-real-time for all major earthquakes. Variability in housing types is reflected across three covariate layers: Shuttle Radar Topography Mission-derived elevation, Visible Infrared Imaging Radiometer Suite (VIIRS) nighttime lights, and distance to Open-Street Map major roads. These data are freely available from the WorldPop Project (WorldPop, 2018) and were selected based on the basis of domain knowledge about social, structural, and physiographic variation in Nepal (Muzzini and Aparicio, 2013;Chaulagain et al., 2015;Gautam and Chaulagain, 2016;Robinson et al., 2018). All Worldpop covariate values are re-sampled onto 1 km 2 grid cells to match the resolution of the USGS Shakemap.
VIIRS nighttime lights data is used as a proxy for urbanization. A majority of engineered structures and other reinforced concrete buildings in Nepal are located in urbanized areas (Gautam, Prajapati, et al., 2016;Gautam and Chaulagain, 2016). Although population density could also be used to similar effect, nighttime lights provides a more direct signal for urbanization of the built environment. Distance to major roads is used as a proxy for remoteness and access to non-local building materials. Large portions of rural Nepal are only accessible on foot and are restricted in the types of building materials available (Muzzini and Aparicio, 2013). Most houses that are far from roads are traditionally constructed with stones and mud mortar. Elevation captures two potentially important variations in building materials. At coarse scales, elevation reflects differences between the predominate building types in the low-lying Terai region (wood and bamboo) versus the mid-hills region (stone with mud/cement mortar). At finer scales, there is a general tendency for reinforced concrete and other engineer structures to be located at lower relative elevations than more traditional construction built into terraced hillslopes. Together, these three covariates capture the dominant features distinguishing housing typologies in Nepal.

Post-Earthquake Damage Estimation from Simulated Survey Clusters
As previously mentioned, in the month following the Gorkha earthquake over 60,000 rapid visual building inspections were performed by volunteer engineers and architects trained by Nepal's National Society for Earthquake Technology (NSET) (Lallemant, Soden, et al., 2017). While these original data are not openly available, a similar type of dataset can be simulated from damage surveys collected further along in the reconstruction process. In this paper, the proposed model is evaluated with simulated datasets containing increasing numbers of surveyed damage clusters (100, 500, and 1000) drawn from a complete set of damage assessments collected for all households across eleven rural districts as part of Nepal's Household Registration for Housing Reconstruction program. While detailed engineering assessments typically assign ordinal damage grades (Lallemant and Kiremidjian, 2015), PDNAs use a less detailed partial-complete damage spectrum (Government of Nepal, 2015). Accordingly, this study focuses only on predicting 'complete damage states', assumed to be those structures assigned the maximum possible damage grade (5 out of 5) in the existing surveys. These structures are priorities for emergency response and often require full reconstruction as opposed to repairs. The cluster simulation procedure is described as follows. Cluster centroids are first drawn as random samples of all household locations in the study area. Each cluster is then assigned a fixed buffer of 250 meters and all households within the buffered zone are assumed to be surveyed irrespective of damage level. The decision to use a 250 meter buffer distance is somewhat arbitrary, but is chosen to be small enough such that a cluster could be surveyed by a small team of engineers reasonably quickly. Other sized clusters or even variable size clusters are usable as long as every structure in a cluster is surveyed. 1,000 surveyed clusters at 250 meters each includes approximately 55,000 households-close to the number of rapid damage assessments that were actually collected in a few weeks following the earthquake. 500 and 100 surveyed clusters cover approximately 28,000 and 5,500 households on average, respectively.
An iterative model fitting procedure with fifty instances for each number of survey clusters is implemented using the R-INLA package in R (Lindgren and Rue, 2015). Iteratively fitting the model to different randomly sampled cluster locations averages out the bias associated with predictions at any single cluster location to give a stronger overall picture of average model performance. For each iterations, a cluster set is generated and covariate values are extracted at the centroid locations. These data are used to fit Equation (1) with a SPDE triangular mesh constructed as a convex hull around the study area (see details in Appendix). Penalized complexity priors with the range set to the median distance between grid points are specified on the spatial random field (Simpson et al., 2017;Fuglstad et al., 2019) and weakly informative normal priors (∼ N (0, 1e5) are placed on the covariate coefficients. Predictive performance is evaluated across all iterations for each number of survey clusters with a combination of three metrics: root mean square error, mean absolute error, correlations between observations and predictions at the grid level. Spatial plots for the median posterior estimates, standard deviations, and differences between predicted and observed values are also included. Figure 3A-E shows the median posterior estimates and standard deviations for the percentage of buildings sustaining complete damage in each 1 km 2 grid cell for 250, 500, and 1000 clusters. The corresponding parameter estimates are included in Table 1 in the Appendix. Damage patterns are positively correlated with MMI and distance to roads and negatively correlated with VIIRS nighttime lights, supporting prior assumptions for the covariates. Elevation switches from positive to negative correlations at 500 grid clusters, likely reflecting the lack of damage in the high mountain regions with no surveyed structures. Among included variables, shaking intensity has the strongest effect on modeled damage percentages. The estimates for the spatial range parameter vary between 1.18 and 1.48 (corresponding to an approximate distance of 11.2 to 15.2 kilometers), indicating strong spatial dependence in damage levels. These ranges are similar to spatial correlations found in previous work (Jayaram and Baker, 2009;Shome, Jayaram, and Rahnama, 2012).

Results
As expected, increasing the number of survey clusters increases the level of spatial heterogeneity in the modeled estimates and decreases the estimated uncertainties. Moving from 100 to 500 clusters has a larger impact on model estimates than the corresponding increase from 500 to 1000 clusters. This is predominately seen as reduction in uncertainty rather than a significant change among estimated damage levels. Uncertainties remain high in the Himalayas where few if any households are located. Compared to the areal damage model (Figure 1), all three cluster models more accurately capture the lower percentages of complete damage at moderate latitudes. However, it is worth noting that the model predications and uncertainties from any single scenario will be affected to some degree by the survey cluster locations, with higher uncertainties further away from surveyed clusters. Optimizing cluster placement for interpolation accuracy is beyond the scope of this study, but could be a worthwhile direction in future work.
The mean absolute errors for the 100, 500, and 1000 cluster models are . 17,.13,and .11,respectively (.23,.18,.15 RMSE). Correlations between predicated and observed damages (see Figure 6 in Appendix) are .71, .82 and .85. The MAE values imply an average difference between observed and predicted values at the grid cell level of 11-17%-a 25%+ improvement over the areal model with improved spatial resolution. As mentioned previously, model performance improves with additional survey clusters but reasonably accurate predictions can be obtained with a modest number of surveyed clusters. The average number of households included amongst the simulated 100 cluster models is less than 10% of the actual number of households surveyed following the Gorkha Earthquake. Figure 4 shows the percent differences between predicted and observed values at the grid cell level. On the whole, all three models do a reasonable job of reproducing the observed damage levels without any major biases. However, the 100 cluster model shows more correlated errors compared to the 500 and 1000 cluster models. Including additional geospatial covariates or altering the included covariates could potentially improve model fit, although errors may also stem from ground motion uncertainties. USGS Shakemaps are themselves a modeled data product and often interpolate across large areas with no strong motion observations. Figure 4 also shows that all thee models contain cells where estimates are off by a significant margin (50-80%). These prediction errors are concentrated in cells that neighbor survey observations but have large differences in damage levels. The spatial random field favors smooth transitions between observations and is not well suited for capturing rapid changes between neighboring pixels. Similar issues may also occur at sharp boundaries in the covariate layers. Although major prediction errors only exist in a small fraction of total model cells, they should be taken into account when interpreting results. Pixel to pixel differences are generally reliable, but the proposed model should not be used for targeting individual cells.

Discussion & Conclusion
Providing accurate models of earthquake damages on short-term timelines is a critical element for both disaster response and recovery planning (Lallemant, Soden, et al., 2017;Monfort, Negulescu, and Belvaux, 2019). Requests for international  aid following major earthquakes rely on the best damage assessments available approximately one month post-event. It is common for rapid damage assessments to be collected at varying locations across impacted areas to inform a post disaster needs assessment, but these data have not historically been used in damage models. In line with recent research focused on integrating disparate post-disaster datasets (Loos et al., 2018), this paper evaluated the a spatial interpolation model using the INLA-SPDE approach for estimating complete damage states from geolocated cluster level-data. The proposed model links surveyed damage states to gridded covariates that correlate with ground motions and building fragilities. Results show improved performance over more traditional fragility based approaches with as few as 100 clusters or approximately 5,500 surveyed structures. However, the true value of a spatial interpolation framework becomes apparent with closer to 500 or 1000 clusters (28,000 -55,000 surveyed structures) where predictions capture detailed spatial heterogeneities and correlations between observed and predicted damage maps reach 85%.
One of the largest advantages of the particular spatial interpolation approach proposed in this paper is the accessibility of the input data sources. The model only requires geolocated survey clusters, a USGS Shakemap, and globally available gridded covariates-no additional damage estimates are needed. Using observed cluster-level damage states as the sole starting point circumvents issues with poorly constrained fragility curves by making fewer assumptions on the a-priori relationship between ground motions and specific housing types. Rather, a more general statistical relationship is derived from observed damage levels across different local geographies with a spatial random field term capturing spatially correlated errors. In essence, the spatial interpolation model trades off an unknown assumption on fragilities for a known assumption on the distribution of different structural types across landscape characteristics. This particular approach works well in places like Nepal with strong socio-environmental variation among housing types. Although previous studies (e.g. Chaulagain et al., 2016;Robinson et al., 2018) use anywhere from 4-7 housing categories, the differences in associated fragilities between many of the included typologies is quite small. Hence, focusing on distinguishing between reinforced concrete, wood, and various stone/mud structures captures most of the fragility variability in Nepal at a one kilometer pixel resolution. Similar assumptions may not be reasonable in regions with more heterogeneous housing types, less spatial consistency among housing types, or weaker correlations between housing types and physiographic variables. In these situations, or cases where other damage proxies (remote-sensing, etc.) are readily available, the Loos et al. (2018) framework may be more useful for synthesizing datasets into a single damage estimate.
Another limitation in this paper, like similar studies, is the use of comprehensive damage assessment data collected after a major earthquake for validation. Other countries may not have comparable datasets for performing similar simulation studies. While this data is not strictly necessary to adapt and use the proposed model for future events, appropriate caution should be placed on results from areas where the model has not been validated. In a similar vein, the results from this study assume that survey clusters are randomly spread throughout the affected area. The decision to use randomly sampled locations was based off the assumption that systematically designing a survey sampling scheme is beyond the scope of the first few weeks of disaster response activities. In reality, there is a strong chance that rapid damage assessments might be preferentially sampled in areas with easy access, places where qualified engineers are already present, or locations with prior knowledge of severe impacts. Previous research has found the location of field surveys to impact model bias (Loos et al., 2018) and model performance in any of these scenarios may be worse than presented here. Adapting the model to handle preferential sampling is a clear next step for this research.
Despite several limitations, this study offers a new use case for the INLA-SPDE methodology and contributes to a growing set of next-generation damage modeling approaches. There are several promising directions for improving upon this particular approach. Most notably, this study focused on a binomial response of a single damage state. A more complex version of the model could incorporate joint likelihoods on several ordinal damage states to provide more comprehensive impact estimates. This improvement would require more detailed rapid damage assessments, but is likely to be useful even along a 'no damage-partial damage-complete damage' spectrum. The existing model could also be extended to fatality modeling. Fatalities are generally derived as the percentage of population in collapsed buildings which is simply further subset of complete damage states (Kircher, Whitman, and Holmes, 2006;Robinson et al., 2018). A procedure for selecting fatality rates at the grid level would need to be developed, but otherwise the approach remains the same. Finally, model extensions that directly incorporate data on housing types from censuses or other sources could be considered. Including percentages of different housing types as an areal covariate similar to C. Utazi et al. (2018) is a potential first step in this direction.

Damage Assessment Data
The damage assessment data used in this study comes from Nepal's Household Registration for Housing Reconstruction Program. In an effort spearheaded by Kathmandu Living Labs and Nepal's National Planning Commission, door to door structural damage assessments were collected at every residential household across the eleven most affected rural districts. The observed percentage of structures assessed at damage grade five are shown below in Figure 6. An open access version of this dataset (and the survey questions) is available at: https://eq2015.npc.gov.np/

Mesh Construction
The SPDE approach requires constructing a triangulated mesh to model the spatial random field. The triangle knots serve as integration points and the values of any point lying within a triangle are interpolated from the knot estimates. Thus, a finer mesh produces more accurate predictions but increases computational time. The maximum triangle edge length for this study was set at .05 degrees or approximately 2% of the prediction grid region, similar to discretizations in previous work C. Utazi et al. (2018).