The Standardization Problem

Consider comparing cereal overconsumption between rural and urban India. Rural households are disproportionately located in rice-dominant states and have different demographic profiles (more children, different social group composition). A naive rural–urban comparison would mix three effects: the actual sector effect on excess consumption, geographic composition differences, and demographic composition differences.

The prediction grid solves this by constructing a counterfactual population where all non-focal variables are held at the same national distribution.

Grid Construction Logic

The grid is constructed differently depending on the type of comparison. In all cases, the approach follows the same principle: hold the focal variable at its actual values, and standardize everything else to a common national distribution.

Grid Construction by Comparison Type:

STATE comparison [c("nss", "state_code")]
   → Keep: actual geographic structure (regions, rural/urban) within each state
   → Standardize: demographics (social, religion, child, female_headed)
   → Result: isolates state-level excess consumption differences net of demographics

SECTOR comparison [c("nss", "sector")]
   → Standardize: state-region distribution across rural/urban
   → Standardize: demographics nationally
   → Result: isolates PURE rural-urban effect

NSS REGION comparison [c("nss", "nss_region")]
   → Standardize: sector distribution across regions
   → Standardize: demographics nationally
   → Result: isolates regional effects net of urban/rural mix and demographics

DEMOGRAPHIC comparison [c("nss", <focal_var>)]
   → Standardize: full geography (state, region, sector)
   → Standardize: other demographics nationally
   → Result: isolates effect of focal demographic variable

Geographic standardization. When comparing sectors or regions, the grid standardizes the geographic distribution so that, for example, rural and urban get the same state–region composition. This isolates the “pure” sector effect from geographic confounding. Without standardization, sector differences would partly reflect the fact that rural populations are concentrated in poorer states.

Grid Construction Steps (Detailed)

For each comparison type, the grid is built through the following process:

Step 1: Define the core structure. The focal variable × expenditure bin combinations define the “backbone” of the grid. For each cell, the mean MPCE within the bin provides the income value at which predictions are evaluated.

Step 2: Compute geographic weights. For state comparisons, each state keeps its actual nss_region × sector structure. For demographic comparisons, a national geographic distribution is computed:

\[w_{\text{geo}}(s, r, u) = \frac{\sum_{i \in (s,r,u)} w_i}{\sum_i w_i}\]

where \(s\) is state, \(r\) is NSS region, and \(u\) is sector.

Step 3: Compute demographic weights. The national demographic distribution is computed marginally:

\[w_{\text{demo}}(\mathbf{d}) = \frac{\sum_{i: \mathbf{d}_i = \mathbf{d}} w_i}{\sum_i w_i}\]

where \(\mathbf{d}\) is the vector of non-focal demographic variables (social group, religion, child presence, female headship).

Step 4: Cross and combine. The geographic and demographic grids are crossed (Cartesian product) and their weights multiplied:

\[w_{\text{grid}} = w_{\text{geo}} \times w_{\text{demo}}\]

This produces the full prediction grid: every combination of the focal variable, expenditure bin, geographic cell, and demographic cell, with a weight reflecting its representation in the standardized population.

Step 5: Assign grid IDs. Each unique combination of the grouping variable × expenditure bin receives a grid_id, used to aggregate predictions.

Step 6: Build “Overall” grid. A separate grid is constructed using the overall (population-weighted) mean MPCE for each group, rather than decile-specific means. This provides the marginal group-level summary.

Expenditure Binning

Households are assigned to expenditure deciles using group-specific cutpoints from the MPCE distribution:

The cutpoints (\(q_{10}, q_{20}, \ldots, q_{90}\)) are computed within each group (e.g., separately for rural and urban) to ensure that each bin contains approximately 10% of the group’s weighted population. This means the same income level (say ₹2,000/month) may fall in different deciles for different groups.

Overview of the Prediction Pipeline

Step 1: Load pre-fitted models (logit + quantity)
Step 2: Construct standardized prediction grids (from previous section)
Step 3: Draw coefficient vectors from posterior (MVN)
Step 4: Draw dispersion parameters (chi-squared)
Step 5: Compute household-level predictions (chunked for memory)
Step 6: Aggregate to group × decile level (weighted)
Step 7: Compute survey standard errors (svydesign)
Step 8: Inject sampling uncertainty into model draws
Step 9: Summarize (mean, median, 95% credible interval)

Step 1: Covariance Matrix Validation

Before drawing coefficients, the posterior covariance matrix \(\mathbf{V}_p\) must be positive definite. Numerical issues in large GAM fits can occasionally produce matrices with near-zero or negative eigenvalues. Two repair strategies are available:

