uvarpro: Unsupervised Variable Selection and Dependence Graphs using...

In varPro: Model-Independent Variable Selection via the Rule-Based Variable Priority

Unsupervised Variable Selection and Dependence Graphs using UVarPro

Description

Performs unsupervised variable selection using the Variable Priority framework. uvarpro constructs forest rules and region-release contrasts from unlabeled data. get.beta.entropy uses local lasso regressions to summarize which variables distinguish these contrasts. sdependent displays the resulting importance matrix as a variable-dependence graph and identifies candidate signal variables.

Usage

uvarpro(data,
        method = c("auto", "unsupv", "rnd"),
        ntree = 200, nodesize = NULL,
        max.rules.tree = 20, max.tree = 200,
        verbose = FALSE, seed = NULL,
        ...)

get.beta.entropy(o,
                 second.stage = FALSE,
                 pre.filter = TRUE,
                 papply = mclapply,
                 vimp.min = 0,
                 nfolds = 10,
                 maxit = 2500,
                 thresh = 1e-3,
                 parallel = FALSE,
                 use.cv = TRUE,
                 lambda.sel = c("lambda.1se", "lambda.min"),
                 nlambda = NULL,
                 lambda.min.ratio = NULL,
                 nlambda.lasso = nlambda,
                 lambda.min.ratio.lasso = lambda.min.ratio)

sdependent(I,
           threshold = .25,
           layout = "grid",
           q.signal = .75,
           directed = TRUE,
           min.degree = NULL,
           title = "s-Dependent Variable Detection",
           plot = TRUE)

Arguments

data	Data frame containing the features to be analyzed, without an outcome or class-label column. Rows with missing values are omitted, unused factor levels are dropped, and categorical features are hot-encoded. See Details.
method	Forest construction method. "auto" uses a multivariate forest autoencoder, "unsupv" uses unsupervised splitting, and "rnd" uses random splitting.
ntree	Number of trees to grow.
nodesize	Minimum terminal node size. If NULL, the requested size is n / 10 for n < 100, and max(n / 10, 20) otherwise, where n is the number of complete observations.
max.rules.tree	Maximum number of rules sampled per tree for the region-release analysis.
max.tree	Maximum number of trees used to extract rules.
verbose	Currently unused.
seed	Seed argument for uvarpro. Seed the lasso cross-validation separately; see Note.
...	Additional named forest arguments recognized by rfsrc, such as mtry or ytry, and a custom scoring function supplied as entropy; see Details. Supply lasso controls to get.beta.entropy.
o	An object returned by uvarpro, containing the processed feature data in x and the region membership lists in entropy.
second.stage	If FALSE, return the mean absolute logistic-lasso coefficient matrix. If TRUE, perform an additional local linear lasso using the predictors selected by the logistic stage, and return its mean absolute coefficients divided by the standard deviation of the released variable. See Details.
pre.filter	Should preliminary importance from uvarpro restrict the lasso analysis to variables with get.vimp(o, pretty = FALSE) > vimp.min? The filter applies to both predictors and release variables. It is skipped when preliminary importance is unavailable.
papply	An lapply-like function for processing released variables. Use lapply for serial execution. The default mclapply is resolved through the package parallel backend using getOption("mc.cores", 1L). The parallel argument controls parallelism across cross-validation folds.
vimp.min	Strict lower cutoff for preliminary importance when pre.filter = TRUE. The default retains values greater than zero. Ignored when pre.filter = FALSE.
nfolds	Number of cross-validation folds for each local logistic-lasso analysis. At least three folds are required. Ignored when use.cv = FALSE.
maxit	Maximum number of iterations passed to glmnet in both lasso stages.
thresh	Convergence tolerance passed to glmnet.
parallel	Should cv.glmnet process its folds in parallel? A suitable foreach backend must be registered. Ignored when use.cv = FALSE. Account for available resources when enabling both levels of parallelism.
use.cv	Should the logistic-lasso penalty be chosen by cross-validation? If FALSE, coefficients are taken at the last, smallest penalty on the computed glmnet path.
lambda.sel	Cross-validation penalty choice for the logistic stage. "lambda.1se" selects the largest penalty whose cross-validation error is within one standard error of the minimum; "lambda.min" selects the penalty minimizing that error. Ignored when use.cv = FALSE.
nlambda	Number of penalty values requested for the logistic stage. NULL leaves the glmnet default unchanged. The computed path can end before all requested values are reached.
lambda.min.ratio	Smallest requested logistic-stage penalty relative to the largest penalty on the path. NULL uses the glmnet default.
nlambda.lasso, lambda.min.ratio.lasso	Corresponding path controls for the optional second-stage linear lasso. They default to the logistic-stage settings and are used only when second.stage = TRUE. The second stage uses the last computed penalty.
I	A finite, nonnegative numeric importance matrix, usually obtained from get.beta.entropy, with released variables as rows and predictors as columns. Supply unique, nonempty row and column names. Row names must belong to the column names and may occur in any order. Rows are aligned by name; missing release rows are filled with zeros.
threshold	Positive edge cutoff on the scale of I. An off-diagonal entry greater than or equal to this value creates an edge. Increasing the cutoff removes edges; decreasing it can recover weaker connections.
layout	Graph layout. Supported names are "fr", "dh", "gem", "kk", "lgl", "mds", "sugiyama", "graphopt", "nicely", "random", "sphere", "grid", "circle", "star", "tree", and "bipartite"; unambiguous abbreviations are accepted. Layout-specific requirements of igraph still apply. Alternatively, supply a numeric coordinate matrix with one row per nonisolated vertex, in graph vertex order.
q.signal	Quantile of the global importance scores among nonisolated vertices used for signal designation. Values range from zero to one; the default is the upper quartile. A signal must meet both this cutoff and the degree condition.
directed	Should the graph be directed? In the directed graph, an edge goes from the released variable to the predictor used to distinguish its region-release contrasts. For an undirected display, supply a symmetric importance matrix; see Note.
min.degree	Minimum degree for candidate signal designation. With directed graphs this is the out-degree. The default is 1 for directed graphs and 2 for undirected graphs.
title	Title for the graph.
plot	Should the graph be plotted? If FALSE, graph construction and signal selection are still performed. The igraph package is required in either case.

