diff --git a/docs/SAE-CATALOGUE.md b/docs/SAE-CATALOGUE.md index 0cfd5c6..7e18371 100644 --- a/docs/SAE-CATALOGUE.md +++ b/docs/SAE-CATALOGUE.md @@ -15,6 +15,7 @@ Use these tokens consistently across all R and Stata templates: | `{{TARGET_VAR}}` | Name of the target (outcome) variable | | `{{AREA_ID}}` | Name of the small area identifier variable | | `{{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` | | `{{WEIGHT_VAR}}` | Sampling weight variable name | | `{{DIRECT_EST_VAR}}` | Pre-computed direct estimate column | diff --git a/src/catalogue/bhf-eblup.ts b/src/catalogue/bhf-eblup.ts index 129c5fc..3cc1840 100644 --- a/src/catalogue/bhf-eblup.ts +++ b/src/catalogue/bhf-eblup.ts @@ -61,7 +61,7 @@ census_data <- read.csv("{{CENSUS_DATA}}") # census_data must have: {{AREA_ID}}, population size (Ni), and means of {{AUX_VARS_R}} # Area population means matrix (columns = auxiliary variable means, rows = areas) -meanxpop <- census_data[, c("{{AREA_ID}}", {{AUX_VARS_R_QUOTED}})] +meanxpop <- census_data[, c("{{AREA_ID}}", {{AUX_VARS_R_VEC}})] popnsize <- data.frame(area = census_data\${{AREA_ID}}, Ni = census_data$Ni) diff --git a/src/catalogue/hb-unit.ts b/src/catalogue/hb-unit.ts index 9522ad4..906237f 100644 --- a/src/catalogue/hb-unit.ts +++ b/src/catalogue/hb-unit.ts @@ -71,7 +71,7 @@ hb_result <- hbNormal( formula = {{TARGET_VAR}} ~ {{AUX_VARS_R}}, dom = survey_data\${{AREA_ID}}, Ni = census_data$Ni, - Xmean = census_data[, c({{AUX_VARS_R_QUOTED}})], + Xmean = census_data[, c({{AUX_VARS_R_VEC}})], data = survey_data, iter.update = 3, iter.mcmc = {{N_SIM}}, diff --git a/src/catalogue/reblup.ts b/src/catalogue/reblup.ts index ee43ff3..7a56672 100644 --- a/src/catalogue/reblup.ts +++ b/src/catalogue/reblup.ts @@ -64,7 +64,7 @@ reblup_result <- reblup( formula = {{TARGET_VAR}} ~ {{AUX_VARS_R}}, domain = survey_data\${{AREA_ID}}, sizedom = census_data$Ni, - meanxpop = census_data[, c({{AUX_VARS_R_QUOTED}})], + meanxpop = census_data[, c({{AUX_VARS_R_VEC}})], data = survey_data ) diff --git a/src/engine/codegen.test.ts b/src/engine/codegen.test.ts new file mode 100644 index 0000000..1014804 --- /dev/null +++ b/src/engine/codegen.test.ts @@ -0,0 +1,293 @@ +import { describe, it, expect } from 'vitest' +import { generateCode } from './codegen.js' +import type { UserInputs } from './codegen.js' +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' + +// ── Helpers ──────────────────────────────────────────────────────────────────── + +function assertNoUnreplacedTokens(script: string, label: string): void { + const match = script.match(/\{\{[A-Z_]+\}\}/) + expect(match, `${label} still contains unreplaced token: ${match?.[0]}`).toBeNull() +} + +function assertAllVarsPresent(script: string, vars: string[], label: string): void { + for (const v of vars) { + expect(script, `${label} missing variable "${v}"`).toContain(v) + } +} + +// ── Shared variable sets ─────────────────────────────────────────────────────── + +const FH_INPUTS: UserInputs = { + targetVar: 'poverty_rate', + areaIdVar: 'district', + auxiliaryVars: ['edu_rate', 'urban_pct'], + directEstVar: 'dir_est', + directVarVar: 'dir_var', + surveyDataPath: 'survey.csv', + censusDataPath: 'area_data.csv', + stataVersion: 14, + nSimulations: 200, + mseMethod: 'analytic', +} + +const BHF_INPUTS: UserInputs = { + targetVar: 'income', + areaIdVar: 'region', + auxiliaryVars: ['x1', 'x2', 'x3'], + surveyDataPath: 'survey.csv', + censusDataPath: 'census.csv', + stataVersion: 14, + nSimulations: 200, +} + +const EBP_INPUTS: UserInputs = { + targetVar: 'cons_pc', + areaIdVar: 'ea_code', + auxiliaryVars: ['hhsize', 'edu_head', 'urban'], + surveyDataPath: 'survey.csv', + censusDataPath: 'census.csv', + stataVersion: 14, // below stataMinVersion:17 → triggers fallback + nSimulations: 200, +} + +const GLMM_INPUTS: UserInputs = { + targetVar: 'employed', + areaIdVar: 'province', + auxiliaryVars: ['age', 'education'], + surveyDataPath: 'survey.csv', + censusDataPath: 'census.csv', + stataVersion: 14, + nSimulations: 200, +} + +// ── Fay–Herriot EBLUP ───────────────────────────────────────────────────────── +describe('generateCode — fh-eblup (area-level, standard)', () => { + const code = generateCode(fhEblup, FH_INPUTS) + + it('R script contains library(sae)', () => { + expect(code.r).toContain('library(sae)') + }) + + it('R script contains library(', () => { + expect(code.r).toContain('library(') + }) + + it('Stata script contains the fhsae command', () => { + expect(code.stata).toContain('fhsae') + }) + + it('R script contains all FH-relevant variable names', () => { + // FH is area-level: uses direct estimate + variance columns, not raw targetVar + assertAllVarsPresent(code.r, ['district', 'edu_rate', 'urban_pct', 'dir_est', 'dir_var'], 'R') + }) + + it('Stata script contains all FH-relevant variable names', () => { + // fhsae does not take an area-id argument — district is absent from Stata template + assertAllVarsPresent(code.stata, ['edu_rate', 'urban_pct', 'dir_est', 'dir_var'], 'Stata') + }) + + it('R script has no unreplaced {{ tokens', () => assertNoUnreplacedTokens(code.r, 'R')) + it('Stata script has no unreplaced {{ tokens', () => assertNoUnreplacedTokens(code.stata, 'Stata')) + + it('does not use fallback (stataMinVersion 14)', () => { + expect(code.usedFallback).toBe(false) + expect(code.fallbackNote).toBeUndefined() + }) + + it('R script contains MSE footer', () => { + expect(code.r).toContain('MSE / RMSE INTERPRETATION') + }) + + it('Stata script contains MSE footer', () => { + expect(code.stata).toContain('MSE / RMSE INTERPRETATION') + }) + + it('AUX_VARS_R produces formula form (+ separated)', () => { + expect(code.r).toContain('edu_rate + urban_pct') + }) + + it('AUX_VARS_STATA produces space-separated form', () => { + expect(code.stata).toContain('edu_rate urban_pct') + }) +}) + +// ── BHF-EBLUP (unit-level) ──────────────────────────────────────────────────── +describe('generateCode — bhf-eblup (unit-level, standard)', () => { + const code = generateCode(bhfEblup, BHF_INPUTS) + + it('R script contains library(sae)', () => { + expect(code.r).toContain('library(sae)') + }) + + it('Stata script contains the mixed command', () => { + expect(code.stata).toContain('mixed') + }) + + it('R script contains all user-supplied variable names', () => { + assertAllVarsPresent(code.r, ['income', 'region', 'x1', 'x2', 'x3'], 'R') + }) + + it('Stata script contains all user-supplied variable names', () => { + assertAllVarsPresent(code.stata, ['income', 'region', 'x1', 'x2', 'x3'], 'Stata') + }) + + it('R script has no unreplaced {{ tokens', () => assertNoUnreplacedTokens(code.r, 'R')) + it('Stata script has no unreplaced {{ tokens', () => assertNoUnreplacedTokens(code.stata, 'Stata')) + + it('does not use fallback (stataMinVersion 14)', () => { + expect(code.usedFallback).toBe(false) + }) + + it('AUX_VARS_R_VEC produces quoted column-selection form', () => { + // bhf-eblup template uses {{AUX_VARS_R_VEC}} inside c(...) for census column selection + expect(code.r).toContain('"x1", "x2", "x3"') + }) + + it('census data path is substituted', () => { + expect(code.r).toContain('census.csv') + expect(code.stata).toContain('survey.csv') + }) +}) + +// ── EBP / CensusEB — poverty mapping with Stata v14 fallback ───────────────── +describe('generateCode — ebp-censuseb (poverty mapping, stataVersion:14)', () => { + const code = generateCode(ebpCensuseb, EBP_INPUTS) + + it('R script contains library(povmap)', () => { + expect(code.r).toContain('library(povmap)') + }) + + it('usedFallback is true (stataVersion 14 < stataMinVersion 17)', () => { + expect(code.usedFallback).toBe(true) + }) + + it('fallbackNote is a non-empty string', () => { + expect(typeof code.fallbackNote).toBe('string') + expect(code.fallbackNote!.length).toBeGreaterThan(0) + }) + + it('Stata output uses the mixed command (fallback)', () => { + expect(code.stata).toContain('mixed') + }) + + it('Stata output does NOT use the sae model command (full impl not available on v14)', () => { + expect(code.stata).not.toContain('sae model') + }) + + it('R script contains all user-supplied variable names', () => { + assertAllVarsPresent(code.r, ['cons_pc', 'ea_code', 'hhsize', 'edu_head', 'urban'], 'R') + }) + + it('Stata script contains all user-supplied variable names', () => { + assertAllVarsPresent(code.stata, ['cons_pc', 'ea_code', 'hhsize', 'edu_head', 'urban'], 'Stata') + }) + + it('R script has no unreplaced {{ tokens', () => assertNoUnreplacedTokens(code.r, 'R')) + it('Stata script has no unreplaced {{ tokens', () => assertNoUnreplacedTokens(code.stata, 'Stata')) + + it('N_SIM is substituted in R script', () => { + expect(code.r).toContain('200') + }) +}) + +// ── GLMM-binary ─────────────────────────────────────────────────────────────── +describe('generateCode — glmm-binary (binary target)', () => { + const code = generateCode(glmmBinary, GLMM_INPUTS) + + it('R script contains library(lme4)', () => { + expect(code.r).toContain('library(lme4)') + }) + + it('Stata script contains the meglm command', () => { + expect(code.stata).toContain('meglm') + }) + + it('R script contains all user-supplied variable names', () => { + assertAllVarsPresent(code.r, ['employed', 'province', 'age', 'education'], 'R') + }) + + it('Stata script contains all user-supplied variable names', () => { + assertAllVarsPresent(code.stata, ['employed', 'province', 'age', 'education'], 'Stata') + }) + + it('R script has no unreplaced {{ tokens', () => assertNoUnreplacedTokens(code.r, 'R')) + it('Stata script has no unreplaced {{ tokens', () => assertNoUnreplacedTokens(code.stata, 'Stata')) + + it('does not use fallback (stataMinVersion 14)', () => { + expect(code.usedFallback).toBe(false) + }) + + it('R script uses logit random-intercept structure', () => { + expect(code.r).toContain('binomial(link = "logit")') + }) + + it('Stata script uses bernoulli logit family', () => { + expect(code.stata).toContain('family(bernoulli)') + }) +}) + +// ── v14 fallback — general property tests ───────────────────────────────────── +describe('generateCode — Stata v14 fallback logic', () => { + it('usedFallback true when stataVersion < stataMinVersion and fallback exists', () => { + const code = generateCode(ebpCensuseb, { ...EBP_INPUTS, stataVersion: 16 }) + expect(code.usedFallback).toBe(true) + }) + + it('usedFallback false when stataVersion >= stataMinVersion', () => { + const code = generateCode(ebpCensuseb, { ...EBP_INPUTS, stataVersion: 17 }) + expect(code.usedFallback).toBe(false) + expect(code.fallbackNote).toBeUndefined() + }) + + it('Stata v17 output uses the full sae model command (not mixed)', () => { + const code = generateCode(ebpCensuseb, { ...EBP_INPUTS, stataVersion: 17 }) + expect(code.stata).toContain('sae model') + }) + + it('no unreplaced tokens in fallback Stata output', () => { + const code = generateCode(ebpCensuseb, { ...EBP_INPUTS, stataVersion: 14 }) + assertNoUnreplacedTokens(code.stata, 'Stata fallback') + }) +}) + +// ── Token substitution edge cases ────────────────────────────────────────────── +describe('generateCode — token substitution', () => { + it('single auxiliary variable: AUX_VARS_R produces plain name (no +)', () => { + const code = generateCode(fhEblup, { ...FH_INPUTS, auxiliaryVars: ['x1'] }) + expect(code.r).toContain('x1') + expect(code.r).not.toContain('x1 + ') + }) + + it('three auxiliary variables: AUX_VARS_R produces x1 + x2 + x3', () => { + const code = generateCode(fhEblup, { ...FH_INPUTS, auxiliaryVars: ['x1', 'x2', 'x3'] }) + expect(code.r).toContain('x1 + x2 + x3') + }) + + it('AUX_VARS_R_VEC with two vars produces "x1", "x2"', () => { + const code = generateCode(bhfEblup, { ...BHF_INPUTS, auxiliaryVars: ['x1', 'x2'] }) + expect(code.r).toContain('"x1", "x2"') + }) + + it('N_SIM default is 200 when nSimulations not provided', () => { + const code = generateCode(ebpCensuseb, { ...EBP_INPUTS, stataVersion: 17, nSimulations: undefined }) + expect(code.r).toContain('200') + }) + + it('custom nSimulations is substituted', () => { + const code = generateCode(ebpCensuseb, { ...EBP_INPUTS, stataVersion: 17, nSimulations: 500 }) + expect(code.r).toContain('500') + }) + + it('generateCode is pure: two calls with same inputs produce identical output (modulo date)', () => { + const a = generateCode(fhEblup, FH_INPUTS) + const b = generateCode(fhEblup, FH_INPUTS) + // Strip the date line before comparing (dates may differ by a second in slow CI) + const stripDate = (s: string) => s.replace(/\d{4}-\d{2}-\d{2}/, 'DATE') + expect(stripDate(a.r)).toBe(stripDate(b.r)) + expect(stripDate(a.stata)).toBe(stripDate(b.stata)) + }) +}) diff --git a/src/engine/codegen.ts b/src/engine/codegen.ts new file mode 100644 index 0000000..70d3ae6 --- /dev/null +++ b/src/engine/codegen.ts @@ -0,0 +1,121 @@ +import type { CatalogueEntry } from '../types/index.js' + +export interface UserInputs { + targetVar: string + areaIdVar: string + auxiliaryVars: string[] + weightVar?: string + directEstVar?: string + directVarVar?: string + coordinateLat?: string + coordinateLon?: string + contiguityMatrixPath?: string + censusDataPath?: string + surveyDataPath?: string + stataVersion: number + nSimulations?: number + mseMethod?: 'analytic' | 'bootstrap' +} + +export interface GeneratedCode { + r: string + stata: string + usedFallback: boolean + fallbackNote?: string +} + +// ── Token substitution ───────────────────────────────────────────────────────── + +function buildTokenMap(inputs: UserInputs): Record { + const date = new Date().toISOString().slice(0, 10) + const auxR = inputs.auxiliaryVars.join(' + ') + const auxRVec = inputs.auxiliaryVars.map(v => `"${v}"`).join(', ') + const auxStata = inputs.auxiliaryVars.join(' ') + const survey = inputs.surveyDataPath ?? 'survey.csv' + const census = inputs.censusDataPath ?? 'census.csv' + + return { + DATE: date, + TARGET_VAR: inputs.targetVar, + AREA_ID: inputs.areaIdVar, + AUX_VARS_R: auxR, + AUX_VARS_R_VEC: auxRVec, + AUX_VARS_STATA: auxStata, + WEIGHT_VAR: inputs.weightVar ?? '', + DIRECT_EST_VAR: inputs.directEstVar ?? '', + DIRECT_VAR_VAR: inputs.directVarVar ?? '', + SURVEY_DATA: survey, + CENSUS_DATA: census, + AREA_DATA: census, // area-level aggregates default to the same path as census + CONTIG_MATRIX: inputs.contiguityMatrixPath ?? 'contiguity_matrix.csv', + N_SIM: String(inputs.nSimulations ?? 200), + MSE_METHOD: inputs.mseMethod ?? 'bootstrap', + } +} + +function substituteTokens(template: string, tokens: Record): string { + let result = template + for (const [key, value] of Object.entries(tokens)) { + result = result.replace(new RegExp(`\\{\\{${key}\\}\\}`, 'g'), value) + } + return result +} + +// ── Footer ───────────────────────────────────────────────────────────────────── + +const R_FOOTER = ` +# ============================================================ +# MSE / RMSE INTERPRETATION +# ============================================================ +# MSE — Mean Squared Error; measures average squared deviation from the +# true area value. Estimated by the method noted above. +# RMSE — sqrt(MSE); on the same scale as the target variable. +# CV — Coefficient of Variation = RMSE / estimate. +# Areas with CV > 0.25 have high relative uncertainty; +# interpret model-based estimates cautiously for these areas. +# +# NEXT STEPS +# - Inspect the model diagnostics before publishing any estimates. +# - Compare model-based estimates against direct estimates as a sanity check. +# - See docs/adding-a-method.md to extend or adapt this script. +` + +const STATA_FOOTER = ` +* ============================================================ +* MSE / RMSE INTERPRETATION +* ============================================================ +* MSE - Mean Squared Error; measures average squared deviation from the +* true area value. Estimated by the method noted above. +* RMSE - sqrt(MSE); on the same scale as the target variable. +* CV - Coefficient of Variation = RMSE / estimate. +* Areas with CV > 0.25 have high relative uncertainty; +* interpret model-based estimates cautiously for these areas. +* +* NEXT STEPS +* - Inspect model diagnostics before publishing any estimates. +* - Compare model-based estimates against direct estimates as a sanity check. +* - See docs/adding-a-method.md to extend or adapt this script. +` + +// ── Public API ───────────────────────────────────────────────────────────────── + +export function generateCode(entry: CatalogueEntry, inputs: UserInputs): GeneratedCode { + const usedFallback = + inputs.stataVersion < entry.stataMinVersion && entry.stataV14Fallback !== null + + const tokens = buildTokenMap(inputs) + + const r = substituteTokens(entry.rTemplate, tokens) + R_FOOTER + + const stataRaw = usedFallback ? entry.stataV14Fallback! : entry.stataTemplate + const stata = substituteTokens(stataRaw, tokens) + STATA_FOOTER + + const fallbackNote = usedFallback + ? `Stata ${inputs.stataVersion} is below the minimum required version ` + + `(v${entry.stataMinVersion}) for ${entry.displayName}. ` + + `The generated .do file uses the base fallback command instead. ` + + `Use the R script for the full implementation.` + : undefined + + return { r, stata, usedFallback, fallbackNote } +}