Factorial Regression: Models & Methods

• Factorial regression is a framework that uses structured experimental contrasts and latent low-rank decompositions to capture main effects and interactions.
• It leverages classical design-based methods and modern machine learning approaches, including sparsity, reduced-rank models, and neural network architectures.
• The methodology balances bias and variance through model specification choices such as saturated versus unsaturated designs, aiding both causal inference and prediction.

Factorial regression denotes a family of regression frameworks in which structure is imposed through either explicit factorial contrasts among multiple treatment factors or latent factor decompositions of high-dimensional predictors. In the design-based literature, the central object is a $2^K$ or more general multi-factor experiment whose main effects and interactions are encoded as orthogonal contrasts and estimated by regression; in high-dimensional statistics and machine learning, related usage refers to regression models that learn low-rank latent factors, factorized interaction coefficients, or factor-extraction layers as part of the predictor itself (Zhao et al., 2021, Guo et al., 16 Feb 2025, Kharratzadeh et al., 2015). The term therefore spans classical ANOVA-style models, causal estimands defined under randomization, sparse and reduced-rank multivariate models, Bayesian arbitrary-order interaction models, and neural architectures that interleave factor extraction with nonlinear transformations.

1. Formal definitions and representational foundations

In a canonical $2^K$ factorial design, there are $K$ binary factors, each coded either as $\{-1,+1\}$ or $\{0,1\}$, and $Q=2^K$ treatment combinations $z=(z_1,\dots,z_K)$. Under the Neyman–Rubin framework, each unit $i$ has potential outcomes $Y_i(z)$ for all $z\in\mathcal T=\{-1,1\}^K$. Factorial effects are defined by contrast vectors $2^K$0 indexed by subsets $2^K$1, with $2^K$2, $2^K$3, and $2^K$4 for $2^K$5. Stacking these vectors yields an orthogonal contrast matrix $2^K$6 satisfying $2^K$7, and the finite-population factorial effects are

$2^K$8

where $2^K$9 is the vector of finite-population cell means (Shi et al., 2023).

The regression representation of the same structure uses all main effects and interactions as columns of a design matrix. With $K$0 coding, the model matrix $K$1 has columns corresponding to the intercept, main-effect contrasts, and all higher-order interactions obtained by elementwise products. This orthogonality is the algebraic basis for the standard decomposition of factorial effects into regression coefficients and for the coincidence of regression and ANOVA decompositions in balanced designs (Lu, 2016).

A distinct but related usage arises in high-dimensional regression with latent low-dimensional structure. In that setting one observes $K$2 and $K$3, posits $K$4 latent factors $K$5, and seeks a predictor $K$6 minimizing empirical squared loss while exploiting low-rank structure and nonlinearity. The Generalized Factor Neural Network (GFNN) formulation makes this explicit by alternating PCA or Soft PCA layers with feed-forward nonlinear layers (Guo et al., 16 Feb 2025). Sparse Multivariate Factor Regression (SMFR) states the contrast with classical factorial regression directly: classical factorial regression assumes known categorical factors encoded by dummies and interactions, whereas SMFR learns continuous latent factors from the data and uses regularization to induce sparsity in both factor construction and factor effects (Kharratzadeh et al., 2015).

2. Regression estimators in $K$7 factorial experiments

With $K$8 two-level factors coded $K$9, the observed outcomes can be written as

$\{-1,+1\}$0

In matrix form, $\{-1,+1\}$1, where each row of $\{-1,+1\}$2 contains the intercept, the factor codings, and all interaction monomials up to order $\{-1,+1\}$3 (Lu, 2016).

Under $\{-1,+1\}$4 coding, each coefficient is a contrast of the cell means. If $\{-1,+1\}$5 is the mean observed outcome in treatment cell $\{-1,+1\}$6, then the OLS estimator for any factorial coefficient satisfies

$\{-1,+1\}$7

