--- title: "mmbcv: Bias-corrected sandwich variance for clustered multistate Cox models" output: rmarkdown::html_vignette: css: vignette.css bibliography: references.bib link-citations: true vignette: > %\VignetteIndexEntry{mmbcv: Bias-corrected sandwich variance for clustered multistate Cox models} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE ) has_survival <- requireNamespace("survival", quietly = TRUE) ``` # Overview The **mmbcv** package computes cluster-robust and bias-corrected sandwich variance estimators for multistate Cox models fitted in counting-process form. Wang, Turner, and Li [@wang2023] developed a class of finite-sample bias corrections for clustered marginal Cox models; **mmbcv** extends those corrections to the multistate setting. The main user-facing function is **`MMBCV()`**. It takes - a fitted multistate Cox model object (`fit`), - the original counting-process dataset (`data`), - the start and stop time variables, - the cluster identifier, - the subject identifier, - the event/state identifier, and returns a collection of robust and bias-corrected variance estimators for the fitted regression coefficients. The package also provides three standalone procedures for testing patterns in a set of correlated coefficient estimates: `heterogeneity_test()` for any departure from equality, `linear_trend_test()` for a score-based linear trend, and `order_restricted_test()` for a monotone (increasing or decreasing) alternative. Each procedure accepts a coefficient vector together with its full variance-covariance matrix, including off-diagonal covariances. For background on fitting multistate Cox models with `survival::coxph()` using list-formulas, see the survival vignette *Multi-state models and competing risks*: - `vignette("compete", package = "survival")` - # Cluster-robust and bias-corrected variance estimators ## Technical details ### Notation Suppose the data arise from a clustered multistate survival process. Let \(\beta\) denote the full regression coefficient vector, of length \(p\), and let \(\hat\beta\) be its partial-likelihood estimate. For clusters \(i = 1,\ldots,n\), let \(U_i(\hat\beta)\) denote the cluster-level score contribution (the sum of the score residuals of the cluster's members), and let \(V_m\) denote the model-based variance estimator (the naive, non-robust covariance of \(\hat\beta\)). \(V_m\) plays the role of the "bread" and appears on both sides of every sandwich below. Following Wang et al. [@wang2023], the baseline cluster-robust sandwich variance estimator is \[ \widehat{\mathrm{Var}}_{\mathrm{robust}}(\hat\beta) = V_m \left( \sum_{i=1}^n U_i(\hat\beta)\, U_i(\hat\beta)^\top \right) V_m . \] In **mmbcv** this estimator is returned as `robust`. When the number of clusters is small, this sandwich estimator tends to underestimate the true variance, which motivates the finite-sample corrections that follow. ### A common sandwich form Many of the corrected estimators share a single template that differs only in the cluster score and/or in a cluster-specific correction matrix: \[ \widehat{\mathrm{Var}}(\hat\beta) = V_m \left( \sum_{i=1}^n C_i \, U_i^\star(\hat\beta)\, U_i^\star(\hat\beta)^\top C_i^\top \right) V_m , \] where \(U_i^\star(\hat\beta)\) is either the original cluster score \(U_i(\hat\beta)\) or the martingale-residual-corrected score \(U_i^{\mathrm{MR}}(\hat\beta)\) introduced below, and \(C_i\) is a cluster-specific correction matrix. The correction matrices are built from the *cluster leverage* \(\Omega_i V_m\), where \(\Omega_i\) is a cluster-specific gradient matrix: the negative derivative of cluster \(i\)'s martingale score with respect to \(\beta\). It is the multistate analogue of the cluster gradient matrix in Wang et al. [@wang2023], obtained by summing the negative score-derivative across the transitions of the multistate model. Because it derives from the martingale-based cluster score rather than the ordinary partial-likelihood score, it is a cluster-level gradient rather than a per-cluster information matrix, and its explicit form which also depends on the tie-handling method [@wang2023]. The product \(\Omega_i V_m\), which plays a role analogous to a cluster leverage [@meng2023], enters each correction matrix \(C_i\). Throughout, \(I_p\) denotes the \(p \times p\) identity matrix. ### Martingale-residual (MR) correction The martingale-residual correction replaces the raw cluster score \(U_i(\hat\beta)\) with a bias-adjusted score \(U_i^{\mathrm{MR}}(\hat\beta)\) constructed from the cluster's martingale residuals; the adjustment is designed to reduce the finite-sample bias of the score cross-products. Its exact construction is given in Wang et al. [@wang2023]. The resulting estimator is \[ \widehat{\mathrm{Var}}_{\mathrm{MR}}(\hat\beta) = V_m \left( \sum_{i=1}^n U_i^{\mathrm{MR}}(\hat\beta)\, U_i^{\mathrm{MR}}(\hat\beta)^\top \right) V_m , \] returned by **mmbcv** as `varMR`. ### Multiplicative corrections: KC, FG, and MD Wang et al. [@wang2023] summarize three multiplicative small-sample corrections that use the common sandwich form with different choices of \(C_i\). The **Kauermann-Carroll (KC)** correction [@kauermann2001; @wang2023] uses \[ C^{KC}_i = \left(I_p - \Omega_i V_m \right)^{-1/2}, \] giving \[ \widehat{\mathrm{Var}}_{\mathrm{KC}}(\hat\beta) = V_m \left( \sum_{i=1}^n C^{KC}_i\, U_i(\hat\beta)\, U_i(\hat\beta)^\top C^{KC \top}_i \right) V_m , \] returned as `varKC`. The **Fay-Graubard (FG)** correction [@fay2001; @wang2023] uses \[ C^{FG}_i = \mathrm{diag} \left\{ \left(1 - \min\{\xi,\ [\Omega_i V_m]_{jj}\}\right)^{-1/2} \right\}_{j=1,\ldots,p}, \] where \([\Omega_i V_m]_{jj}\) is the \(j\)th diagonal element of the cluster leverage and \(\xi = 0.75\) is a truncation constant that keeps the correction finite when a leverage value approaches one. The resulting estimator is \[ \widehat{\mathrm{Var}}_{\mathrm{FG}}(\hat\beta) = V_m \left( \sum_{i=1}^n C^{FG}_i\, U_i(\hat\beta)\, U_i(\hat\beta)^\top C^{FG \top}_i \right) V_m , \] returned as `varFG`. The **Mancl-DeRouen (MD)** correction [@mancl2001; @wang2023] uses \[ C^{MD}_i = \left(I_p - \Omega_i V_m \right)^{-1}, \] giving \[ \widehat{\mathrm{Var}}_{\mathrm{MD}}(\hat\beta) = V_m \left( \sum_{i=1}^n C^{MD}_i\, U_i(\hat\beta)\, U_i(\hat\beta)^\top C^{MD \top}_i \right) V_m , \] returned as `varMD`. The KC and MD corrections apply the same leverage term \(\Omega_i V_m\) with exponents \(-1/2\) and \(-1\), respectively, so MD applies the stronger inflation of the two. ### Additive correction: MBN The **Morel-Bokossa-Neerchal (MBN)** correction [@morel2003; @wang2023] takes a different form: it is a linear combination of the robust sandwich variance and the model-based variance, which guards against underestimation when the number of clusters is small, \[ \widehat{\mathrm{Var}}_{\mathrm{MBN}}(\hat\beta) = c_1\, \widehat{\mathrm{Var}}_{\mathrm{robust}}(\hat\beta) + c_2\, \hat{\phi}\, V_m , \] where the scalars \(c_1\), \(c_2\), and \(\hat{\phi}\) are defined as in Wang et al. [@wang2023]. This is returned as `varMBN`. ### Hybrid MR estimators Wang et al. [@wang2023] further combine the MR correction with each of the four corrections above. The KC, FG, and MD hybrids replace \(U_i(\hat\beta)\) with \(U_i^{\mathrm{MR}}(\hat\beta)\) in the corresponding sandwich, and the MBN hybrid replaces the robust variance with the MR variance (and \(\hat\phi\) with its MR analogue \(\hat\phi_{\mathrm{MR}}\)): \[ \widehat{\mathrm{Var}}_{\mathrm{KCMR}}(\hat\beta) = V_m \left( \sum_{i=1}^n C^{KC}_i\, U_i^{\mathrm{MR}}(\hat\beta)\, U_i^{\mathrm{MR}}(\hat\beta)^\top C_i^{KC \top} \right) V_m , \] \[ \widehat{\mathrm{Var}}_{\mathrm{FGMR}}(\hat\beta) = V_m \left( \sum_{i=1}^n C^{FG}_i\, U_i^{\mathrm{MR}}(\hat\beta)\, U_i^{\mathrm{MR}}(\hat\beta)^\top C_i^{FG \top} \right) V_m , \] \[ \widehat{\mathrm{Var}}_{\mathrm{MDMR}}(\hat\beta) = V_m \left( \sum_{i=1}^n C^{MD}_i\, U_i^{\mathrm{MR}}(\hat\beta)\, U_i^{\mathrm{MR}}(\hat\beta)^\top C_i^{MD \top} \right) V_m , \] and \[ \widehat{\mathrm{Var}}_{\mathrm{MBNMR}}(\hat\beta) = c_1\, \widehat{\mathrm{Var}}_{\mathrm{MR}}(\hat\beta) + c_2\, \hat{\phi}_{\mathrm{MR}}\, V_m . \] These are returned as `varKCMR`, `varFGMR`, `varMDMR`, and `varMBNMR`, respectively. Every quantity returned by `MMBCV()` is a variance-covariance matrix for \(\hat\beta\). Standard errors from any such estimator \(V\) are obtained from its diagonal, \[ \mathrm{SE}(\hat\beta_j) = \sqrt{V_{jj}} . \] ## Main function and arguments Once a multistate Cox model has been fitted with `survival::coxph()`, all of the covariance estimators described above are obtained from a single function, `MMBCV()`. Its call is ```r MMBCV(fit, data, StartTime, StopTime, ClusterID, SubjectID, Event, tie = c("breslow", "efron"), details = FALSE, detail_names = NULL, kc_tol = 1e-10) ``` with arguments summarized below; see `?MMBCV` for the full reference. | Argument | Description | Default | |:--|:--|:--| | `fit` | Fitted multistate Cox model object, with components `cmap`, `rmap`, `states`, and `coefficients`; typically produced by a `survival` workflow. | — | | `data` | Data frame in counting-process form, matching the data used to fit `fit`. | — | | `StartTime`, `StopTime` | Unquoted column names in `data` giving the interval start and stop times. | — | | `ClusterID` | Unquoted column name in `data` identifying the clustering unit. | — | | `SubjectID` | Unquoted column name in `data` identifying the subject. | — | | `Event` | Unquoted column name in `data` giving the event/state indicator used to define transition-specific events. | — | | `tie` | Tie-handling method used in the variance calculation, `"breslow"` or `"efron"`. | `"breslow"` | | `details` | Logical; whether to return additional intermediate quantities beyond the variance-covariance matrices. | `FALSE` | | `detail_names` | Optional character vector selecting which intermediate components to append when `details = TRUE`. | `NULL` | | `kc_tol` | Numeric tolerance for deciding whether the cluster-specific KC correction matrix is non-singular enough to admit a real inverse square root. | `1e-10` | `fit` is the fitted multistate Cox model object. `MMBCV()` uses its multistate mapping components (`cmap`, `rmap`, `states`) and the fitted coefficient vector to reconstruct the transition-specific design and to aggregate the cluster-level score contributions; for this reason `data` must be the same counting-process data set used in `coxph()`, containing the interval start/stop times, the subject and cluster identifiers, and the event/state variable. `StartTime` and `StopTime` give the interval endpoints; `ClusterID` and `SubjectID` identify the clustering unit and the subject; and `Event` gives the state occupied at `StopTime`, from which the transition-specific event indicators are formed. Internally, `MMBCV()` partitions the data by transition, builds the corresponding event indicators, and computes the requested sandwich estimators by aggregating the resulting cluster-level scores. `tie` selects the tie-handling approximation (Breslow or Efron) used in the variance calculation and should match the tie method used to fit the model. `details` controls whether intermediate quantities—such as the model-based covariance matrix and the cluster-specific correction objects—are returned alongside the variance-covariance matrices; when `details = TRUE`, `detail_names` selects specific components, and otherwise a default subset is returned. `kc_tol` sets the tolerance for the eigenvalue check that decides whether the KC correction matrix \(C^{KC}_i\) admits a stable real inverse square root. `MMBCV()` returns a named list of variance-covariance matrices for the full coefficient vector: `robust` for the cluster-robust sandwich estimator, together with the bias-corrected estimators `varMR`, `varMD`, `varMDMR`, `varFG`, `varFGMR`, `varKC`, `varKCMR`, `varMBN`, and `varMBNMR`. Their individual meanings, and how to convert them into standard errors, are given in *Interpreting the output* below. ## Worked example ### Example data The package ships with a small simulated clustered multistate dataset, `msdat3`. ```{r} library(mmbcv) data("msdat3") dim(msdat3) length(unique(msdat3$id)) length(unique(msdat3$clus_id)) table(msdat3$event) ``` `msdat3` is in long (counting-process) format, with one row per subject-interval: - **`id`**: subject id - **`clus_id`**: cluster id - **`Z`**: cluster-level treatment assignment (0/1) - **`X`**: individual-level baseline covariate - **`from`, `to`**: start and end states for the interval - **`Tstart`, `Tstop`**: interval start and stop times - **`event`**: end state at `Tstop` (e.g., `censor`, `S1`, `S2`, `S3`, `D`) The state labels are: - **`(s0)`**: initial state - **`S1`**: first intermediate state - **`S2`**: second intermediate state - **`S3`**: third intermediate state - **`D`**: absorbing state - **`censor`**: right-censoring at the end of the interval A quick look at the observed transitions: ```{r} with(msdat3, table(from, to)) ``` ### Fitting a multistate Cox model with survival `MMBCV()` takes a fitted model object as input. In a typical workflow this model is fit with `survival::coxph()` using list-formulas. ```{r} library(survival) fit <- coxph( list( Surv(Tstart, Tstop, event) ~ 1, state("(s0)"):state("S1") + state("S1"):state("S2") + state("S2"):state("S3") ~ Z + X, state("(s0)"):state("D") + state("S1"):state("D") + state("S2"):state("D") + state("S3"):state("D") ~ (Z + X)/common ), data = msdat3, id = id, ties = "breslow", timefix = FALSE ) fit ``` If the **survival** package is not installed, install it in an interactive R session with `install.packages("survival")`. ### Computing variance estimators Once the multistate Cox model is fitted, pass it (with the original dataset) to `MMBCV()`: ```{r} out <- MMBCV( fit, msdat3, StartTime = Tstart, StopTime = Tstop, ClusterID = clus_id, SubjectID = id, Event = event, tie = "breslow", details = FALSE ) names(out) ``` ### Interpreting the output `MMBCV()` returns a list of variance-covariance matrices for the fitted regression coefficients. Each matrix can be used to compute standard errors, Wald tests, and confidence intervals. The elements are: - `robust`: cluster-robust sandwich variance estimator - `varMR`: martingale-residual (MR) bias-corrected estimator - `varMD`: Mancl-DeRouen (MD) bias-corrected estimator - `varMDMR`: hybrid MD + MR estimator - `varKC`: Kauermann-Carroll (KC) bias-corrected estimator - `varKCMR`: hybrid KC + MR estimator - `varFG`: Fay-Graubard (FG) bias-corrected estimator - `varFGMR`: hybrid FG + MR estimator - `varMBN`: Morel-Bokossa-Neerchal (MBN) bias-corrected estimator - `varMBNMR`: hybrid MBN + MR estimator Standard errors are obtained from any of these matrices with `sqrt(diag(...))`: ```{r} Vlist <- out[c("robust","varMR","varMD","varMDMR","varFG","varFGMR","varKC","varKCMR","varMBN","varMBNMR")] SE <- sapply(Vlist, function(V) sqrt(diag(V))) SE ``` ### Efron ties If the model was fitted in `survival::coxph()` with Efron's tie handling, use Efron ties in the variance calculation as well, so that the two steps are consistent: ```{r} out_efron <- MMBCV( fit, msdat3, StartTime = Tstart, StopTime = Tstop, ClusterID = clus_id, SubjectID = id, Event = event, tie = "efron", details = FALSE ) sqrt(diag(out_efron$robust)) ``` # Tests for correlated coefficient patterns ## 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: ```r heterogeneity_test(beta, vcov, index = NULL) ``` ```r linear_trend_test(beta, vcov, index = NULL, scores = NULL, alternative = c("greater", "less", "two.sided")) ``` ```r 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) ``` ### Common inputs and output Let \[ \widehat{\boldsymbol\beta} = (\widehat\beta_1,\ldots,\widehat\beta_R)^\top \] denote the selected coefficient estimates, and let \(\widehat V\) denote their variance-covariance matrix. In every test, the order of the rows and columns of \(\widehat V\) must match the order of the estimates, and the off-diagonal covariances must be retained, because all three statistics depend on the joint covariance structure. The optional `index` argument selects coefficients from a larger `beta` vector; the same selection is applied automatically to `vcov`. Integer, logical, and character indices are all accepted. For the linear-trend and order-restricted tests, the resolved order defines the sequence being tested. The heterogeneity statistic itself does not require a natural ordering, although the order determines how its returned adjacent differences are displayed. Each function returns an object of class `"htest"`. Printing the object shows the primary statistic and p-value; additional quantities—such as the selected coefficients, the covariance matrix, the contrasts, the fitted trend, the cone projection, and the chi-bar-square weights—can be accessed with `$`. The functions do not take a significance level; a decision at a user-selected level \(\alpha\) is obtained by comparing `result$p.value < alpha`. ### 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 [@wald1943; @coxHinkley1974; @lehmannRomano2022]. 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<\cdots0\); - `"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 [@robertson1988; @silvapulle2005]. 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 [@kudo1963; @shapiro1985; @silvapulle2005]. 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. ```{r test-model-inputs} 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] ``` ```{r test-model-examples} 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 L_Z O_Z ``` The p-values and decisions at a chosen level can be collected explicitly: ```{r test-model-decisions} test_p_values <- c( heterogeneity = H_Z$p.value, linear_trend = L_Z$p.value, order_restricted = O_Z$p.value ) test_p_values test_p_values < 0.05 ``` 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: ```{r test-seven-coefficients} 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 L7 O7 ``` Additional returned quantities are available beyond the standard printed `"htest"` output: ```{r test-seven-components} H7$contrast_matrix H7$differences L7$slope L7$stderr L7$scores O7$weights O7$weight_mcse O7$p.value.mcse ``` ### 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: ```{r test-subset} 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 ) ``` 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. # Practical remarks ## Choosing a variance estimator The bias-corrected estimators are most useful when the number of clusters is not large. Wang et al. (2023) report that, in their simulations for clustered time-to-event outcomes, the MD estimator performed well when cluster-size variation was modest, while KCMR performed best under larger cluster-size variation; they also note that MR alone could be unstable in some scenarios. It is therefore useful to compare `robust`, `varMD`, and `varKCMR` in practice. ## Choosing among the coefficient-pattern tests The three tests answer related but different questions: - use `heterogeneity_test()` when any departure from equality is scientifically relevant and no direction or shape should be imposed; - use `linear_trend_test()` when change with respect to a meaningful numeric score is expected to be approximately linear; - use `order_restricted_test()` when the direction is prespecified and a monotone, possibly nonlinear, pattern is scientifically plausible. The tests are not interchangeable. The choice should follow from the scientific hypothesis, not from which function produces the smallest p-value. ## Reporting For any of the three coefficient-pattern tests, the covariance estimator is part of the analysis and can materially affect the statistic and the p-value, so the chosen estimator should be stated. When the order-restricted weights are estimated by simulation, set a reproducible `seed`, use an `nsim` appropriate for the desired precision, and examine `p.value.mcse` when interpreting a p-value close to the chosen significance threshold. # References