Quantifying the dynamics of migration after Hurricane Maria in Puerto Rico

ACS.

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.

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.

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.

Aviation Records.

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.

Estimating Population Size.

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 Nt for time t using the following:Nt=N0×mtm0,

where N0=3,325,001, the ACS population estimate for 1 July 2017; mt is a source-specific measurement for time t; and m0 is a source-specific baseline. For the APT data, baseline was defined as 1 July 2017, mt corresponds to the sum of N0 and the cumulative net passenger movement for month t, and m0 is the sum of N0 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:m0=N0+(passengersout0−passengersin0)mt=N0+∑i=0t(passengersouti−passengersini),

where i represents the ith month after baseline and passengerouti and passengerini 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; mt and m0 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, mt corresponds to the Facebook population at time t, and m0 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.

Imputation.

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 jth geotag had a corresponding value, cj(new), 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 cj(new)≤ τ, 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 jth geotag has a corresponding value, cj(home), 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.