From 4a66ee5fab8c2acab43f53bb49b20167a7bd64a5 Mon Sep 17 00:00:00 2001 From: bakodramane Date: Wed, 17 Jun 2026 23:28:26 +0200 Subject: [PATCH] feat: add Fay-Herriot measurement-error method for sample-based auxiliaries (fh-me) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Trimmed Phase 7 (v1.1.0). Adds one new area-level method, fh-me (Ybarra & Lohr, 2008), for auxiliary variables that come from a sample rather than a full census, plus the inputs, recommender logic, code generation, wizard changes, and docs to support it. No existing method behaviour changes. - types: add auxiliaryFromSample/hasAuxiliaryVariances to DataAvailability, requiresAuxiliaryVariances + mseMethod jackknife to CatalogueEntry, auxiliaryVarianceVars to UserInputs, and an auxiliary-variance VariableRole - catalogue: new src/catalogue/fh-me.ts (17 entries); schema test accepts jackknife - recommender: rank fh-me first for sample-based continuous/proportion targets; sampling-error caveat on fh-eblup/spatial-fh/robust-fh; blocking caveat on fh-me when variances are missing; hide fh-me for register auxiliaries - codegen: {{AUX_VAR_VARIANCES_R}} and {{CI_ARRAY_BUILDER_R}} tokens + emdi fh(method=me) - wizard: Step 3 auxiliary-source question + follow-up; Step 2 auxiliary-variance pairing; Step 4 badge + prominent amber caveat; Step 5 variance mapping in summary - tests: 3 recommender scenarios, fh-me codegen snapshot, Playwright sample-aux flow - docs: SAE-CATALOGUE.md §7 + tokens, adding-a-method.md, README feature Co-Authored-By: Claude Opus 4.8 --- README.md | 5 +- docs/SAE-CATALOGUE.md | 128 +++++++++++++++++++++++++++++++ docs/adding-a-method.md | 5 +- e2e/wizard.spec.ts | 47 ++++++++++++ src/catalogue/catalogue.test.ts | 7 +- src/catalogue/fh-me.ts | 126 ++++++++++++++++++++++++++++++ src/catalogue/index.ts | 2 + src/engine/codegen.test.ts | 53 +++++++++++++ src/engine/codegen.ts | 33 ++++++++ src/engine/recommender.test.ts | 73 ++++++++++++++++++ src/engine/recommender.ts | 48 ++++++++++++ src/types/index.ts | 12 ++- src/wizard/Step2Roles.tsx | 23 +++++- src/wizard/Step3Availability.tsx | 79 +++++++++++++++++++ src/wizard/Step4Methods.tsx | 48 +++++++++--- src/wizard/Step5Generate.tsx | 20 +++++ src/wizard/WizardContext.tsx | 2 + 17 files changed, 694 insertions(+), 17 deletions(-) create mode 100644 src/catalogue/fh-me.ts diff --git a/README.md b/README.md index 5de79f7..66ad890 100644 --- a/README.md +++ b/README.md @@ -14,9 +14,12 @@ Small area estimation bridges the gap between nationally representative surveys for reliable estimates at district, county, or municipality level. This tool: 1. Lets you describe your data (variable types, what auxiliary data you have, Stata version). -2. Recommends the most appropriate SAE method from a catalogue of 16 methods. +2. Recommends the most appropriate SAE method from a catalogue of 17 methods. 3. Generates a complete, commented R script and Stata `.do` file — ready to run with your real variable names filled in. +4. Supports auxiliary variables that come from a sample (an agricultural census or other large + survey) rather than a full census or register, via the measurement-error Fay–Herriot model + (Ybarra & Lohr, 2008), which corrects for the sampling error those auxiliaries carry. Target users: statisticians in national statistical offices and development organisations, including those working in countries where Stata 14 is the installed standard. diff --git a/docs/SAE-CATALOGUE.md b/docs/SAE-CATALOGUE.md index 7e18371..565d964 100644 --- a/docs/SAE-CATALOGUE.md +++ b/docs/SAE-CATALOGUE.md @@ -17,6 +17,8 @@ Use these tokens consistently across all R and Stata templates: | `{{AUX_VARS_R}}` | Auxiliary variables as R formula terms, e.g. `x1 + x2` | | `{{AUX_VARS_R_VEC}}` | Auxiliary variable names as a quoted R character vector, e.g. `"x1", "x2"` — for column-selection contexts such as `c({{AUX_VARS_R_VEC}})` | | `{{AUX_VARS_STATA}}` | Auxiliary variables space-separated, e.g. `x1 x2` | +| `{{AUX_VAR_VARIANCES_R}}` | Auxiliary sampling-variance column names as a quoted R character vector, e.g. `"var_x1", "var_x2"` | +| `{{CI_ARRAY_BUILDER_R}}` | Generated R code assembling the per-domain measurement-error variance–covariance array `Ci` from the variance columns | | `{{WEIGHT_VAR}}` | Sampling weight variable name | | `{{DIRECT_EST_VAR}}` | Pre-computed direct estimate column | | `{{DIRECT_VAR_VAR}}` | Sampling variance of the direct estimate | @@ -1311,3 +1313,129 @@ Rao, J.N.K. & Molina, I. (2015). *Small Area Estimation*, 2nd ed. Wiley. Sinha, S.K. & Rao, J.N.K. (2009). Robust small area estimation. *Canadian Journal of Statistics* 37, 381–399. Asian Development Bank (2020). *Introduction to Small Area Estimation Techniques: A Practical Guide for National Statistical Offices.* https://www.adb.org/publications/introduction-small-area-estimation-techniques + +--- + +## §7 Survey-based auxiliary variables (added v1.1.0) + +Some users have auxiliary variables that come from a large *sample* (such as an agricultural +census run as a sample) rather than a full census or register. Those auxiliaries carry sampling +error. The standard Fay–Herriot model assumes auxiliaries are known exactly; using it here biases +the model, understates uncertainty, and can perform worse than the direct estimate. This section +adds the measurement-error Fay–Herriot model (Ybarra & Lohr, 2008). + +When the user declares sample-based auxiliaries, the recommender ranks `fh-me` first for +area-level continuous or proportion targets and attaches a prominent caveat to every standard +known-covariate area-level method (`fh-eblup`, `spatial-fh`, `robust-fh`). If the sampling +variances of the auxiliary estimates are not available, `fh-me` carries a blocking caveat because +the correction cannot be applied without them. + +### 7.01 Fay–Herriot with Measurement Error (Ybarra–Lohr) + +``` +id: fh-me +displayName: Fay–Herriot with Measurement Error (Area-Level, Sample-Based Auxiliaries) +level: area +inferenceType: frequentist +targetTypes: [continuous, proportion] +requiredInputs: + microdata: false + areaAggregates: true + censusAuxiliaries: area + weights: false + contiguityMatrix: false + coordinates: false +requiresAuxiliaryVariances: true +spatial: false +robust: false +mseMethod: jackknife +rPackage: emdi (or saeME) +rFunction: fh(method = "me") / saeME::eblupME +stataPackage: base +stataCommand: (no standard Stata command; use R) +stataMinVersion: 14 +stataV14Fallback: | + * There is no standard Stata command for the measurement-error Fay-Herriot model. + * Run the R script for this method (emdi or saeME package). + * If you must stay in Stata, the standard fhsae command ignores auxiliary + * sampling error and may bias the estimates — interpret with caution: + * fhsae {{DIRECT_EST_VAR}} {{AUX_VARS_STATA}}, vardir({{DIRECT_VAR_VAR}}) method(reml) +plainDescription: | + An extension of the Fay–Herriot model for when your auxiliary variables come from a + sample (for example, an agricultural census run as a large sample) rather than a full + census, so the auxiliaries carry their own sampling error. The model adjusts for that + error, leaning more on the direct survey estimate in areas where the auxiliaries are + noisier. It needs the sampling variance of each auxiliary estimate for each area. +whyChooseThis: | + Choose this when your area-level auxiliary totals or means are themselves survey + estimates rather than known population values. Using a standard Fay–Herriot model in + that situation can bias the results and overstate their precision, and can even do worse + than the direct estimate. This method is designed to avoid that. +assumptions: + - Sampling variances of the auxiliary estimates are available for each area. + - The auxiliary measurement error is of the classical type (observed = true + noise). + - All other Fay–Herriot assumptions hold (known direct-estimate variances; correct + linking model; enough areas for variance estimation). +references: + - Ybarra, L.M.R. & Lohr, S.L. (2008). Biometrika 95(4), 919–931. + https://doi.org/10.1093/biomet/asn048 + - Harmening, S. et al. (2023). The R Journal 15(1), RJ-2023-039. + https://journal.r-project.org/articles/RJ-2023-039/ + - saeME package: https://cran.r-project.org/package=saeME +caveats: + - Requires the sampling variances of the auxiliary estimates; without them the + correction cannot be applied. + - MSE is estimated by jackknife only. + - No Stata equivalent; R is required. +``` + +**R template:** +```r +# ============================================================ +# Fay–Herriot with Measurement Error (Ybarra–Lohr) +# For auxiliary variables that come from a sample (e.g. an agricultural +# census run as a large sample), and therefore carry sampling error. +# Generated by SAE Syntax Generator on {{DATE}} +# Reference: Ybarra & Lohr (2008); Harmening et al. (2023) +# R package: emdi +# Area-level data: {{AREA_DATA}} +# ============================================================ + +if (!requireNamespace("emdi", quietly = TRUE)) install.packages("emdi") +library(emdi) + +area_data <- read.csv("{{AREA_DATA}}") +# Required columns: +# {{DIRECT_EST_VAR}} direct estimate of the target, per area +# {{DIRECT_VAR_VAR}} sampling variance of the direct estimate, per area +# {{AUX_VARS_R}} auxiliary estimates (from the large-sample census) +# {{AUX_VAR_VARIANCES_R}} sampling variance of each auxiliary estimate, per area +# {{AREA_ID}} area identifier + +# Build the per-domain measurement-error variance–covariance array Ci. +# emdi expects an array of (p+1) x (p+1) x m, where p = number of auxiliaries, +# the leading row/column is the intercept (zero variance), and off-diagonal +# covariances between auxiliaries are assumed zero unless supplied. +{{CI_ARRAY_BUILDER_R}} + +# Fit the measurement-error Fay–Herriot model +fh_me <- fh( + fixed = {{DIRECT_EST_VAR}} ~ {{AUX_VARS_R}}, + vardir = "{{DIRECT_VAR_VAR}}", + combined_data = area_data, + domains = "{{AREA_ID}}", + method = "me", + Ci = Ci, + MSE = TRUE, + mse_type = "jackknife" +) + +summary(fh_me) +estimators(fh_me, MSE = TRUE, CV = TRUE) + +# NOTE: the modified shrinkage factor leans more on the direct estimate in +# areas where the auxiliary variances are large. Compare against the direct +# estimate (the benchmark) to confirm the model is helping rather than harming. +``` + +**Stata template:** use the `stataV14Fallback` note above (R-only method). diff --git a/docs/adding-a-method.md b/docs/adding-a-method.md index 423fb93..9a86405 100644 --- a/docs/adding-a-method.md +++ b/docs/adding-a-method.md @@ -72,7 +72,8 @@ Open `src/catalogue/my-new-method.ts` and edit each field. The full schema is de |-------|------|-------------| | `spatial` | `boolean` | Does the method explicitly model spatial dependence? | | `robust` | `boolean` | Is the method robust to outliers (M-estimation or similar)? | -| `mseMethod` | `'prasad-rao' \| 'bootstrap' \| 'both' \| 'posterior'` | How mean squared error is estimated | +| `requiresAuxiliaryVariances` | `boolean` (optional) | Set `true` for measurement-error methods that need the sampling variances of the auxiliary estimates (e.g. `fh-me`). Treated as `false` when absent. When `true`, the recommender only offers the method if the user has declared sample-based auxiliaries, and it expects the variance columns to be supplied. | +| `mseMethod` | `'prasad-rao' \| 'bootstrap' \| 'both' \| 'posterior' \| 'jackknife'` | How mean squared error is estimated. Use `'jackknife'` for the measurement-error Fay–Herriot model, whose MSE is estimated by jackknife only. | ### Software @@ -112,6 +113,8 @@ be substituted from user input. | `{{AUX_VARS_R}}` | Auxiliary variables as `var1 + var2 + ...` (R formula syntax) | | `{{AUX_VARS_STATA}}` | Auxiliary variables as `var1 var2 ...` (space-separated) | | `{{AUX_VARS_R_VEC}}` | Auxiliary variables as `c("var1", "var2", ...)` (R vector) | +| `{{AUX_VAR_VARIANCES_R}}` | Auxiliary sampling-variance column names as a quoted R character vector, e.g. `"var_x1", "var_x2"` (measurement-error methods) | +| `{{CI_ARRAY_BUILDER_R}}` | Generated R code that assembles the per-domain measurement-error variance–covariance array `Ci` from the variance columns (used by `fh-me`) | | `{{SURVEY_DATA}}` | Path to the survey CSV file | | `{{AREA_DATA}}` | Path to the area-level CSV file | | `{{CENSUS_DATA}}` | Path to the census CSV file | diff --git a/e2e/wizard.spec.ts b/e2e/wizard.spec.ts index eca4cf2..9128393 100644 --- a/e2e/wizard.spec.ts +++ b/e2e/wizard.spec.ts @@ -93,4 +93,51 @@ test.describe('Wizard happy path', () => { const dl = await downloadPromise expect(dl.suggestedFilename()).toContain('.R') }) + + test('sample-based auxiliaries: FH-ME appears and standard FH shows the sampling-error caveat', async ({ page }) => { + // ── Step 1: Upload codebook ───────────────────────────────────────────── + const fileInput = page.locator('[data-testid="file-input"]') + const fixturePath = path.resolve(__dirname, 'fixtures/test-codebook.csv') + await fileInput.setInputFiles(fixturePath) + await page.getByRole('button', { name: /next/i }).click() + + // ── Step 2: Assign roles ──────────────────────────────────────────────── + await expect(page.getByText('Variable roles')).toBeVisible() + await setRole(page, 'income', 'Target') + await setRole(page, 'edu_rate', 'Auxiliary') + await setRole(page, 'urban_pct', 'Auxiliary') + await setRole(page, 'dir_est', 'Direct estimate') + await setRole(page, 'dir_var', 'Sampling variance') + await page.getByRole('button', { name: /next/i }).click() + + // ── Step 3: Data availability ─────────────────────────────────────────── + await expect(page.getByText('Data availability')).toBeVisible() + + // Continuous target + await page.getByText('Continuous (e.g. income', { exact: false }).first().click() + + // Area aggregates + census auxiliaries so FH-EBLUP and FH-ME are both eligible + await page.getByText('area-level direct estimates', { exact: false }).first().click() + await page.getByText('population-level auxiliary variables', { exact: false }).first().click() + + // Declare sample-based auxiliaries, with variances available + await page.locator('[data-testid="aux-source-sample"]').click() + await page.locator('[data-testid="aux-has-variances"]').check() + + await page.getByRole('button', { name: /next/i }).click() + + // ── Step 4: Methods ───────────────────────────────────────────────────── + await expect(page.getByText('Recommended methods')).toBeVisible() + + // FH-ME card is present and selectable + const fhMeRadio = page.locator('[data-testid="select-fh-me"]') + await expect(fhMeRadio).toBeVisible() + + // Standard FH-EBLUP card shows the prominent sampling-error caveat + await expect(page.locator('[data-testid="sampling-error-note-fh-eblup"]')).toBeVisible() + + // FH-ME is selectable + await fhMeRadio.click() + await expect(fhMeRadio).toBeChecked() + }) }) diff --git a/src/catalogue/catalogue.test.ts b/src/catalogue/catalogue.test.ts index 72724cd..7ae2a87 100644 --- a/src/catalogue/catalogue.test.ts +++ b/src/catalogue/catalogue.test.ts @@ -19,11 +19,12 @@ const EXPECTED_IDS = [ 'glmm-count', 'two-part-zinfl', 'hb-unit', + 'fh-me', ] describe('Catalogue index', () => { - it('exports exactly 16 entries', () => { - expect(catalogue).toHaveLength(16) + it('exports exactly 17 entries', () => { + expect(catalogue).toHaveLength(17) }) it('contains all expected method IDs', () => { @@ -128,7 +129,7 @@ describe('Catalogue schema validation', () => { }) it('has valid mseMethod', () => { - expect(['prasad-rao', 'bootstrap', 'both', 'posterior']).toContain( + expect(['prasad-rao', 'bootstrap', 'both', 'posterior', 'jackknife']).toContain( entry.mseMethod, ) }) diff --git a/src/catalogue/fh-me.ts b/src/catalogue/fh-me.ts new file mode 100644 index 0000000..c2e56a4 --- /dev/null +++ b/src/catalogue/fh-me.ts @@ -0,0 +1,126 @@ +import type { CatalogueEntry } from '../types/index.js' + +const entry: CatalogueEntry = { + id: 'fh-me', + displayName: 'Fay–Herriot with Measurement Error (Area-Level, Sample-Based Auxiliaries)', + level: 'area', + inferenceType: 'frequentist', + targetTypes: ['continuous', 'proportion'], + requiredInputs: { + microdata: false, + areaAggregates: true, + censusAuxiliaries: 'area', + weights: false, + contiguityMatrix: false, + coordinates: false, + }, + requiresAuxiliaryVariances: true, + spatial: false, + robust: false, + mseMethod: 'jackknife', + rPackage: 'emdi (or saeME)', + rFunction: 'fh(method = "me") / saeME::eblupME', + stataPackage: 'base', + stataCommand: '(no standard Stata command; use R)', + stataMinVersion: 14, + stataV14Fallback: + '* There is no standard Stata command for the measurement-error Fay-Herriot model.\n' + + '* Run the R script for this method (emdi or saeME package).\n' + + '* If you must stay in Stata, the standard fhsae command ignores auxiliary\n' + + '* sampling error and may bias the estimates — interpret with caution:\n' + + '* fhsae {{DIRECT_EST_VAR}} {{AUX_VARS_STATA}}, vardir({{DIRECT_VAR_VAR}}) method(reml)', + plainDescription: + 'An extension of the Fay–Herriot model for when your auxiliary variables come from a ' + + 'sample (for example, an agricultural census run as a large sample) rather than a full ' + + 'census, so the auxiliaries carry their own sampling error. The model adjusts for that ' + + 'error, leaning more on the direct survey estimate in areas where the auxiliaries are ' + + 'noisier. It needs the sampling variance of each auxiliary estimate for each area.', + whyChooseThis: + 'Choose this when your area-level auxiliary totals or means are themselves survey ' + + 'estimates rather than known population values. Using a standard Fay–Herriot model in ' + + 'that situation can bias the results and overstate their precision, and can even do worse ' + + 'than the direct estimate. This method is designed to avoid that.', + assumptions: [ + 'Sampling variances of the auxiliary estimates are available for each area.', + 'The auxiliary measurement error is of the classical type (observed = true + noise).', + 'All other Fay–Herriot assumptions hold (known direct-estimate variances; correct ' + + 'linking model; enough areas for variance estimation).', + ], + references: [ + 'Ybarra, L.M.R. & Lohr, S.L. (2008). Biometrika 95(4), 919–931. https://doi.org/10.1093/biomet/asn048', + 'Harmening, S. et al. (2023). The R Journal 15(1), RJ-2023-039. https://journal.r-project.org/articles/RJ-2023-039/', + 'saeME package: https://cran.r-project.org/package=saeME', + ], + caveats: [ + 'Requires the sampling variances of the auxiliary estimates; without them the ' + + 'correction cannot be applied.', + 'MSE is estimated by jackknife only.', + 'No Stata equivalent; R is required.', + ], + rTemplate: `# ============================================================ +# Fay–Herriot with Measurement Error (Ybarra–Lohr) +# For auxiliary variables that come from a sample (e.g. an agricultural +# census run as a large sample), and therefore carry sampling error. +# Generated by SAE Syntax Generator on {{DATE}} +# Reference: Ybarra & Lohr (2008); Harmening et al. (2023) +# R package: emdi +# Area-level data: {{AREA_DATA}} +# ============================================================ + +if (!requireNamespace("emdi", quietly = TRUE)) install.packages("emdi") +library(emdi) + +area_data <- read.csv("{{AREA_DATA}}") +# Required columns: +# {{DIRECT_EST_VAR}} direct estimate of the target, per area +# {{DIRECT_VAR_VAR}} sampling variance of the direct estimate, per area +# {{AUX_VARS_R}} auxiliary estimates (from the large-sample census) +# {{AUX_VAR_VARIANCES_R}} sampling variance of each auxiliary estimate, per area +# {{AREA_ID}} area identifier + +# Build the per-domain measurement-error variance–covariance array Ci. +# emdi expects an array of (p+1) x (p+1) x m, where p = number of auxiliaries, +# the leading row/column is the intercept (zero variance), and off-diagonal +# covariances between auxiliaries are assumed zero unless supplied. +{{CI_ARRAY_BUILDER_R}} + +# Fit the measurement-error Fay–Herriot model +fh_me <- fh( + fixed = {{DIRECT_EST_VAR}} ~ {{AUX_VARS_R}}, + vardir = "{{DIRECT_VAR_VAR}}", + combined_data = area_data, + domains = "{{AREA_ID}}", + method = "me", + Ci = Ci, + MSE = TRUE, + mse_type = "jackknife" +) + +summary(fh_me) +estimators(fh_me, MSE = TRUE, CV = TRUE) + +# NOTE: the modified shrinkage factor leans more on the direct estimate in +# areas where the auxiliary variances are large. Compare against the direct +# estimate (the benchmark) to confirm the model is helping rather than harming. +`, + stataTemplate: `* ============================================================ +* Fay–Herriot with Measurement Error — No Standard Stata Equivalent +* Generated by SAE Syntax Generator on {{DATE}} +* Please use the R script (emdi or saeME package) for this method. +* Area-level data: {{AREA_DATA}} +* ============================================================ + +* There is no standard Stata command for the measurement-error Fay–Herriot model. +* The correction for sampling error in the auxiliary variables requires the R +* implementation (emdi::fh(method = "me") or saeME::eblupME). +* +* If you must stay in Stata, the standard fhsae command ignores auxiliary +* sampling error and may bias the estimates — interpret with caution: +* +* net install fhsae, from("https://raw.github.com/jpazvd/fhsae/master/") +* use "{{AREA_DATA}}", clear +* fhsae {{DIRECT_EST_VAR}} {{AUX_VARS_STATA}}, vardir({{DIRECT_VAR_VAR}}) method(reml) +`, +} + +export default entry diff --git a/src/catalogue/index.ts b/src/catalogue/index.ts index f10c486..aafb0ed 100644 --- a/src/catalogue/index.ts +++ b/src/catalogue/index.ts @@ -15,6 +15,7 @@ import glmmBinary from './glmm-binary.js' import glmmCount from './glmm-count.js' import twoPartZinfl from './two-part-zinfl.js' import hbUnit from './hb-unit.js' +import fhMe from './fh-me.js' export const catalogue: CatalogueEntry[] = [ direct, @@ -33,6 +34,7 @@ export const catalogue: CatalogueEntry[] = [ glmmCount, twoPartZinfl, hbUnit, + fhMe, ] export default catalogue diff --git a/src/engine/codegen.test.ts b/src/engine/codegen.test.ts index 24be9bd..0e57a5a 100644 --- a/src/engine/codegen.test.ts +++ b/src/engine/codegen.test.ts @@ -5,6 +5,7 @@ import fhEblup from '../catalogue/fh-eblup.js' import bhfEblup from '../catalogue/bhf-eblup.js' import ebpCensuseb from '../catalogue/ebp-censuseb.js' import glmmBinary from '../catalogue/glmm-binary.js' +import fhMe from '../catalogue/fh-me.js' // ── Helpers ──────────────────────────────────────────────────────────────────── @@ -230,6 +231,58 @@ describe('generateCode — glmm-binary (binary target)', () => { }) }) +// ── FH-ME (measurement error, sample-based auxiliaries) ──────────────────────── +describe('generateCode — fh-me (measurement-error Fay–Herriot)', () => { + const FH_ME_INPUTS: UserInputs = { + targetVar: 'mean_yield', + areaIdVar: 'district', + auxiliaryVars: ['agcensus_area', 'fertiliser_use'], + auxiliaryVarianceVars: ['var_agcensus_area', 'var_fertiliser_use'], + directEstVar: 'dir_est', + directVarVar: 'dir_var', + areaDataPath: 'area_data.csv', + stataVersion: 14, + } + const code = generateCode(fhMe, FH_ME_INPUTS) + + it('R script loads the emdi package', () => { + expect(code.r).toContain('library(emdi)') + }) + + it('R script calls fh(... method = "me" ...)', () => { + expect(code.r).toContain('fh(') + expect(code.r).toContain('method = "me"') + }) + + it('R script references all auxiliary variance columns', () => { + assertAllVarsPresent(code.r, ['var_agcensus_area', 'var_fertiliser_use'], 'R') + }) + + it('R script references all auxiliary columns and direct estimate columns', () => { + assertAllVarsPresent(code.r, ['agcensus_area', 'fertiliser_use', 'dir_est', 'dir_var'], 'R') + }) + + it('R script builds the Ci array', () => { + expect(code.r).toContain('Ci <- array(0') + expect(code.r).toContain('aux_var_cols <- c("var_agcensus_area", "var_fertiliser_use")') + }) + + it('R script documents the zero-covariance assumption', () => { + expect(code.r.toLowerCase()).toContain('covariances between auxiliaries are assumed') + }) + + it('R script has no unreplaced {{ tokens', () => assertNoUnreplacedTokens(code.r, 'R')) + it('Stata script has no unreplaced {{ tokens', () => assertNoUnreplacedTokens(code.stata, 'Stata')) + + it('Stata output is an R-only explanatory note', () => { + expect(code.stata.toLowerCase()).toContain('use the r script') + }) + + it('does not use fallback at Stata 14 (stataMinVersion 14)', () => { + expect(code.usedFallback).toBe(false) + }) +}) + // ── v14 fallback — general property tests ───────────────────────────────────── describe('generateCode — Stata v14 fallback logic', () => { it('usedFallback true when stataVersion < stataMinVersion and fallback exists', () => { diff --git a/src/engine/codegen.ts b/src/engine/codegen.ts index 9e45bd3..b725571 100644 --- a/src/engine/codegen.ts +++ b/src/engine/codegen.ts @@ -4,6 +4,9 @@ export interface UserInputs { targetVar: string areaIdVar: string auxiliaryVars: string[] + // Column names holding the sampling variance of each auxiliary estimate, in the + // same order as auxiliaryVars. Used by the measurement-error Fay–Herriot model. + auxiliaryVarianceVars?: string[] weightVar?: string directEstVar?: string directVarVar?: string @@ -25,6 +28,32 @@ export interface GeneratedCode { fallbackNote?: string } +// ── Ci array builder (measurement-error Fay–Herriot) ───────────────────────────── +// Generates the R code that assembles the per-domain measurement-error +// variance–covariance array `Ci` that emdi::fh(method = "me") requires. Ci has +// dimension (p+1) x (p+1) x m, where p is the number of auxiliary variables and m +// is the number of areas. The leading row/column is the intercept (zero variance), +// and off-diagonal covariances between auxiliaries are assumed zero unless +// covariance columns are supplied. +function buildCiArrayBuilderR(auxVarVarsRVec: string): string { + return `# Build the per-domain measurement-error variance–covariance array Ci. +# Ci has dimension (p+1) x (p+1) x m, where p = number of auxiliary variables and +# m = number of areas. Position 1 of each matrix is the intercept and carries zero +# variance; positions 2..(p+1) hold the sampling variance of each auxiliary. +# Off-diagonal covariances between auxiliaries are assumed to be zero, because no +# covariance columns were supplied. If your auxiliary estimates are correlated, +# fill in the relevant off-diagonal entries of Ci by hand. +aux_var_cols <- c(${auxVarVarsRVec}) +m <- nrow(area_data) +p <- length(aux_var_cols) +Ci <- array(0, dim = c(p + 1, p + 1, m)) +for (i in seq_len(m)) { + for (j in seq_len(p)) { + Ci[j + 1, j + 1, i] <- area_data[i, aux_var_cols[j]] + } +}` +} + // ── Token substitution ───────────────────────────────────────────────────────── function buildTokenMap(inputs: UserInputs): Record { @@ -32,6 +61,8 @@ function buildTokenMap(inputs: UserInputs): Record { const auxR = inputs.auxiliaryVars.join(' + ') const auxRVec = inputs.auxiliaryVars.map(v => `"${v}"`).join(', ') const auxStata = inputs.auxiliaryVars.join(' ') + const auxVarVars = inputs.auxiliaryVarianceVars ?? [] + const auxVarVarsRVec = auxVarVars.map(v => `"${v}"`).join(', ') const survey = inputs.surveyDataPath ?? 'survey.csv' const census = inputs.censusDataPath ?? 'census.csv' const area = inputs.areaDataPath ?? 'area_data.csv' @@ -43,6 +74,8 @@ function buildTokenMap(inputs: UserInputs): Record { AUX_VARS_R: auxR, AUX_VARS_R_VEC: auxRVec, AUX_VARS_STATA: auxStata, + AUX_VAR_VARIANCES_R: auxVarVarsRVec, + CI_ARRAY_BUILDER_R: buildCiArrayBuilderR(auxVarVarsRVec), WEIGHT_VAR: inputs.weightVar ?? '', DIRECT_EST_VAR: inputs.directEstVar ?? '', DIRECT_VAR_VAR: inputs.directVarVar ?? '', diff --git a/src/engine/recommender.test.ts b/src/engine/recommender.test.ts index 53819a4..e0decc6 100644 --- a/src/engine/recommender.test.ts +++ b/src/engine/recommender.test.ts @@ -14,6 +14,8 @@ const base: DataAvailability = { targetType: 'continuous', likelyOutliers: false, likelySpatialCorrelation: false, + auxiliaryFromSample: false, + hasAuxiliaryVariances: false, } function ids(recs: ReturnType): string[] { @@ -307,3 +309,74 @@ describe('Scenario 15 — purity and determinism', () => { expect(a.map(r => r.rank)).toEqual(b.map(r => r.rank)) }) }) + +// ── Scenario 16: Sample-based auxiliaries with variances → FH-ME first ───────── +describe('Scenario 16 — sample auxiliaries, variances available', () => { + const avail: DataAvailability = { + ...base, + hasAreaAggregates: true, + hasCensusAuxiliaries: true, + targetType: 'continuous', + auxiliaryFromSample: true, + hasAuxiliaryVariances: true, + } + const recs = recommend(avail) + + it('FH-ME is present', () => expect(ids(recs)).toContain('fh-me')) + it('FH-ME ranks first', () => expect(recs[0].entry.id).toBe('fh-me')) + it('FH-ME ranks above standard FH-EBLUP', () => + expect(rankOf(recs, 'fh-me')).toBeLessThan(rankOf(recs, 'fh-eblup'))) + it('standard FH-EBLUP carries the sampling-error caveat', () => { + const fh = recs.find(r => r.entry.id === 'fh-eblup')! + expect(fh.caveats.some(c => c.includes('come from a sample'))).toBe(true) + }) + it('FH-ME does NOT carry the blocking caveat (variances available)', () => { + const fhMe = recs.find(r => r.entry.id === 'fh-me')! + expect(fhMe.caveats.some(c => c.includes('provide the variance columns'))).toBe(false) + }) + it('direct estimator ranks last', () => expect(recs[recs.length - 1].entry.id).toBe('direct')) +}) + +// ── Scenario 17: Sample-based auxiliaries, variances missing → blocking caveat ── +describe('Scenario 17 — sample auxiliaries, variances missing', () => { + const avail: DataAvailability = { + ...base, + hasAreaAggregates: true, + hasCensusAuxiliaries: true, + targetType: 'proportion', + auxiliaryFromSample: true, + hasAuxiliaryVariances: false, + } + const recs = recommend(avail) + + it('FH-ME is present', () => expect(ids(recs)).toContain('fh-me')) + it('FH-ME carries the blocking caveat', () => { + const fhMe = recs.find(r => r.entry.id === 'fh-me')! + expect(fhMe.caveats.some(c => c.includes('provide the variance columns'))).toBe(true) + }) + it('standard FH-EBLUP still carries the sampling-error caveat', () => { + const fh = recs.find(r => r.entry.id === 'fh-eblup')! + expect(fh.caveats.some(c => c.includes('come from a sample'))).toBe(true) + }) +}) + +// ── Scenario 18: Register auxiliaries → no change, FH-ME hidden ───────────────── +describe('Scenario 18 — register auxiliaries (known exactly)', () => { + const avail: DataAvailability = { + ...base, + hasAreaAggregates: true, + hasCensusAuxiliaries: true, + targetType: 'continuous', + auxiliaryFromSample: false, + hasAuxiliaryVariances: false, + } + const recs = recommend(avail) + + it('FH-ME is NOT present', () => expect(ids(recs)).not.toContain('fh-me')) + it('FH-EBLUP ranks first (unchanged behaviour)', () => expect(recs[0].entry.id).toBe('fh-eblup')) + it('FH-EBLUP carries no sampling-error caveat', () => { + const fh = recs.find(r => r.entry.id === 'fh-eblup')! + expect(fh.caveats.some(c => c.includes('come from a sample'))).toBe(false) + }) + it('direct estimator ranks last', () => expect(recs[recs.length - 1].entry.id).toBe('direct')) +}) diff --git a/src/engine/recommender.ts b/src/engine/recommender.ts index 05b67b6..120250f 100644 --- a/src/engine/recommender.ts +++ b/src/engine/recommender.ts @@ -10,6 +10,22 @@ export interface Recommendation { } const DIRECT_ID = 'direct' +const FH_ME_ID = 'fh-me' + +// Standard area-level methods that assume the auxiliary variables are known exactly. +// When the user declares sample-based auxiliaries, these carry a sampling-error caveat. +const KNOWN_COVARIATE_AREA_METHODS = ['fh-eblup', 'spatial-fh', 'robust-fh'] + +const SAMPLE_AUX_CAVEAT = + 'Your auxiliary variables come from a sample and carry sampling error. This method ' + + 'assumes they are known exactly, which can bias the estimates and understate their ' + + 'uncertainty. Prefer the measurement-error model (FH-ME) unless the auxiliary sample ' + + 'is very large.' + +const FH_ME_BLOCKING_CAVEAT = + 'The measurement-error correction needs the sampling variances of your auxiliary ' + + 'estimates. Without them this method cannot be applied — provide the variance columns, ' + + 'or treat the auxiliaries as approximate and interpret results with caution.' // ── Eligibility ──────────────────────────────────────────────────────────────── @@ -34,6 +50,10 @@ function isEligible(entry: CatalogueEntry, availability: DataAvailability): bool if (req.coordinates && !availability.hasCoordinates) return false if (req.censusAuxiliaries !== 'none' && !availability.hasCensusAuxiliaries) return false + // Measurement-error methods only make sense when auxiliaries come from a sample. + // Hide them entirely when the user has register/census auxiliaries known exactly. + if (entry.requiresAuxiliaryVariances && !availability.auxiliaryFromSample) return false + return true } @@ -57,6 +77,17 @@ function computeScore(entry: CatalogueEntry, availability: DataAvailability): nu if (id === DIRECT_ID) return 9999 + // ── Sample-based auxiliaries ────────────────────────────────────────────────── + // When the auxiliaries carry sampling error, the measurement-error Fay–Herriot + // model ranks first for area-level continuous or proportion targets. + if ( + availability.auxiliaryFromSample && + id === FH_ME_ID && + (targetType === 'continuous' || targetType === 'proportion' || targetType === 'unknown') + ) { + return 5 + } + // ── Poverty target ────────────────────────────────────────────────────────── if (targetType === 'poverty') { if (id === 'ebp-censuseb') return 10 @@ -190,6 +221,12 @@ function buildWhyApplicable(entry: CatalogueEntry, availability: DataAvailabilit return 'Use as a benchmark to validate model-based estimates against the raw survey data.' } + if (id === FH_ME_ID) { + return availability.auxiliaryFromSample + ? 'Ranked first because your auxiliary variables come from a sample: this measurement-error Fay–Herriot model corrects for the sampling error they carry, leaning more on the direct estimate where the auxiliaries are noisier.' + : 'Measurement-error Fay–Herriot model for auxiliary variables that come from a sample rather than a full census or register.' + } + // Poverty-specific explanations if (targetType === 'poverty') { if (id === 'ebp-censuseb') { @@ -351,6 +388,17 @@ function buildCaveats(entry: CatalogueEntry, availability: DataAvailability): st const caveats: string[] = [...(entry.caveats ?? [])] const { likelyOutliers, likelySpatialCorrelation, hasOutOfSampleAreas } = availability + // Sample-based auxiliaries: warn the standard known-covariate area-level methods, + // and block FH-ME itself if the required auxiliary variances are not available. + if (availability.auxiliaryFromSample) { + if (KNOWN_COVARIATE_AREA_METHODS.includes(entry.id)) { + caveats.push(SAMPLE_AUX_CAVEAT) + } + if (entry.id === FH_ME_ID && !availability.hasAuxiliaryVariances) { + caveats.push(FH_ME_BLOCKING_CAVEAT) + } + } + // Outlier warning for non-robust methods if (likelyOutliers && !entry.robust && entry.id !== DIRECT_ID) { caveats.push( diff --git a/src/types/index.ts b/src/types/index.ts index 4adaef9..440d26d 100644 --- a/src/types/index.ts +++ b/src/types/index.ts @@ -3,6 +3,7 @@ export type VariableRole = | 'area-id' | 'weight' | 'auxiliary' + | 'auxiliary-variance' | 'coordinate' | 'direct-est' | 'direct-var' @@ -37,6 +38,12 @@ export interface DataAvailability { targetType: 'continuous' | 'binary' | 'proportion' | 'count' | 'poverty' | 'unknown' likelyOutliers: boolean likelySpatialCorrelation: boolean + // True when auxiliary variables are estimated from a survey or large sample + // (e.g. an agricultural census run as a sample) and therefore carry sampling error. + auxiliaryFromSample: boolean + // True when the user can supply the sampling variances of the auxiliary estimates, + // area by area. Required by the measurement-error Fay–Herriot model (fh-me). + hasAuxiliaryVariances: boolean } export interface CatalogueEntry { @@ -55,7 +62,10 @@ export interface CatalogueEntry { } spatial: boolean robust: boolean - mseMethod: 'prasad-rao' | 'bootstrap' | 'both' | 'posterior' + // True when the method needs the sampling variances of the auxiliary estimates + // (the measurement-error Fay–Herriot model). Treated as false when absent. + requiresAuxiliaryVariances?: boolean + mseMethod: 'prasad-rao' | 'bootstrap' | 'both' | 'posterior' | 'jackknife' rPackage: string rFunction: string rTemplate: string diff --git a/src/wizard/Step2Roles.tsx b/src/wizard/Step2Roles.tsx index 3fda446..a6f1fa7 100644 --- a/src/wizard/Step2Roles.tsx +++ b/src/wizard/Step2Roles.tsx @@ -7,6 +7,7 @@ const ROLES: { value: VariableRole; label: string; tip: string }[] = [ { value: 'target', label: 'Target', tip: 'The outcome variable you want to estimate for small areas (e.g. poverty rate, income).' }, { value: 'area-id', label: 'Area identifier', tip: 'The column that identifies which small area each record belongs to (e.g. district code).' }, { value: 'auxiliary', label: 'Auxiliary', tip: 'A covariate available in both the survey and the population data that is correlated with the target.' }, + { value: 'auxiliary-variance', label: 'Auxiliary variance', tip: 'The sampling variance of a sample-based auxiliary estimate, area by area. Pair it with its auxiliary, in the same order. Required by the measurement-error Fay–Herriot model.' }, { value: 'weight', label: 'Survey weight', tip: 'The sampling weight for each unit, reflecting the probability of selection.' }, { value: 'direct-est', label: 'Direct estimate', tip: 'A pre-computed area-level estimate (e.g. from the direct estimator). Used only by area-level methods.' }, { value: 'direct-var', label: 'Sampling variance', tip: 'The estimated sampling variance of the direct estimate. Required by Fay–Herriot models.' }, @@ -49,14 +50,25 @@ export function Step2Roles() { dispatch({ type: 'SET_VARIABLES', payload: updated }) } + const fromSample = state.availability.auxiliaryFromSample + // The auxiliary-variance role is only relevant for sample-based auxiliaries. + const availableRoles = fromSample ? ROLES : ROLES.filter(r => r.value !== 'auxiliary-variance') + const targets = state.variables.filter(v => v.role === 'target') const areaIds = state.variables.filter(v => v.role === 'area-id') const auxVars = state.variables.filter(v => v.role === 'auxiliary') + const auxVarVars = state.variables.filter(v => v.role === 'auxiliary-variance') const errors: string[] = [] if (targets.length !== 1) errors.push(`Exactly one variable must be assigned the "Target" role (currently ${targets.length}).`) if (areaIds.length < 1) errors.push('At least one variable must be assigned the "Area identifier" role.') if (auxVars.length < 1) errors.push('At least one variable must be assigned the "Auxiliary" role.') + if (fromSample && auxVarVars.length > 0 && auxVarVars.length !== auxVars.length) { + errors.push( + `Pair each auxiliary with its variance column: ${auxVars.length} auxiliary variable(s) ` + + `but ${auxVarVars.length} "Auxiliary variance" column(s). They must match, in order.`, + ) + } return (
@@ -67,6 +79,15 @@ export function Step2Roles() {

+ {fromSample && ( +
+ Your auxiliary variables come from a sample. If you can supply their sampling variances, + assign each variance column the "Auxiliary variance"{' '} + role, in the same order as the auxiliaries they correspond to. The measurement-error + Fay–Herriot model uses these to correct for sampling error. +
+ )} +
@@ -96,7 +117,7 @@ export function Step2Roles() { onChange={e => setRole(i, e.target.value as VariableRole)} className="border border-gray-300 rounded px-2 py-1 text-xs w-full focus:ring-2 focus:ring-indigo-500 focus:outline-none" > - {ROLES.map(r => ( + {availableRoles.map(r => ( ))} diff --git a/src/wizard/Step3Availability.tsx b/src/wizard/Step3Availability.tsx index b01cfe1..15468dd 100644 --- a/src/wizard/Step3Availability.tsx +++ b/src/wizard/Step3Availability.tsx @@ -80,6 +80,25 @@ export function Step3Availability() { dispatch({ type: 'SET_AVAILABILITY', payload: { ...availability, targetType: v } }) } + function setAuxiliaryFromSample(fromSample: boolean) { + dispatch({ + type: 'SET_AVAILABILITY', + payload: { + ...availability, + auxiliaryFromSample: fromSample, + // Reset the follow-up answer when switching back to a register source. + hasAuxiliaryVariances: fromSample ? availability.hasAuxiliaryVariances : false, + }, + }) + } + + function setHasAuxiliaryVariances(has: boolean) { + dispatch({ + type: 'SET_AVAILABILITY', + payload: { ...availability, hasAuxiliaryVariances: has }, + }) + } + return (
@@ -140,6 +159,66 @@ export function Step3Availability() {
+ {/* Auxiliary variable source */} +
+ + Where do your auxiliary variables come from?{' '} + + + + +
+ + +
+ + {/* Follow-up: only shown when auxiliaries come from a sample */} + {availability.auxiliaryFromSample && ( +
+ +
+ )} +
+ +
+ {/* Stata version */}
- {/* Caveats */} - {rec.caveats.length > 0 && ( -
    - {rec.caveats.map((c, i) => ( -
  • - •{c} -
    • - ))} -
- )} + {/* Prominent sampling-error caveats — shown as a visible amber note */} + {(() => { + const isProminent = (c: string) => + c.includes('come from a sample') || c.includes('provide the variance columns') + const prominent = rec.caveats.filter(isProminent) + const rest = rec.caveats.filter(c => !isProminent(c)) + return ( + <> + {prominent.length > 0 && ( +
+ {prominent.map((c, i) => ( +

+ ⚠{c} +

+ ))} +
+ )} + {rest.length > 0 && ( +
    + {rest.map((c, i) => ( +
  • + •{c} +
    • + ))} +
+ )} + + ) + })()} {/* Why this method — expandable */}
diff --git a/src/wizard/Step5Generate.tsx b/src/wizard/Step5Generate.tsx index ae81c7a..f589b56 100644 --- a/src/wizard/Step5Generate.tsx +++ b/src/wizard/Step5Generate.tsx @@ -22,6 +22,9 @@ function buildInputs(state: ReturnType['state']): UserInputs { const targetVar = state.variables.find(v => v.role === 'target')?.name ?? '' const areaIdVar = state.variables.find(v => v.role === 'area-id')?.name ?? '' const auxiliaryVars = state.variables.filter(v => v.role === 'auxiliary').map(v => v.name) + const auxiliaryVarianceVars = state.variables + .filter(v => v.role === 'auxiliary-variance') + .map(v => v.name) const weightVar = state.variables.find(v => v.role === 'weight')?.name const directEstVar = state.variables.find(v => v.role === 'direct-est')?.name const directVarVar = state.variables.find(v => v.role === 'direct-var')?.name @@ -30,6 +33,7 @@ function buildInputs(state: ReturnType['state']): UserInputs { targetVar, areaIdVar, auxiliaryVars, + auxiliaryVarianceVars: auxiliaryVarianceVars.length > 0 ? auxiliaryVarianceVars : undefined, weightVar, directEstVar, directVarVar, @@ -188,6 +192,8 @@ export function Step5Generate() { const targetVar = variables.find(v => v.role === 'target') const areaIdVar = variables.find(v => v.role === 'area-id') const auxVars = variables.filter(v => v.role === 'auxiliary') + const auxVarVars = variables.filter(v => v.role === 'auxiliary-variance') + const showAuxVariances = selectedMethodId === 'fh-me' || auxVarVars.length > 0 return (
@@ -205,6 +211,20 @@ export function Step5Generate() {

Auxiliary variables:{' '} {auxVars.map(v => {v.name})}

+ {showAuxVariances && ( +

+ Auxiliary variance columns:{' '} + {auxVarVars.length > 0 + ? auxVarVars.map((v, i) => ( + + {v.name} + ↔ {auxVars[i]?.name ?? '—'} + {i < auxVarVars.length - 1 ? , : null} + + )) + : none assigned — FH-ME needs the sampling variance of each auxiliary} +

+ )}

Stata version: {state.stataVersion}

diff --git a/src/wizard/WizardContext.tsx b/src/wizard/WizardContext.tsx index 02baafb..28c2013 100644 --- a/src/wizard/WizardContext.tsx +++ b/src/wizard/WizardContext.tsx @@ -47,6 +47,8 @@ const defaultAvailability: DataAvailability = { targetType: 'unknown', likelyOutliers: false, likelySpatialCorrelation: false, + auxiliaryFromSample: false, + hasAuxiliaryVariances: false, } const initialState: WizardState = {