What is a partial chi-square statistic according to Frank Harrell?

Question

In his RMS course (section 4.1.1), Frank Harrell mentions the use of a partial chi square statistic for measuring the strength of association between a predictor and an outcome. See below for a screenshot of this passage:

Given a model, the differences between partial chi-squares and their respective degrees of freedom can be plotted using plot(anova(mod). Here is an example

From what test are these chi-squares obtained? How might I compute these quantities in base R?

Answer 1

For other than ordinary least squares (OLS) regression, the anova() function in Harrell's rms package performs Wald tests on individual coefficients and sets of related coefficients; Wald tests are an option for OLS models. The Wald $\chi^2$ statistic used in the test for a coefficient or a set of coefficients is the "partial $\chi^2$ statistic."

The code is in the ava() function defined at the start of rms:::anova.rms. For a vector of coefficient estimates coefs and the corresponding subset of the covariance matrix vcov(coefs), it's just the quadratic form combining the coefs with the inverse of vcov(coefs). See this answer for a simple implementation in base R. Subtracting the number of degrees of freedom corrects for the mean $\chi^2$ under the null hypothesis.

For OLS, plot.anova.rms() by default multiplies similar partial F-statistics by the (numerator) degrees of freedom to get $\chi^2$ values. It can, upon request, display corresponding partial $R^2$ values.

Answer 2

This is an appendix to @EdM answer (+1). I look at the implementation (in base R) to make sure I understand the partial $\chi^2$ statistic. The R code borrows heavily from the rms::anova. This is for illustration only, so there are restrictions; most importantly, it's assumed the predictors have unique names.

library("rms")

getHdata(nhgh)

g <- function(x) 0.09 - x^-(1 / 1.75)
formula <- g(gh) ~ rcs(age, 4) + re + sex + rcs(bmi, 4)
plot(anova(
  ols(formula, data = nhgh)
))

# Fit model in base R
model <- lm(formula, data = nhgh)

# `age` is transformed into a restricted cubic spline with 4 knots,
# so there are 3 components.
associated_terms(model, "age")
#> [1] "rcs(age, 4)age"   "rcs(age, 4)age'"  "rcs(age, 4)age''"

The partial $\chi^2$ statistic is $\hat{\beta}_{S}^\top\widehat{\Sigma}_{S}^{-1}\hat{\beta}_S$ where $S$ is the set of terms associated with the predictor (linear, nonlinear, interactions), $\hat{\beta}_S$ is the corresponding subset of coefficient estimates and $\widehat{\Sigma}_S$ is their covariance matrix.

rbind(
  partial_chisq(model, "sex"),
  partial_chisq(model, "re"),
  partial_chisq(model, "bmi"),
  partial_chisq(model, "age")
)
#>   predictor       chi2 df             P
#> 1       sex   15.17428  1  9.802941e-05
#> 2        re  172.89056  4  2.506276e-36
#> 3       bmi  332.38234  3  9.730894e-72
#> 4       age 1324.07373  3 8.795631e-287

---

The R implementation in full.

library("rms")

# Find model terms which are a function of the given predictor,
# including (linear and nonlinear) main effects and interactions.
#
# @param model: A fitted `lm` or `glm` model.
# @param predictor: The name of a single predictor.
#
# Caution!
# This function assumes predictors have unique names.
associated_terms <- function(model, predictor) {
  terms <- names(coef(model))
  terms[grepl(predictor, terms, perl = TRUE)]
}

# Compute t(x) @ inv(V) @ x
# @param V: a square n-by-n matrix.
# @param x: a n-dimensional vector.
compute_quadratic <- function(V, x) {
  x %*% solvet(V, x, tol = 1e-9)
}

# Compute the Wald chi squared statistic for a subset of model terms.
partial_chisq <- function(model, predictor) {
  terms <- associated_terms(model, predictor)

  b <- coef(model)
  V <- vcov(model)

  idx <- names(b) %in% terms

  chi2 <- compute_quadratic(V[idx, idx], b[idx])
  df <- sum(idx)

  data.frame(
    predictor, chi2, df,
    P = pchisq(chi2, df, lower.tail = FALSE)
  )
}


# BBR, Section 4.3.5
getHdata(nhgh)

g <- function(x) 0.09 - x^-(1 / 1.75)
ginverse <- function(y) (0.09 - y)^-1.75

formula <- g(gh) ~ rcs(age, 4) + re + sex + rcs(bmi, 4)
plot(anova(
  ols(formula, data = nhgh)
))

# Fit model in base R
model <- lm(formula, data = nhgh)

# `age` is transformed into a restricted cubic spline with 4 knots,
# so there are 3 components.
associated_terms(model, "age")

rbind(
  partial_chisq(model, "sex"),
  partial_chisq(model, "re"),
  partial_chisq(model, "bmi"),
  partial_chisq(model, "age")
)

# BBR, Section 19.8
getHdata(acath)
acath <- subset(acath, !is.na(choleste))

formula <- sigdz ~ sex * rcs(age, 5)
plot(anova(
  lrm(formula, data = acath)
))

# Fit model in base R
model <- glm(formula, family = binomial, data = acath)

# The "contribution" of `age` includes a restricted cubuc spline *and*
# the interaction with `sex`.
associated_terms(model, "age")

rbind(
  partial_chisq(model, "age"),
  partial_chisq(model, "sex")
)