COmpositional Nonlinear-debiased Inference, Fay–Herriot with Extremal conformal Regions
Robust. Compositional. Confident.

⚠️ Preview Release — CONIFER is under active development. APIs may change before the first stable 1.0. This README describes the current v0.3 release; earlier versions and their results are in the git history.

Design-aware small-area estimation of forest structure as distributions, not just totals. CONIFER estimates the diameter distribution — stem density split across DBH classes — for small forest areas where the field sample is too thin for a reliable direct estimate, and attaches an honest, checked statement of uncertainty.

In a nutshell. You have a cruise — a few stands measured well, many with just a plot or two, some not at all. CONIFER returns the stems per acre in each diameter class for every stand, including the thin ones, each with an uncertainty range you can defend to whoever paid for the inventory. It borrows strength from your LiDAR/aerial metrics and from comparable stands, and it tells you per stand how much of the answer came from that stand's own plots versus the model. No statistics background needed: pip install "conifer-sae[app]" then conifer-studio, or press Load the demo cruise and click through. Full‑code and no‑code paths give the same estimator.

▶ Watch the 60-second, zero-code tour — from pip install to a checked stand report, no code.

For a stand i the target is the whole stand table s_i = N_i · p_i: a total density N_i (stems per acre) times a composition p_i — the shares across DBH classes, a point on the simplex. That is three hard problems at once, and no prior small-area method handles them together:

• the response is compositional (shares that sum to one, with a nonlinear link to canopy structure),
• the large-diameter tail is structurally zero (many stands genuinely tally nothing above 15–20″), and
• the stands you care about are data-poor (a handful of plots, sometimes none).

Under the hood CONIFER is a compositional area-level Fay–Herriot estimator with a cross-fitted, one-step-debiased machine-learned mean, a zero-robust hurdle for the empty tail, an analytic design-based sampling covariance, design-aware conformal prediction sets on the simplex, and — new in v0.3 — a support-aware reduce-to-direct deferral gate that hands accuracy back to the design-direct exactly where it becomes reliable. It reduces exactly to classical Fay–Herriot when the mean is linear (verified ratio 1.0002), so the familiar estimator is a special case, not a competitor.

pip install conifer-sae

Each release fixed one honest, measured failure mode. The method never changed shape; it got more correct in the regime where it was weakest.

You have a tree list and, ideally, stand-level remote-sensing metrics. That is all CONIFER needs.

import conifer

inv = conifer.from_treelist(
    trees,                       # one row per tallied tree
    stand_col="STAND", plot_col="PLOT", dbh_col="DBH_IN",
    plot_area=0.2,               # fixed-area cruise; use baf=20 for a prism cruise
    aux=stand_metrics,           # LiDAR / spectral / terrain, one row per stand
    aux_stand_col="STAND",
)
print(inv.describe())            # what CONIFER sees
print(inv.issue_table())         # problems, in language you can act on

est = inv.fit()

# prediction intervals, calibrated from your own plots — no known truth needed
lo, hi = conifer.calibration.calibrated_intervals(est, alpha=0.10)
print(est.interval_report_.summary)

cov = conifer.coverage_check(est)                 # and check that they actually hold up
print(cov.summary)
# "Across 10 independent splits, the 90% interval for a given diameter class contained
#  the held-out field value 91% of the time. That meets the stated level."

conifer.report.to_html(est, "stand_report.html", coverage=cov)
conifer.report.to_excel(est, "results.xlsx")

from_treelist bins DBH, derives effective area, and wires the correct sampling covariance for your plot design — a prism cruise gets a design-based covariance built from the plot replicates, because the multinomial assumption behind the analytic one does not hold when each tree carries a DBH-dependent expansion factor. Stand identifiers are carried through every table, figure and export, so a mis-sorted input can no longer silently corrupt the estimate.

Metric units, FIA 1″/2″ class breaks, stand polygons in and out (read_stands, attach_estimates), and a min_dbh merchantability threshold are all supported.

pip install "conifer-sae[app]"
conifer-studio