Thus the regression coefficients are not merely convenient summary parameters; they are explicit linear contrasts of the observed cell means (Lu, 2016).

From the design-based perspective, the same coefficients estimate finite-population causal quantities. Zhao and Ding give a unified definition of unit-level and population factorial effects for arbitrary subsets $\{-1,+1\}$8: $\{-1,+1\}$9 With classical $\{0,1\}$0 coding, the saturated factor-based regression recovers these effects directly, and a location-shift strategy $\{0,1\}$1 extends the interpretation to general product weighting schemes $\{0,1\}$2 (Zhao et al., 2021).

A central result is the equivalence between randomization-based and regression-based inference in $\{0,1\}$3 designs. Lu shows that the OLS estimator of the factorial-effect vector is exactly the randomization-based estimator and that the amended Huber–White robust covariance coincides algebraically with the Neyman-style randomization variance estimator. In balanced designs, even the usual homoskedastic OLS variance yields correct element-wise variances for the $\{0,1\}$4 coefficients (Lu, 2016). This resolves a long-standing methodological concern: regression-based factorial analysis is not only algebraically convenient but also justified from a finite-population randomization perspective.

3. Model specification, covariate adjustment, and post-selection inference

The principal specification question in factorial regression is whether to fit a saturated model containing all $\{0,1\}$5 nontrivial factorial terms or an unsaturated model that omits some higher-order interactions. In centered regressors $\{0,1\}$6, the saturated model includes all monomials $\{0,1\}$7 for nonempty $\{0,1\}$8, whereas unsaturated models retain only a subset, such as main effects only or main effects plus two-way interactions (Zhao et al., 2021).

The trade-off is explicit. If the omitted higher-order interactions are nonzero, the unsaturated estimator is biased for the target factorial effects, and the bias does not vanish asymptotically. If those nuisance effects are truly absent, the unsaturated model can reduce variance. Zhao and Ding formalize this through a design-based Gauss–Markov theorem: when omitted interactions do not exist, the parsimonious specification gains efficiency; when they do exist, saturation preserves unbiasedness for the full set of effects (Zhao et al., 2021).

Covariate adjustment introduces an additional layer of specification. In factorial experiments with baseline covariates $\{0,1\}$9, Zhao and Ding distinguish additive adjustment,

$Q=2^K$0

from fully interacted adjustment,

$Q=2^K$1

The fully interacted estimator is asymptotically most efficient, and both the additive and unadjusted estimators are asymptotically dominated except in special cases such as equal-correlation or zero-correlation regimes. Yet the fully interacted fit may be unstable in moderate samples when $Q=2^K$2 interaction terms must be estimated (Zhao et al., 2021).

Restricted Least Squares (RLS) addresses this finite-sample instability by imposing linear restrictions $Q=2^K$3 on the fully interacted model. Additive and unadjusted fits are special cases corresponding to constraints on the treatment–covariate interaction coefficients. The design-based theory gives three guarantees: asymptotic efficiency when the restriction is correctly specified, consistency under mild misspecification provided the restriction on treatment coefficients is correctly specified and separate from that on treatment–covariate interactions, and possible finite-sample gains even when the restriction is only moderately misspecified (Zhao et al., 2021).

Model selection in large factorial designs is treated by Shi et al. through a forward selection procedure governed by sparsity, hierarchy, and heredity. Starting from $Q=2^K$4, the algorithm expands to order $Q=2^K$5, prunes by a weak or strong heredity operator $Q=2^K$6, fits WLS on the corresponding columns of $Q=2^K$7, and applies Bonferroni-adjusted marginal $Q=2^K$8-tests to the newly introduced effects. Under nearly uniform designs, bounded outcomes and covariances, heredity, and signal conditions such as $Q=2^K$9, the selected model is consistent. When consistency fails at higher interaction orders, the theory distinguishes under-selection, which can retain validity if the target contrast is orthogonal to omitted interactions, from over-selection, which preserves consistency and asymptotic normality at the cost of increased variance (Shi et al., 2023).