Details

Overview and forest construction

UVarPro studies dependence among features in unlabeled data. The population framework of Zhou et al. (2026) seeks a signal set that accounts for dependence among the remaining variables. Redundant or interchangeable variables can give rise to multiple signal sets.

The analysis has three steps: uvarpro constructs local region-release comparisons, get.beta.entropy estimates a lasso coefficient importance matrix, and sdependent displays its dependence graph.

For method = "auto", copies of the processed features serve as the multivariate response and the original features serve as predictors in rfsrc. Regressing the feature vector on itself produces a forest autoencoder that partitions feature space.

For method = "unsupv", the forest uses unsupervised splitting with internally selected pseudo-responses. The default ytry is min(ceiling(sqrt(p)), p - 1), where p is the number of processed features. For method = "rnd", an independent Gaussian response is generated and splitrule = "random" is used. All three methods use the same subsequent region-release and lasso calculations.

Complete observations are hot-encoded and saved as o$x. Membership indices refer to these observations; lasso rows and columns and graph vertices refer to the processed features. Each encoded column is analyzed separately.

Rules, release regions, and local classification

A tree rule defines a region R through restrictions on feature values. Releasing variable s removes its restrictions while retaining the others. Write the enlarged region as R^{(s)} and its complementary, or near-miss, region as

C_s = R^{(s)} \setminus R.

The complementary region contains observations admitted by releasing s that still satisfy the other rule conditions.

uvarpro extracts original-region out-of-bag membership and complementary-region membership. Comparisons are retained when both sets are nonempty. The default o$entropy stores these index pairs by released variable, with comp first and oob second.

