We have found some issues when running the preprocessing steps of the code including:

• sections that are commented out which are needed only once as a set up
• difficulty in understanding how the code is meant to be run - no defined pipeline

03_HH_Model_Workflow_2024

Some Key Statistical terms

Gamma distribution: The gamma distribution is a continuous probability distribution commonly used to model waiting times until a certain number of events occur. It is well suited for data that are positively skewed. The distribution is described by two parameters: a shape parameter (α) and a rate parameter (β).

Fixed effects: Fixed effects are variables in a model that estimate the relationship between specific predictors and the outcome, while accounting for the presence of other predictors in the model.

Random effects: Random effects account for differences between groups or individuals that are not explained by the fixed effects. They represent unexplained variation across groups or individuals.

z-scores: A z-score measures how far a value is from the mean in terms of standard deviations. A positive z-score indicates the value is above the mean, while a negative z-score indicates it is below the mean.

Matérn covariance function: The Matérn covariance function is widely used in spatial statistics to model spatial correlation. It reflects the idea that observations closer together in space tend to be more similar than those farther apart.

Packages & Files

This script loads several packages used for spatial analysis, statistical modelling, data wrangling, and data visualization. It uses R‑INLA, which performs fast approximate Bayesian inference for latent Gaussian models. The inlabru package is also used, providing a higher‑level and more user‑friendly interface for building spatial latent Gaussian models based on the Integrated Nested Laplace Approximation method implemented in R‑INLA. In addition, the script uses the sf package to manage and analyse spatial vector data; gstat for geostatistical tasks such as variogram estimation, spatial interpolation, and spatio‑temporal modelling; and spdep for analysing spatial dependence and spatial autocorrelation using neighbour structures and spatial weights.

After loading the required packages, the script imports the datasets, including a CSV file with household and building information and a shapefile defining the spatial boundaries of enumeration areas. A unique numeric identifier is then created for each district, followed by the creation of a categorical rural–urban classification variable. The script then produces a summary of household counts.

In the next step, a new variable—household density—is calculated as the ratio of the observed number of households to the number of detected buildings (google_v2_5). Records with missing (NA) or infinite values of household density are removed, and summary statistics for this variable are produced. Finally, the script generates a set of exploratory plots to examine the distributions of household density, household counts, and building counts across enumeration areas. Boxplots, histograms, and density plots are used to assess central tendencies, variability, and the presence of outliers.

Preparation of Data for Modelling

The next section of the script prepares the covariates for modelling. First, all predictor variables with names beginning with “x” are selected from the dataset. Summary statistics are calculated for each covariate, after which a standardization function is defined and applied to convert all covariates to z‑scores. This ensures that all predictors are on a comparable scale. The standardized covariates are then combined with the response variable, hh_density, to create the modelling dataset (see lines 211–224). Any columns containing only missing values are removed to keep the dataset clean

The script then performs covariate selection to identify the most relevant predictors of household density. First, a Gaussian generalized linear model (GLM) is fitted using all available covariates, and MASS::stepAIC() is applied to select a more parsimonious model (see lines 233–243). The Akaike Information Criterion (AIC) balances model fit and complexity, with lower values indicating a better model. AIC selects the model that best explains the data while using as few parameters as possible, with lower values indicating a better trade‑off between model fit and complexity.

Next, variance inflation factors (VIFs) are examined to assess multicollinearity among predictors. An iterative procedure removes predictors with high VIF values to reduce collinearity (see lines 265–268). The script then applies an additional backward elimination step based on statistical significance. Using a manually specified set of predictors, Gaussian GLMs are repeatedly fitted, and the predictor with the highest p‑value is removed if it exceeds 0.05. This process continues until all remaining predictors are statistically significant (see lines 287–320). Finally, the summary of the final model and its VIF values are reviewed.

Fixed effects selected:

Covariate labels​ | Descriptive names | Covariate names​

• x13 | Google (BCB_gl): Building metrics are calculated based on the pixel in which their respective centroids are located (Building Centroid Based – BCB method) | mean.MOS_MLW_buildings_mean_area_BCB_gl_100m_v1_1
• x44​ | Distance (in metres) from each 100m pixel to the nearest edge of flooded freshwater forest (ESA CCI class 160), derived from the 2022 land cover map, for Malawi — WorldPop covariate version 1 | mean.MOS_MLW_esalc_160_dst_2022_100m_v1
• x47 | Same as above for Water bodies (210) | mean.MOS_MLW_esalc_210_dst_2022_100m_v1
• x45 | Same as above but for Urban / built-up (190) | mean.MOS_MLW_esalc_190_dst_2022_100m_v1
• x39 | Terrain elevation (metres above sea level) at 100m resolution, sourced from the MERIT DEM, tile/product code 103 | mean.MOS_MLW_elevation_merit103_100m_v1
• x49 | Distance (in metres) from each pixel to the nearest major road/highway, derived from OpenStreetMap data as of 2023 | mean.MOS_MLW_highway_dist_osm_2023_100m_v1
• x56 | Distance (in metres) from each pixel to the nearest road intersection, derived from OpenStreetMap 2023, at 100 m resolution | mean.MOS_MLW_rd_intrs_dist_osm_2023_100m_v1
• x61 | "Distance (in metres) from each pixel to the nearest water body (lake, river, pond, etc.), derived from OpenStreetMap 2023, at 100m resolution | mean.MOS_MLW_waterbodies_dist_osm_2023_100m_v1

Model 1 - Fixed Effect + Urban_Rural_Random_Effect (see lines 387-454)

