Generalized Engression Models

• A GEM uses a neural generator, a data-type-specific link, and learned-scale perturbation, while keeping the generator, loss, and fitting algorithm shared across outcome types.The links include identity for continuous coordinates, sign for binary ones, largest-coordinate indicators for categorical ones, rounding for ordinal ones, and rank vectors for rankings.Source: §1 · pp. 2–3
• Under positivity and continuity conditions, a continuous generator can represent any conditional distribution of an outcome type and be approximated by a neural network.Without those conditions, the representation result still holds for a measurable generator.Source: §1 · p. 3
• Perturbation makes discontinuous discrete links trainable: averaging over it yields a differentiable expected loss whose gradient can be estimated by Monte Carlo.The method uses the energy score as its loss.Source: §1 · p. 3
• 0.015 TV error: GEM matches all six ranking probabilities where Plackett–Luce models with linear or spline scores cannot.Across the figure’s panels, GEM also recovers binary-label co-occurrence and the conditional laws of a mixed continuous-categorical outcome from one model.Source: §1 · p. 3, Figure 1 · p. 4
• On 242-species vegetation data, GEM matches or exceeds purpose-built joint models on joint prediction and is better calibrated for species richness than models fit one species at a time.One GEM fit takes seconds on one GPU, compared with minutes to hours for the other joint models.Source: §1 · pp. 4–5
• GEM improves on stacked-marginal, Gaussian-copula, and tree-based alternatives on the health survey’s joint distribution on every split, while coming close on the ordinal margin.The survey includes seven continuous measurements, five binary indicators, and an ordinal self-rated health score from 27,983 participants.Source: §1 · p. 5

2 Background

• The section situates the proposal by reviewing type-specific statistical regression, generative models, and engression.Source: §2 · p. 5

2.1 Regression models for outcomes of a given type

• Classical regression models use type-specific links on latent predictors; multivariate dependence often comes from Gaussian noise or low-rank factors, while likelihood fitting can require integration or simulation.Flexible marginal laws can also be paired with parametric copulas, while distributional random forests reweight observed outcome vectors.Source: §2.1 · pp. 5–6

2.2 Generative models

• Generative approaches differ in their generator and fitting criterion, including adversarial, variational, likelihood-based, denoising, and scoring-rule methods; discrete and mixed outcomes often use type-specific approaches.Source: §2.2 · p. 6

2.3 Engression and the energy score

• Engression represents the conditional law through generator draws and fits it by minimizing the proper energy score, using stochastic gradient descent with fresh noise draws.The score is strictly proper for distributions with finite moment of order β, and the paper uses β = 1.Source: §2.3 · p. 6
• Engression’s continuous-valued draws are invalid for discrete outcomes, although its energy score remains proper for discrete and mixed distributions; GEM constrains generator outputs through links and derives gradients for discontinuous links.Source: §2.3 · p. 7

3.1 Model and links

• GEM represents any-type or mixed conditional outcomes with a neural generator, a fixed type-specific link, and a learned-scale perturbation; ε generates distributional variation while η makes discontinuous links tractable.The recommended default uses a common scale across coordinates; coordinate-wise scales can suit mixed blocks on different scales. The ordinal link rounds to the nearest level, while rankings use rank vectors.Source: §3.1 · p. 7, Table 1 · p. 8

3.2 Objective and gradient

• When the model class contains the true conditional law, the population energy-score objective is minimized exactly by models that represent that law.Averaging over perturbations makes the loss continuously differentiable even when the link is discontinuous.Source: §3.2 · p. 8

3.3 Algorithm

• For 242 labels, the binary-coordinate control variate reduced gradient-estimator variance by two orders of magnitude; continuous coordinates instead use pathwise derivatives.The coordinate-wise baseline leaves estimator means unchanged, and the variance reduction was reported on the species data.Source: §3.3 · pp. 9–10

3.4 Other kernel scores

• Replacing the energy score with a conditionally negative-definite kernel score fits squared discrepancy between conditional laws; on finite outcome spaces, the indicator kernel gives half the Brier score.The Gaussian kernel is another example; strict propriety on finite outcome spaces holds for the examples described.Source: §3.4 · p. 11

4 Universality theory