For each comparison, get.beta.entropy combines the complementary observations (class 0) and original observations (class 1). It predicts this membership label using the eligible features other than the released variable. Excluding the released variable prevents the classifier from simply recovering the rule restriction that separates the groups. The remaining predictors identify features that distinguish the local contrast. Repeating this analysis across rules allows different predictors to contribute in different parts of feature space.

Preliminary dispersion importance

The default entropy function computes a dispersion ratio. Let d(M) be the average column standard deviation of a feature matrix M. For the original and complementary feature matrices X_O and X_C, the stored rule score is

H = \frac{d(X_C) + d(X_O)}{2 d(X_C \cup X_O)}.

Here the union denotes row concatenation. The ratio compares average within-group dispersion with pooled dispersion. Before scoring, features are divided by their full-data standard deviations without centering; near-zero standard deviations are replaced by one. All features, including the released variable, enter this calculation.

importance(o) summarizes the rule scores. By default, get.beta.entropy uses preliminary importance to screen predictors and release variables before the local regressions. Set pre.filter = FALSE to include all processed features.

Logistic lasso and the importance matrix

Each local classification problem uses binomial lasso logistic regression from glmnet, with penalty parameter \alpha = 1. With use.cv = TRUE, cv.glmnet selects a penalty separately for each comparison using binomial deviance; the default choice is lambda.1se. With use.cv = FALSE, the last computed penalty on the regularization path is used.

For released variable s and retained comparison r, let \widehat\beta_{srj} denote the coefficient of predictor j. The entries of the importance matrix are mean absolute coefficients across retained comparisons. The intercept is discarded; unselected predictors and the released coordinate receive zero. Writing \mathcal K_s for the retained comparisons, the default matrix has entries

I_{sj} = \frac{1}{|\mathcal K_s|} \sum_{r \in \mathcal K_s} |\widehat\beta_{srj}|, \qquad j \ne s, \qquad I_{ss} = 0.

Rows identify released variables and columns identify predictors. Thus I[s, j] measures predictor j's contribution when variable s is released. The matrix is generally asymmetric. Column sums aggregate each predictor's contributions across release tasks; column means give the same ranking.

A comparison is retained when fitting and coefficient extraction succeed and at least one finite, nonzero predictor coefficient remains. Failed fits and intercept-only solutions are excluded. Retained comparisons receive equal weight, including zeros for unselected predictors. A release row is omitted when it has no retained comparisons, so predictor columns can outnumber rows. The result is NULL when no rows remain.

Coefficient scaling and the optional second stage

Before each logistic fit, scale(X, center = FALSE) divides each predictor column by \sqrt{\sum_i X_{ij}^2/(m-1)}, where m is the local sample size. glmnet also applies its default internal standardization and returns coefficients on the scale of the matrix supplied to it.

With second.stage = TRUE, predictors selected by the logistic lasso enter a local Gaussian lasso with the released variable as its response. This regression uses the combined complementary and original observations, the same predictor scaling, and the response from o$x. Coefficients are taken at the last computed penalty on the linear-lasso path.

Absolute second-stage coefficients are averaged over the comparisons retained by the logistic stage, with zeros for unselected predictors. A second-stage fit that fails or selects no predictors contributes zeros to this average. Row s is then divided by the full-data standard deviation of the released variable:

I^{(2)}_{sj} = \frac{1}{\widehat\sigma_s |\mathcal K_s|} \sum_{r \in \mathcal K_s} |\widehat\gamma_{srj}|.

Here \widehat\gamma_{srj} is a second-stage coefficient. A zero or nonfinite response standard deviation produces NA in that row.

get.beta.entropy returns this matrix when second.stage = TRUE. Choose the graph cutoff on the coefficient scale of the selected stage.

Adjacency matrix and interpretation of the graph

In UVarPro, s-dependence concerns which variables remain informative about the local contrast induced by releasing s. The graph summarizes their estimated contributions.

The theoretical motivation distinguishes releasing a noise variable, whose contrast is explained by associated signal variables, from releasing a signal variable, whose contrast can also involve associated noisy variables. This asymmetry motivates aggregating contributions across release tasks.

