A Promising Tool to Assess Long Term Public Health Effects of Natural Disasters: Combining Routine Health Survey Data and Geographic Information Systems to Assess Stunting after the 2001 Earthquake in Peru

Background Research on long-term health effects of earthquakes is scarce, especially in low- and middle-income countries, which are disproportionately affected by disasters. To date, progress in this area has been hampered by the lack of tools to accurately measure these effects. Here, we explored whether long-term public health effects of earthquakes can be assessed using a combination of readily available data sources on public health and geographic distribution of seismic activity. Methods We used childhood stunting as a proxy for public health effects. Data on stunting were attained from Demographic and Health Surveys. Earthquake data were obtained from U.S. Geological Survey’s ShakeMaps, geographic information system-based maps that divide earthquake affected areas into different shaking intensity zones. We combined these two data sources to categorize the surveyed children into different earthquake exposure groups, based on how much their area of residence was affected by the earthquake. We assessed the feasibility of the approach using a real earthquake case – an 8.4 magnitude earthquake that hit southern Peru in 2001. Results and conclusions Our results indicate that the combination of health survey data and disaster data may offer a readily accessible and accurate method for determining the long-term public health consequences of a natural disaster. Our work allowed us to make pre- and post- earthquake comparisons of stunting, an important indicator of the well-being of a society, as well as comparisons between populations with different levels of exposure to the earthquake. Furthermore, the detailed GIS based data provided a precise and objective definition of earthquake exposure. Our approach should be considered in future public health and disaster research exploring the long-term effects of earthquakes and potentially other natural disasters.


Introduction
Earthquakes are natural disasters that impact health and society. During the past two decades alone, earthquakes have caused nearly one million deaths globally [1]. Other severe consequences of earthquakes include injuries, loss of assets and livelihood, as well as disruption of infrastructure and societal functions [2].
Research into the health effects of earthquakes has overwhelmingly focused on immediate consequences, such as on injury and mortality. Furthermore, most studies concentrate on high-income settings [3], although the risk of death and damage from earthquakes is significantly higher in low and middle-income countries [4]. It has also been reported that in these lower income settings, sudden shocks have greater long-term negative impacts on society and development than in higher income regions [5,6]. However, such prolonged effects are difficult to define and quantify and, thus, have not been well studied.
One way to quantify long-term effects of earthquakes could be to monitor changes in important public health indicators. Such indicators often reflect living conditions, socioeconomic status and other factors that affect the general well-being of a population, and that can be affected by an earthquake. Furthermore, such indicators are often routinely collected in national health surveys. For example, Zolala et al considered maternal mortality a useful indicator to measure long-term health effects of the 2003 Bam earthquake, since maternal deaths were likely to be affected by the earthquake's impact on women's living conditions and access to health care [7]. Childhood stunting is another important public health indicator. Stunting correlates with the degree of food insecurity, poverty and disease, and could thus indicate long term effects on public health [8,9]. Previous studies have found high levels of stunting in areas affected by earthquakes [10]. Data on childhood stunting is regularly collected in the Demographic and Health Surveys (DHS), which is an important source of health information in lowand middle-income countries. In most recent DHS surveys, the surveyed clusters are geo-referenced, which may allow levels of stunting to be tracked according to the geographical distribution of a disaster such as an earthquake [11].
Data on the distribution of earthquakes and other natural disasters have also recently become more accessible due to the increasing use of geographic Information Systems (GIS) [12]. For example, the U.S Geological Survey's (USGS) ShakeMap atlas is a rich source of GIS data for earthquakes [13]. ShakeMaps divide earthquake-affected areas into different zones of seismic intensity. Seismic intensity represents the level of ground shaking in a given location and is closely related to the level of human impact and destruction [14]. Combined with data on population distribution, the population affected by the disaster can be estimated with high precision [15].
In this study we explored whether long-term public health effects of earthquakes can be assessed using a combination of readily available data sources on public health and geographic distribution of seismic activity.