If your shell reports conifer-studio: command not found (common on Windows, where pip's scripts folder is often not on PATH), this does the same thing and never depends on PATH:

python -m conifer.studio

Upload a tree list — or press Load the demo cruise — and get stand tables, maps, an Excel workbook and a printable stand report, with a per-stand "% from your own plots" readout so you can see how far each estimate was borrowed. Everything runs on the machine you start it on, so proprietary inventory never leaves your network. See apps/forester/.

CONIFER was not tuned on one dataset. The full pipeline was developed on the Idaho origin region and then independently replicated, refit from scratch, on three Southeastern ownerships — a deliberate test of whether it generalizes rather than overfits where it was born. A design-based Monte-Carlo simulation with a known truth backs the empirical ordering and hardened the estimator.

Beyond the "% better than direct" summary, each Southern region was run as a fair head-to-head against published area-level SAE competitors — the multivariate Fay–Herriot of msae, the compositional sae.prop / Esteban family, and a tri-compositional FH (TriCompFH) — on identical inputs (same subsample, covariates, truth), scored as log-density RMSE against the full-cruise held-out truth. Two regimes tell different stories. The origin region, Idaho, is deliberately excluded from this head-to-head — scoring competitors on the region the pipeline was developed on would not be a fair test — so its result is reported separately above (+18% versus direct, with the coverage study), and the three Southern plantations are the clean out-of-region comparison.

In the data-poor regime (2–4 plots per stand) CONIFER beats the design-direct by a wide margin — the small-area use case:

In the data-rich regime the design-direct and the mature FH become reliable, and v0.2 over-shrank against them. The v0.3 support-aware gate closes that gap to the frontier without touching the data-poor win:

(log-density RMSE vs held-out truth; lower is better.) The gate is a single global rule, w_ik = [k_i/(k_i+8)] · [n_ik/(n_ik+1)] — a data-adequacy factor (defer more as plot support k grows) times an add-one support prior (defer a class only where the direct has real tally). South Carolina is the discriminating case: its structural-zero tail makes the raw direct unreliable even when rich (0.794), so the support prior keeps the hurdle there, where a plot-count-only gate would over-defer. The gate is validated no-harm on both the held-out simulation (it improves at every plot count) and all three real regions, and it needs plot-level input — where plots are absent it reduces exactly to v0.2.

St. Joe, Idaho — the origin study. ~50,000 trees on 375 industrial stands, wall-to-wall 3D-NAIP canopy metrics, benchmarked against design-based FIA. In the sparse regime (1–2 plots per stand) CONIFER cuts log-density RMSE from the direct estimator's 1.517 to 1.306 (+14%, gated v0.3.8 identity), and for never-sampled stands to 1.222 — 40% below the global-mean stand table (2.022). On merchantable stems (≥ 5″) the quadratic-mean-diameter EBLUP beats direct by +20%. It leads a broad competitor slate — kNN, Weibull, a random forest, BART-FH, a Dirichlet-multinomial, a spatial Fay–Herriot, KBAABB, and a multivariate FH on raw densities — and against that raw multivariate FH the compositional treatment is decisive (1.229 vs 1.475, p < 0.001): the simplex geometry, not just "more covariates," is the source of the gain. The estimate stays coherent with the design-based FIA totals it rolls up to (461 vs 524 stems/acre, ≥ 2″ composition Aitchison distance 0.36). An external transfer to Moscow Mountain — a different ownership 50 miles away — is reported honestly as moderate (Hellinger 0.30), the limit a canopy-surface sensor should have.

The Southeastern regions — a different forest, refit from scratch. Arkansas (Bradley–Drew, 96% loblolly), Mississippi (Meridian, mature + young-clearcut loblolly) and South Carolina (Greenwood, young loblolly with a near-empty large-diameter tail) are FIA-benchmarked loblolly-pine plantations — ecologically the opposite of Inland-NW conifer. Refit on each, CONIFER wins or ties the real-competitor shootout, delivers valid conformal coverage everywhere (0.90–0.94), and in South Carolina simply kept running on the degenerate tail that made BART-FH, hierarchical BART-FH and a Dirichlet-multinomial crash — "the crashes are the finding." Mississippi is included precisely because it is the case where the direct estimate is competitive: a small, plot-rich, homogeneous population where small-area borrowing has little to add, and CONIFER does not manufacture an advantage the sample size doesn't support — and where v0.3 now closes the data-rich gap.

What is intrinsic vs what you recalibrate. Across all three Southern regions the single load-bearing capability is the zero-robust compositional hurdle — removing it inflates error by +117% (Mississippi) / +160% (Arkansas) / +308% (South Carolina), by far the largest effect of any component, and it is species-agnostic: exactly what a degenerate large-diameter tail needs, whether the tail is Idaho conifer or Carolina pine. What must be calibrated per region, never transferred, is the plot-density convention (audited every time). The earlier per-region adequacy-gate threshold (τ* = 1 / 3 / 5) is superseded in v0.3 by the single global support-aware gate above, which needs no per-region tuning. Bolting on spatial or regeneration structure did not help (neutral to harmful); CONIFER borrows strength through covariates and the gate, not through region-specific machinery. That is the honest reason the method travels: the transferable core is geometry, not an Idaho-shaped default model.

The simulation (v2, registered design). Thirty-four populations — a 3×5×2 factorial of link form (linear / hinge-saturation / sine-interaction) × stand structure × latent signal, two plasmode populations carrying the real NAIP covariate matrices, and two chronosequence mechanism cells — sampled under an operational mixed cruise (25% reference stands at 25 plots, targets at 1–25 plots, 10% never sampled), 200 replicates, every competitor the published implementation (dbarts BART, sae::mseSFH spatial FH with a genuine geographic W, DirichletReg, sae.prop, kNN, Weibull, global mean). Headlines: 52% lower log-density RMSE than the design-direct at one plot per stand (no harm when plots are dense); the flattest risk profile in the slate (worst cell 1.70× the per-cell best, against 2.6–4.7× for every non-compositional competitor); best-in-slate on the real-covariate plasmodes at sparse k; and conformal joint coverage 0.910 / 0.960 (90/95% nominal) essentially constant across all 34 populations and every sampling intensity — worst cell in the study 0.902 — where the analytic Gaussian joint interval collapses to 0.08–0.39. Reference pooling wins aggregate RMSE only in weak-signal homogeneous cells and scores exactly zero on stand discrimination (it predicts one composition for every stand); under real covariates or real structure it degrades 1.3–1.4. The staged draws, per-stand predictions from every method, and scores are released so any estimator can be benchmarked on the identical Monte-Carlo samples (Zenodo DOI in the paper). The earlier eleven-population v1 study, on which the defer gate was originally tuned, is superseded by this design.

Two things here are worth stating plainly, because both are easy to get wrong.

Calibration without a known truth. Conformal prediction needs a calibration set where the truth is known. A real inventory never has one. Calibrating against the design-direct estimate of the same plots looks reasonable and is invalid: CONIFER's estimate is a shrinkage of that estimate, so the residual is mechanically smaller than the true error — measured against a known truth it under-covers (0.68 for a nominal 0.90 set; and a naive Gaussian box under-covers jointly at 0.475 against nominal 0.95). conformalize_holdout() instead splits each stand's plots in two, fits on one half and calibrates against the other. Independent by construction, and it errs wide rather than narrow.

Per-class by default, joint on request. For any one diameter class you name, the reported interval contains the truth at the stated rate — the question a forester actually asks, and several times narrower than a set guaranteed to contain all classes at once. Pass joint=True for the simultaneous claim; every table and figure states which of the two it is reporting.

The set is checked against a known truth, and it holds. Measured on the St. Joe study, the design-aware joint $L_\infty$ / max-score conformal set (the featured, robustly valid band) achieves joint five-class coverage 0.958 at a nominal 0.95 — versus 0.475 for a Gaussian box on the same data — and an independent audit reproduced it at 0.944. (A minimum-volume ellipsoid is available as a data-rich variant; with an in-sample shape it under-covers, so the $L_\infty$ band is the reported set.) It is 20% tighter than that Gaussian box at matched validity, and the zero-robust log-ratio shrinks it by a further ~80% at unchanged coverage. Coverage holds conditionally too (SD 0.032 across canopy strata), not just on average. Across the three Southeastern regions the joint conformal coverage is 0.897–0.935 (target 0.90), while the analytic Gaussian box on a common joint basis collapses to 0.33–0.60 — so the conformal set is necessary, not merely preferable.

Why conformal, not the analytic MSE. The analytic (plug-in or double-bootstrap) MSE is honestly anti-conservative — it recovers only about half of the empirical error, most severely for regeneration. So CONIFER reports the conformal set as its operational uncertainty statement, and coverage_check() measures realised coverage on held-out plots and hands it back in a sentence you can put in front of a client. On the classes anyone acts on, the interval runs about ±90% of the estimate; read the near-empty tail classes in stems per acre, not percent, where the denominator is near zero and a percentage reads alarmingly and means little.

CONIFER carries the University of Idaho / Intermountain Forestry Cooperative provenance, and it is fair to ask whether it is therefore locked to Pacific/Inland-Northwest conifers. It is not — but the distinction matters:

• It is a method you refit on your own cruise, not a frozen pretrained model. Every region above was fit from scratch on its own field data. There are no Idaho-baked coefficients shipped as a default.
• The validated envelope is Inland-Northwest mixed conifer and Southeastern loblolly-pine plantations — two very different systems. The capability that carries the method (the zero-robust compositional geometry) is species-agnostic.
• What you recalibrate locally is the plot-density convention (audited every time) plus the fit itself against your plots. The v0.3 deferral gate is global — no per-region threshold to set.
• What is not yet validated: hardwood and mixed-hardwood / bottomland systems. All Southeastern validation is pine plantation; no hardwood field data has been run. Treat hardwood as unproven — not claimed to work, not claimed to fail.

In short: the honest warning is "audit your density convention per region," not "this only works on Douglas-fir."

• Debiased-ML mean. The synthetic mean is a cross-fitted ensemble (random-feature ridge + a gradient-boosted residual correction) with a one-step (Riesz) debiasing applied to the out-of-fold mean — the discrete analogue of Neyman orthogonality. It shrinks toward the out-of-fold mean, which stops a flexible learner from quietly leaking the field estimate it is meant to improve on (the debiasing restores first-order coverage to near its nominal level). Honest caveat: the cross-fitted mean buys a first-order orthogonal MSE, not a genuine second-order expansion.
• Total × composition. The total N_i (a univariate log-scale FH) and the shares p_i (a compositional FH in additive-log-ratio space) are estimated as two coupled models and recombined by the delta method, with a structural-zero hurdle for class presence and additive benchmarking to the design total.
• Design-aware conformal. Split / Mondrian (group-conditional) conformal on isometric-log-ratio residuals, scored by the design-and-model covariance; basis-equivariant on the simplex. The featured joint set is the $L_\infty$ / max-score band (robustly valid and tightest-valid; the engine default), with a minimum-volume ellipsoid (Braun et al. 2025) offered as a data-rich variant that under-covers on an in-sample shape. Informative sampling is handled by a weighted conformal product weight (covariate-shift ratio × design weight).
• Support-aware deferral (v0.3). When plot support is supplied, the deployed estimate is blended toward the design-direct density by w_ik = [k_i/(k_i+c)] · [n_ik/(n_ik+a)] (c = 8, a = 1), so it reduces to the direct as plots accumulate in classes the direct supports, while the hurdle keeps structural-zero classes. See docs/deferral-v0.3.md.
• Auxiliary-information adequacy gates (v0.3). fit_gated makes covariate use no-harm: a learner-adequacy gate picks the mean (linear ridge vs boosted-tree) by an out-of-fold contest, and a covariate-adequacy gate caps the covariate mean by its out-of-fold trust rho (defers fully below a floor). Never worse than the covariate-free estimate; gains kept where the auxiliary signal is real.

The v0.1 matrix interface is unchanged, for when the data is already aligned:

import conifer
est = conifer.DiameterDistribution(seed=0).fit(counts, area_eff, X)
est.conformalize(s_truth_cal, cal_idx, joint=True, alpha=0.10)
lo, hi = est.predict_interval(joint=True)
est.benchmark(totals, var_totals=var_totals)      # coherence with FIA class totals

# opt in to the v0.3 data-rich deferral by passing per-stand plot support:
est = conifer.DiameterDistribution(seed=0).fit(
    counts, area_eff, X, total_logN=logN, var_logN=vlogN,
    plots=plots,          # list of per-stand plot arrays (supplies k and n)
    direct_dens=D,        # the design-direct class density to defer toward
)

# v0.3.4+: the audited gated identity — reference-subset contest, reliability-corrected
# ramp trust, deployed learner, and the v0.3.8 unsampled-stand rule
gated = conifer.fit_gated(
    counts, area_eff, X, total_logN=logN, var_logN=vlogN,
    adequacy_target=ref_direct,        # reference stands' direct estimate (the trust target)
    adequacy_idx=ref_rows,             # which stands carry the reference cruise
    adequacy_reliability=lam,          # split-half reliability of that target (disattenuates rho)
    debias=True,
    fit_kwargs=dict(plots=plots, direct_dens=D),   # keeps the defer gate live in the gated fit
)
gated.s_hat_    # gated.rho / gated.rho_eff / gated.learner / gated.deployed_mode_ explain the decision
                # stands with zero counts receive the deployed covariate mean directly (v0.3.8)

Run the minimal worked example:

python examples/quickstart.py

An executed walkthrough of the path above — the tree list you have, from_treelist and its data checks, the fit, how much of each estimate came from that stand's own plots, calibrated intervals with a coverage check, the tables, your cruise beside CONIFER, and the plain-language narrative. Every number and figure is produced by the code above it, on a synthetic but silviculturally realistic stocked cruise.

• Notebook (renders on GitHub): docs/vignettes/conifer-getting-started.ipynb
• Rendered HTML: open the executed vignette

pip install conifer-sae
jupyter notebook docs/vignettes/conifer-getting-started.ipynb

CONIFER is a package; you build estimators on a shared engine.

A new region is a run, not a name — DiameterDistribution(spatial=True, regen_aware=True) toggles capabilities; you don't fork the package per state.

Small-area estimation and conformal prediction are both mature — but not together, not on the simplex, and not for forestry. CONIFER is, to our knowledge, the first estimator to unify in one engine: (i) a compositional area-level Fay–Herriot model for the full DBH-class stem-density vector, (ii) a cross-fitted, one-step-debiased ML mean that provably nests the linear FH mean, (iii) a zero-robust hurdle for the empty tail, (iv) design-aware distribution-free conformal sets on the simplex (a joint $L_\infty$ / max-score band, with a minimum-volume-ellipsoid variant), and (v) a support-aware reduce-to-direct deferral gate that hands back to the design-direct exactly where it is reliable — validated across four regions and a known-truth simulation.

• emdi / sae (R) give area-level FH with parametric MSE, but no compositional target and no distribution-free prediction sets.
• MAPIE / crepes (Python) give conformal prediction, but no small-area borrowing and no simplex geometry.
• rFIA / FIESTA give design-based forest estimates, but do not model or borrow strength — CONIFER consumes these as its benchmark rather than competing with them.
• The nearest prior art — Esteban et al. (2020) compositional FH, Georgakis et al. (2025) multivariate FH for volume/basal-area/height, Amaral et al. (2025) Dirichlet-HDR conformal, White et al. (2025) KBAABB — each supplies one ingredient; none targets the diameter distribution with valid, finite-sample, set-valued uncertainty.

On a variable-radius (prism/BAF) cruise, expect CONIFER's estimate in the smallest diameter class to sit well below the field-only estimate. That is the method working, not a miss: a prism selects trees with probability proportional to basal area, so a single tallied small tree carries a very large trees-per-acre expansion, and only a minority of stands catch one at all. The field-only small-tree estimate is therefore high-variance and inflated in exactly those stands, and CONIFER shrinks it toward what comparable stands and the covariates support. The % from this stand's own plots column tells you how far that shrinkage went.

Point it at three aligned CSV matrices:

conifer fit --counts counts.csv --area area.csv --aux aux.csv --out s_hat.csv

This project states its limits rather than burying them.

• CONIFER defers where it should. The gain is in genuinely thin samples; as plots accumulate the v0.3 gate hands accuracy back to the design-direct — which is what small-area estimation is for.
• v0.3 reaches the frontier, not past it. In the data-rich regime CONIFER matches the direct and mature FH (Arkansas 0.363 vs TriCompFH's 0.326), the correct SAE target when the direct is reliable — no claim of dominance there.
• Surface sensors miss the understory. 3D-NAIP canopy metrics don't see sub-canopy regeneration, driving the 0–2″ errors and moderate cross-ownership transfer; canopy-penetrating LiDAR or local plots are the remedy.
• Geography adds little here. A spatial Fay–Herriot ties within Monte-Carlo noise and residuals are spatially white — a spatial term isn't warranted on these data.
• A set-efficiency gap, largely closed. A Dirichlet-HDR set (Amaral et al. 2025) is tighter in principle; the zero-robust log-ratio closed ~80% of the gap — a modeling choice, not a defect.
• Where it loses: the direct estimate wins the 1-Wasserstein shape metric (CONIFER optimizes simplex geometry, not diameter-axis transport), and hardwood systems are unvalidated.
• The bundled demo is synthetic. Sparse shows the thin-sample gain, rich shows convergence to direct — but the numbers are simulated; the accuracy claims rest on the four real studies.

A methodology manuscript describing the estimator is in preparation; this README will be updated with the citation once it is available. Until then, cite the software:

@software{poolakkal_conifer,
  author  = {Poolakkal, Jaslam},
  title   = {{CONIFER}: COmpositional Nonlinear-debiased Inference, Fay–Herriot with
             Extremal conformal Regions},
  url     = {https://github.com/IFC-UIDAHO/conifer},
  note    = {Methodology manuscript in preparation}
}

CITATION.cff in this repository carries the machine-readable version, which GitHub renders as a Cite this repository button.

Supported by the NCASI Foundation through the Partnership for Small Area Estimation, under Robust small-area estimation strategies for developing accurate stand-level diameter distributions (PI: Jaslam Poolakkal, University of Idaho), funded by the USDA Forest Service, Rocky Mountain Research Station.

MIT. See LICENSE.