The validation checks all eigenvalues: \(\lambda_{\min}(\mathbf{V}_p) > 10^{-10}\).

Step 2: Draw Coefficient Vectors

Let \(\hat{\boldsymbol{\beta}}\) denote the estimated coefficients and \(\mathbf{V}_p\) the Bayesian posterior covariance matrix¹. For each simulation \(m = 1, \ldots, M\):

\[\boldsymbol{\beta}^{(m)} \sim \mathcal{N}\!\left(\hat{\boldsymbol{\beta}},\; \mathbf{V}_p\right)\]

drawn via MASS::mvrnorm(). Separate draws are made for the logit model (\(\boldsymbol{\alpha}^{(m)}\)) and the quantity model (\(\boldsymbol{\beta}^{(m)}\)).

Step 3: Draw Dispersion Parameters

For the Log-Normal model, the residual variance is drawn from its posterior:

\[\sigma^{(m)} = \sqrt{\frac{\hat{\sigma}^2 \cdot \nu}{\chi^2_\nu}}\]

where \(\nu\) is the residual degrees of freedom and \(\chi^2_\nu\) is a chi-squared draw with \(\nu\) degrees of freedom. A bias-correction term is precomputed: \(\delta^{(m)} = \tfrac{1}{2}(\sigma^{(m)})^2\).

For the Gamma model, the dispersion parameter \(\phi\) (inverse of shape) is drawn analogously:

\[\phi^{(m)} = \frac{\hat{\phi} \cdot \nu}{\chi^2_\nu}\]

yielding shape \(= 1 / \phi^{(m)}\) for each draw.

Step 4: Compute Household-Level Predictions (Chunked)

Predictions are computed in memory-efficient blocks (default: 2,000 rows per block) to handle the large prediction grids. For each block of rows in the grid:

4a. Linear predictors via the lpmatrix:

\[\boldsymbol{\eta}_L = \mathbf{X}_L \cdot \mathbf{B}_L^\top \qquad \boldsymbol{\eta}_Q = \mathbf{X}_Q \cdot \mathbf{B}_Q^\top\]

where \(\mathbf{X}_L\) and \(\mathbf{X}_Q\) are the design matrices from predict(..., type = "lpmatrix"), and \(\mathbf{B}_L\), \(\mathbf{B}_Q\) are the \(M \times p\) matrices of coefficient draws. The matrix cross-product tcrossprod() yields an \(n_{\text{block}} \times M\) matrix of predictions in a single operation.

4b. Transform to response scale:

Probability: \(\hat{p}_i^{(m)} = \text{logit}^{-1}(\eta_{L,i}^{(m)})\)

Conditional mean (log-normal): \(\hat{\mu}_i^{(m)} = \exp(\eta_{Q,i}^{(m)} + \delta^{(m)})\)

Conditional mean (Gamma): \(\hat{\mu}_i^{(m)} = \exp(\eta_{Q,i}^{(m)})\)

4c. Unconditional excess quantity:

\[\widehat{EY}_i^{(m)} = \hat{p}_i^{(m)} \times \hat{\mu}_i^{(m)}\]

4d. Posterior predictive draws (full distributional draws, not just means):

The unconditional posterior predictive combines a Bernoulli draw with the conditional draw:

\[Z_i^{(m)} \sim \text{Bernoulli}(\hat{p}_i^{(m)}), \qquad Y_{\text{uncond},i}^{(m)} = Z_i^{(m)} \times Y_{\text{pos},i}^{(m)}\]

where \(Y_{\text{pos},i}^{(m)}\) is drawn from the appropriate conditional distribution (log-normal or Gamma) using the draw-specific parameters.

Step 5: Weighted Aggregation to Group × Decile Level

Within each block, predictions are accumulated into group-level weighted sums. For each group \(g\) (defined by the grouping variable × expenditure bin):

\[\bar{p}_{g}^{(m)} = \sum_{i \in g} \tilde{w}_i \cdot \hat{p}_i^{(m)}\]

where \(\tilde{w}_i = w_i / \sum_{j \in g} w_j\) are the normalized within-group weights. The same weighted averaging is applied to \(\hat{\mu}_i^{(m)}\), \(\widehat{EY}_i^{(m)}\), and \(Y_{\text{uncond},i}^{(m)}\).

The result is a set of \(M\) draws for each group × bin combination, stored as list-columns:

Step 6: Survey Standard Errors

To quantify sampling uncertainty, plug-in predictions (\(\hat{p}_i\), \(\hat{\mu}_i\), \(\widehat{EY}_i\) at coefficient point estimates) are computed for every household in the microdata. These are then summarized using the survey package with the full complex survey design (PSU clustering, stratification, survey weights):

