Design-based test for nested survey-weighted GLMs

Compare two nested survey::svyglm() models with a single design-based hypothesis test. The function reports whether the extra terms in the larger ("full") model jointly improve on the smaller ("reduced") model, accounting for the survey design.

The motivating use case is checking for nonlinearity. Comparing a model that enters a predictor linearly against a model that enters the same predictor through a spline basis such as splines::ns(). A small p-value is evidence that the linear form is inadequate.

Usage

svyglm_nested_test(
  reduced,
  full,
  design,
  family = quasibinomial(),
  method = c("LRT", "Wald"),
  tol = 1e-07
)

Arguments

reduced: (formula)
The smaller model. Must include the response on the left-hand side.
full: (formula)
The larger model, which must contain reduced as a nested special case. If supplied one-sided, the response is taken from reduced.
design: (survey.design or svyrep.design)
A survey design object created by survey::svydesign(), survey::svrepdesign(), or a calibrated or post-stratified variant.
family: (family object: default quasibinomial())
The error family passed to survey::svyglm(). The default suits logistic regression of a binary outcome and avoids the spurious non-integer-count warning produced by binomial().
method: (Scalar character: default "LRT")
The design-based test to apply. "LRT" gives the Rao-Scott working likelihood ratio test. "Wald" gives the design-based Wald test.
tol: (Positive scalar numeric: default 1e-7)
Relative tolerance used to verify nesting and to confirm that the extra directions are cleanly separated from the reduced model.

Value

A list of class svyglm_nested_test with components:

test: the test object returned by stats::anova() or survey::regTermTest().
method: the test method used, "LRT" or "Wald".
df: the number of extra directions, equal to the test degrees of freedom.
n: the number of observations used by both fitted models.
fit_full: the fitted augmented full model, a svyglm object. Coefficients on the .nl* columns are on an arbitrary orthonormal scale and are not interpretable as effect sizes. The joint test is invariant to this scaling.
fit_reduced: the fitted reduced model, a svyglm object.

The print method shows the method, degrees of freedom, sample size, and the underlying test.

Details

Why not call `anova()` or `regTermTest()` directly

Two obstacles make a direct comparison not possible.

When the linear and spline models are fitted by separate survey::svyglm() calls, each call silently drops its own complete cases. If both formulas use the same variables the row sets happen to be identical, but that coincidence does not hold when either formula introduces additional predictors or interactions. Calling anova() on two separately-fitted objects can therefore compare fits on different samples, making the test ambiguous. This function removes the ambiguity: it pools all variables from both formulas, builds one common complete-case set, and fits both models on that exact set.
A linear term and a spline term are different parameterizations of the same predictor. The spline basis columns do not split into a "linear part" and a "nonlinear part", so regTermTest() cannot be pointed at the nonlinear piece of the spline model directly.

How the function works

The reduced model is nested within the full model whenever the column space of the reduced design matrix is contained in the column space of the full design matrix. This holds for the linear-versus-spline comparison because a natural spline basis spans the linear function of the predictor.

The function refits the full model in an equivalent parameterization:

Build the design matrices of both models on a common complete-case set.
Confirm the reduced model is nested within the full model.
Compute the part of the full design matrix orthogonal to the reduced design matrix. The number of these "extra" directions equals \(\text{rank}(\text{full}) - \text{rank}(\text{reduced})\). An orthonormal basis is extracted by singular value decomposition, which is numerically stable even though the residualized columns that belong to the reduced model are only zero up to roundoff.
Fit an augmented model that contains the original reduced terms plus the isolated extra directions as plain numeric columns.

The augmented full model has the same column space, and therefore the same fit, as the original full model. Because it contains no spline call, the design-based tests apply cleanly.

The test

Survey-weighted GLMs are fitted by pseudo-likelihood and have no genuine likelihood, so the test is a design-based analogue of a likelihood ratio test.

method = "LRT" uses stats::anova() on the two fitted models. This is the Rao-Scott working likelihood ratio test. The working \(2\log\text{LR}\) statistic is referred to a weighted sum of \(\chi^2_1\) variables whose scale factors are estimated from the design.
method = "Wald" uses survey::regTermTest() to test that the coefficients of the extra directions are jointly zero, using the design-based variance estimator.

The working LRT is generally preferred because it is less sensitive to the parameterization and to small-sample instability of the variance estimator, but the two should agree closely in large samples.

The test always has \(\text{rank}(\text{full}) - \text{rank}(\text{reduced})\) degrees of freedom and is a single joint test of all extra directions at once.

Assumptions and scope

The reduced model must be nested within the full model. The function checks this and errors otherwise.
Both models are fitted on the same rows. Rows with missing values in any variable used by either formula are dropped. If the full formula introduces variables, with their own missingness, that are absent from the reduced formula, the function errors rather than comparing models fitted on different samples.
The orthogonalization that isolates the extra directions is unweighted. This is intentional. Equivalence of column spaces is a property of the design matrices alone, so the survey weights still enter the fit and the test through survey::svyglm().

Interpreting the result

The p-value answers one question: do the extra terms in the full model jointly carry signal beyond the reduced model, given the survey design.

Linear versus spline, main effect (y ~ pred + cov against y ~ ns(pred, df) + cov):