• The universality results concern the zero-perturbation limit: perturbation is needed for fitting, not representation, and the paper leaves quantitative guarantees based on derivative norms to future work.For discrete outcomes, the stated convergence to the unperturbed representation is in total variation.Source: §4 · p. 11

4.1 Measurable map g

• For outcomes in standard Borel spaces, a measurable generator with independent Gaussian noise can represent the conditional law through a measurable link with a measurable right inverse.A surjective Borel link is one sufficient condition; even univariate noise suffices, though the resulting generator may be highly discontinuous.Source: §4.1 · p. 12

4.2 Continuous map g for continuous outcomes

• Under a positive continuous joint density and finite continuous marginal densities, a continuous generator represents Euclidean outcomes using sequential conditional quantiles.The construction generates coordinates in turn from their conditional quantiles, extending the univariate quantile map through the inverse Rosenblatt transformation.Source: §4.2 · p. 13

4.3 Continuous map g for discrete outcomes

• For discrete outcomes, a regular positive latent lift through the link yields a continuous generator; the stated condition requires interior link cells and positive continuous class subdensities.The discrete representation theorem applies under Condition 2, which includes these cell and subdensity requirements.Source: §4.3 · pp. 13–14

4.4 Continuous map g for mixed-type outcomes

• For mixed continuous and finite-valued outcomes, a continuous generator exists under Condition 3 by generating the continuous block first and the discrete block conditional on it.This ordering preserves dependence between the two outcome blocks; the conditions include regularity of the continuous density and positive continuous subdensities for each discrete value.Source: §4.4 · p. 15

5 Simulation study

• The simulations compare GEM with type-specific baselines and an oracle across three outcome types, assessing both joint distributions and marginal fit over 50 replications.The reported metric curves use bands of one standard error.Source: §5 · p. 15

5.1 Multi-label outcome

• GEM stays close to the oracle on multi-label joint metrics as dependence ρ grows, while independent baselines increasingly miss co-occurrence and the conditional count distribution.On marginal metrics, GEM is slightly behind the flexible NN-logistic baseline; Figure 3 reports joint energy score, mean per-label CRPS, cross-entropy, and co-occurrence error, while Figure 4 shows the conditional count distribution.Source: §5.1 · p. 16, Figure 3 · p. 16, Figure 4 · p. 16

5.2 Mixed-type outcome

• GEM has a lower joint energy score than the stacked-marginals pipeline at every sample size, while its cross-type dependence errors fall to a fraction of the pipeline’s by n = 1,000.Marginal scores are generally comparable, though GEM’s categorical cross-entropy is higher; GEM’s joint score is closer to the oracle despite comparable or marginally less sharp marginals.Source: §5.2 · pp. 17–18, Figure 5 · p. 17

5.3 Ranking outcome

• In the coupled ranking setting, GEM is ahead of Plackett–Luce on all four metrics from n = 2,000, whereas Plackett–Luce remains far from the truth at every sample size.In the independent setting, Plackett–Luce performs better at n = 500, the methods are close at n = 2,000, and GEM leads at n = 10,000.Source: §5.3 · p. 18, Figure 6 · p. 19

6.1 Joint species distribution modeling

• The vegetation benchmark covers 242 plant and lichen species across 1,200 plots, where residual richness variance is about nine times the independence expectation.GEM uses a sign link and an ensemble of five seeded fits; the benchmark compares it with ten independent-draw and five joint baselines.Source: §6.1 · p. 19
• First on six of eight metrics and second on the other two, the GEM ensemble also leads the corrected overall benchmark score.It is second on richness interval error and AUC.Source: Figure 7 · p. 20, §6.1 · pp. 20–21, Table 2 · p. 21
• Independent-draw models with strong marginal scores remain worse on richness than the best joint models, while GEM is best on both.GEM’s collective-dependence performance coexists with weaker pairwise-association scores than the other joint baselines.Source: §6.1 · pp. 20–21, Figure 7 · p. 20
• GEM fits in seconds per fit on one GPU, while joint baselines take minutes to hours on CPU; its five-fit ensemble takes under two minutes.Each further GEM draw is one network pass, and each gradient step scales linearly with the number of plots.Source: §6.1 · pp. 21–22

6.2 Health examination surveys