Selection of earthquake case
We tested the approach using a real case-an earthquake that hit southern Peru on June 23rd 2001. Three criteria were used to select the case. First, DHS data had to be available both before and after the earthquake, with significant time (at least two years) elapsing between the earthquake and the last measuring-point. Second, the ShakeMap atlas had to include GIS-based intensity data for that specific earthquake. Finally, the clusters surveyed in DHS had to be georeferenced. The latter two criteria were deemed necessary to geographically link households to the corresponding earthquake intensity level.

Peru and the earthquake
Peru is a middle income country [16]. Rapid economic growth over recent years has yielded a decline in poverty, but economic disparities within the country remain high [17,18]. The majority, 91 percent, of the population is ethnic Spanish, followed by the indigenous ethnic group Quechua, which constitutes 7 percent of the population [19]. Stunting remains a significant health problem in Peru, especially in poor highland regions. In 2006, the prevalence ranged from 9% in the richest region, to 57% in the poorest [20]. Stunting prevalence is higher among indigenous groups and in rural areas [21].
Peru lies on an area of intense seismic activity and is frequently affected by earthquakes. The earthquake investigated in this study struck outside the southern coast and measured 8.4 on the Richter scale (see Fig 1). Over 300,000 people were affected, and 145 were killed [1]. The earthquake caused devastation in the affected regions, especially to adobe built houses [22]. According to the EM-DAT disaster database, four administrative regions were affected-Arequipa, Moquegua, Tacna, and Ayacucho [1].

Data sources
Two DHS surveys were considered for our study: one standard DHS [23] and one continuous DHS [19,24]. The standard DHS was conducted in 2000 (one year before the earthquake) and sampled 1,414 randomly selected clusters in different parts of Peru. In the continuous DHS, five equally-sized subsamples of these 1,414 clusters were surveyed in five separate rounds over five consecutive years (2004)(2005)(2006)(2007)(2008). For our study, this implied that most clusters surveyed before the earthquake were surveyed again a few years after. Anthropometric indicators (weight and height measurements) were collected for children under five years of age in the surveyed households. In the continuous DHS survey, anthropometric data was collected in the 2005, 2007, and 2008 rounds of the survey. The total response rate was 95% in the 2000 standard survey [23] and 98% in the continuous survey [19]. About 90% of children under five years of age were measured for height in both surveys [19,23,24].
The surveyed clusters were geo-referenced, and data files with the geographical location (GPS) of the clusters were provided in addition to the survey data sets. To assure respondents' confidentiality, the GPS positions were randomly displaced up to two kilometres in urban areas and ten kilometres in rural areas [11].
In USGS ShakeMaps, seismic intensity levels are derived instrumentally, by combining instrumental measurements with information about local seismic conditions and the earthquake source [25]. Seismic intensity is expressed in the Modified Mercalli Intensity Scale (MMI), which ranges from I (shaking not felt, no potential damage) to X (extreme perceived shaking, very heavy potential damage). The 2001 Peru earthquake produced intensities up to VIII MMI (see Fig 1).
We used a GIS software program (ArcGIS version 10.1) to link each surveyed cluster to the corresponding MMI value in the area. The MMI zone map and the cluster GPS-point map were overlaid and joined using the software's spatial join tool. This created a data set that contained the unique cluster numbers and the corresponding MMI values. This data set was subsequently merged with the full DHS data set, using the cluster number as the common variable.

Ethical consideration
This study was based on secondary analysis, using data from publically available databases (DHS and USGS Shakemap archive). The DHS surveys were approved by the the Peruvian National Institute for Statistics (INEI) and the ICF Macro at Calverton in the USA. Data were collected with informed consent and anonymity of respondents was assured. Geographical data did not disclose the location of individual households [11].

