Wald Tests for Fixed-Effect Contrasts Using Predicted Values

Performs Wald tests on linear contrasts of predicted values obtained from asreml::predict.asreml(). Working directly with predict() output means contrasts are specified through meaningful factor-level labels rather than raw coefficient indices, and the prediction error variance–covariance matrix (pred$vcov) is used directly – no access to model internals required.

Two test types are available for each element of cc:

"con" – contrasts

Test \(H_0: \bm{c}^\top\hat{\bm{\tau}} = 0\) for a specified contrast vector \(\bm{c}\). The Wald statistic is \(W = (\bm{c}^\top\hat{\bm{\tau}})^2 / \bm{c}^\top\text{Vcov}\,\bm{c} \sim \chi^2_1\) (or \(F_{1,\nu}\) with test = "F"). Set comp = "pairwise" to automatically generate all \(\binom{k}{2}\) pairwise contrasts for \(k\) specified levels.

"zero" – joint zero equality

Test \(H_0: \bm{Z}\hat{\bm{\tau}} = \bm{0}\) for a set of \(q\) coefficients simultaneously. The statistic is \(W = \hat{\bm{\tau}}^\top\bm{Z}^\top(\bm{Z}\,\text{Vcov}\,\bm{Z}^\top)^{-1} \bm{Z}\hat{\bm{\tau}} \sim \chi^2_q\) (or \(F_{q,\nu}/q\)).

Usage

waldTest(
  pred,
  cc,
  by = NULL,
  test = c("Wald", "F"),
  df_error = NULL,
  adjust = c("none", "bonferroni", "holm", "fdr", "BH", "BY")
)

waldTest.asreml(
  object,
  classify,
  cc,
  by = NULL,
  test = c("Wald", "F"),
  adjust = c("none", "bonferroni", "holm", "fdr", "BH", "BY"),
  ...
)

Arguments

pred

List returned by predict(model, classify = ..., vcov = TRUE). Must contain elements pvals (data frame of predictions) and vcov (prediction error variance–covariance matrix).

cc

Named list of test specifications. Each element is itself a list with fields:

coef

Integer indices or character labels identifying which rows of pred$pvals are involved. When by is used, indices/labels refer to rows within the group. Character labels are matched against the pasted factor-level combination (e.g. "T0:G01" for a Treatment:Genotype prediction).

type

"con" (contrast) or "zero" (joint zero test).

comp

For type = "con" only. One of: a numeric contrast vector (length = length(coef)); a numeric contrast matrix (rows = contrasts, cols = length(coef)); "pairwise" to generate all pairwise differences automatically; NULL or "helmert" for the default last Helmert contrast (a warning is issued).

group

Optional named list list(left = ..., right = ...) to override the auto-generated "A vs B" row labels.

by

Character string or character vector naming one or more factor columns in pred$pvals to split by. When a single name is supplied the rows are split by that column's levels. When a vector of names is supplied (e.g. by = c("Site", "Year")) the levels of those columns are pasted together to form an interaction grouping variable, and the output grouping column is named with the same pasted string (e.g. "Site:Year"). The same cc specification is applied independently within each group, with coef referring to within-group row positions. Default NULL (no splitting).

test

"Wald" (default) for asymptotic \(\chi^2\) tests, or "F" for \(F\)-tests. "F" requires df_error.

df_error

Numeric error degrees of freedom for F-tests (e.g. model$nedf). Ignored when test = "Wald".

adjust

P-value adjustment method passed to p.adjust(). Applied across all contrasts within each group. Options:

"none"

No adjustment (default).

"bonferroni"

Bonferroni correction – controls family-wise error rate (FWER) conservatively.

"holm"

Holm step-down – controls FWER less conservatively than Bonferroni.

"fdr" / "BH"

Benjamini–Hochberg False Discovery Rate. Controls the expected proportion of false discoveries among rejected hypotheses. Best choice when many contrasts are tested and some true effects are expected.

"BY"

Benjamini–Yekutieli FDR. Valid under arbitrary (including positive) dependence between test statistics; more conservative than "BH" but provides stronger guarantees.

object

An ASReml-R V4 model object.

classify

Character string passed to asreml::predict.asreml() as the classify argument.

...

Additional arguments forwarded to asreml::predict.asreml().

Value

A named list with elements:

Contrasts

Data frame of contrast test results (or NULL if no "con" elements were supplied). Columns: grouping variable (if by used), Comparison, Estimate, Std.Error, Wald.Statistic or F.Statistic, df (numerator), P.Value.

Zero

Data frame of joint zero test results (or NULL). Columns: grouping variable (if by used), Test, Wald.Statistic or F.Statistic, df, P.Value.

test

The test argument used.

adjust

The adjust argument used.

Functions

• waldTest.asreml(): Convenience method that calls asreml::predict.asreml() internally. classify is required; test = "F" uses object$nedf automatically.

See also

asreml::predict.asreml()

Examples

if (FALSE) { # \dontrun{
## --- predict first, then test ---
pred <- predict(model, classify = "Treatment", vcov = TRUE)

## Pairwise contrasts among all three treatment levels
res <- waldTest(pred,
                cc = list(list(coef = c("N0","N1","N2"),
                               type = "con",
                               comp = "pairwise")))

## Manual contrast: N2 vs N0
res2 <- waldTest(pred,
                 cc = list(list(coef = c("N0","N2"),
                                type = "con",
                                comp = c(-1, 1),
                                group = list(left="N2", right="N0"))))

## Joint zero test: are N1 and N2 effects simultaneously zero?
res3 <- waldTest(pred,
                 cc = list(list(coef = c("N1","N2"),
                                type = "zero",
                                group = "N1:N2 joint")))

## F-tests with denominator df from model
res4 <- waldTest(pred, cc = list(...),
                 test = "F", df_error = model$nedf)

## Run within each site (by-group, single factor)
pred_site <- predict(model, classify = "Treatment:Site", vcov = TRUE)
res5 <- waldTest(pred_site,
                 cc   = list(list(coef = c("N0","N1","N2"),
                                  type = "con", comp = "pairwise")),
                 by   = "Site",
                 test = "F", df_error = model$nedf,
                 adjust = "fdr")

## Run within each Site x Year combination (multi-factor by)
pred_sy <- predict(model, classify = "Treatment:Site:Year", vcov = TRUE)
res5b <- waldTest(pred_sy,
                  cc  = list(list(coef = c("N0","N1","N2"),
                                  type = "con", comp = "pairwise")),
                  by  = c("Site", "Year"),
                  test = "F", df_error = model$nedf)

## Convenience: pass the model directly
res6 <- waldTest(model, classify = "Treatment",
                 cc = list(list(coef=c("N0","N2"), type="con", comp=c(-1,1))),
                 test = "F")
} # }