High-dimensional projection heteroscedasticity test

Source: R/modern_diagnostics.R · Help: ?performHighDimensionalTest

Applies principal component projections to the design matrix and evaluates a Breusch-Pagan style statistic on the reduced representation. The procedure is designed for settings where the number of predictors rivals or exceeds the sample size.

Usage

performHighDimensionalTest(
  model,
  data = NULL,
  variance_threshold = 0.9,
  max_components = NULL
)

Arguments

model

A fitted stats::lm object describing the mean structure whose residual variance is to be assessed.

data

The data frame model was fitted to, or NULL, the default. The response, design, weights and offset are read from the fit itself, so a model that transforms its variables, has an offset() term, or takes its weights or offset from a column of the data needs no data, and the model frame of model is accepted as data (since 0.16.0; before, such models needed data). Data that are given are matched to the observations of the fit: rows with row names of their own by name, in any order; rows that only number themselves (a tibble, a data.table, a data frame whose row names were reset) by the row numbers of the fit, or in the order of the fit when they hold one row per observation. The rows the fit dropped may be there or not. Data whose rows cannot be shown to hold the observations of the fit are refused; see test_results.

variance_threshold

Proportion of predictor variance that the selected principal components should explain. Defaults to 0.9.

max_components

Optional upper bound on the number of principal components to retain. Defaults to min(10, n - 5).

Value

An htest object summarising the chi-squared statistic of the auxiliary regression on principal component scores. Since 0.18.0 the test refuses components that fit the squared residuals exactly, by the rule and with the class of performKoenkerTest() (see its Exact fit), where the statistic was the number of observations: on a paired design of 30 subjects of two observations with variance_threshold = 1, which keeps every component, 60 on 30. The default keeps 26 there, which do not fit them exactly, and gives a value. Only a fit the design forces is refused; one the data make exact gives a value, as up to 0.17.0.

Supported models

It accepts an lm and a glm with the Gaussian family and the identity link, unweighted or with positive weights, a fitWLS() fit and, since 0.16.0, a rank-deficient fit, which it tests on the estimable part of its design: it gives the value of the same model fitted without the aliased columns. It refuses a glm of another family or link, a survey::svyglm() fit, a fit with a zero weight, a coefficient that is not finite for another reason, and data that lack rows or variables of the fit. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Fan, J., & Lv, J. (2008). Sure independence screening for ultrahigh dimensional feature space. Journal of the Royal Statistical Society: Series B, 70(5), 849–911.

Examples

set.seed(123)
X <- matrix(rnorm(80 * 20), nrow = 80)
beta <- c(rep(1, 5), rep(0, 15))
y <- X %*% beta + rnorm(80, sd = 0.5 + 0.2 * scale(X[, 1]))
df <- as.data.frame(cbind(y = as.numeric(y), X))
model <- lm(y ~ ., data = df)
performHighDimensionalTest(model, df)