4. Latent-factor and reduced-rank formulations

A separate line of work studies regression when the predictors themselves possess low-rank structure. Kneip and Sarda decompose each predictor vector as

$z=(z_1,\dots,z_K)$0

where $z=(z_1,\dots,z_K)$1 is a common factor part and $z=(z_1,\dots,z_K)$2 is an idiosyncratic part with mutually uncorrelated components. The response is modeled by an augmented regression

$z=(z_1,\dots,z_K)$3

where $z=(z_1,\dots,z_K)$4 are normalized factor scores derived from the leading eigenvectors of the common covariance. Estimation proceeds by computing empirical principal components, projecting $z=(z_1,\dots,z_K)$5 orthogonally to those components, and fitting a Lasso on the combined design $z=(z_1,\dots,z_K)$6. The paper derives finite-sample inequalities for eigenvalue and eigenvector estimation, establishes a restricted-eigenvalue condition for the projected design, and obtains oracle-style bounds for the coefficient and prediction errors (Kneip et al., 2012).

SMFR extends latent-factor regression to multivariate responses $z=(z_1,\dots,z_K)$7 by factorizing the coefficient matrix as $z=(z_1,\dots,z_K)$8 with

$z=(z_1,\dots,z_K)$9

where $i$0 maps predictors to latent factors and $i$1 maps factors to responses. The optimization penalizes $i$2 by an elastic net and $i$3 by an $i$4 norm: $i$5 Because the objective is biconvex but nonconvex overall, SMFR uses alternating minimization, with a prox-linear variant reported as best in practice. The method also estimates the number of factors by starting from an upper bound $i$6 and decrementing until the fitted $i$7 and $i$8 both have full rank $i$9. Empirically, SMFR attained 10–40% lower test-MSE than Lasso, group-Lasso, sparse reduced-rank approaches, RemMap, sparse-PLS, trace-norm, and Ridge, and it recovered the number of latent factors more reliably than cross-validated SRRR/SPLS (Kharratzadeh et al., 2015).

GOFAR generalizes factor regression to mixed response types within a vector generalized linear model. With continuous, binary, or count outcomes $Y_i(z)$0, it models the integrated natural parameter matrix as

$Y_i(z)$1

where $Y_i(z)$2 and $Y_i(z)$3 are sparse left and right singular vectors and $Y_i(z)$4 contains nonnegative singular values. Co-sparsity is enforced by an elastic-net-style penalty on each rank-1 component $Y_i(z)$5, and estimation is reduced to sequential or parallel unit-rank subproblems solved by alternating majorization–minimization. Theoretical results include monotone descent of the objective and local consistency under sub-exponential errors, positive-definite covariance limits, distinct singular values, bounded second derivatives, and $Y_i(z)$6-consistent initialization for the parallel algorithm (Mishra et al., 2020). This suggests that latent-factor regression can be interpreted not only as dimension reduction but also as structured estimation of the natural-parameter surface in mixed-outcome models.

5. Arbitrary-order interaction models and factorized coefficients

In settings where factorial structure is understood as interactions among predictors rather than treatment assignments, the central difficulty is combinatorial. The fully general interaction regression takes the form

$Y_i(z)$7

and the number of possible interactions is $Y_i(z)$8. Ordinary penalties such as LASSO become computationally infeasible when arbitrary higher-order interactions are allowed (Yurochkin et al., 2017).

MiFM addresses this by representing only $Y_i(z)$9 active interactions $z\in\mathcal T=\{-1,1\}^K$0 and factorizing the corresponding coefficients. The regression function is

$z\in\mathcal T=\{-1,1\}^K$1

