Quantifying the dynamics of migration after Hurricane Maria in Puerto Rico

We compared four independent datasets to estimate population changes over the course of a year after the hurricane. First, we obtained the intercensal yearly population estimates provided by the American Community Survey (ACS) for 2010 to 2018, a yearly survey conducted by the US Census Bureau. The ACS estimates are for 1 July of each year. We considered this estimate to be our gold standard.

Second, we extracted data from Disaster Maps provided by Facebook’s Data for Good team. Specifically, for a group of people determined to be living in Puerto Rico the week before Hurricane Maria, weekly estimates from 21 August 2017 to 30 July 2018 on the number of these users residing in 77 of the 78 municipalities on the island were available. Note that these data are a closed cohort of known individuals, so no new users appear in the data, and there is a constant rate of attrition of users expected. However, this is also the only dataset on a municipality, rather than an island-wide, level.

Third, Teralytics provided island-wide daily population proportion estimates from 31 May 2017 to 30 April 2018 based on cell phone usage patterns, analyzed in partnership with an unnamed mobile operator on the island. Since we do not know which operator the data are from, it was impossible to assess the geographic bias in ownership in this group. The proportions were calculated relative to a baseline population determined by Teralytics using a method unavailable to the authors. Due to the unreliability of data for the 4 wk after Hurricane Maria, presumably due primarily to low connectivity, proportion estimates were not provided for this time period.

Finally, we obtained Airline Passenger Traffic (APT) data from the US Bureau of Transportation Statistics (BTS) through the Puerto Rico Institute of Statistics. The data are composed of monthly counts of passengers who arrived and left the island per month from January 2010 to February 2018. These data are unbiased, but they have a coarse temporal and geographic scale and cannot account for the same individuals on repeated trips or distinguish Puerto Rican residents from short-term aid workers and other visitors.

Details on how each of these datasets was constructed are included in Methods, and a side-by-side comparison is found in SI Appendix, Table S1.

We found agreement across all data sources of a consistent loss of population from 1 July 2017 to 1 July 2018; however, the dynamics and magnitude of the loss differ between data types (Fig. 1). The ACS predicts a population loss of 129,848, a 4% decrease; the APT data predict a population loss of 168,295, a 5% decrease; the Teralytics data predict a population loss of 235,375, an 8% decrease; and the Facebook data predict a 17% decrease, which equates to a total estimated population loss of 475,779 on the island. In all cases, we observed a sharp drop after the hurricane. Both the APT data and Teralytics data show a rebound in population and stabilization after 31 December 2017. The Facebook data do not show stabilization until April 2018. We note that Facebook data are based on a closed cohort that is not able to measure population increases due to immigration to Puerto Rico. However, it does represent the cohort that likely experienced Hurricane Maria since the selection criteria used by Facebook’s algorithm ensure that transient populations, such as tourists, are not included in the closed cohort.

Facebook provided municipality-level data that allowed us to assess within-island migration. For each municipality, the data included weekly counts of individuals still in the municipality and new to the municipality. For privacy protection, the new to municipality data were not included if the number of users was below a prespecified threshold (Methods has details on the data imputation techniques used, as well as a sensitivity analysis comparing different imputation approaches). We found that the share of users belonging to each municipality was highly correlated with the baseline population size of each municipality (SI Appendix, Fig. S1). On average, municipalities with a smaller population size lost a bigger proportion of their population during the study period (Fig. 2). The municipality of Toa Baja was a large outlier (SI Appendix, Fig. S2) due to a reported surge of new individuals moving there at the end of the study period. Possible explanations for Toa Baja being an outlier are provided in Discussion. Apart from this outlier, San Juan, which is the capital of Puerto Rico, is the only municipality with more individuals at the end of the study period relative to baseline.

All municipalities experienced loss of their baseline resident populations. Much of this loss is explained by off-island immigration (Fig. 3). However, the locality gaining the most new residents was San Juan (Fig. 3). In fact, for urban areas, the baseline resident population loss was compensated by in migration from other municipalities. For example, San Juan experienced a 19% loss of its baseline resident population but ended the study period with a 7% increase in total population, suggesting an in migration of 26% by the end of the study period. Another appealing aspect of these data is that they provide information on the destination of those displaced (Fig. 3B). Note that the top two destinations for Facebook users in our cohort were Miami, FL, and Jacksonville, FL, respectively. These results taken together suggest a migration pipeline from rural to urban municipalities and likely off island.