\[\text{SE}(\bar{p}_g) = \sqrt{\text{Var}_{\text{svy}}\!\left(\frac{\sum_{i \in g} w_i \hat{p}_i}{\sum_{i \in g} w_i}\right)}\]

computed via survey::svyby() with survey::svymean().

Lonely PSU adjustment. The option survey.lonely.psu = "adjust" is set to handle strata with a single PSU, ensuring that variance estimation does not fail.

Step 7: Injecting Sampling Uncertainty

The survey standard errors are injected into the model draws using scale-appropriate transformations. This combines both sources of uncertainty into a single set of draws.

Step 8: Final Summary Statistics

The \(M\) survey-adjusted draws for each group × bin are summarized into:

Engel Curve Plots

The primary visualization is a pair of Engel curves — the relationship between household expenditure and excess food consumption — plotted for each food item across demographic groups:

Top panel: Proportion consuming in excess (%)
   → Y-axis: Estimated share of people consuming above the daily requirement
   → X-axis: Monthly per capita expenditure (₹, 2011-12 prices, log scale)
   → Faceted by group level within grouping variable

Bottom panel: Average excess burden (grams)
   → Y-axis: Estimated daily excess consumption per AFE
   → X-axis: Same as above
   → Annotation: % of daily requirement at Overall mark

Key visual elements:

• Points and lines: Posterior means at each decile bin, connected by smoothed curves
• Ribbons: 95% credible intervals (shaded bands)
• Color: Survey round (2011–12 vs. 2023–24) — enables visual comparison of changes over time
• Annotation: Mean values (overall average) marked with dashed crosshairs and labels showing the exact excess quantity and its percentage relative to the daily requirement

Cross-Sectional Plots (Single Survey Round)

For state-level comparisons within HCES 2023–24, an alternative visualization displays:

• All states within a geographic region (Northern, Central, Eastern, Northeast, Western, Southern) on a single panel
• Smoothed Engel curves with labeled mean values
• Choice of excess quantity or proportion in excess on the y-axis

These are combined into multi-panel grids (3 × 2 for six regions) using the magick package for image composition, with an annotated title bar.

Suite of Analyses

The batch processing generates figures across all combinations of:

For state-level analyses, figures are further organized by geographic region.

Parallel Model Estimation

Each of the 2 food items requires two models (logit + quantity) for each of two survey rounds, yielding 8 model fits. These are estimated in parallel:

Architecture:
   → doParallel/foreach with 4 worker processes
   → Each worker: fits models for one food item × round
   → Models saved to item-specific directories

Prediction Data Generation Pipeline

The figure data generation is the most computationally intensive step. It is orchestrated by a batch processing function that runs across all item × grouping variable combinations in parallel.

Batch Processing Architecture:

   Input: 2 food items × 6 grouping variables = 12 jobs
          (each job = 2 survey rounds)

   For each job (item × grouping_var):
      1. Load pre-fitted models for the item
      2. Load MPCE decile cutpoints for the grouping variable
      3. Construct standardized prediction grid
      4. Draw M=1000 coefficient vectors (MVN) and dispersion parameters
      5. Compute predictions in blocks of 2000 rows
      6. Compute survey standard errors (svydesign)
      7. Inject sampling uncertainty
      8. Summarize to mean + 95% credible interval
      9. Save to item-specific .RData file

   Parallelization:
      → future_lapply() with auto-detected workers
      → Worker count: min(cores-1, memory/2GB per job, n_jobs)
      → Each job: ~2GB peak memory, ~3-5 min per item × group × round

System Resource Management

The pipeline includes adaptive resource management to prevent memory exhaustion:

Output File Organization

Each food item × grouping variable combination produces a prediction file:

~/data/bam_models/excess_food/
   ├── cereal_millets/
   │   ├── cereal_millets_model.RData
   │   ├── data_prop_quantity_cereal_millets_sector.RData
   │   ├── data_prop_quantity_cereal_millets_rel.RData
   │   ├── data_prop_quantity_cereal_millets_social.RData
   │   ├── data_prop_quantity_cereal_millets_child.RData
   │   ├── data_prop_quantity_cereal_millets_female_headed.RData
   │   └── data_prop_quantity_state_cereal_millets_sector.RData
   ├── fats_oils/
   │   └── ... (same structure)
   └── all_excess_food_models.RData

All results are also compiled into a single Excel workbook (data_prop_quantity_excess.xlsx) with two sheets: the full dataset and a variable descriptions sheet.