The same matrix $z\in\mathcal T=\{-1,1\}^K$2 therefore encodes interaction coefficients across orders. Interaction structure is controlled by a prior on a binary incidence matrix $z\in\mathcal T=\{-1,1\}^K$3, using either a finite-feature Beta–Bernoulli prior or its generalized version $z\in\mathcal T=\{-1,1\}^K$4, which modulates a “rich-gets-richer” versus “poor-gets-richer” tendency in interaction depth (Yurochkin et al., 2017).

Posterior inference is performed by Gibbs sampling. Conjugate updates handle $z\in\mathcal T=\{-1,1\}^K$5, while each $z\in\mathcal T=\{-1,1\}^K$6 and $z\in\mathcal T=\{-1,1\}^K$7 has a Gaussian full conditional because the regression is affine in any single coordinate once the others are held fixed. The incidence matrix $z\in\mathcal T=\{-1,1\}^K$8 is updated from a closed-form prior full conditional combined with a Gaussian likelihood ratio (Yurochkin et al., 2017).

The paper’s main theoretical result is posterior consistency of the regression function. Informally, if the truth is generated by a finite set of interactions $z\in\mathcal T=\{-1,1\}^K$9 with nonzero coefficients and $2^K$00 is chosen at least as large as a threshold $2^K$01, then the MiFM posterior concentrates around the true regression function. The argument relies on a PARAFAC approximation lemma showing that finite collections of interaction coefficients can be represented exactly once $2^K$02 is large enough (Yurochkin et al., 2017).

Empirically, MiFM handled interaction orders from 2 up to 8 or higher, maintained low RMSE as the true interaction depth increased, achieved exact recovery of the presence or absence of each true interaction up to approximately 80% in challenging high-depth scenarios, recovered a hidden 5-way interaction in at least 95% of posterior samples in a genetics example, and identified 62–63 meaningful 3-way and 4-way interactions in retail demand forecasting while remaining competitive with a production factorization machine on AMAPE (Yurochkin et al., 2017). A plausible implication is that “factorial regression” in this usage refers not to experimental factors but to factorized parameterization of combinatorially many interaction effects.

6. Neural and operator-level extensions

GFNN integrates non-parametric regression, factor models, and neural networks for high-dimensional regression. Its defining architectural move is to alternate factor-extraction layers with nonlinear blocks. A PCA layer computes

$2^K$03

while feed-forward layers apply

$2^K$04

Soft PCA replaces explicit eigendecomposition with a differentiable loss contribution, either through a nuclear-norm reconstruction objective

$2^K$05

or through a surrogate combining variance maximization and orthogonality penalties. Because these operations are differentiable, the Soft PCA layer is trained end-to-end by backpropagation (Guo et al., 16 Feb 2025).

The overall training criterion combines prediction loss, factor regularization, and weight decay: $2^K$06 Architecturally, the model may interleave PCA, ReLU blocks, further PCA layers, and additive subnetworks. In GFADNN variants, the factor output $2^K$07 is partitioned into $2^K$08 blocks, each passed through a separate subnetwork and summed, enforcing $2^K$09. The stated motivation is to peel off linear structure, model residual nonlinearity, and then re-factorize in later stages; this is described as especially effective for hierarchical compositional data (Guo et al., 16 Feb 2025).

The empirical summary is unusually explicit. In simulations under both linear factor designs $2^K$10 and nonlinear designs $2^K$11, GFADNN with Soft PCA $2^K$12 NN $2^K$13 Soft PCA $2^K$14 Additive cut test MSE by 30–60% relative to FAR-NN and by an order of magnitude over vanilla NN or PCR. In equity-ETF forecasting with moving-average returns of 11 sector ETFs ($2^K$15) predicting next-day S&P 500 returns, GFADNN achieved Sharpe ratios above 1.2, maximum drawdown under 10%, and directional accuracy around 54%, outperforming Lasso, PCR, vanilla NN, and FAR-NN. In macro nowcasting on FRED-MD across 127 monthly series ($2^K$16), GFANN beat Lasso in 77 cases and attained the best average rank, $2^K$17 of 5, whereas FAST-NN beat Lasso in 30 cases (Guo et al., 16 Feb 2025).