sdependent aligns rows to column names, adds zero rows for missing release variables, and then clears the diagonal. Column order determines vertex order before isolated vertices are removed. It calculates the column scores G_j = \sum_s I_{sj} and constructs the binary adjacency matrix

A_{sj} = \mathbf{1}\{I_{sj} \ge t\}, \quad s \ne j, \qquad A_{ss} = 0,

where t is threshold. Each retained connection is one edge, with its plotted width determined by the corresponding importance entry.

With directed = TRUE, A_{sj} = 1 creates the arrow s \longrightarrow j: predictor j contributes to distinguishing the comparisons formed by releasing s. Arrows describe predictive contributions. Reciprocal arrows and cycles are allowed.

The out-degree of s counts thresholded entries in its row, giving the number of predictors connected to its release tasks. The in-degree of j counts thresholded entries in its column, giving the number of release variables connected to that predictor. Global importance G_j sums the unthresholded magnitudes.

Isolated vertices are removed. Global scores are calculated before thresholding and removal, so they include subthreshold contributions and contributions from rows whose vertices are later removed. The signal-designation quantile uses the remaining vertices.

Signal designation and graphical appearance

A retained vertex enters signal.vars when its degree is at least min.degree and its global importance is at least the q.signal quantile of retained scores. For directed graphs, these conditions use out-degree and column-sum importance. The defaults require at least one outgoing edge and importance at or above the upper quartile of nonisolated vertices, including ties.

Signal vertices are blue; other retained vertices are gray. A directed edge is blue when its source is designated as signal and gray otherwise. Vertex size is 6 + 2 * log1p(imp.score). An edge with importance w has width 1 + 3w/w_{\max}, where w_{\max} is the largest plotted edge value. Widths are rescaled for each plot.

threshold controls the displayed connections; q.signal and min.degree control vertex highlighting. Changing the edge cutoff can also change the quantile reference set through vertex removal. Signal designations are exploratory.

Latent-variable example

The simulation in Examples contains four independent standard normal variables w, x, y, and z, each accompanied by five noisy copies with independent Gaussian errors of variance 0.1. Two further features are h_1 = w + x + \epsilon_{21} and h_2 = y + z + \epsilon_{22}, with independent Gaussian errors of variance 0.4. All errors are independent of the four source variables. All 26 features, including the sources, are supplied to uvarpro.

The signal set is \{w, x, y, z\}: conditioning on these four variables makes the other 22 features mutually independent. Releasing any of wi.1, ..., wi.5 has w as its s-dependent signal; the xi, yi, and zi blocks similarly point to x, y, and z. Releasing h1 has s-dependent signals w and x, while releasing h2 has y and z.

When a source such as w is released, its noisy copies and h1 can help predict local membership; x can also contribute by accounting for its component of h1. The graph can therefore contain connections in both directions. Repeated contributions across the noisy-feature releases motivate the column-sum ranking of the four signals.

The example uses pre.filter = FALSE to keep all 26 features eligible as predictors and release variables. It compares the estimated importance entries and graph connections with this known dependence structure.

Custom scoring functions

A custom entropy function accepts xC, xO, and .... The first two arguments contain the scaled complementary and original features, with the released feature first. Additional information includes compMembership, oobMembership, the processed data, and xvar.names.

Return a numeric rule score, or a list containing the score first and auxiliary information second. For subsequent use of get.beta.entropy, supply list(comp = compMembership, oob = oobMembership) as the second element, preserving that order. See Examples.

Value

uvarpro returns an object of class "uvarpro". Its principal components are:

rf

The forest used to generate rules.

results

Rule-level information with columns tree, branch, variable, n.oob, and imp. The variable index refers to xvar.names. Unevaluated comparisons have NA importance.

x

The complete-case, hot-encoded feature data, before the scaling used in scoring or local regressions.

xvar.names, xvar.org.names

Processed and original feature names, respectively.

entropy

With the default scoring function, membership pairs grouped by released variable. Each pair contains comp and oob observation indices into x. May be NULL when no usable memberships are available. Custom callbacks can change the auxiliary information stored here.