• The NHANES application models seven continuous measures, five binary indicators, and one five-level ordinal health outcome given demographic, socioeconomic, and behavioral covariates.Comparators include separately fitted marginals, a Gaussian copula, and distributional random forests that resample observed outcome vectors.Source: §6.2 · p. 22
• Positive on every split against every baseline, GEM’s joint energy-score differences show its strongest health-survey performance is on the joint distribution.Its self-rated-health ranked probability score is 0.4980 versus 0.4954 for the softmax-based baselines; against DRF, HbA1c is the noted exception to the otherwise consistent per-outcome ordering.Source: §6.2 · pp. 22–23, Figure 8 · p. 23

7 Discussion

• GEM generates draws to estimate the joint conditional distribution, capturing both coordinate-level uncertainty and dependence.Source: §7 · p. 23
• The framework is nonparametric in covariates and response distribution, which is not restricted to a particular family.Source: §7 · p. 23
• GEM covers continuous, binary, categorical, ordinal, ranking, and mixed-type outcomes.Source: §7 · p. 23
• The same loss function and algorithm apply across outcome types, unlike GLMs, which require a likelihood specified for each family.Source: §7 · p. 23
• GEM uses a single sample-based loss without an adversary or explicit likelihood, and its one-step sampler supports Monte Carlo estimates of functionals.Source: §7 · p. 23
• In two applications, GEM matches or improves on specialized joint species models and improves health-outcome joint prediction over several baselines, with little or no loss on individual outcomes.The health outcome mixes continuous, binary, and ordinal coordinates.Source: §7 · pp. 23–24
• Type-specific methods can occasionally be more accurate on an individual coordinate, including softmax for categorical outcomes and Plackett–Luce when correctly specified with small samples.Source: §7 · p. 24
• The representation results do not establish estimation rates, meaningful finite-sample results remain future work, and extrapolation for discrete outcomes remains open.Source: §7 · p. 24

A.1 Proof of Proposition 1

• The smoothed loss is continuously differentiable in its location and positive scale parameters when the perturbation density meets the stated conditions and the function has at most linear growth.The proof applies this result to a Gaussian perturbation and the loss built from two independent perturbed draws.Source: §A.1 · pp. 29–30
• Changing variables moves parameter dependence from the potentially nondifferentiable function into the density, so derivatives can be taken without differentiating that function.The density derivatives and dominated convergence justify differentiating the expectation.Source: §A.1 · pp. 29–30
• The proposition’s loss equals an expectation over perturbations, and its derivatives follow by applying the smoothing lemma and the chain rule.The proof uses two independent perturbations with the same conditional distribution.Source: §A.1 · p. 30

A.2 Proof of Proposition 2

• A measurable right inverse of the link outside a probability-zero set lets the proof lift the response and apply the representation result to the lifted outcome.Applying the link afterward recovers the original response almost surely.Source: §A.2 · pp. 30–31

A.3 Proof of Lemma 3

• Positive continuous conditional densities yield continuous conditional CDFs and unique finite quantiles for every probability in (0, 1).The proof also establishes continuity of the quantile functions in their conditioning arguments and probability level.Source: §A.3 · p. 31

A.4 Proof of Theorem 1

• Sequential conditional quantiles driven by independent uniform variables construct a continuous generator whose conditional distributions match those of the outcome.The uniforms are obtained by applying the standard normal CDF to independent Gaussian noise.Source: §A.4 · pp. 31–32

A.5 Proof of Theorem 2

• A continuous generator for a lifted continuous outcome yields a generator for the original outcome after applying the link, without requiring the link itself to be continuous.The proof uses the defining property of the lift to recover the target distribution.Source: §A.5 · p. 32

A.6 Proof of Theorem 3

• For finite discrete outcomes, the proof constructs a regular positive latent lift through the link and applies the lifted-outcome representation result.This gives a continuous generator whose linked output has the required joint distribution with the covariates.Source: §A.6 · pp. 32–33

A.7 Proof of Theorem 4

• For mixed outcomes, the proof first generates the continuous coordinates, then uses them as inputs to a discrete generator, preserving the joint distribution.The combined map is continuous because both component generators are continuous.Source: §A.7 · pp. 33–34

B.1 Implementation of GEM