Disaster Maps provided information about the proportion of people whose location is unknown for a particular week. Unknown locations may be the result of 1) people stopping use of Facebook, 2) people having changed their Facebook behavior, or 3) loss of electricity and communication infrastructure. In the first full week of data collection immediately after the hurricane, approximately half of the cohort members did not register a location (Fig. 4A), likely due to the widespread loss of infrastructure. During this week, the proportion of users with reported unknown locations was substantially higher in the rural areas compared with in urban areas (Fig. 4B and SI Appendix, Fig. S3). We focused our analysis to only include the week immediately following the impact of Hurricane Maria to ensure that any loss of individuals in a region was not likely to be due to population migration. This is especially important as there was widespread damage to roads and loss of transportation infrastructure during this time (12).

In SI Appendix, Table S2, we show the top 15 municipalities that had a greater proportion of cohort members with their location unknown compared with the island-wide average. While rurality is directly linked with availability of infrastructure, the Facebook data highlight the heterogeneity in loss of access to resources in the immediate aftermath of a disaster. Some municipalities continued to have missing data relative to the island-wide average for many months.

We obtained yearly population estimates from 2010 to 2018 from the intercensal yearly population estimates provided by the ACS, a yearly survey conducted by the US Census Bureau. Specifically, we have the vintage 2018 estimates that were published last year.

Facebook is a social media company with over 2 billion monthly active users. In 2017, the Data for Good team within Facebook launched Disaster Maps with the goal of aiding response organizations with information vital to optimal resource allocation in postdisaster settings. These maps are built using privacy-preserving techniques including aggregation and deidentification that protect individual privacy. The Harvard team accessed these data free of charge through Data for Good’s standard license agreement, which allows partners working in humanitarian operations and research to improve their work through the use of Disaster Maps.

For this study, we used the original version of Facebook’s Displacement Maps, a specific product within Disaster Maps (27), which has since been updated to improve its data sources and methodology. However, work is needed to assess the functionality of these improvements in the context described here. Displacement Maps are generated by first defining a geographic bounding box along with an index date defining the disaster event of interest. In this case, the geographic bounding box consisted of the entire island of Puerto Rico, and the index date was 20 September 2017, the day that Hurricane Maria made landfall on the island. In these data, Puerto Rico is divided into 188 nonoverlapping geographical tags (geotags), and users are assigned a home location defined as the geotags where they had the most interactions with the Facebook platform through a browser during the 5 wk prior to the index date. In the original Displacement Maps methodology, the location of these interactions was determined from associated internet protocol addresses; the new methodology instead utilizes location-based data from cell phones. Displacement Maps generate a closed population consisting of people using Facebook satisfying the following two conditions: 1) they registered an interaction with Facebook services from in the geographic bounding box during the 5 wk preceding the index date, and 2) they were present in their home location during the week before the disaster (28).

They then followed this cohort through July 2018 and calculated the most commonly occurring geotag each week, aggregating total numbers of cohort members per geotag. Geotags with fewer than 100 people at baseline were excluded to protect privacy. Then, for each of the 49 wk after the index date, the dataset includes the number of crisis-affected people in their home location, the number of new users, and the number of unknown users for each geotag. Cohort members were defined as having an unknown location if they did not register an interaction during the week for which data were aggregated. For our results, we assumed these people with unknown locations were in their home location. We then combined geotag counts into Puerto Rico’s 78 municipalities. Due to low counts, the municipality of Las María had no data. Below, we describe how we imputed other missing data.

Teralytics is a tech company that works with governments and private clients to assess human movement by partnering with mobile network operators. Specifically, the partnership allows the company to access and analyze the data that cell towers receive from mobile devices. From them, we obtained island-wide daily population proportion estimates relative to an undisclosed baseline from 31 May 2017 to 30 April 2018, based on all subscribers of a major undisclosed telecom company that created events in Puerto Rico. Events were defined as signal exchanges between a cell phone and the nearest cell phone tower. These signal exchanges occurred, for example, when a phone call was made or a text message was sent. The data were filtered by the provider, and only subscribers with activity all over the Teralytics analysis period were considered (31 May 31 2017 to 30 April 2018). Activity was defined as an event on at least 10 d/mo for all 12 mo. Due to the unreliability of data, proportion estimates for 4 wk after Hurricane Maria were not provided. For every day, a distinct number of subscribers in Puerto Rico is computed by considering events generated from different mobile devices. It is important to note that Teralytics is a commercial company that operates by competitively sourcing, cleaning, and preprocessing these data. Therefore, much of this analytic pipeline is proprietary and a black box to researchers. We have taken measures in our analysis to more readily index and compare this source with others; however, as noted in SI Appendix, Table S1, all of these sources have their benefits and drawbacks.

