1.1 Motivation

In health systems, patient populations are highly heterogeneous. A patient hospitalized for heart failure has very different clinical risk factors from one admitted for a traumatic injury. Fitting a single pooled model ignores this heterogeneity, while fitting fully separate models per DRG is statistically infeasible when subgroup sample sizes are small and the number of predictors is large.

hierNest resolves this tension by decomposing each covariate’s effect into three interpretable, hierarchically nested components:

\[\beta_d = \mu + \eta^{M(d)} + \delta_d\]

where:

• \(\mu\) is the overall common effect shared across all subgroups,
• \(\eta^{M(d)}\) is the MDC-specific deviation for the Major Diagnostic Category \(M(d)\) that DRG \(d\) belongs to,
• \(\delta_d\) is the DRG-specific deviation unique to subgroup \(d\).

Pairing this reparameterization with structured regularization (lasso or overlapping group lasso) allows the model to borrow strength across related subgroups while remaining flexible enough to capture genuine heterogeneity.

1.2 Two penalization strategies

The package provides two penalization methods:

1. hierNest-Lasso — applies a standard lasso penalty to the reparameterized coefficients, with penalty weights adjusted for subgroup sample size. This is implemented via a modified design matrix that can be passed directly to glmnet. It is fast and serves as a useful baseline.

2. hierNest-OGLasso — applies an overlapping group lasso penalty that additionally enforces effect hierarchy: an MDC-specific effect can be non-zero only if the corresponding overall effect is non-zero, and a DRG-specific effect can be non-zero only if the corresponding MDC effect is non-zero. This is the recommended method when the hierarchical structure is plausible.

1.3 Package scope

This vignette covers:

• The structure of the required hier_info input.
• Simulating or loading data suitable for hierNest.
• Fitting a model with hierNest() at a fixed penalty.
• Selecting tuning parameters with cv.hierNest().
• Extracting and interpreting coefficients.
• Generating predictions on new data with predict_hierNest().
• Visualizing results with plot().

5.1 hierNest() — fit at a single penalty level

hierNest() fits the full regularization path for a single pair of hyperparameters \((\alpha_1, \alpha_2)\). It is useful for exploration or when tuning parameters have already been selected.

library(hierNest)
fit <- hierNest(
  x         = example_data$X,
  y         = example_data$Y,
  hier_info = example_data$hier_info,
  family    = "binomial",
  asparse1  = 1,    # weight for MDC-level group penalty
  asparse2  = 1     # weight for DRG-level subgroup penalty
)

# The fit contains a lambda sequence and corresponding coefficient paths
length(fit$lambda)
dim(fit$beta)   # rows = reparameterized coefficients, cols = lambda values

Key arguments:

5.2 cv.hierNest() — cross-validated tuning parameter selection

In practice, the penalty hyperparameters \(\lambda\), \(\alpha_1\), and \(\alpha_2\) need to be chosen by cross-validation. cv.hierNest() wraps the fitting procedure with cross-validation and returns the optimal parameter combination.

cv.fit <- cv.hierNest(
  x           = example_data$X,
  y           = example_data$Y,
  method      = "overlapping",
  hier_info   = example_data$hier_info,
  family      = "binomial",
  partition   = "subgroup",    # stratify CV folds within subgroups
  cvmethod    = "grid_search",
  asparse1    = c(0.5, 20),    # search alpha_1 over [0.5, 20]
  asparse2    = c(0.05, 0.20), # search alpha_2 over [0.05, 0.20]
  asparse1_num = 3,            # 3 x 3 = 9 grid points
  asparse2_num = 3,
  nlambda     = 50,            # lambda values per grid point
  pred.loss   = "ROC"          # maximize AUROC
)

cv.fit.simple <- cv.hierNest(
  x         = example_data$X,
  y         = example_data$Y,
  method    = "overlapping",
  hier_info = example_data$hier_info,
  family    = "binomial",
  partition = "subgroup",
  cvmethod  = "general",
  asparse1  = 1,
  asparse2  = 0.1,
  nlambda   = 100,
  pred.loss = "ROC"
)

cv.fit.user <- cv.hierNest(
  x         = example_data$X,
  y         = example_data$Y,
  method    = "overlapping",
  hier_info = example_data$hier_info,
  family    = "binomial",
  partition = "subgroup",
  cvmethod  = "user_supply",
  asparse1  = c(0.5, 1, 5, 10),   # explicit (alpha1, alpha2) pairs
  asparse2  = c(0.05, 0.1, 0.1, 0.2),
  nlambda   = 50
)

7.1 Coefficient heatmap

Plotting with type = "coefficients" shows the estimated overall, MDC-specific, and DRG-specific components for each selected predictor, at the cross-validation optimal \(\lambda\).

# Returns a plotly interactive heatmap
p_coef <- plot(cv.fit, type = "coefficients")
p_coef

Each row of the heatmap corresponds to one level in the hierarchy (overall mean, then each MDC, then each DRG within that MDC). Each column is a selected predictor. Blue cells indicate a positive contribution; red cells indicate negative. Rows with all zeros are suppressed.

7.2 Subgroup effect heatmap

Plotting with type = "Subgroup effects" reconstructs the total effect for each DRG subgroup: \(\hat{\beta}_d = \hat{\mu} + \hat{\eta}^{M(d)} + \hat{\delta}_d\).

p_sub <- plot(cv.fit, type = "Subgroup effects")
p_sub

This visualization is useful for comparing risk profiles across subgroups — subgroups with similar colors for a given predictor share similar risk processes for that covariate.

8.1 Hierarchical reparameterization

For participant \(i\) with DRG \(D_i = d\), the logistic regression model is:

\[\text{logit}(\Pr(Y_i = 1 \mid X_i, D_i = d)) = X_i^\top \beta_d\]

The hierNest reparameterization decomposes \(\beta_d\) as:

\[\beta_d = \mu + \eta^{M(d)} + \delta_d\]

This is encoded via a modified design matrix \(X_H\) constructed as the Khatri–Rao (column-wise Kronecker) product of \(X\) and the hierarchy indicator matrix \(H = [\mathbf{1}_n \;;\; H_M \;;\; H_D]\).

8.2 Overlapping group lasso penalty

The overlapping group lasso penalty on the grouped coefficient vector \(\theta_j = \{\mu_j, \{\eta_{Mj}\}_{M}, \{\delta_{dj}\}_d\}\) for the \(j\)-th predictor is:

\[P_{\text{og}}(\theta) = \lambda \sum_{j=1}^{p} \left[ |\theta_j|_2 + \sum_{M \in \mathcal{M}} \left(\alpha_1 |\theta_{Mj}|_2 + \alpha_2 |\eta_{Mj}| \right) \right]\]

This penalty enforces a hierarchical zero pattern: \(\mu_j = 0\) implies all \(\eta_{Mj} = 0\) and \(\delta_{dj} = 0\); \(\eta_{Mj} = 0\) implies all \(\delta_{dj} = 0\) for \(d \in M\).

8.3 Computation

For method = "overlapping", the majorization-minimization (MM) algorithm described in Algorithm 1 of the paper is used, with sequential strong rules to skip inactive predictor groups and accelerate computation. The design matrix is stored as a sparse matrix to minimise memory requirements.