• Table 3 records GEM settings across five experiments, with d and k denoting covariate and coded-outcome dimensions.The common settings are described in the accompanying text.Source: Table 3 · p. 34
• GEM uses a four-layer StoNet generator with fresh Gaussian noise at each hidden layer and a separate two-layer scale network of width 100.The generator has three hidden ReLU layers followed by a linear output layer; the scale network applies softplus to a ReLU network.Source: §B.1 · p. 35
• The objective is estimated with two draws per observation and optimized by full-batch Adam; outside species data, the validation energy score selects the iterate, but GEM is otherwise untuned.The shared default hyperparameters are used across simulation settings, sample sizes, and NHANES data.Source: §B.1 · p. 35

B.2 Multi-label outcome

• At ρ = 1, the joint energy score is 1.35 versus 1.84 for NN logistic and 1.91 for linear logistic, while co-occurrence error is 0.08 versus 0.39 and 0.40.At n = 1,000 GEM has the lowest joint energy score for ρ ≥ 0.3 and co-occurrence error for ρ ≥ 0.2; at n = 100 those thresholds are ρ ≥ 0.5 and ρ ≥ 0.3.Source: Figure 9 · p. 36, §B.2 · p. 36

B.3 Mixed-type outcome

• The mixed-type experiment measures cross-type dependence through residual correlation ratios; independently drawn stacked marginals have errors near 0.48 for binary labels and 0.43 for the class.The true ratios range from 0.39 to 0.54 for binary labels and equal 0.43 for the categorical class.Source: §B.3 · pp. 36–37

B.4 Ranking outcome

• Spline rather than linear ranking scores change total variation by +0.014, +0.005, and −0.002 at n = 500, 2,000, and 10,000 in the independent setting, and by −0.003 to −0.004 in the coupled setting.Both Plackett–Luce models draw rankings by adding independent standard Gumbel noise to fitted scores.Source: §B.4 · p. 37

C.1 Data and splits

• 1,200 benchmark plots are split into five random halves, each with 600 training and 600 test plots.The benchmark reduces six environmental variables to four principal components carrying 88% of their variance.Source: §C.1 · p. 38
• 290–316 shared plots between splits yield per-plot score-difference correlations of about 0.7–0.9, so results are reported as split wins rather than paired t-statistics.Residual presence at a plot correlates only 0.13 on average with its mean residual at the eight nearest plots.Source: §C.1 · p. 38
• 148 observed richness variance contrasts with 17 expected under independent species with the fitted probabilities.The variance estimates use test plots from the benchmark split and species-wise logistic regressions.Source: §C.1 · p. 38

C.2 Models and settings

• The comparison includes 15 baselines, with HMSC using spatial latent factors run only on the benchmark split and taking 3.9 hours per fit.Each baseline is fitted once per split and generates 1,000 draws per test plot; nine models come from the benchmark.Source: §C.2 · pp. 38–39, Table 4 · p. 39
• GEM uses four standardized environmental components, 4,000 iterations, and a learning rate annealed from 10^-3 to zero; the selected settings are unchanged across splits.Tuning used held-out training plots, selected 4,000 iterations at 10^-3, and each fit produced 1,000 test-plot draws.Source: §C.2 · pp. 39–40
• A single GEM fit beats 13 of 15 baselines on energy score and 14 on richness CRPS on every split, although its marginals trail those of the best single-species models.The ensemble pools five fits.Source: §C.2 · p. 40
• Across five seeds, a single fit’s energy score has within-split standard deviation 0.007, motivating an ensemble of five fits.Each baseline is a single fit at benchmark or package settings, although some tune internally.Source: §C.2 · p. 40
• 242–606-fold variance reduction from the control variate leaves the gradient estimator unbiased within 0.98–1.08 standard errors across six study cells.Without the control variate, the tuned network needs 16 times as many iterations and reaches energy score 2.916 against 2.795 with it.Source: §C.2 · p. 40, Table 5 · p. 40

C.3 Metrics

