Multivariate mean inference and linear contrasts

The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

This vignette uses synthetic data: a trivariate vector autoregression with cross-dependence. It illustrates joint inference on the mean vector, the increment-hull confidence region, simultaneous intervals for linear contrasts, and the distinction between simultaneous and marginal intervals.

library(aersn)
set.seed(11)
n <- 250; q <- 3
A <- matrix(c(0.4, 0.1, 0, -0.2, 0.3, 0.1, 0, 0.1, 0.5), 3, 3, byrow = TRUE)
e <- matrix(rnorm(n * q), n, q) %*% chol(matrix(c(1, .5, .2, .5, 1, .3, .2, .3, 1), 3))
Y <- e
for (t in 2:n) Y[t, ] <- Y[t - 1, ] %*% t(A) + e[t, ]
Y <- sweep(Y, 2, c(0, 0.2, -0.1), `+`)
fit <- aersn_mean(Y, names = c("m1", "m2", "m3"))
fit
#> Affine-equivariant adjusted-range self-normalization (increment hull)
#>   model: sample mean (psi_t = Y_t - Ybar) 
#>   n = 250 observations; q = 3 parameters
#>   centering: calendar time, tau(r) = r 
#>   estimate:
#>        m1        m2        m3 
#> -0.007784  0.216400  0.054080 
#>   hull diagnostics: numerical rank 3 of 3  ; condition 5.042 ; min projected range 0.9794 ; spread 3.069 
#> Use aersn_test(), confint(), aersn_contrast(), aersn_region(), plot().

The hull diagnostics report the numerical rank of the centered path (it must equal q), its condition number, and the minimum projected adjusted range over the coordinate axes and 1,000 prespecified directions.

Joint test and confidence region

ref <- aersn_reference(fit, draws = 4000, seed = 1)   # matched grid, q = 3, n = 250
ref
#> Reference distribution (aersn_reference): increment-hull gauge
#>   matched-grid Monte Carlo law for increment-hull gauge: q = 3, uniform grid with n = 250 intervals, 4000 draws, seed 1
#>   quantiles (type 8):
#>      90.0%: 2.7628  (Monte Carlo s.e. 0.025)
#>      95.0%: 3.1539  (Monte Carlo s.e. 0.029)
#>      99.0%: 4.1278  (Monte Carlo s.e. 0.067)
aersn_test(fit, null = c(0, 0, 0), reference = ref)
#> 
#>  Adjusted-range increment hull test, q = 3
#> 
#> data:  fit
#> T = 3.173, q = 3, n = 250
#> alternative hypothesis: true parameter is not equal to the null value
#> null value: m1 = 0, m2 = 0, m3 = 0 
#> estimate:   m1 = -0.007784, m2 = 0.2164, m3 = 0.05408 
#> critical value at level 0.95: 3.154  (reject)
#> p-value = 0.0495 (Monte Carlo s.e. 0.0034, resolution 0.00025)
#> reference: matched-grid Monte Carlo law for increment-hull gauge: q = 3, uniform grid with n = 250 intervals, 4000 draws, seed 1
reg <- aersn_region(fit, reference = ref)
reg
#> Adjusted-range increment hull confidence region, level 0.95
#>   q = 3  n = 250  critical value 3.154  scale c/sqrt(n) = 0.1995 
#>   center (estimate):
#>        m1        m2        m3 
#> -0.007784  0.216400  0.054080 
#>   coordinate projections (simultaneous):
#>       lower  upper
#> m1 -0.41970 0.4041
#> m2 -0.07076 0.5036
#> m3 -0.40060 0.5088
#>   reference: matched-grid Monte Carlo law for increment-hull gauge: q = 3, uniform grid with n = 250 intervals, 4000 draws, seed 1

The region is \(\hat\theta_n + (c/\sqrt n)K(\hat G_n)\), a translated and scaled copy of the increment hull. It is convex and centrally symmetric but not an ellipsoid. Membership of candidate vectors is tested through the gauge:

aersn_contains(reg, rbind(c(0, 0, 0), c(0, 0.2, -0.1)))
#> [1] FALSE  TRUE
#> attr(,"gauge")
#> [1] 1.0059492 0.5963886

For three or more parameters, plot() shows two-dimensional views: the default type = "projection" draws the exact projections of the joint region onto coordinate planes (simultaneous coverage), while type = "slice" draws cross-sections through the estimate.

plot(reg, null = c(0, 0, 0))

Simultaneous intervals for contrasts

Each contrast \(a^\top\theta\) has the projection interval \(a^\top\hat\theta_n \pm (c_q/\sqrt n)\{\max_k a^\top\hat G_{n,k} - \min_k a^\top\hat G_{n,k}\}\), which uses the \(q\)-dimensional critical value and therefore holds simultaneously for every linear contrast.

contrasts <- rbind("m2 - m1" = c(-1, 1, 0),
                   "m3 - m2" = c(0, -1, 1),
                   "average" = c(1, 1, 1) / 3)
aersn_contrast(fit, contrasts, reference = ref)
#> Simultaneous (joint-region projection) Adjusted-range increment hull intervals for linear contrasts, level 0.95
#>         estimate   lower  upper
#> m2 - m1  0.22420 -0.1222 0.5706
#> m3 - m2 -0.16230 -0.5318 0.2071
#> average  0.08757 -0.2518 0.4270
#> critical value 3.154 (reference dimension 3); matched-grid Monte Carlo law for increment-hull gauge: q = 3, uniform grid with n = 250 intervals, 4000 draws, seed 1

Three interval types are distinguished:

type = "simultaneous" (default): projections of the full joint region, critical value \(c_q\);
type = "joint": projections of the joint region of the specified contrasts only (dimension \(m\)), critical value \(c_m\); and
type = "marginal": separately constructed scalar intervals with \(c_1\), which cover each contrast individually and are not simultaneous.

confint(fit, reference = ref)
#> Simultaneous (joint-region projection) Adjusted-range increment hull confidence intervals, level 0.95
#>       lower  upper
#> m1 -0.41970 0.4041
#> m2 -0.07076 0.5036
#> m3 -0.40060 0.5088
#> critical value 3.154 (reference dimension 3); matched-grid Monte Carlo law for increment-hull gauge: q = 3, uniform grid with n = 250 intervals, 4000 draws, seed 1
confint(fit, type = "marginal", draws = 20000, seed = 1)
#> Marginal (separately constructed) Adjusted-range increment hull confidence intervals, level 0.95
#>       lower  upper
#> m1 -0.25290 0.2374
#> m2  0.04551 0.3873
#> m3 -0.21650 0.3247
#> critical value 1.877 (reference dimension 1); matched-grid Monte Carlo law for increment-hull gauge: q = 1, uniform grid with n = 250 intervals, 20000 draws, seed 1
#> Note: marginal intervals do not provide simultaneous coverage.

A vector can lie inside every marginal interval and still be excluded from the joint region, because some linear combination of the coordinates violates the joint restriction (manuscript Section 3).

Affine equivariance

Rescaling, rotating or shearing the data leaves the test statistic unchanged and maps the region accordingly:

H <- matrix(c(2, 0.5, 0, -0.3, 1, 0.2, 0, 0, 0.7), 3, 3)
fitH <- aersn_mean(Y %*% t(H))
c(original = aersn_gauge(fit, c(0, 0, 0)), transformed = aersn_gauge(fitH, c(0, 0, 0)))
#>    original transformed 
#>     3.17269     3.17269

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.