After selecting the influential variables as fixed effects, Model 1 fits a Bayesian gamma regression with household density as the outcome variable. The model includes the selected fixed effects and a random intercept for rural–urban category, allowing rural and urban areas to have their own unobserved deviations from the overall regression line. This random effect captures differences between rural and urban areas that are not explained by the geospatial predictors. In this way, the model estimates household density based on the observed geospatial variables while also accounting for variation between rural and urban areas. The output provides posterior summaries for both the fixed effects and the rural–urban random effects.

In‑sample predictions are then produced for each observation in EA_pop using the estimated fixed effects and the corresponding rural–urban random effect. Since these predictions are produced on the linear predictor scale, they are back‑transformed using the exponential function to obtain predicted household density on the original response scale. Finally, model performance is evaluated by comparing observed and predicted household densities and household counts using summary accuracy measures.

Model 2 - Fixed Effect + Urban_Rural_Random_Effect + Dist_Random_Effect (see lines 459-530)

Model 2 extends Model 1 by adding an additional independent (iid) random intercept at the district level, while keeping the same fixed effects and the rural–urban random effect. This district‑level random effect allows household density to vary across districts due to unobserved district‑specific differences that are not explained by the geospatial predictors or by the rural–urban grouping. That is to say, the model captures variation at two levels: overall differences between rural and urban areas, and additional variation between districts. As in Model 1, posterior summaries of the model parameters are examined, in‑sample predictions are produced and back‑transformed to the response scale, and model performance is evaluated by comparing observed and predicted household densities and household counts.

Model 3 - Fixed Effect + Urban_Rural_Random_Effect + Dist_Random_Effect + EA Random_Effect (see lines 532-603)

Model 3 further extends the hierarchical structure by adding an independent (iid) random intercept at the enumeration area (EA) level, while retaining the same fixed effects and the rural–urban and district random effects. This additional random effect captures unobserved, EA‑specific variation in household density that is not explained by the selected predictors or by higher‑level grouping structures. That is to say, the model accounts for variation at three levels: rural–urban category, district, and enumeration area. As in the previous models, posterior summaries are examined, in‑sample predictions are produced and back‑transformed to the response scale, and model performance is assessed by comparing observed and predicted household densities and household counts.

Geostatistical Model (see lines 610-673)

The next code block fits a geostatistical model, extending the earlier hierarchical models by adding a spatial random effect using the INLA–SPDE framework. This allows the model to capture spatial dependence, reflecting the idea that enumeration areas located closer together are more likely to have similar household densities than those farther apart.To construct the spatial component, centroid coordinates of the enumeration areas are used to define the spatial domain. A non‑convex hull and a triangular mesh are then created over the study area. This mesh provides the basis for estimating a Matérn spatial random field, which models residual spatial variation as a smooth surface across space. This way, the model accounts not only for variation between rural–urban categories, districts, and enumeration areas, but also for spatial autocorrelation that is not explained by the observed covariates. After fitting the model, posterior summaries are examined, in‑sample predictions are produced and back‑transformed to the response scale, and model performance is evaluated using the same household‑density and household‑count accuracy measures as in the previous models. Finally, model comparison statistics are compiled to compare all four model specifications.

Cross validation (see lines 728-1663)

This section evaluates the predictive performance of Models 3 and 4 using two cross-validation approaches. In both cases, the data are divided into ten folds, and the models are repeatedly fitted and evaluated across different training–test splits. The aim is to assess how well the models generalise to new data, rather than relying solely on in‑sample fit. Model performance is examined for both household density and the derived household count.

The first approach applies standard 10‑fold cross‑validation to both Model 3 and Model 4. For each fold, the model is fitted using the training data and then used to generate predictions for both the training and test sets. Predicted household density values are back‑transformed to the original response scale and used to estimate household counts from building counts. For each fold, several accuracy metrics are computed, including root mean squared error (RMSE), correlation, mean absolute error (MAE), and bias. Both the performance statistics and observation‑level predictions are stored for later summarisation.

The second approach uses INLA’s grouped cross‑validation procedure. In this case, the training and test observations are combined into a single dataset, but the outcome values for the test observations are set to missing before model fitting. This allows INLA to generate predictions for the held‑out observations directly within the modelling framework. This procedure is applied to both Model 3 and Model 4, with Model 4 additionally including a spatial random effect through the SPDE specification. As in the standard cross‑validation approach, predictions are back‑transformed to the response scale, converted to household counts, and evaluated using the same set of performance metrics.

Finally, the code summarises the average training and test performance across all four cross‑validation analyses: standard cross‑validation for Model 3, INLA grouped cross‑validation for Model 3, standard cross‑validation for Model 4, and INLA grouped cross‑validation for Model 4.

Make Predictions for HH Count Models (see lines 1670 - 2218)

This section generates household count predictions for all four models using a common prediction dataset and a consistent workflow. For each model, household density is predicted for every grid cell and then converted into household counts using building‑count information. The results are summarised at the national, district, and pixel levels, with uncertainty intervals reported throughout to reflect variability in the predicted values. The key difference between the four models lies in the type of geographic variation they capture, ranging from broad rural–urban and district‑level effects to more localised and spatially structured variation.

The final section (see lines 2238–2280) exports the predictions from the selected best‑performing model for rasterisation. The pixel‑level predictions are first converted to a spatial format, after which the main prediction summaries are rasterised to match the reference grid of the study area. These raster outputs are intended for mapping and further spatial analysis.