• Energy score evaluates the full species vector, while CRPS evaluates scalar richness and Sørensen dissimilarity across test plots or plot pairs.The energy score uses all species, coded 0/1; richness is the number of present species, and dissimilarity is computed for 300 sampled pairs.Source: §C.3 · p. 41
• Marginal occurrence evaluation includes log score and plain AUC, while randomized probability-integral-transform interval error measures calibration of richness and pairwise dissimilarity.The log score, AUC, interval errors, and selected overall measures use the same fixed subset of 100 draws for every model.Source: §C.3 · p. 41
• Dependence gain is the mean log-ratio, in nats, between the model’s probability for observed common-species pairs and the product of its own marginals.For models with draws independent across species, dependence gain is zero up to Monte Carlo error; co-occurrence correlation compares predicted and observed residual co-occurrence.Source: §C.3 · p. 42
• The overall score averages nine standardized, clipped, direction-aligned measures, while the benchmark’s draw-based accuracy and rank-correlation measures can reward overly narrow predictions.On richness, GEM scores 12.72 against 9.56 for the random forest with coordinates under the benchmark coding, but leads 11.44 against 11.80 on predictive-mean RMSE on all five splits.Source: §C.3 · p. 42, Table 6 · p. 43
• Figure 10 plots paired per-split differences for nine baselines against the GEM ensemble: each dot is a split, and positive values favor GEM.Its vertical bar is the mean across splits; HMSC with spatial latent factors contributes one split.Source: Figure 10 · p. 43

C.4 Further results

• 0.0104 is GEM’s mean energy-score margin over the random forest with coordinates, with paired-difference standard deviation 0.012 across splits.Across models, energy score varies by about 4%, with split standard deviations of 0.028–0.039; paired differences remove much of the shared split variation.Source: §C.4 · p. 43
• 0.746 versus 0.717: HMSC without associations has higher AUC than GEM on species below 5% training prevalence; GEM ranks eighth on rare and third on common species.On common species, AUC is 0.734 for HMSC without associations and 0.755 for GEM.Source: §C.4 · p. 44, Table 7 · p. 44
• 0.0078 nats per pair is GEM’s dependence gain, below the five joint baselines’ 0.0088–0.0142 range on every split.Predicted versus observed residual co-occurrence correlation is 0.81 for GEM, 0.80 for gllvm, and 0.89–0.92 for HMSC and the chains; latent-factor models recover more pairwise association at the expense of marginals.Source: §C.4 · p. 44
• Adding coordinates improves the random forest’s energy score from 2.856 to 2.792 and AUC from 0.709 to 0.734; a spatial smooth raises GAM AUC by 0.012.A spatial basis retains most of the classifier chain’s dependence gain, 0.0096 against 0.0123.Source: §C.4 · p. 44

D Details for the health examination survey application

• The 20 splits each contain 16,789 training, 5,596 validation, and 5,598 test participants.Age is top-coded at 80; missing education, marital, or smoking status is coded by setting that variable's indicators to zero.Source: D Details for the health examination… · p. 45
• GEM is fitted to a 13-dimensional outcome coding, then its draws are converted to the 17-dimensional coding used to score all methods.Baselines include independent stacked marginal fits, their Gaussian-copula coupling, and DRF, which resamples outcome vectors using forest weights.Source: D Details for the health examination… · p. 45
• All methods are scored using 500 draws per test point in the 17-dimensional coding, with joint energy score, marginal scores, and fit time reported.Fit times were measured on shared hardware, with GEM and stacked marginals on one GPU and DRF on CPU threads.Source: D Details for the health examination… · p. 45, Table 8 · p. 46
• Because splits re-partition the same participants and 36 comparisons invite chance stars, the authors emphasize consistent signs over individual entries.The re-partitioning makes standard errors somewhat understate uncertainty; paired-difference standard errors are smaller because split-to-split variation is shared.Source: D Details for the health examination… · p. 46
• GEM has lower CRPS than every baseline on six continuous outcomes, while DRF is better on HbA1c by 2.5% of the score.Against the copula and stacked marginals, GEM's gains are 0.1–0.4% of the score across all seven outcomes.Source: D Details for the health examination… · p. 46, Table 9 · p. 46
• GEM has lower log-loss than the copula and stacked marginals on all five binary labels, and beats DRF on four; cholesterol medication slightly favors DRF.Medication indicators and diagnosed diabetes show the largest marginal-baseline improvements, at 0.8–1.1% of the score.Source: D Details for the health examination… · pp. 46–47, Table 10 · p. 47