Finally, through the Puerto Rico Institute of Statistics we obtained APT data from the US BTS. The data are composed of monthly counts of passengers who arrived and left the island per month from January 2010 to February 2018. The per-month difference between these two numbers will be referred to as net migration. We added the monthly net migrations to the vintage 2017 population estimates corresponding to the same date and interpolated using a linear model between these data points. This resulted in daily population estimates from July 2010 to July 2018 that account for flight passenger movement (29).

SI Appendix, Table S1 has a side-by-side comparison of all of the data sources. This research proposal was reviewed by the Harvard T. H. Chan School of Public Health Institutional Review Board and was deemed exempt as nonhuman subjects research.

We estimated island-wide population sizes using each of the four data sources. For the ACS data, we simply interpolated the points for each year (SI Appendix, Fig. S3). For the other three sources, we defined population size estimate $N_t$ for time t using the following:

where $N_0=\mathrm{3,325,001}$, the ACS population estimate for 1 July 2017; $m_t$ is a source-specific measurement for time $t$; and $m_0$ is a source-specific baseline. For the APT data, baseline was defined as 1 July 2017, $m_t$ corresponds to the sum of $N_0$ and the cumulative net passenger movement for month $t$, and $m_0$ is the sum of $N_0$ and the cumulative net passenger movement at baseline. The formula above, therefore, is simply the cumulative sum, starting on 1 July 2017, of passengers arriving and passengers leaving:

where $i$ represents the $i\text{th}$ month after baseline and $passengerou{t}_i$ and $passengeri{n}_i$ are the passengers leaving and entering Puerto Rico in month $i$, respectively. For the cell phone data, the baseline was also defined as 1 July 2017; $m_t$ and $m_0$ represent the proportions provided by Teralytics for time $t$ and the proportion corresponding to baseline, respectively. Finally, for the Facebook data baseline was defined as 21 August 2017, which corresponds to the first observation in the dataset. Here, $m_t$ corresponds to the Facebook population at time $t$, and $m_0$ represents the Facebook population at baseline.

For the municipality-level data, we aggregated the city-specific data from Displacement Maps by municipality in Puerto Rico while maintaining the temporal resolution at weeks. For all analyses, we assumed that cohort members who were defined as unknown for the week were present in their home town locations and were either unable to or chose not to interact with Facebook services using a browser in that period.

To evaluate the utility of Displacement Maps, a tool to target municipalities for resource allocation, we evaluated the proportion of the population with unknown locations every week compared with baseline. Our primary assumption here is that in the immediate aftermath of a disaster, any discontinuation of interaction with Facebook services defined at baseline would primarily be caused by loss of access to infrastructure, death, or other factors related to the event. Therefore, a higher proportion of municipality-specific cohort members whose location was unknown would be a proxy for impacts of the disaster in that region.

In Facebook’s dataset, Puerto Rico is divided into 188 nonoverlapping geotags. For each geotag, the dataset contains the number of crisis-affected people in their home location at time t, denoted here as home; the number of new users at time t, denoted here as new; and the number of users whose location is unknown at time t, denoted here as unknown. For privacy-preserving reasons, geotags were excluded from the data if they had fewer than 100 Facebook cohort members at time t. We, therefore, applied an imputation approach along with a sensitivity analysis to provide intervals of possible values. First, we computed the distribution of the maximum number of consecutive new missing values for each geotag. Thus, the $j\text{th}$ geotag had a corresponding value, $c_j^{\left(new\right)},$ that represents the maximum number of consecutive weeks when the new variable was missing. Second, we denote a threshold τ that corresponds to the maximum number of successive weeks we are willing to accept and a default imputation value δ to be used later. Then, we split the data into the geotags that comply with the threshold and those that did not [i.e., if $c_j^{\left(new\right)}\le$ τ, then geotag j complies with the threshold]. For the geotags that complied with the threshold, we imputed the missing values using a linear interpolation of the observed values. For the geotags that did not comply with the threshold, we imputed the missing values with δ. Finally, we computed the distribution of the maximum number of consecutive home missing values for each tag. Hence, similar to above, the $j\text{th}$ geotag has a corresponding value, $c_j^{\left(home\right)}$, that corresponds to the maximum number of consecutive weeks when the home variable was missing. Then, we excluded all geotags where the newly computed statistic was greater than τ (SI Appendix, Fig. S6). For our results, we used τ = 8 and δ = 99 (SI Appendix, Fig. S7). Note that δ = 99 is the maximum number that the missing values can take. We conducted a sensitivity analysis where δ = 0, the minimum number that the missing values can take (SI Appendix, Fig. S8). Our sensitivity analysis shows that the results do not change much.