max.rules.tree, max.tree

Rule extraction limits.

family

The value "unsupv", irrespective of the forest construction method.

get.beta.entropy returns a named numeric matrix, with released variables as rows and eligible predictors as columns. Entries are mean absolute logistic-lasso coefficients when second.stage = FALSE, and response-scaled mean absolute linear-lasso coefficients when second.stage = TRUE. Rows with no retained comparisons are omitted. Returns NULL if no usable release rows remain.

sdependent invisibly returns a list with components:

signal.vars

Names of vertices meeting both the degree and importance-quantile conditions. Can be empty even when the graph is nonempty.

imp.score

Named global column-sum scores for nonisolated vertices, sorted in increasing order.

degree

Named degrees for nonisolated vertices, in graph vertex order. These are out-degrees when directed = TRUE.

If no edges remain, sdependent returns a character diagnostic. See Note.

Note

For reproducible serial lasso calculations, call set.seed immediately before get.beta.entropy(o, papply = lapply) to control its randomized cross-validation folds. Seed forest construction separately. Small local classes, degenerate predictors, and excessive pre-filtering can leave no usable comparisons. Individual lasso errors are caught and warnings suppressed.

For an undirected display, supply a symmetric importance matrix. For square input with matching row and column order, pmax(I, t(I)) retains either direction. Symmetrization can change global scores and signal designations. Undirected edge colors use the first endpoint in the internal edge list.

To admit weaker positive edges in an empty graph, decrease threshold. Keep the cutoff positive, since >= admits zero-valued entries at a zero cutoff. Check for NULL lasso output before graph construction.

Author(s)

Min Lu and Hemant Ishwaran

References

Tang F. and Ishwaran H. (2017). Random forest missing data algorithms. Statistical Analysis and Data Mining, 10:363-377.

Zhou L., Lu M. and Ishwaran H. (2026). Variable priority for unsupervised variable selection. Pattern Recognition, 172:112727. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.patcog.2025.112727")}.

See Also

varpro, importance, cv.glmnet, glmnet, scale, graph_from_adjacency_matrix

Examples

## Small forest for a quick check
set.seed(1)
o <- uvarpro(mtcars, ntree = 1)


## ------------------------------------------------------------
## Latent-variable model: known signals and s-dependence
## ------------------------------------------------------------
set.seed(21)
n <- 1000
w <- rnorm(n)
x <- rnorm(n)
y <- rnorm(n)
z <- rnorm(n)
ei <- matrix(rnorm(n * 20, sd = sqrt(.1)), ncol = 20)
e21 <- rnorm(n, sd = sqrt(.4))
e22 <- rnorm(n, sd = sqrt(.4))
wi <- w + ei[, 1:5]
xi <- x + ei[, 6:10]
yi <- y + ei[, 11:15]
zi <- z + ei[, 16:20]
h1 <- w + x + e21
h2 <- y + z + e22
dta <- data.frame(w = w, wi = wi, x = x, xi = xi,
                  y = y, yi = yi, z = z, zi = zi, h1 = h1, h2 = h2)
signal <- c("w", "x", "y", "z")

## All 26 columns are observed features. The remaining 22 are
## mutually independent conditional on the four signal variables.
o <- uvarpro(dta)
print(importance(o))

## Retain every feature for the local lasso analysis, including
## the noisy-copy and mixed-variable releases that reveal signals.
set.seed(22)
beta <- get.beta.entropy(o, pre.filter = FALSE,
                         papply = lapply, nfolds = 5)

