Package {pecme}

pecme: Penalized ECME Estimation for Censored Linear Mixed Models

Description

Fits Gaussian random-intercept linear mixed models with a left-, right-, and/or interval-censored (and optionally missing) response via the ECME algorithm, with simultaneous penalized variable selection (Lasso, Adaptive Lasso, SCAD, MCP, Elastic Net, Ridge, or no penalty). See [pecme()] for a single fit and [grid_pecme()] for penalty-parameter selection.

Model

Y^*_{ij} = x_{ij}'\beta + b_j + e_{ij}, \quad e_{ij}\sim N(0,\sigma^2), \quad b_j \sim N(0,\sigma_b^2)

with Y^* observed exactly, left-censored, right-censored, interval-censored, or completely missing, depending on the row.

Algorithm

Liu and Rubin's (1994) ECME algorithm: an E-step (Gauss-Hermite quadrature for the random intercept – a numerical approximation to the continuous integral, not an exact evaluation of it; see [pecme_quadrature_check()] – with closed-form truncated-normal moments for censored/missing rows) alternates with two conditional-maximization steps – a penalized coordinate-descent CM-step for \beta that maximizes the expected complete-data objective, and an ECME CM-step for (\sigma,\sigma_b) that maximizes the actual observed-data marginal likelihood directly.

Author(s)

Maintainer: Ali Asghar Haeri-Mehrizi haeri.stat@gmail.com

Authors:

• Ali Asghar Haeri-Mehrizi haeri.stat@gmail.com

• Arshia Haeri-Mehrizi haeri.arshia@gmail.com

Other contributors:

• Adel Mohammadpour adel.mohammadpour@ucalgary.ca [reviewer]

Stable log(Phi(b) - Phi(a)) for possibly-vectorized a < b

Description

Stable log(Phi(b) - Phi(a)) for possibly-vectorized a < b

Usage

.pecme_interval_logprob(a, b)

log(exp(log_a) - exp(log_b)) for log_a > log_b, stable near 0

Description

log(exp(log_a) - exp(log_b)) for log_a > log_b, stable near 0

Usage

.pecme_logspace_sub(log_a, log_b)

Gauss-Hermite nodes and weights for N(0, 1) integration

Description

Returns nodes z and weights w such that, for a smooth function f, E[f(Z)] \approx \sum_k w_k f(z_k) when Z \sim N(0,1). Results are cached per quadrature order because the Golub-Welsch eigendecomposition is by far the most expensive part of the E-step if it is repeated every ECME iteration.

Usage

gauss_hermite_normal(n = 40L)

Arguments

Value

