Technical details
The package provides three standalone tests for scientifically
distinct questions about a set of correlated coefficient estimates:
heterogeneity_test(): a non-directional omnibus test of
equality;linear_trend_test(): a test of whether the coefficients
follow a linear trend with respect to assigned numeric scores;order_restricted_test(): a test of equality against a
monotone increasing or decreasing alternative, without requiring a
linear pattern.
These functions are especially useful for transition-specific
treatment effects, but they apply to any coefficient vector accompanied
by a full variance-covariance matrix. The call of the three test
functions are shown below here:
heterogeneity_test(beta, vcov, index = NULL)
linear_trend_test(beta, vcov, index = NULL, scores = NULL, alternative = c("greater", "less", "two.sided"))
order_restricted_test(beta, vcov, index = NULL, alternative = c("increasing", "decreasing"),
weight_method = c("auto", "simulation"), nsim = 100000L, seed = NULL, weights = NULL, tol = 1e-10)
Omnibus heterogeneity test
The heterogeneity test evaluates
\[
H_0:\beta_1=\cdots=\beta_R
\qquad\text{against}\qquad
H_A:\text{not all }\beta_r\text{ are equal.}
\]
This is a multivariate Wald test of a set of linear equality
constraints (Wald
1943; Cox and Hinkley 1974; Lehmann and Romano 2022). Let
\(D\) be the \((R-1)\times R\) adjacent-contrast matrix,
so that
\[
D\widehat{\boldsymbol\beta}
=
(\widehat\beta_2-\widehat\beta_1,\ldots,
\widehat\beta_R-\widehat\beta_{R-1})^\top .
\]
The omnibus Wald statistic is
\[
W=
(D\widehat{\boldsymbol\beta})^\top
(D\widehat V D^\top)^{-1}
(D\widehat{\boldsymbol\beta}) .
\]
Under \(H_0\), \(W\) has an asymptotic chi-square
distribution with \(R-1\) degrees of
freedom. The test is non-directional: positive and negative departures
from equality are treated symmetrically. Using the upper tail of the
chi-square distribution for the p-value does not make it a one-sided
test. With two coefficients it is equivalent to a two-sided Wald \(Z\) test, since \(W=Z^2\).
Generalized least-squares linear-trend test
The linear-trend test models the selected coefficients as
\[
\beta_r=\eta+\gamma s_r,\qquad r=1,\ldots,R,
\]
where \(s_1<\cdots<s_R\) are
numeric scores, \(\eta\) is an
intercept, and \(\gamma\) is the trend
slope. The linear model specifies the relationship being tested;
generalized least squares (GLS) is the estimation method, which accounts
for the unequal variances and correlations in \(\widehat V\). Score-based tests for linear
trend have a long history, including the work of Cochran and Armitage
(Cochran 1954; Armitage 1955).
With design matrix \(A=(\boldsymbol 1,\
\boldsymbol s)\), the GLS estimator treats \(\widehat{\boldsymbol\beta}\) as the
response with known covariance \(\widehat
V\):
\[
\widehat{\boldsymbol\theta}
=
\begin{pmatrix}\widehat\eta\\ \widehat\gamma\end{pmatrix}
=
(A^\top\widehat V^{-1}A)^{-1}
A^\top\widehat V^{-1}\widehat{\boldsymbol\beta},
\qquad
\widehat{\mathrm{Var}}(\widehat{\boldsymbol\theta})
=
(A^\top\widehat V^{-1}A)^{-1} .
\]
This covariance-weighted estimator follows the generalized
least-squares principle associated with Aitken’s result and its later
regression and GLS treatments (Aitken 1935; Seber and Lee 2003; Kariya and Kurata 2004). Because
the covariance \(\widehat V\) is
treated as known, the slope is tested with a standard-normal reference
rather than a \(t\) reference:
\[
Z=\frac{\widehat\gamma}{\mathrm{SE}(\widehat\gamma)},
\qquad
\mathrm{SE}(\widehat\gamma)=\sqrt{\big[(A^\top\widehat
V^{-1}A)^{-1}\big]_{22}} ,
\]
which is asymptotically standard normal at \(\gamma=0\). The available options for
alternative are:
"greater": \(H_0:\gamma\leq0\) versus \(H_A:\gamma>0\);"less": \(H_0:\gamma\geq0\) versus \(H_A:\gamma<0\);"two.sided": \(H_0:\gamma=0\) versus \(H_A:\gamma\neq0\).
If scores is omitted, the equally spaced values
0, 1, ..., R - 1 are used. User-supplied scores must be
finite and strictly increasing. They may have length
length(beta), in which case the selected entries are
extracted with index, or length equal to the number of
selected coefficients. Scores must be supplied when index
selects nonconsecutive positions, so that the intended spacing is
explicit.
A positive slope means the coefficients tend to increase as the
scores increase, and a negative slope means they tend to decrease. The
observed estimates need not be monotone. A significant slope is evidence
of a linear trend; it does not by itself establish that a linear model
describes the pattern adequately.
Order-restricted test
For alternative = "increasing", the order-restricted
test evaluates
\[
H_0:\beta_1=\cdots=\beta_R
\qquad\text{against}\qquad
H_A:\beta_1\leq\cdots\leq\beta_R,
\]
with at least one strict inequality. For
alternative = "decreasing", the inequalities are reversed.
Unlike the linear-trend test, this procedure does not assume constant or
linearly changing adjacent effects. The general theory and methods of
order-restricted inference are described by Robertson, Wright, and
Dykstra and by Silvapulle and Sen (Robertson, Wright, and Dykstra 1988; Silvapulle and Sen 2005).
For an increasing alternative, define the adjacent differences and
their covariance,
\[
\widehat{\boldsymbol\delta}
=D\widehat{\boldsymbol\beta}
=(\widehat\beta_2-\widehat\beta_1,\ldots,
\widehat\beta_R-\widehat\beta_{R-1})^\top,
\qquad
\widehat\Omega=D\widehat V D^\top .
\]
For a decreasing alternative, the sign of \(D\) is reversed, so that positive contrasts
always point in the direction of the specified alternative. Let \(\widetilde{\boldsymbol\delta}\) be the
covariance-weighted projection of \(\widehat{\boldsymbol\delta}\) onto the
nonnegative orthant:
\[
\widetilde{\boldsymbol\delta}
=
\underset{\boldsymbol\theta\geq0}{\operatorname{argmin}}\;
(\widehat{\boldsymbol\delta}-\boldsymbol\theta)^\top
\widehat\Omega^{-1}
(\widehat{\boldsymbol\delta}-\boldsymbol\theta) .
\]
The order-restricted statistic is
\[
T=
\widehat{\boldsymbol\delta}^\top
\widehat\Omega^{-1}\widehat{\boldsymbol\delta}
-
(\widehat{\boldsymbol\delta}-\widetilde{\boldsymbol\delta})^\top
\widehat\Omega^{-1}
(\widehat{\boldsymbol\delta}-\widetilde{\boldsymbol\delta}) .
\]
It measures the improvement in covariance-weighted fit obtained by
moving from the equality-constrained null point to the closest point in
the monotone cone. By the projection identity for a convex cone, \(T\) also equals \(\widetilde{\boldsymbol\delta}^\top\widehat\Omega^{-1}\widetilde{\boldsymbol\delta}\),
so it is always non-negative. Under \(H_0\), \(T\) has an asymptotic chi-bar-square
distribution,
\[
\overline\chi^2
=
\sum_{k=0}^{R-1}w_k\,\chi_k^2,
\qquad
w_k\geq0,
\quad
\sum_{k=0}^{R-1}w_k=1,
\]
a mixture of chi-square distributions with degrees of freedom \(0\) through \(R-1\). This distribution arises in
multivariate one-sided and inequality-constrained testing (Kudo 1963; Shapiro 1985; Silvapulle and Sen 2005).
The mixture weights depend on the contrast covariance matrix and on
the geometry of the cone. With one or two order constraints (two or
three selected coefficients), weight_method = "auto" uses
analytic weights; for larger problems it uses projection simulation. The
nsim and seed arguments control this
simulation, and the estimated weights, their Monte Carlo standard
errors, and the p-value Monte Carlo standard error are returned as
weights, weight_mcse, and
p.value.mcse. Simulation currently supports at most 12
order constraints; for larger problems, precomputed weights can be
supplied through weights.
Only the coefficient ordering matters for this test; the numeric
score spacing is not used. If the observed differences already lie in
the monotone cone, \(T\) equals the
omnibus heterogeneity statistic \(W\),
but the two reference distributions—and hence the p-values—still
differ.
Worked examples
Applying the tests to an MMBCV() analysis
The following code selects the three transition-specific
Z coefficients for the sequential intermediate-state
transitions and uses the MD covariance estimator. Because these log
hazard-ratio coefficients are negative for the earlier transitions, an
“increasing” trend corresponds to effects moving toward zero.
recurrent_transitions <- c("1:2", "2:3", "3:4")
z_index <- unname(fit$cmap["Z", recurrent_transitions])
beta_Z <- fit$coefficients[z_index]
V_Z_MD <- out$varMD[z_index, z_index, drop = FALSE]
H_Z <- heterogeneity_test(beta_Z, V_Z_MD)
L_Z <- linear_trend_test(
beta_Z,
V_Z_MD,
scores = 0:2,
alternative = "greater"
)
O_Z <- order_restricted_test(
beta_Z,
V_Z_MD,
alternative = "increasing"
)
H_Z
#>
#> Omnibus Wald test of equality of transition-specific treatment effects
#>
#> data:
#> W = 2.5408, df = 2, p-value = 0.2807
#> alternative hypothesis: not all selected transition-specific treatment effects are equal
L_Z
#>
#> Generalized least-squares linear trend test
#>
#> data:
#> Z = 0.28447, p-value = 0.388
#> alternative hypothesis: true trend slope is greater than 0
O_Z
#>
#> Order-restricted test with a chi-bar-square reference distribution
#>
#> data:
#> T = 1.6863, p-value = 0.1821
#> alternative hypothesis: selected coefficients are monotonically increasing, with at least one strict inequality
The p-values and decisions at a chosen level can be collected
explicitly:
test_p_values <- c(
heterogeneity = H_Z$p.value,
linear_trend = L_Z$p.value,
order_restricted = O_Z$p.value
)
test_p_values
#> heterogeneity linear_trend order_restricted
#> 0.2807254 0.3880257 0.1821005
test_p_values < 0.05
#> heterogeneity linear_trend order_restricted
#> FALSE FALSE FALSE
The same coefficient vector can be paired with robust,
varMR, varMDMR, varFG, or any
other covariance estimator returned by MMBCV(). The
selected rows and columns must correspond to beta_Z in the
same order.
Standalone example with seven coefficients
The test functions do not require an MMBCV() model
object. The following example constructs seven correlated coefficient
estimates and a positive-definite covariance matrix with an
autoregressive (AR(1)) correlation pattern:
beta7 <- c(
beta1 = -0.31,
beta2 = -0.21,
beta3 = -0.12,
beta4 = -0.02,
beta5 = 0.09,
beta6 = 0.21,
beta7 = 0.34
)
se7 <- c(0.12, 0.11, 0.10, 0.10, 0.11, 0.12, 0.13)
position <- seq_along(beta7)
cor7 <- 0.35 ^ abs(outer(position, position, "-"))
V7 <- diag(se7) %*% cor7 %*% diag(se7)
dimnames(V7) <- list(names(beta7), names(beta7))
H7 <- heterogeneity_test(beta7, V7)
L7 <- linear_trend_test(
beta7,
V7,
scores = 0:6,
alternative = "greater"
)
O7 <- order_restricted_test(
beta7,
V7,
alternative = "increasing",
nsim = 5000,
seed = 20260713
)
H7
#>
#> Omnibus Wald test of equality of transition-specific treatment effects
#>
#> data:
#> W = 15.998, df = 6, p-value = 0.01376
#> alternative hypothesis: not all selected transition-specific treatment effects are equal
L7
#>
#> Generalized least-squares linear trend test
#>
#> data:
#> Z = 3.9882, p-value = 3.329e-05
#> alternative hypothesis: true trend slope is greater than 0
O7
#>
#> Order-restricted test with a chi-bar-square reference distribution
#>
#> data:
#> T = 15.998, p-value = 0.0007513
#> alternative hypothesis: selected coefficients are monotonically increasing, with at least one strict inequality
Additional returned quantities are available beyond the standard
printed "htest" output:
H7$contrast_matrix
#> beta1 beta2 beta3 beta4 beta5 beta6 beta7
#> beta2 - beta1 -1 1 0 0 0 0 0
#> beta3 - beta2 0 -1 1 0 0 0 0
#> beta4 - beta3 0 0 -1 1 0 0 0
#> beta5 - beta4 0 0 0 -1 1 0 0
#> beta6 - beta5 0 0 0 0 -1 1 0
#> beta7 - beta6 0 0 0 0 0 -1 1
H7$differences
#> beta2 - beta1 beta3 - beta2 beta4 - beta3 beta5 - beta4 beta6 - beta5
#> 0.10 0.09 0.10 0.11 0.12
#> beta7 - beta6
#> 0.13
L7$slope
#> [1] 0.1070206
L7$stderr
#> [1] 0.02683445
L7$scores
#> [1] 0 1 2 3 4 5 6
O7$weights
#> df0 df1 df2 df3 df4 df5 df6
#> 0.0894 0.2678 0.3330 0.2104 0.0798 0.0184 0.0012
O7$weight_mcse
#> df0 df1 df2 df3 df4 df5
#> 0.0040350375 0.0062623184 0.0066649981 0.0057642318 0.0038322829 0.0019006020
#> df6
#> 0.0004896039
O7$p.value.mcse
#> [1] 1.769086e-05
Selecting a subset and specifying unequal positions
Suppose only coefficients 1, 3, 4, and 5 are of interest. The same
index selects the coefficients and the covariance submatrix
for all three functions. For the linear-trend test, the scores preserve
the intended positions of the selected coefficients:
keep <- c(1, 3, 4, 5)
selected_scores <- c(0, 2, 3, 4)
H_subset <- heterogeneity_test(
beta7,
V7,
index = keep
)
L_subset <- linear_trend_test(
beta7,
V7,
index = keep,
scores = selected_scores,
alternative = "greater"
)
O_subset <- order_restricted_test(
beta7,
V7,
index = keep,
alternative = "increasing",
nsim = 2000,
seed = 20260714
)
c(
heterogeneity = H_subset$p.value,
linear_trend = L_subset$p.value,
order_restricted = O_subset$p.value
)
#> heterogeneity linear_trend order_restricted
#> 0.097634752 0.006052228 0.024381584
For linear_trend_test(), supplying the full score vector
0:6 would be equivalent, because the function extracts
positions 1, 3, 4, and 5 using index. In this example the
heterogeneity null is \(\beta_1=\beta_3=\beta_4=\beta_5\), and the
increasing order-restricted alternative is \(\beta_1\leq\beta_3\leq\beta_4\leq\beta_5\),
with at least one strict inequality. Neither procedure makes any
statement about the omitted coefficients, and the numeric scores affect
only the linear-trend test.