Reconstructing subdistrict-level population denominators in Yemen after six years of armed conflict and forced displacement

2.2.1 Population estimates

Yemen conducted its last census in 2004 [10]. The UN World Population Prospects [11] provide country-wide yearly projections from this baseline, reflecting assumed natural growth. The WorldPop project redistributes these projections across space, with 100 m² resolution, by applying a geospatial model that predicts population density using a variety of remotely sensed climate, topography, illumination, transport network, urbanisation and other land use variables (see https://www.worldpop.org/methods, Stevens et al. [12] and Linard et al. [13]). WorldPop annual estimates were used as the baseline (June 2014) and to quantify subsequent yearly relative changes by subdistrict due to migration (Table 3). We speculated that WorldPop estimates might not accurately capture forced displacement, since many Yemeni IDPs live in rented accommodation or communal buildings [14] that would not appear changed in remotely sensed observations.

2.2.2 Internal displacement data

Available displacement datasets covered both IDPs and ‘returnees’ (IDPs who have returned to their communities of origin). Returnee data are acknowledged to feature underestimation [15]. Data were also classifiable as ‘prevalent’ (i.e. information on IDPs or returnees present at a specific time in a given location) or ‘incident’ (novel displacements or returns). Prevalent data sources consisted of baseline or repeat site assessments carried out by the UNHCR-Population Movement Tracking (PMT; 2015-2016) and the International Organisation for Migration (IOM)’s Displacement Tracking Matrix (DTM; 2016-2018) projects. Assessments were sometimes done in-person by agency staff, but largely relied on a network of key informants working with standard templates [16]. The November 2018 assessment round achieved the highest geographical coverage (Supplementary file). During 2018, the IOM also collected incident data on displacements due to insecurity in Al Hudeydah governorate. Since 2019, only incident data have been published.

Displacement datasets were available as unprotected Microsoft Excel worksheets on the IOM DTM site (https://displacement.iom.int/yemen); variable sets, names and formats changed repeatedly over time. We retained the following variables, as available: date of displacement or return; date of assessment (prevalent data only); governorate, district, subdistrict, locality name, geographic codes (hereafter, geocodes) of these administrative levels, and coordinates of the location of arrival / refuge; governorate, district and subdistrict of origin, or of last displacement for returnees, with their geocodes; and number of households. Locality geocodes had an 11-digit structure: digits 1-2 identify the governorate, 1-4 the district and 1-6 the subdistrict. We validated recorded geocodes against the OCHA and CSO gazetteers [6] of all place names (for localities); we also used available locality geocodes to work out missing governorate, district or subdistrict geocodes.

We appended all displacement datasets into one. After applying range and consistency checks, and deleting duplicate records (these were only identifiable for prevalent data, based on identical dates of assessment, displacement/return and locality geocode), the appended dataset consisted of 222,069 records, of which 206,109 (92.8%) concerned IDPs and 15.960 (7.2%) returnees; 195,477 (88.0%) were prevalent-type data. Only IDP data were carried into further analysis, as we assumed returnee data were too incomplete. After removing 2784 (1.3%) records with missing year or month of displacement and 13,188 (6.4%) with missing district of origin or arrival, we retained 190,137 IDP records.

2.2.3 Predictors of displacement

We used multivariate predictive models to impute missing subdistrict data and quantify IDP movements after displacement (see below). While some predictors were built from the population and displacement data themselves, we searched for additional candidate predictor datasets available at month-year and subdistrict resolution. These included (i) a CSO geospatial dataset of Yemen’s road network [17], which we transformed into road density (Km per Km² area); (ii) a crowd-sourced dataset of health facilities [18], which we combined with WorldPop data to estimate health facility density per 100,000 inhabitants; and (iii) the Armed Conflict Location and Event Data Project (ACLED) as a source of georeferenced insecurity event information [19]. Since 2015, ACLED has carried out particularly intensive data collection on Yemen through media monitoring and networks of in-country civil society sources [20]. We mapped each insecurity event to subdistricts based on the event’s coordinates (87/62,629 or 0.1% of records did not map to a subdistrict and were excluded), and aggregated data by subdistrict-month.

2.3.1 Standardising place names

Place names in the displacement dataset were a mixture of Arabic and inconsistent Latin character transliterations, and only a fraction of data had unique geocodes, with some geocodes mapping to places that differed from those recorded. We therefore came up with equivalences, down to subdistrict level, between the displacement dataset and the OCHA gazetteer. The dataset featured 5833 unique instances (‘sets’) of missing or non-missing governorate, district, subdistrict names and geocodes. We identified OCHA gazetteer matches for each such set based on the following sequentially applied criteria: (1) the recorded geocode sets (e.g. ‘142024’ for a set down to subdistrict level; ‘1115’ for a set down to district level) also existed in the OCHA gazetteer, and the OCHA place names they mapped to (e.g. governorate 14, Al Bayda, district 1420, Al Malajim and subdistrict 142024, Dhi Khirah) were the same within two characters as the names recorded on the dataset; (2) geocodes were missing, but the place name sets matched a place name set in the OCHA gazetteer within two characters, after applying eight approximate character string matching techniques (stringdist package [21]) and choosing the most common match; (3) the recorded geocode sets also existed in the OCHA gazetteer, and at least the recorded governorate and district names matched with the OCHA names corresponding to the same geocodes (a less stringent version of criterion 1); and (4) a combination of approximate string matching and manual searches applied to any remaining unmatched sets, with matches established at governorate, then district, then subdistrict level so as to restrict each successive search to place names within the same higher-level administrative unit. Criteria 1, 2, 3 and 4 led to a match for 69.2% (4038/5833), 17.3% (1012/5833), 6.5% (380/5833) and 6.9% (403/5833) of unique instances, respectively. All instances were successfully matched.

2.3.2 Imputing missing subdistrict data

After applying the above equivalence, the subdistrict was still missing for 25.6% (48,642/190,137) locations of arrival (31.5%, 30.0%, 1.2%, 49.4%, 0.0% and 100.0% for data collection years 2016 to 2021, respectively) and 79.5% (151,144/190,137) locations of origin (97.8%, 98.7%, 19.1%, 100.0%, 100.0% and 100.0%). Missingness was higher for incident (45.2% and 100.0% for locations of arrival and origin, respectively) than prevalent (22.6% and 76.3%) data. We identified most missing subdistricts through four sequential steps, applied to the individual records (Figure 1):

• Download figure
• Open in new tab

Figure 1.

Reductions in subdistrict missingness achieved after each successive data management step. Only eligible IDP records are included in the denominator.

1. Many (n = 117) districts consisted of only one subdistrict;

2. For records that featured longitude and latitude, we identified the subdistrict based on the OCHA/CSO administrative boundaries that the coordinates fell within;

3. We used the above string matching techniques to match the locality name (below subdistrict level), if recorded, to the locality name of other records within the displacement dataset, restricting the matching search to the same governorate and district; if any of the matching localities had a non-missing subdistrict, the latter was applied to matching records with missing subdistrict;

4. We also matched the recorded locality name to locality names at administrative levels below subdistrict (city, neighbourhood, harrah, village, sub-village) in the CSO gazetteer, again restricting the matching search to within the same governorate and district, adopting the closest match across the above administrative levels and looking up the subdistrict the CSO match fell within.

We resolved additional missingness through machine learning models. We first predicted how many subdistricts IDPs came from or went to, out of all possible subdistricts in each ‘parent’ district of origin/arrival, by month-year of displacement (Model 1); we then predicted which specific subdistricts IDPs came from or to (Model 2), and the relative share of all IDPs that came from/went to each such subdistrict, out of the parent district total (Model 3). All models were trained on DTM data from November 2018 (n = 41,375), which had 89.8% subdistrict of origin and 100.0% subdistrict of arrival completeness. Training data were aggregated by month-year of displacement, subdistrict of arrival and subdistrict of origin, and augmented to feature all other subdistricts of origin/arrival within the same district, with outcome = 1 attributed to the subdistricts of origin/arrival that any IDPs did come from, and 0 otherwise. For Model 1, training data were further aggregated by district of arrival or origin. For Model 3, training data excluded subdistricts that IDPs did not come from/move to as well as districts with a single subdistrict of origin/arrival.

For each model, random forest algorithms were grown using the ranger R package [22] (weighted for class imbalance and tuned to 500 trees, up to 3 variables to split each node on and maximum tree depth of 20), using the following candidate predictors (at district level for model 1; at subdistrict level otherwise): total population, distance between the geodesic centroids of the subdistricts of origin and arrival, cumulative incidence per capita of insecurity events and fatalities during the current and previous month, health facilities per capita, road density per surface area, surface area, number of candidate subdistricts of origin/arrival within the parent district and the natural log of the number of IDP households from/to the parent district. We evaluated models’ performance out-of-sample using ten-fold cross-validation.

For subdistricts of origin, model 1 yielded fair predictions for the single subdistrict category, but was downward-biased for multi-subdistrict observations (Table 1). As nearly all (97.3%) instances in the training data had only one subdistrict of origin, we applied a simplifying assumption that, for any district of origin - time - subdistrict of arrival combination, all IDPs came from a single subdistrict. For subdistricts of arrival, model 1’s performance was reasonable (Table 1), and we applied the corresponding predictions. Model 2, constrained to the number of subdistricts predicted by model 1, correctly guessed ∼80% of the true subdistricts of origin and ∼63% of true subdistricts of arrival, while incorrectly classifying ∼4% of subdistricts that IDPs did not come from and ∼3% of subdistricts that IDPs did not move to. For model 3, we categorised the outcome (percent share of IDPs) into five classes from 1-19% to 80-99%, with the predicted value at the mid-point of each class. As shown in Table 2, the model had low predictive performance for both origin and arrival subdistricts, though it performed better than a random guess. This low performance was inconsequential for subdistricts of origin (since we assumed no multi-subdistrict instances) and of minor influence for subdistricts of arrival, since model 1 predicted multi-subdistrict instances for only 12.2% of the dataset. Further, we scaled model 3’s predictions to ensure a denominator of 100% for each district.

Overall, we imputed missing subdistricts for all but a small minority of displacement records (Figure 1), which were excluded from further analysis.

2.4.1 General equations

Let N be a two-dimensional matrix with dimensions i ∈ (1, 2, 3 … / substricts of arrival) and t ∈ (1, 2, 3 … T, with 1 = June 2014 and T = September 2021, and each increment = one month) where n_{i t}is the population of subdistrict i at the start of month t. Let F be a three-dimensional matrix of forced displacement, with dimensions t, i and j ∈ (1, 2, 3 … ] subdistricts of origin), where f_{i jt}is the net flow of IDPs from j to i during month t (IDPs can be displaced within their subdistrict, in which case i = j). We further define B, D and M as matrices with dimensions i and t where b_itand d_itare birth and death rates per capita during each subdistrict-month, and m_itis the percent change in subdistrict i’s population during month t resulting from migration other than forced displacement. Lastly, φ is the proportion of forced displacement that is already taken into account by M, i.e. the extent to which data on migration also capture forced displacement (if φ= 0, the two data sources have no overlap; cp = 1 implies complete overlap). It follows that Otherwise put, the following month’s population is this month’s starting population multiplied by the net rate of natural growth and migration, plus any net change in IDPs that isn’t already captured by migration estimates.

2.4.2 Estimating displacement flows

To populate matrix F, we needed to combine prevalent and incident data. We aggregated all data to identify unique IDP groups that moved from a given subdistrict of origin to a given subdistrict of refuge during a given month (we call these ‘instances’, denoting discrete waves of primary displacement). Each such instance was subject to one or more longitudinal observations (DTM site assessments). For a minority (30.2% or 6258/20,744) of multi-observation instances, later observations in the dataset featured a higher number of IDPs than at previous assessment points, which is theoretically impossible beyond marginal increases due to natural growth. Most such cases were moderate and occurred when small numbers of IDP households were involved. We assumed these were due to clerical error, data collection problems or the discovery of previously undetected IDP households. We converted these problematic instances into constant or monotonically decreasing series based on alternative assumptions (Table 3).

As depicted in Figure 2, prevalent observations may underestimate past displacement, since during the period before assessments all or some of the IDPs may have returned home or moved to another subdistrict. This bias is dependent on the rate of return or onward movement. Equally, how IDP populations evolve after the last assessment is unknown.

• Download figure
• Open in new tab

Figure 2.

Illustration of four hypothetical scenarios (A to D) in which a group of IDPs arrives to a subdistrict from another subdistrict. In scenario A, nearly all IDPs remain in the subdistrict of refuge throughout the period of interest. In scenario B, only a fraction are left by the second assessment round. In scenario C, all IDPs have left the subdistrict (either returned to their subdistrict of origin, or moved elsewhere) by the second assessment, and in scenario D IDPs have left even before the first assessment, thereby potentially being missed altogether by the displacement tracking system.

To quantify this evolution and thereby predict IDP populations before and after the timeframe of available observations, we used the gamlss framework [23] to fit a generalised additive mixed growth model to all unique instances, as defined above, for which at least two prevalent assessment observations existed, excluding records for which the subdistrict was imputed. The model predicted the number of IDP households as a monotonic penalised-spline smoothed function of time since displacement, with IDP instance as a random effect, and, as predictors, distance between subdistricts of origin and arrival, health facility density (arrival), road coverage (origin and arrival) and a monotonic spline of insecurity event incidence in the subdistrict of origin since the previous assessment. We assumed a quasi-Poisson distribution for the data, as this provided reasonable model diagnostics. Figure 3 suggests that, on average, about 75% of IDPs left the subdistrict of arrival within two years of first displacement; thereafter, departures reduced.

• Download figure
• Open in new tab

Figure 3.

Predictions of a generalised additive mixed growth model of IDP population as a function of time. Panel A shows the model’s predicted percent change in IDP group size by month since displacement: the thick line is the main analysis prediction, and the shaded area indicates the 95%CI. Predictions using reasonable-high and -low adjustments for non-monotonic series are also shown. Panel B shows the model’s predictions (blue dots and 95% confidence bands) and observed total prevalent number of IDPs assessed (orange squares) during any given month after displacement, for the main analysis only.

We used the model to predict the evolution of IDP household counts across time for both the model-training data instances and all other (i.e. single-assessment and incident) instances: for the latter, we made predictions using only the fixed-effects model coefficients, and scaled these to the single recorded IDP household count. We made a simplifying assumption that IDPs either stayed in the subdistrict of first arrival or returned to their subdistrict of origin (i.e. zero secondary displacement).

2.4.3 Input values

Input values for all equation parameters are detailed in Table 3 and Figure 4.

• Download figure
• Open in new tab

Figure 4.

Countrywide average values of crude birth rate and crude death rate assumed, by scenario. The grey dotted line shows secular trends based on UN projections; the latter are represented by squares and centred at the mid-point of their period of reference. The shaded area indicates the crisis period.