A small p-value is evidence that the relationship between pred and the outcome departs from linearity on the link (here, logit) scale. The spline model is preferred.
A large p-value means there is no evidence of nonlinearity in the available data. This supports keeping the simpler linear term, but it is absence of evidence, not proof that the effect is exactly linear. Power is limited by sample size, design effects, and the chosen spline degrees of freedom.

Linear versus spline, with an interaction (y ~ pred * cov against y ~ ns(pred, df) * cov):

This is a joint test. It tests nonlinearity in the main effect of pred and nonlinearity in the pred-by-cov interaction simultaneously.
A small p-value means the linear-by-linear specification is inadequate somewhere, but does not say where. To localize the source, also run the main-effect comparison above. A nonlinear main effect with a non-significant joint-minus-main difference points to curvature that does not vary across cov.

General nested comparison:

For any nested pair, the result tests whether the additional terms in the full model are jointly zero. A small p-value favors the full model. A large p-value indicates the reduced model is adequate for the data at hand.

Cautions:

The degrees of freedom equal the number of extra directions. A spline with df = 4 against a linear term contributes 3 nonlinear degrees of freedom. The same comparison crossed with a continuous covariate contributes 6. Larger df spreads power across more directions.
Statistical significance is not the same as practical importance. Inspect the fitted spline curve, for example with survey::svycontrast() or predicted probabilities, before deciding on a final model.
If many such tests are run, adjust for multiple comparisons.

Examples

#----------------------------------------------------------------------------
# svyglm_nested_test() examples
#----------------------------------------------------------------------------
library(bkstat)
library(survey)
#> Loading required package: grid
#> Loading required package: Matrix
#> 
#> Attaching package: ‘survey’
#> The following object is masked from ‘package:rms’:
#> 
#>     calibrate
#> The following object is masked from ‘package:Hmisc’:
#> 
#>     deff
#> The following object is masked from ‘package:graphics’:
#> 
#>     dotchart

# Stratified sample of California schools, with a binary outcome.
data(api)
apistrat$sw <- as.integer(apistrat$sch.wide == "Yes")
design <- svydesign(
  id = ~1,
  strata = ~stype,
  weights = ~pw,
  fpc = ~fpc,
  data = apistrat
)

#----------------------------------------------------------------------------
# 1. Is the effect of `ell` nonlinear, adjusting for school type?
#    Linear term versus a natural spline with 4 degrees of freedom.
#----------------------------------------------------------------------------
fit1 <- svyglm_nested_test(
  reduced = sw ~ ell + stype,
  full = sw ~ splines::ns(ell, 4) + stype,
  design = design
)
fit1
#> Survey nested-model test (LRT), 3 df, n = 200
#> 
#> Working (Rao-Scott+F) LRT for .nl1 .nl2 .nl3
#>  in svyglm(formula = f_full, design = des, family = family)
#> Working 2logLR =  8.670567 p= 0.03706 
#> (scale factors:  1.2 0.95 0.85 );  denominator df= 191
# A 3 df test. A small p-value favors the spline; a large p-value
# indicates no evidence against a linear `ell` effect.

#----------------------------------------------------------------------------
# 2. Joint nonlinearity test when `ell` interacts with school type.
#----------------------------------------------------------------------------
fit2 <- svyglm_nested_test(
  reduced = sw ~ ell * stype,
  full = sw ~ splines::ns(ell, 4) * stype,
  design = design
)
fit2
#> Survey nested-model test (LRT), 9 df, n = 200
#> 
#> Working (Rao-Scott+F) LRT for .nl1 .nl2 .nl3 .nl4 .nl5 .nl6 .nl7 .nl8 .nl9
#>  in svyglm(formula = f_full, design = des, family = family)
#> Working 2logLR =  23.03396 p= 0.02386 
#> (scale factors:  2.8 1.7 1.1 0.87 0.74 0.58 0.54 0.43 0.3 );  denominator df= 183
# A joint test of nonlinearity in the main effect and in the interaction.

#----------------------------------------------------------------------------
# 3. Sanity check: with no splines the result must match anova().
#----------------------------------------------------------------------------
fit3 <- svyglm_nested_test(
  reduced = sw ~ ell + meals,
  full = sw ~ ell * meals,
  design = design
)
fit3
#> Survey nested-model test (LRT), 1 df, n = 200
#> 
#> Working (Rao-Scott+F) LRT for .nl1
#>  in svyglm(formula = f_full, design = des, family = family)
#> Working 2logLR =  0.7034813 p= 0.40262 
#> df=1;  denominator df= 194

mod_r <- svyglm(sw ~ ell + meals, design = design, family = quasibinomial())
mod_f <- svyglm(sw ~ ell * meals, design = design, family = quasibinomial())
anova(mod_r, mod_f)
#> Working (Rao-Scott+F) LRT for ell:meals
#>  in svyglm(formula = sw ~ ell * meals, design = design, family = quasibinomial())
#> Working 2logLR =  0.7034813 p= 0.40262 
#> df=1;  denominator df= 194

#----------------------------------------------------------------------------
# 4. The design-based Wald test instead of the working LRT.
#----------------------------------------------------------------------------
svyglm_nested_test(
  reduced = sw ~ ell + stype,
  full = sw ~ splines::ns(ell, 4) + stype,
  design = design,
  method = "Wald"
)
#> Survey nested-model test (Wald), 3 df, n = 200
#> 
#> Wald test for .nl1 .nl2 .nl3
#>  in svyglm(formula = f_full, design = des, family = family)
#> F =  2.944538  on  3  and  191  df: p= 0.034203