Model comparison for nested regression models that are not symbolically nested

Question

Say I have two nested models and I want to compare them, but they are parameterized differently so that no simple constraint (i.e., setting coefficients to 0) on the larger model corresponds exactly to the smaller model (e.g., because they use different basis functions). We might say they are nested but not symbolically nested.

For example, suppose model 1 is a regression of the outcome on a 4-df natural cubic spline of a predictor, and model 2 is a regression of the outcome on just the predictor as a linear term. These two models are nested, and yet there is no immediately obvious constraint on the spline model that would yield the same parameterization as the linear model. Usually, this isn't a problem, because we can compare nested models using an ANOVA or likelihood ratio (LR) test. But in my case, the model doesn't use the usual OLS or MLE variance matrix, so these tests would not be expected to be accurate. Normally, I would use a Wald test to test the constraints of the larger model that yield a model equivalent to the smaller, but that only works when the two models have coefficients that meaningfully correspond.

Essentially, I'm hoping to find a way to take in the output of the larger model and the output of the smaller model and compute the Wald test statistic assuming the models are nested without needing to specify the constraints that are usually required for a Wald test, but correctly incorporate custom coefficient variances that otherwise cannot be used with the ANOVA or LR test.

Below is an example in R to help fix ideas:

data("lalonde", package = "MatchIt")

fit1 <- lm(re78 ~ splines::ns(age, df = 4), data = lalonde)
fit2 <- lm(re78 ~ age, data = lalonde)

#Works for basic comparison
anova(fit2, fit1)
#> Analysis of Variance Table
#> 
#> Model 1: re78 ~ age
#> Model 2: re78 ~ splines::ns(age, df = 4)
#>   Res.Df        RSS Df Sum of Sq      F  Pr(>F)   
#> 1    612 3.3826e+10                               
#> 2    609 3.3007e+10  3 818987331 5.0369 0.00187 **
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Doesn't work for sandwich SEs
lmtest::waldtest(fit2, fit1, vcov = sandwich::sandwich)
#> Error in modelCompare(objects[[i - 1]], objects[[i]], vfun = vcov0): nesting of models cannot be determined

^{Created on 2024-08-09 with reprex v2.1.1}

If one knew how to parameterize the larger model so that the two were symbolically nested, one could perform a Wald test without difficulty, but this isn't always straightforward, and I'm wondering if there is a general solution. The internals of survey:::anova.svyglm() provide some clue for comparing such models, but it seems like the method they use is specific to the type of model that was fit (i.e., wouldn't generalize to arbitrary models like a Wald test on symbolically nested models or LR test would otherwise).

Answer 1

For a linear model, say that the full model is $\mathbf{Y} = \mathbf{X} \boldsymbol\beta + \mathbf{e}_f$ and the reduced model is $\mathbf{Y} = \mathbf{R} \boldsymbol\alpha + \mathbf{e}_r$. $\mathbf{X}$ is $N \times f$ and $\mathbf{R}$ is $N \times r$ where $r < f$. If the models are nested, then there must be a $(f - r) \times f$ matrix $\mathbf{C}$ such that $\mathbf{X} \boldsymbol\beta = \mathbf{R} \boldsymbol\alpha$ whenever $\mathbf{C}\boldsymbol\beta = \mathbf{0}$. The question is how to find $\mathbf{C}$.

Regressing $\mathbf{X}$ on $\mathbf{R}$ will give an $N \times f$ matrix of residuals that is orthogonal to $\mathbf{R}$. Let $\boldsymbol\gamma = \left(\mathbf{R}'\mathbf{R}\right)^{-1} \mathbf{R}' \mathbf{X}$ so that $$ \mathbf{X} = \mathbf{R} \boldsymbol\gamma + \mathbf{\tilde{X}}, $$ where $\mathbf{\tilde{X}} = \mathbf{X} - \mathbf{R}\boldsymbol\gamma$. $\mathbf{\tilde{X}}$ is $N \times f$ but has column rank of $f - r$. Take the singular value decomposition of $\mathbf{\tilde{X}}$ so $\mathbf{\tilde{X}} = \mathbf{U}\mathbf{D} \mathbf{V}'$ where $\mathbf{U}$ is $N \times (f - r)$ and orthogonal, $\mathbf{D}$ is $(f - r) \times (f - r)$ and diagonal, and $\mathbf{V}$ is $f \times (f - r)$ and orthogonal. Then $$ \mathbf{X} \boldsymbol\beta = \mathbf{R} \boldsymbol\gamma \boldsymbol\beta + \mathbf{U}\mathbf{D}\mathbf{V}'\boldsymbol\beta. $$ If $\mathbf{V}' \boldsymbol\beta = 0$ then $\mathbf{X} \boldsymbol\beta = \mathbf{R} \boldsymbol\gamma \boldsymbol\beta$. So we can use $\mathbf{C} = \mathbf{V}'$, and $\boldsymbol\alpha = \boldsymbol\gamma \boldsymbol\beta$.

In R:

data("lalonde", package = "MatchIt")
fit1 <- lm(re78 ~ splines::ns(age, df = 4), data = lalonde)
fit2 <- lm(re78 ~ age, data = lalonde)

X <- model.matrix(fit1)
U <- model.matrix(fit2)
UX_fit <- lm.fit(U, X)
gamma <- coefficients(UX_fit)
X_tilde <- residuals(UX_fit)
X_svd <- svd(X_tilde)
keep <- X_svd$d > 1e-7
Cmat <- t(X_svd$v[,keep,drop=FALSE])

car::linearHypothesis(
  model = fit1, 
  hypothesis.matrix = Cmat,
  vcov. = sandwich::sandwich
)