A different extension uses factorial regression as an analysis tool for machine-learning architectures rather than as a predictor form. Jiao et al. decompose each molecular MPNN layer into three operator families—message-seed initialization, node-edge fusion, and node update—with 4, 7, and 3 levels respectively, yielding $2^K$18 architectures under an orthogonal factorial design. For regression tasks they analyze standardized prediction errors by a blocked factorial ANOVA with dataset as a blocking factor: $2^K$19 On MoleculeNet regression datasets, initialization had $2^K$20, $2^K$21, partial $2^K$22; fusion had $2^K$23, $2^K$24, partial $2^K$25; initialization–fusion interaction had $2^K$26, $2^K$27, partial $2^K$28; fusion–update interaction had $2^K$29, $2^K$30, partial $2^K$31; and update had a weak, non-significant Friedman test result of approximately $2^K$32 (Jiao et al., 28 May 2026).

Their mechanistic comparison between Hadamard gating and Concat4 further links the ANOVA findings to representation geometry. Concat4 computes

$2^K$33

whereas Hadamard gating uses

$2^K$34

Pairwise cosine distances and MDS visualizations on Quinethazone showed that Concat4 preserved larger inter-group distances among chemically distinct atoms and resisted oversmoothing better than Hadamard gating (Jiao et al., 28 May 2026). This suggests that factorial regression can function both as a causal or predictive model and as a structured inferential framework for decomposing performance variation across modular design choices.

7. Conceptual distinctions and recurrent methodological issues

A recurrent source of confusion is the conflation of factorial regression with generic factor regression. The distinction is explicit in the SMFR summary: classical factorial regression uses prespecified categorical factors and interaction terms, whereas SMFR learns continuous latent factors from numeric predictors (Kharratzadeh et al., 2015). The former is primarily about treatment coding and causal contrasts; the latter is about low-rank representation, sparse loading structure, and shared latent dimensions.

A second issue concerns the relation between regression-based and randomization-based inference. The design-based literature rejects the idea that ordinary regression inference is merely model-based approximation in factorial experiments. Lu proves that, for $2^K$35 designs, OLS point estimates and amended Huber–White robust covariance estimates coincide with the corresponding randomization-based estimators and Neyman-style standard errors (Lu, 2016). Zhao and Ding extend this logic to more general factor-based specifications, clarifying the design-based properties of both saturated and unsaturated regressions and recommending robust standard errors for Wald-type inference (Zhao et al., 2021).

A third issue is the bias–variance frontier induced by parsimony. Forward selection, additive adjustment, and RLS are all motivated by the exponential growth of factorial terms with $2^K$36. The modern theory does not deny the efficiency benefits of parsimonious models; rather, it makes them conditional on structural assumptions such as heredity, vanishing nuisance effects, or correct linear restrictions. When those assumptions fail, omitted interaction bias is non-negligible and does not wash out asymptotically (Shi et al., 2023, Zhao et al., 2021).

A fourth issue concerns interpretability versus flexibility in high-dimensional settings. MiFM, SMFR, GOFAR, and GFNN all use factorization or low-rank structure to control complexity, but they do so at different levels: MiFM factorizes interaction coefficients, SMFR factorizes the multivariate coefficient matrix, GOFAR factorizes the natural-parameter matrix under mixed outcomes, and GFNN factorizes intermediate representations inside a deep network (Yurochkin et al., 2017, Mishra et al., 2020, Guo et al., 16 Feb 2025). This suggests that the modern landscape of factorial regression is best understood as a spectrum of structured regressions in which orthogonal contrasts, sparse active sets, low-rank factors, or operator-level decompositions are used to make otherwise intractable interaction or dependence structure estimable.