Performs stability selection for Structural Equation Models. The underlying arrow selection algorithm (e.g. regularised Structural Equation Modelling) is run with different combinations of parameters controlling the sparsity (e.g. penalty parameter) and thresholds in selection proportions. These two hyper-parameters are jointly calibrated by maximisation of the stability score.
StructuralModel(
xdata,
adjacency,
residual_covariance = NULL,
Lambda = NULL,
pi_list = seq(0.01, 0.99, by = 0.01),
K = 100,
tau = 0.5,
seed = 1,
n_cat = NULL,
implementation = PenalisedLinearSystem,
resampling = "subsampling",
cpss = FALSE,
PFER_method = "MB",
PFER_thr = Inf,
FDP_thr = Inf,
Lambda_cardinal = 100,
optimisation = c("grid_search", "nloptr"),
n_cores = 1,
output_data = FALSE,
verbose = TRUE,
...
)An object of class variable_selection. A list with:
a matrix of the best stability scores for different parameters controlling the level of sparsity in the underlying algorithm.
a matrix of parameters controlling the level of sparsity in the underlying algorithm.
a matrix of the average number of selected features by the underlying algorithm with different parameters controlling the level of sparsity.
a matrix of the calibrated number of stably selected features with different parameters controlling the level of sparsity.
a matrix of calibrated thresholds in selection proportions for different parameters controlling the level of sparsity in the underlying algorithm.
a matrix of upper-bounds in PFER of calibrated stability selection models with different parameters controlling the level of sparsity.
a matrix of upper-bounds in FDP of calibrated stability selection models with different parameters controlling the level of sparsity.
a matrix of stability scores obtained with different combinations of parameters. Columns correspond to different thresholds in selection proportions.
a matrix of upper-bounds in FDP obtained with different combinations of parameters. Columns correspond to different thresholds in selection proportions.
a matrix of upper-bounds in PFER obtained with different combinations of parameters. Columns correspond to different thresholds in selection proportions.
a matrix of selection proportions.
Columns correspond to predictors from xdata.
an array
of model coefficients. Columns correspond to predictors from xdata.
Indices along the third dimension correspond to different resampling
iterations. With multivariate outcomes, indices along the fourth dimension
correspond to outcome-specific coefficients.
a list with
type="variable_selection" and values used for arguments
implementation, family, resampling, cpss and
PFER_method.
a list with values used for arguments
K, pi_list, tau, n_cat, pk, n
(number of observations), PFER_thr, FDP_thr and seed.
The datasets xdata and ydata are also included if
output_data=TRUE.
For all matrices and arrays returned, the rows are ordered in the same way and correspond to parameter values stored in
Lambda.
matrix with observations as rows and variables as columns.
Column names must be defined and in line with the row and column names of
adjacency.
binary adjacency matrix of the Directed Acyclic Graph (transpose of the asymmetric matrix A in Reticular Action Model notation). The row and column names of this matrix must be defined.
binary and symmetric matrix encoding the nonzero entries in the residual covariance matrix (symmetric matrix S in Reticular Action Model notation). By default, this is the identity matrix (no residual covariance).
matrix of parameters controlling the level of sparsity in the
underlying feature selection algorithm specified in implementation.
vector of thresholds in selection proportions. If
n_cat=NULL or n_cat=2, these values must be >0 and
<1. If n_cat=3, these values must be >0.5 and
<1.
number of resampling iterations.
subsample size. Only used if resampling="subsampling" and
cpss=FALSE.
value of the seed to initialise the random number generator and
ensure reproducibility of the results (see set.seed).
computation options for the stability score. Default is
NULL to use the score based on a z test. Other possible values are 2
or 3 to use the score based on the negative log-likelihood.
function to use for variable selection.
resampling approach. Possible values are:
"subsampling" for sampling without replacement of a proportion
tau of the observations, or "bootstrap" for sampling with
replacement generating a resampled dataset with as many observations as in
the full sample. Alternatively, this argument can be a function to use for
resampling. This function must use arguments named data and
tau and return the IDs of observations to be included in the
resampled dataset.
logical indicating if complementary pair stability selection
should be done. For this, the algorithm is applied on two non-overlapping
subsets of half of the observations. A feature is considered as selected if
it is selected for both subsamples. With this method, the data is split
K/2 times (K models are fitted). Only used if
PFER_method="MB".
method used to compute the upper-bound of the expected
number of False Positives (or Per Family Error Rate, PFER). If
PFER_method="MB", the method proposed by Meinshausen and Bühlmann
(2010) is used. If PFER_method="SS", the method proposed by Shah and
Samworth (2013) under the assumption of unimodality is used.
threshold in PFER for constrained calibration by error
control. If PFER_thr=Inf and FDP_thr=Inf, unconstrained
calibration is used (the default).
threshold in the expected proportion of falsely selected
features (or False Discovery Proportion) for constrained calibration by
error control. If PFER_thr=Inf and FDP_thr=Inf, unconstrained
calibration is used (the default).
number of values in the grid of parameters controlling
the level of sparsity in the underlying algorithm. Only used if
Lambda=NULL.
character string indicating the type of optimisation
method. With optimisation="grid_search" (the default), all values in
Lambda are visited. Alternatively, optimisation algorithms
implemented in nloptr can be used with
optimisation="nloptr". By default, we use
"algorithm"="NLOPT_GN_DIRECT_L", "xtol_abs"=0.1,
"ftol_abs"=0.1 and "maxeval"=Lambda_cardinal. These values
can be changed by providing the argument opts (see
nloptr). For stability selection using penalised
regression, optimisation="grid_search" may be faster as it allows
for warm start.
number of cores to use for parallel computing (see argument
workers in multisession). Using
n_cores>1 is only supported with optimisation="grid_search".
logical indicating if the input datasets xdata and
ydata should be included in the output.
logical indicating if a loading bar and messages should be printed.
additional parameters passed to the functions provided in
implementation or resampling.
In stability selection, a feature selection algorithm is fitted on
K subsamples (or bootstrap samples) of the data with different
parameters controlling the sparsity (Lambda). For a given (set of)
sparsity parameter(s), the proportion out of the K models in which
each feature is selected is calculated. Features with selection proportions
above a threshold pi are considered stably selected. The stability
selection model is controlled by the sparsity parameter(s) for the
underlying algorithm, and the threshold in selection proportion:
\(V_{\lambda, \pi} = \{ j: p_{\lambda}(j) \ge \pi \} \)
In Structural Equation Modelling, "feature" refers to an arrow in the corresponding Directed Acyclic Graph.
These parameters can be calibrated by maximisation of a stability score
(see ConsensusScore if n_cat=NULL or
StabilityScore otherwise) calculated under the null
hypothesis of equiprobability of selection.
It is strongly recommended to examine the calibration plot carefully to
check that the grids of parameters Lambda and pi_list do not
restrict the calibration to a region that would not include the global
maximum (see CalibrationPlot). In particular, the grid
Lambda may need to be extended when the maximum stability is
observed on the left or right edges of the calibration heatmap. In some
instances, multiple peaks of stability score can be observed. Simulation
studies suggest that the peak corresponding to the largest number of
selected features tend to give better selection performances. This is not
necessarily the highest peak (which is automatically retained by the
functions in this package). The user can decide to manually choose another
peak.
To control the expected number of False Positives (Per Family Error Rate)
in the results, a threshold PFER_thr can be specified. The
optimisation problem is then constrained to sets of parameters that
generate models with an upper-bound in PFER below PFER_thr (see
Meinshausen and Bühlmann (2010) and Shah and Samworth (2013)).
Possible resampling procedures include defining (i) K subsamples of
a proportion tau of the observations, (ii) K bootstrap
samples with the full sample size (obtained with replacement), and (iii)
K/2 splits of the data in half for complementary pair stability
selection (see arguments resampling and cpss). In
complementary pair stability selection, a feature is considered selected at
a given resampling iteration if it is selected in the two complementary
subsamples.
To ensure reproducibility of the results, the starting number of the random
number generator is set to seed.
For parallelisation, stability selection with different sets of parameters
can be run on n_cores cores. Using n_cores > 1 creates a
multisession. Alternatively,
the function can be run manually with different seeds and all other
parameters equal. The results can then be combined using
Combine.
ourstabilityselectionsharp
stabilityselectionMBsharp
stabilityselectionSSsharp
RegSEMsharp
SelectionAlgo,
Resample, StabilityScore
Other stability functions:
BiSelection(),
Clustering(),
GraphicalModel(),
VariableSelection()
oldpar <- par(no.readonly = TRUE)
par(mar = rep(7, 4))
# \donttest{
# Data simulation
set.seed(1)
pk <- c(3, 2, 3)
simul <- SimulateStructural(
n = 500,
pk = pk,
nu_between = 0.5,
v_between = 1,
v_sign = 1
)
# Stability selection (using glmnet)
dag <- LayeredDAG(layers = pk)
stab <- StructuralModel(
xdata = simul$data,
adjacency = dag
)
CalibrationPlot(stab)
LinearSystemMatrix(vect = Stable(stab), adjacency = dag)
# Stability selection (using OpenMx)
if (requireNamespace("OpenMx", quietly = TRUE)) {
stab <- StructuralModel(
xdata = simul$data,
implementation = PenalisedOpenMx,
Lambda = seq(50, 500, by = 50),
adjacency = dag
)
CalibrationPlot(stab)
OpenMxMatrix(SelectedVariables(stab), adjacency = dag)
}
# }
if (FALSE) {
# Data simulation with latent variables
set.seed(1)
pk <- c(3, 2, 3)
simul <- SimulateStructural(
n = 500,
pk = pk,
nu_between = 0.5,
v_sign = 1,
v_between = 1,
n_manifest = 3,
ev_manifest = 0.95
)
# Stability selection (using OpenMx)
if (requireNamespace("OpenMx", quietly = TRUE)) {
dag <- LayeredDAG(layers = pk, n_manifest = 3)
penalised <- dag
penalised[, seq_len(ncol(simul$data))] <- 0
stab <- StructuralModel(
xdata = simul$data,
implementation = PenalisedOpenMx,
adjacency = dag,
penalised = penalised,
Lambda = seq(10, 100, by = 20),
K = 10 # to increase for real use
)
CalibrationPlot(stab)
ids_latent <- grep("f", colnames(dag))
OpenMxMatrix(SelectedVariables(stab),
adjacency = dag
)[ids_latent, ids_latent]
}
}
par(oldpar)
Run the code above in your browser using DataLab