Variables
Stunting was defined as a height-for-age more than two standard deviations below the median of the World Health Organization (WHO) growth reference population [26]. Biologically implausible z-scores (z-scores below -6 or above 6, according to WHO Child Growths Standard recommendations [26]) were recoded as missing data.
The MMI values were categorized into three intensity level zones. MMI values I through IV, were not portrayed on the intensity map and were therefore collapsed into the low intensity zone. The remaining observations were divided into medium and high intensity zone by the median MMI value. As a result, the medium intensity zone included MMI value V and the high exposure zone values VI through VIII.
Children were categorized into pre-and post-earthquake groups according to their year of birth. This categorization minimized the risk of including prevalent cases of stunting in the post-earthquake group, i.e. children who were already stunted before the earthquake. Children born before the earthquake (1995-2000) comprised the pre-earthquake reference group. The two post-earthquake groups consisted of children born shortly after (2001)(2002)(2003), and a few years after (2004-2008) the earthquake respectively. The cut-offs between these two groups was defined by the median birth year.
Other covariates collected from the DHS survey data sets included age of the child (0, 1, 2, 3 and 4 years of age), motherhe child (0,al level (no education, primary, secondary, higher), and availability of household electricity (available, not available). The latter two were considered proxies for socioeconomic status.

Statistical analysis
We first conducted descriptive analysis by comparing the prevalence of stunting between the birth year groups and, in turn, comparing these differences between the intensity zones. The descriptive analysis was restricted to children over two years of age, as few incident cases occur after this age [27]. Multilevel logistic regression analysis was subsequently performed to account for potential cluster effects (within clusters and within families). We defined three separate regression models-one for each intensity zone-with birth year group as independent variable and stunting as dependent variable. Model 1 adjusted for age of the child whereas Model 2, in addition, adjusted for the level of the motherof for istic regavailability of household electricity, both proxy indicators for socioeconomic status. All statistical analyses were performed in STATA, version 12.0. P-values less than 0.05 were considered statistically significant.

Results
Using the results of the spatial join between the ShakeMap and the DHS clusters, we developed a GIS-based map that portrays the earthquake affected southern part of Peru divided into the three intensity zones, with the triangle shaped points marking the locations of the surveyed clusters (Fig 2). This map revealed that the intensity zones cut across administrative borders.
For example, the Ayacucho region was exposed to all three levels of intensity. The high and medium intensity zones cover a relatively small area by the coast.
The median and high intensity zones included just over 1,000 children each (Table 1). Over 11,000 of the surveyed children lived in the low intensity zone, corresponding to more than 80% of the total sample. The majority of the children in the sample were born before the earthquake, between years 1995 and 2000.
We determined that the high intensity zone was mainly urban and the majority of the population was of Spanish ethnicity ( Table 2). This zone also had significantly better living standards, higher educational level and better maternal nutritional status (height) than the low and medium intensity zones.
Due to these large pre-earthquake differences in socioeconomic indicators and other important risk factors of stunting, we stratified our analyses by intensity zone. We found that in the high intensity zone 14% of the children born before the earthquake were stunted (Fig 3). The corresponding proportions were 42% in the low intensity zone, and 50% in the medium intensity zone. In the both low and medium intensity zones, stunting prevalence was lower among children born shortly after the earthquake (2001-2003) and even lower in the last birth year group (children born 2004-2008). In the high intensity zone, however, the prevalence increased to 20% for children born shortly after the earthquake. In the last birth year group, the prevalence decreased but remained higher than before the earthquake. When we investigated the association between birth year group and stunting in the multilevel analysis and adjusted for age of the child, the odds of stunting were not significantly different between any of the birth year groups. The exception was the low intensity zone, in which the odds of stunting among children born several year after the earthquake (2004-2008) was 20% lower than in the pre-earthquake reference group (Table 3, Model 1). However, in the second model, in which we also adjusted for the socioeconomic indicators, children born in the high intensity zone shortly after the earthquake had significantly higher odds of stunting than children born in the same zone before the earthquake (OR 2.01 95% CI 1.03-3.90). For children in the third birth year group, stunting decreased in comparison to the second birth year groups, and its rate was not significantly different from that in the pre-earthquake reference group, suggesting that by this time, sufficient economic and societal recovery had occurred to support healthy development.

Usefulness of the approach
Our results indicate that long-term public health effects of earthquakes, more specifically effects on childhood stunting, may be evaluated using a combination of ShakeMap and DHS data. Both data sources are readily available online, and may be combined using the geographical coordinates of the surveyed clusters. The repeated measuring-points of DHS surveys allow comparing stunting levels before and after an earthquake, a prerequisite to draw conclusions about its effects. Follow-up data was available for up to seven years after the earthquake, which gives a unique opportunity to study long-term effects. Most previous information about public health effects of earthquakes come from studies with short-term follow up, and/or  without baseline data [28,29]. DHS data also included several socio-economic indicators, which allowed us to control for decreasing trends in stunting driven by socioeconomic improvements. Our work shows that the level of stunting in Peru was elevated among children born after the 2001 earthquake only in areas exposed to high seismic intensity levels. These results are consistent with previous findings indicating that the societal impact of earthquakes swells with increasing levels of intensity [6]. The fact that the level of stunting is lower among children born several years after the earthquake than among those born immediately after the quake, correlates with sufficient time for recovery, both economic and societal, following the disaster.
Our results underline the importance of GIS-based disaster data for assessing effects of a disaster on a given population with high precision. As the intensity levels of the earthquake cut across administrative borders, a division based on administrative boundaries would have resulted in misclassification of earthquake exposure. Our results therefore confirm previous findings that GIS data provide a more precise estimation of disaster affected populations than when administrative borders are used to demarcate disaster zones [15].
We identified a few factors that limit the usefulness of the approach. Firstly, DHS samples are not designated to be broken down into intensity zones. As the sample sizes in the medium and high intensity zones were relatively small in our Peruvian model, they yielded decreased statistical power, especially in the high intensity zone, where stunting prevalence was low compared to the other zones. The same issue has been discussed previously in a Tsunami Recovery Impact Assessment and Monitoring System (TRIAMS) workshop, where the use of national household survey data was recommended to follow the recovery process of the 2004 Indian Ocean tsunami [30]. A suggested solution was to oversample the population in disaster affected areas. Another recommendation was to modify routine health surveys to provide representative statistics of the geographical area affected by the disaster.
Another disadvantage was the large difference between the earthquake intensity zones at "baseline" (i.e. before the earthquake hit). The area most affected by the quake was mainly urban, had relatively high socioeconomic status, and significantly lower levels of stunting than the rest of the country. According to Ukoumunne et al., such baseline imbalance might bias results of effect evaluation studies, because changes over time are dependent on the levels at baseline [31]. For example, the level of stunting might be more likely to decrease in a population with high initial levels, than in a population where levels are already low. This could be an alternative explanation of why stunting decreased only in the low and medium intensity zone. It is, however, less conceivable that this would explain the increased level of stunting in the high intensity zone. Going forward, it will be important to study the effect of baseline imbalance between differently affected populations, as well as to develop analysis methods to account for it.
The cross-sectional design of DHS is another potential weakness of the approach. Large migration from the earthquake-affected areas might change population composition and bias prevalence estimates. For example, if well-off families migrate to a larger extent than poorer families, an artificial increase of stunting would occur in the remaining population. To rule out this possibility, the same children would have to be followed over time in a longitudinal study. To assess the validity of results based on cross-sectional data, it will be important to study patterns of population displacement from affected areas.

Study limitations
A limitation of our study is that we did not have a gold standard of comparison to know if the approach generated valid results. The lack of previous research also precludes comparing our results to other studies. The social and economic impacts of earthquakes are well-known, at least in the short-term, as is the link between socioeconomic factors and stunting. However, a direct link between earthquakes and stunting has not been reported in previous literature. Nevertheless, the main aim of this study was to test the feasibility of the approach, not to determine the validity of the results.
Our study only evaluated a single earthquake, and some of our assumptions and results may not be applicable to other locales and earthquakes. In addition, the availability of DHS and Sha-keMap data differs from case to case. For example, not all DHS surveys collect GPS data, which is necessary to link the respondents to the respective seismic intensity zones. Considering the wealth of both DHS and ShakeMap data, however, it is likely that our approach can be applied to other earthquakes and potentially even other natural disasters, using other sources of GIS disaster data. For example, Guha-Sapir et al. [15] reported on several initiatives that have made GIS-based flood data available for disaster research.

Conclusions
Long-term public health effects of earthquakes are understudied in disaster research, which may be attributed to the lack of systematically collected health data in affected regions. Our study suggests that a combination of readily available data sources can be used to increase knowledge about such long-term effects. Although limitations exist, given the wealth, easy access, and good quality of DHS and ShakeMap data, it is an approach worth investigating further and considering in future public health and disaster research.