if (!is.null(beta)) {
  ## Compare the column-sum ranking with the known signal set.
  score <- sort(colSums(beta), decreasing = TRUE)
  print(data.frame(variable = names(score),
                   importance = unname(score),
                   signal = names(score) %in% signal))

  ## Inspect one copy from each block and the two mixed variables.
  ## Population targets: each copy points to its source; h1 points
  ## to w and x, and h2 to y and z. See Details for source releases.
  release <- c("wi.1", "xi.1", "yi.1", "zi.1", "h1", "h2")
  release <- intersect(release, rownames(beta))
  print(round(beta[release, signal, drop = FALSE], 2))

  ## Inspect the corresponding adjacency block.
  ## A[s, j] = 1 gives the arrow s -> j.
  threshold <- 0.25
  A <- 1L * (beta[release, signal, drop = FALSE] >= threshold)
  print(A)

  if (all(is.finite(beta)) &&
      requireNamespace("igraph", quietly = TRUE)) {
    set.seed(23)
    ## Pass the lasso matrix directly; row alignment is automatic.
    graph.info <- sdependent(beta, threshold = threshold, layout = "fr",
                              title = "Latent model: s-dependence")
    if (is.list(graph.info)) {
      print(graph.info$signal.vars)
      print(graph.info$degree)
    } else {
      print(graph.info)
    }

    ## Return graph summaries without drawing the figure:
    ## graph.info <- sdependent(beta, threshold = threshold, plot = FALSE)
  }
}

## Optional second-stage lasso, keeping the same feature set.
## set.seed(22)
## beta.second <- get.beta.entropy(o, pre.filter = FALSE,
##                                  second.stage = TRUE,
##                                  papply = lapply, nfolds = 5)

## ------------------------------------------------------------
## Boston housing: alternative forest and categorical features
## ------------------------------------------------------------
if (requireNamespace("mlbench", quietly = TRUE)) {
  data(BostonHousing, package = "mlbench")

  ## Every column is treated as a feature; no response is specified.
  o.unsupv <- uvarpro(BostonHousing, method = "unsupv")
  print(importance(o.unsupv))

  ## Random splitting uses the same subsequent analysis.
  ## o.random <- uvarpro(BostonHousing, method = "rnd")

  Boston <- BostonHousing
  Boston$zn <- factor(Boston$zn)
  Boston$chas <- factor(Boston$chas)
  Boston$lstat <- factor(round(0.2 * Boston$lstat))
  Boston$nox <- factor(round(20 * Boston$nox))
  Boston$rm <- factor(round(Boston$rm))

  o.factor <- uvarpro(Boston)
  print(importance(o.factor))
  print(get.orgvimp(o.factor))
  print(get.orgvimp(o.factor, pretty = FALSE))
}

## ------------------------------------------------------------
## Iowa housing: larger lasso analysis
## ------------------------------------------------------------
data(housing, package = "randomForestSRC")
iowa <- roughfix(housing)
## Numeric coding below is only a speed-oriented illustration;
## retain factors to use the usual hot-encoding instead.
iowa <- data.frame(data.matrix(iowa))
o.iowa <- uvarpro(iowa, ntree = 50, max.tree = 50,
                  max.rules.tree = 5)
set.seed(22)
beta.iowa <- get.beta.entropy(o.iowa, papply = lapply, nfolds = 5)
if (!is.null(beta.iowa)) {
  print(sort(colMeans(beta.iowa), decreasing = TRUE))
  ## sdependent(beta.iowa)
}

## ------------------------------------------------------------
## Custom scoring callback: reproduces the default dispersion ratio
## ------------------------------------------------------------
my.entropy <- function(xC, xO, ...) {
  mean.sd <- function(x) mean(apply(x, 2, sd, na.rm = TRUE))
  imp <- 0.5 * (mean.sd(xC) + mean.sd(xO)) /
    mean.sd(rbind(xC, xO))
  dots <- list(...)
  list(imp = imp,
       membership = list(comp = dots$compMembership,
                         oob = dots$oobMembership))
}

if (requireNamespace("mlbench", quietly = TRUE)) {
  data(BostonHousing, package = "mlbench")
  o.custom <- uvarpro(BostonHousing, entropy = my.entropy,
                      ntree = 50, max.tree = 50, max.rules.tree = 5)
  print(importance(o.custom))
  ## you can still run the lasso if you want
  ## beta.custom <- get.beta.entropy(o.custom, papply = lapply)
}

varPro documentation built on Sept. 22, 2026, 5:09 p.m.