A list with nodes and weights, both of length n, ordered by increasing node value, calibrated so that they directly integrate against the standard normal density (i.e. the physicists' Hermite rule has already been rescaled).

Grid search over the penalty parameter lambda

Description

Grid search over the penalty parameter lambda

Usage

grid_pecme(
  fixed_formula,
  random_var,
  data,
  censor = "none",
  left = NULL,
  right = NULL,
  penalty = c("lasso", "adaptive", "scad", "mcp", "elastic", "ridge", "none"),
  lambda = NULL,
  nlambda = 50L,
  lambda_min_ratio = NULL,
  alpha = 0.5,
  gamma = NULL,
  adaptive_gamma = 1,
  criterion = c("bic", "ebic", "aic", "gcv", "cv"),
  ebic_gamma = 0.5,
  cv_folds = 5L,
  cv_seed = NULL,
  missing_y = FALSE,
  missing_x = c("error", "cca", "mean", "reg"),
  standardize = TRUE,
  quadrature = 40L,
  workers = 1L,
  parallel_over = c("lambda", "groups"),
  variance_step = c("ecme", "em"),
  refit_at_best = TRUE,
  init_method = c("lmer", "zero", "random"),
  max_iter = 200L,
  tol = 1e-06,
  cd_max_iter = 1000L,
  cd_tol = 1e-08,
  keep_path = TRUE,
  verbose = FALSE
)

Arguments

Value

An object of class '"pecme_grid"': a list with '$best' (the selected '"pecme"' fit, refit per 'refit_at_best'), '$lambda', '$scores' (a data frame of one row per lambda with every criterion, 'df', 'logLik', and convergence status), '$best_index', and '$criterion'.

Numerically stable log(sum(exp(x)))

Description

Numerically stable log(sum(exp(x)))

Usage

logsumexp(x)

Fit a penalized ECME censored/missing-response linear mixed model

Description

Fits Y^*_{ij} = x_{ij}'\beta + b_j + e_{ij}, e_{ij} \sim N(0,\sigma^2), b_j \sim N(0,\sigma_b^2), where the response Y^* may be left-, right-, or interval-censored (and/or missing), by the ECME algorithm of Liu and Rubin (1994). The random intercept is integrated out numerically by fixed-order Gauss-Hermite quadrature (a numerical approximation to the continuous integral, not an exact evaluation of it; see [pecme_quadrature_check()] for a sensitivity diagnostic and 'quadrature' below). \beta is estimated by a penalized CM-step (coordinate descent); (\sigma,\sigma_b) are estimated, by default, with a true ECME CM-step that directly maximizes the observed-data marginal likelihood (not merely its EM lower bound), following the strategy of Vaida and Liu (2009) for censored linear mixed models.

Usage

pecme(
  fixed_formula,
  random_var,
  data,
  censor = "none",
  left = NULL,
  right = NULL,
  penalty = c("none", "lasso", "adaptive", "scad", "mcp", "elastic", "ridge"),
  lambda = 0,
  alpha = 0.5,
  gamma = NULL,
  adaptive_gamma = 1,
  missing_y = FALSE,
  missing_x = c("error", "cca", "mean", "reg"),
  standardize = TRUE,
  quadrature = 40L,
  workers = 1L,
  variance_step = c("ecme", "em"),
  refit = TRUE,
  init_method = c("lmer", "zero", "random"),
  max_iter = 200L,
  tol = 1e-06,
  cd_max_iter = 1000L,
  cd_tol = 1e-08,
  verbose = FALSE
)

Arguments

Value

An object of class '"pecme"'.

References

Liu, C. and Rubin, D. B. (1994). The ECME algorithm. *Biometrika*, 81(4), 633-648.

Vaida, F. and Liu, L. (2009). Fast implementation for normal mixed effects models with censored response. *JCGS*, 18(1), 151-162.

Fan, J. and Li, R. (2001). Variable selection via nonconcave penalized likelihood and its oracle properties. *JASA*, 96(456), 1348-1360.

Zhang, C.-H. (2010). Nearly unbiased variable selection under minimax concave penalty. *Annals of Statistics*, 38(2), 894-942.

See Also

[grid_pecme()], [pecme_mi_combine()]

Cluster (group-level) bootstrap for a pecme fit

Description

Resamples whole groups with replacement, refits the identical model (same formula, 'penalty', 'lambda', 'alpha', 'gamma') to each resampled data set, and summarizes the resulting coefficient distribution – a valid source of standard errors and confidence intervals for a penalized fit, where naive SEs do not apply.

Usage

pecme_bootstrap(
  object,
  R = 200L,
  workers = NULL,
  ci_level = 0.95,
  seed = NULL,
  verbose = FALSE
)

Arguments

Details

Because each replicate is refit at a FIXED penalty/lambda, this captures sampling variability in the coefficient estimates GIVEN that regularization strength, but not uncertainty in the choice of lambda itself, and not formal selection-event-conditional inference for which predictors end up nonzero. A replicate in which a coefficient is shrunk to exactly zero contributes 0 to that coefficient's bootstrap distribution (not 'NA'), so the reported bootstrap SE reflects both estimation and selection variability together for that coefficient.

Value

An object of class '"pecme_bootstrap"': a list with 'table' (a data frame of 'estimate', 'boot_mean', 'boot_se', and percentile 'ci_lower'/'ci_upper' per coefficient), 'replicates' (the raw 'R x p' matrix of bootstrap coefficient estimates), 'R', 'n_failed' (replicates that errored or failed to converge, excluded), and 'ci_level'.

Conditional log-likelihood contribution log P(observation | B = b)

Description

Vectorized over observations; 'mu' may be a matrix ('n x K') for simultaneous evaluation at several quadrature nodes.

Usage

pecme_conditional_loglik(y, mu, sigma, censor, left, right)

Closed-form coordinate update for the penalized quadratic subproblem

Description

Solves, in the single coordinate 'beta_j',

\min_{\beta_j} \frac{v_j \beta_j^2 - 2 z_j \beta_j}{2\sigma^2} + p_\lambda(|\beta_j|)

for each of the supported penalties, where 'z_j' is the partial cross-product of column 'j' with the partial residual and 'v_j' is its sum of squares. 'lambda_star = lambda * sigma^2' rescales the penalty onto the same footing as the (unscaled) least-squares term.

Usage

pecme_coordinate_update(
  zj,
  vj,
  sigma,
  lambda,
  penalty,
  alpha = 0.5,
  gamma = 3.7,
  wj = 1
)

Naive parameter count used for AIC/BIC/EBIC/GCV

Description

Counts nonzero fixed-effects coefficients plus two variance components ('sigma', 'sigma_b'). This is a convenient, commonly used HEURISTIC, not the effective degrees of freedom in the technical sense the AIC/GCV theory assumes. It is a reasonable approximation for Lasso-type penalties under regularity conditions (Zou, Hastie & Tibshirani, 2007, show 'E[#nonzero]' is an unbiased estimate of the Lasso's effective df), but it is NOT generally accurate for Ridge (which keeps essentially every coefficient nonzero, so this count is close to the maximum 'p' regardless of how much shrinkage is actually applied – see [pecme_effective_df()] for a proper, closed-form Ridge alternative) nor for Elastic Net, SCAD, or MCP (whose effective df would need a generalized-df/local-derivative calculation, e.g. Zou et al. 2007's generalized df; not implemented here). AIC/BIC/EBIC/GCV computed from this 'df' should therefore be read as approximate, especially for Ridge, Elastic Net, SCAD, and MCP.

Usage

pecme_df(beta, has_intercept = TRUE, tol = 1e-08)

Sequential (single-core) E-step

Description

Sequential (single-core) E-step

Usage

pecme_e_step(
  y,
  X,
  beta,
  group,
  sigma,
  sigma_b,
  censor,
  left = NULL,
  right = NULL,
  quadrature = 40L
)

Arguments

Value

A list of posterior moments: 'Ey', 'Ey2' (per observation), 'Eb', 'Ebb' (per observation, broadcast from the group), 'Eyb' (per observation), and 'group_Eb', 'group_Ebb' (per group, in the order of 'sort(unique(group))').

E-step dispatcher: sequential or cluster-parallel over groups

Description

With 'workers <= 1', or too few groups to make parallelism worthwhile, this simply calls [pecme_e_step()]. Otherwise the groups are partitioned across a pre-existing PSOCK/FORK cluster (created once per [pecme()] call and reused across iterations).

Usage

pecme_e_step_dispatch(
  y,
  X,
  beta,
  group,
  sigma,
  sigma_b,
  censor,
  left = NULL,
  right = NULL,
  quadrature = 40L,
  workers = 1L,
  cluster = NULL
)

Arguments

Working Ridge-type effective degrees-of-freedom approximation

Description

For Ridge, the naive nonzero-coefficient count from [pecme_df()] is a poor stand-in for effective df. The classical Ridge smoother trace is used as a CONDITIONAL/WORKING approximation on the standardized fixed- effects design. It is not claimed to be the exact effective df of the complete PECME estimator, which also contains censoring, random effects, estimated variance components, quadrature, and possibly adaptive penalization.

Usage

pecme_effective_df(object)

Arguments

Details

No corresponding closed form is implemented here for Elastic Net, SCAD, or MCP (a generalized-df/local-derivative calculation would be needed, e.g. Zou, Hastie & Tibshirani, 2007); for '"none"'/'"lasso"'/ '"adaptive"' the naive nonzero count is already a standard, reasonable approximation and is returned unchanged (plus the 2 variance components).

Value

A working/conditional df approximation, including the intercept if present and the two variance components, for '"ridge"', '"none"', '"lasso"', or '"adaptive"' fits; 'NA_real_' (with an explanatory 'attr(, "note")') for '"elastic"', '"scad"', or '"mcp"', where no closed form is implemented.

References

Hastie, T. and Tibshirani, R. (1990). *Generalized Additive Models*. Chapman and Hall. (Ridge effective df / "trace of the hat matrix".)

Zou, H., Hastie, T. and Tibshirani, R. (2007). On the "degrees of freedom" of the lasso. *Annals of Statistics*, 35(5), 2173-2192.

Recode a response/censoring specification to include missing values

Description

Given raw vectors 'y', 'left', 'right' and a 'censor' label vector ('"none"', '"left"', '"right"', '"interval"'), returns a new 'censor' vector in which any row with 'is.na(y)' (for '"none"'/'"left"'/ '"right"'/'"interval"' rows that have no usable censoring bound either) is relabelled '"missing"'.

Usage

pecme_flag_missing_response(y, left, right, censor)

Handle missing values in the fixed-effects design matrix

Description

Handle missing values in the fixed-effects design matrix

Usage

pecme_handle_missing_x(
  X,
  method = c("error", "cca", "mean", "reg"),
  has_intercept = TRUE,
  max_iter = 10L
)

Arguments

Value

A list with the (possibly modified) 'X', a logical 'keep' vector (all 'TRUE' unless 'method = "cca"'), and 'n_missing'.

Information criteria (AIC/BIC/EBIC/GCV)

Description

AIC/BIC/GCV use 'df', the naive nonzero-coefficient-based parameter count of [pecme_df()] – a heuristic, not a general effective-degrees-of-freedom calculation (see [pecme_df()] and [pecme_effective_df()]). This matters most for Ridge (essentially no coefficients are ever exactly zero, so 'df' sits near the maximum regardless of shrinkage strength) and, to a lesser extent, Elastic Net/SCAD/MCP; treat AIC/BIC/EBIC/GCV from those fits as approximate, and prefer [pecme_effective_df()]-based comparisons for Ridge, or cross-validated criteria ('criterion = "cv"' in [grid_pecme()], [pecme_loo()]) when the choice of 'df' matters for the conclusion.

Usage

pecme_ic(logLik, df, n, p_total, gamma_ebic = 0.5, n_selected = NULL)

Build a log-spaced lambda grid

Description

Build a log-spaced lambda grid

Usage

pecme_lambda_grid(
  lambda_max,
  nlambda = 50L,
  lambda_min_ratio = NULL,
  n = NULL,
  p = NULL
)

Arguments

Value

A numeric vector of length 'nlambda' containing a decreasing log-spaced sequence from 'lambda_max' to 'lambda_min'.

Leave-one-group-out (or leave-one-fold-out) cross-validated log-likelihood

Description

Refits the model with each group (or fold of groups) held out in turn and evaluates the held-out marginal log-likelihood at the refit parameters – a genuine, frequentist out-of-sample predictive criterion for clustered data (leave-one-*observation*-out is not meaningful here since observations within a group share a random intercept).

Usage

pecme_loo(object, method = c("exact", "kfold"), folds = 10L, workers = NULL)

Arguments

Value

A list with 'elpd_loo' (summed held-out log-likelihood, on the same scale as 'logLik'), 'se' (a naive standard error across successful group/fold contributions), 'n_successful_folds', and 'pointwise' (a data frame of one row per group/fold).

Marginal observed-data log-likelihood, random intercept integrated out

Description

\ell(\beta,\sigma,\sigma_b) = \sum_{j=1}^J \log \int \prod_{i \in j} P(y_{ij} \mid x_{ij}'\beta + b, \sigma)\, \phi(b; 0, \sigma_b^2)\, db

Usage

pecme_marginal_loglik(
  y,
  X,
  group,
  beta,
  sigma,
  sigma_b,
  censor,
  left = NULL,
  right = NULL,
  quadrature = 40L
)

Arguments

Details

evaluated by Gauss-Hermite quadrature. This is the objective that the ECME step maximizes (over 'sigma', 'sigma_b') and that the penalized criterion in [pecme()] is built from (over 'beta', given 'sigma', 'sigma_b').

Comprehensive fit-quality metrics for a pecme model

Description

Collects the quantities typically reported for a mixed-model fit: pointwise prediction error (MSE/RMSE/MAE/MAPE/SMAPE), two kinds of pseudo-R-squared, information criteria, and (optionally) a cross-validated predictive log-likelihood.

Usage

pecme_metrics(
  object,
  newdata = NULL,
  newresponse = NULL,
  loo = FALSE,
  loo_method = c("exact", "kfold"),
  loo_folds = 10L,
  workers = NULL,
  ebic_gamma = 0.5
)

Arguments

Details

# Why there is no WAIC WAIC (and DIC) are defined from the variance of the pointwise log-likelihood *across posterior draws* of the parameters. 'pecme' is a frequentist penalized-likelihood estimator: it returns a single point estimate (plus, for an unpenalized fit, an asymptotic Hessian- based covariance), not a posterior sample. Reporting a WAIC number here would therefore not be a real WAIC – rather than approximate it, '$WAIC' is returned as 'NA' with an explanatory 'attr'. For out-of-sample model comparison, use 'AIC'/'BIC'/'EBIC' (fast, from the fitted likelihood) or 'loo = TRUE' (slower, but a genuine held-out predictive log-likelihood, which is the same quantity WAIC and PSIS-LOO both aim to approximate).

# Two R-squared values 'R2_marginal' and 'R2_conditional' follow Nakagawa & Schielzeth (2013): 'R2_marginal = Var(fixed) / (Var(fixed) + sigma_b^2 + sigma^2)‘ uses only the fixed-effects predictor’s variance in the numerator; 'R2_conditional' additionally credits the random intercept, '(Var(fixed) + sigma_b^2) / (Var(fixed) + sigma_b^2 + sigma^2)'. These are the R-squared values normally reported for a mixed model. 'R2_predictive' = '1 - MSE / Var(y)' is the more familiar OLS-style quantity computed directly from predictions; it is provided for intuition/comparability but is not the standard mixed-model R-squared. 'R2_predictive_adj' heuristically applies the usual OLS adjusted-R2 formula with 'object$df'; this adjustment is not a rigorously derived quantity for mixed models and should be reported as a diagnostic, not a citation-grade statistic.

Value

An object of class '"pecme_metrics"' (a list; has a 'print' method).

Combine ‘pecme' fits across multiply-imputed data sets (Rubin’s rules)

Description

A convenience helper for the recommended workflow when fixed-effects covariates (not the response) have missing values: impute 'm' complete data sets externally (e.g. with the 'mice' package), fit the *same* ‘pecme()' model to each, and combine with Rubin’s rules.

Usage

pecme_mi_combine(fits, use_refit = FALSE)

Arguments

Details

Rubin's rules require a (classically) valid standard error for every fit being combined. A penalized fit ('penalty != "none"' and 'lambda > 0') deliberately stores 'se = NA' (see [pecme()]): naive standard errors are not valid for a shrinkage estimate, so silently pooling 'NA's here would produce 'NA' within-imputation variance and a misleading pooled result. Because of this, 'pecme_mi_combine()':

* by default, REQUIRES every element of 'fits' to have valid SEs ('se_valid == TRUE', i.e. 'penalty == "none"' or 'lambda == 0') and errors otherwise, naming which elements fail and why; * with ‘use_refit = TRUE', instead pools each fit’s '$refit' – the automatic unpenalized "relaxed refit" on the selected predictors (see [pecme()]'s 'refit' argument) – which does have valid SEs. This is a practical compromise, not a rigorous solution: (a) the relaxed refit is not valid unconditional post-selection inference to begin with (its SEs ignore the selection step); (b) different imputed data sets can select different variable subsets, so a coefficient absent from one imputation's refit is treated as exactly 0 (with 'se = 0') in that imputation for pooling purposes, which folds selection instability into the between-imputation variance rather than separating it out. Report results obtained this way as approximate/exploratory.

For fully rigorous inference on a penalized, multiply-imputed fit, see [pecme_bootstrap()] (a cluster bootstrap) applied within each imputed data set, or consult a specialist reference on selective/ post-selection inference under multiple imputation.

Value

A data frame with the pooled estimate, within-imputation variance, between-imputation variance, total variance, and pooled standard error for each fixed-effect coefficient.

Penalized observed-data objective: -logLik(beta, sigma, sigma_b) + P_lambda(beta)

Description

Penalized observed-data objective: -logLik(beta, sigma, sigma_b) + P_lambda(beta)

Usage

pecme_objective(
  beta,
  y,
  X,
  sigma,
  sigma_b,
  group,
  censor,
  left = NULL,
  right = NULL,
  penalty,
  lambda,
  alpha = 0.5,
  gamma = 3.7,
  weights = NULL,
  has_intercept = TRUE,
  quadrature = 40L
)

Conditional (given the random intercept) observed-data log-likelihood

Description

Sum over observations of 'log P(observation | mu, sigma)', where 'mu' already includes the random intercept. Used mainly for testing.

Usage

pecme_observed_loglik(y, mu, sigma, censor, left = NULL, right = NULL)

Conditional expected complete-data objective for beta (the CM-step target)

Description

Given E-step moments, the beta-dependent part of the expected complete-data penalized log-likelihood reduces to

Q(\beta) = \frac{1}{2\sigma^2} \sum_i (E[y^*_i] - E[b_{g(i)}] - x_i'\beta)^2 + P_\lambda(\beta) + \text{const}

because E\|y^* - b - X\beta\|^2 = \mathrm{Var}(y^*-b\mid \text{data}) + (E[y^*-b\mid\text{data}] - X\beta)^2 and the variance term does not depend on beta. This is what makes coordinate descent on the adjusted response 'E[y*] - E[b]' exact rather than an approximation.

Usage

pecme_q_beta(
  beta,
  y_star,
  X,
  random_effect = NULL,
  sigma,
  penalty = .pecme_penalties,
  lambda = 0,
  alpha = 0.5,
  gamma = 3.7,
  weights = NULL,
  has_intercept = TRUE
)

Check sensitivity of a fitted model to the Gauss-Hermite quadrature order

Description

Fixed-order Gauss-Hermite quadrature is a numerical approximation to the (otherwise intractable) integral over the random intercept, not an exact evaluation of it. This function reports how much the fit changes as the number of quadrature nodes changes, so 'quadrature' is a checked choice rather than an unexamined default. Approximation error is typically largest when 'sigma_b' is large relative to 'sigma', or with heavy censoring/many missing responses.

Usage

pecme_quadrature_check(
  object,
  quadrature_values = c(15L, 20L, 30L, 50L, 80L, 120L),
  method = c("loglik", "refit")
)

Arguments

Value

An object of class '"pecme_quadrature_check"': a data frame ('quadrature', 'logLik', and, for 'method = "refit"', 'sigma', 'sigma_b', 'max_abs_beta_change' relative to 'object'), with a 'print' method that flags whether the range of 'logLik' values looks small relative to typical evidence thresholds.

Examples

sim <- simulate_pecme_data(n_groups = 25, n_per_group = 8, p = 3, seed = 1)
fit <- pecme(y ~ x1 + x2 + x3, random_var = "group", data = sim$data)
pecme_quadrature_check(fit)

Conditional first/second moments of (possibly censored or missing) Y*

Description

Vectorized over observations and, optionally, over a second dimension (quadrature nodes for the random intercept): 'mu' may be a matrix of shape 'n x K', in which case moments are returned as 'n x K' matrices.

Usage

pecme_truncated_moments(mu, sigma, censor, left, right, y)

Value of the penalty term P_lambda(beta)

Description

Value of the penalty term P_lambda(beta)

Usage

penalty_value(
  beta,
  penalty = .pecme_penalties,
  lambda,
  alpha = 0.5,
  gamma = 3.7,
  weights = NULL,
  has_intercept = TRUE
)

Arguments

Plot the regularization path and/or selection-criterion curve

Description

Plot the regularization path and/or selection-criterion curve

Usage

## S3 method for class 'pecme_grid'
plot(x, which = c("both", "path", "criterion"), ...)

Arguments

Value

Invisibly returns the input 'x', an object of class '"pecme_grid"'.

Predict from a fitted pecme model

Description

Predict from a fitted pecme model

Usage

## S3 method for class 'pecme'
predict(
  object,
  newdata = NULL,
  level = c("population", "group"),
  use_refit = FALSE,
  se_fit = FALSE,
  ...
)

Arguments

Value

If 'se_fit = FALSE', a numeric vector of predictions. If 'se_fit = TRUE', a list with components 'fit' and 'se.fit', containing the predicted values and their estimated standard errors, respectively.

Simulate data from a censored/missing random-intercept linear mixed model

Description

Generates Y^*_{ij} = \beta_0 + x_{ij}'\beta + b_j + e_{ij}, e_{ij}\sim N(0,\sigma^2), b_j\sim N(0,\sigma_b^2), with 'p' standard-normal (optionally correlated) predictors, of which only 'n_active' have a non-zero coefficient – convenient for checking variable-selection performance. Censoring and/or missingness can be injected on top of the latent 'y_star'.

Usage

simulate_pecme_data(
  n_groups = 30L,
  n_per_group = 10L,
  p = 10L,
  n_active = 4L,
  beta_active = 1.5,
  rho = 0,
  sigma = 1,
  sigma_b = 0.7,
  censor_type = c("none", "left", "right", "interval"),
  censor_prob = 0,
  missing_prob = 0,
  seed = NULL
)

Arguments

Value

A list with 'data' (a data frame with columns 'y', 'x1..xp', 'group', 'censor', 'left', 'right'), 'beta_true' (length 'p + 1', intercept first), 'sigma_true', 'sigma_b_true', and 'y_star' (the uncensored, non-missing latent response, for validation).

Examples

sim <- simulate_pecme_data(n_groups = 20, n_per_group = 8, p = 6,
                            n_active = 3, censor_type = "left",
                            censor_prob = 0.2, seed = 1)
head(sim$data)

Soft-thresholding operator

Description

Soft-thresholding operator

Usage

soft_threshold(z, a)