Simulate sparse GLMs and deconvolution problems

Generate blurred spike trains or regression data with Gaussian predictors. Returned data and ground truth can be passed to sparseglm(), or to L0glm.fit() and L0glm.tune() with the returned family and intercept flag.

Usage

simul(
  method = "deconvolution",
  n = 500,
  p = 500,
  s = 0.1,
  sparse = NULL,
  options = list(),
  seed = 1,
  family = NULL,
  n_nonzero = NULL,
  intercept = FALSE,
  m = 1L,
  beta = NULL
)

Arguments

method

Setup: "deconvolution", "normal", "abess", or "bestsubset". Numeric values 1/2 select single-channel Gaussian/Poisson deconvolution, 3/4 are their multichannel versions, and 5 is independent normal-design Gaussian regression.

n, p

Numbers of observations and predictors, excluding any intercept.

s

Fraction of predictor rows selected if n_nonzero is omitted. Numeric deconvolution setups use round(s*n), capped at min(n,p); named setups use round(s*p). Setup 5 uses min(10,p) unless a nonzero count is supplied.

sparse

Return sparse design and coefficient matrices? NULL uses sparse storage for deconvolution and dense storage for Gaussian designs. Sparse storage usually wastes memory for a fully populated Gaussian design.

options

Named list of setup-specific arguments; see Details.

seed

Nonnegative integer seed. The caller's RNG state is restored on exit. NULL uses and advances the current random-number stream.

family

An R GLM family object/function, or a family name. Names include "gaussian", "binomial", "poisson", "Gamma", "inverse.gaussian", "nb", "multinomial" and "cox". Positive-response regression families default to a log link; deconvolution uses identity links. Special-family objects from mgcv are accepted. NULL uses Poisson for named deconvolution, Gaussian for regression, or the numbered setup's family. signal_families objects use their fixed NB1 or Gamma noise parameters; their Gamma detection limits are applied when drawing observations.

n_nonzero

Number of nonzero predictor rows, excluding the intercept. Overrides s; a supplied beta determines its own support. Deconvolution caps this number at min(n,p). See beta.type below for weak sparsity.

intercept

FALSE/NULL omits an intercept; TRUE uses one. A finite numeric scalar or one value per channel/category supplies its value, including zero. Multinomial intercepts are centered across categories. For Cox this multiplies the baseline hazard on the log scale.

m

Number of response channels, or multinomial/ordinal categories. A multinomial family object determines the category count. Numbered setups 3/4 use all informative library mass channels when m is omitted. Explicit m>1 deconvolution uses random sparse, normalized spectra with exactly that many channels.

beta

Optional finite slope vector or p by m matrix, in original predictor units, for Gaussian-design setups. A vector is repeated across channels. Multinomial rows are centered to CLR coordinates. Ordinal models use a single slope vector shared across categories.

Value

A list with X (without an intercept column), y, y_true (conditional means, category probabilities, or Cox relative risks), eta_true, beta_true (including any identifiable intercept row), intercept (logical fitting flag), intercept_true, family, and weights. Unit weights have one shared column, so subsetting response channels does not require subsetting them; binomial weights contain the trial totals. Regression data include covariance (structure and rho), generating dispersion, and NB theta. With nval>0, Xval, yval, yval_true, eta_val and weightsval describe validation data. Deconvolution also returns time and spike grids. Cox responses are right-censored Surv objects.

Details

Deconvolution. Options are arguments of simulate_spike_train(), including mean_beta, sd_logbeta, sd_noise, psf, thresh, theta, clip, spectra and channel.sparsity. Gaussian responses are untruncated by default. Numeric setups 1/3 clip negative responses to zero by default; set options$clip=FALSE for untruncated Gaussian noise. With clipping, y_true is the latent uncensored mean.

Gaussian predictors. "normal" uses independent standard normal columns and random support with uniform beta_range (default c(-1,1)). All regression setups accept rho and correlation ("independent", "ar1", "equicorrelated"), or ABESS cortype=1:3. AR(1) covariance is \(\Sigma_{ij}=\rho^{|i-j|}\); equicorrelation has off-diagonal entries \(\rho\). Generation and population signal-variance calculations avoid allocating a p by p covariance matrix.

ABESS setups. "abess" follows predictor and coefficient distributions in abess::generate.data(): random support and family-specific uniform ranges proportional to \(\sqrt{2\log(p)/n}\). support.size, y.dim, class.num, sigma, snr, weibull.shape and uniform.max are recognized. Non-Gaussian responses are sampled directly from the requested family, without the source's extra latent Gaussian noise or linear-predictor clipping; therefore y_true and beta_true describe the same specified model. "mgaussian" chooses Gaussian multichannel data. "ordinal" generates proportional-odds data with ordered cutpoints (default equally spaced from -1 to 1); L0glm does not fit ordinal models.

Best-subset setups. "bestsubset" implements the bestsubset::sim.xy() coefficient patterns. beta.type=1 places unit coefficients at equally spaced locations; 2 uses the first selected rows; 3 decreases coefficients linearly from 10 to 0.5; 4 uses c(-10,-6,-2,2,6,10) in the first six rows (requires p>=6, ignoring n_nonzero); 5 uses unit coefficients in the first selected rows followed by 0.5^(1:(p-n_nonzero)). Type 5 is weakly sparse: trailing coefficients are also nonzero until numerical underflow. Nonzero rho defaults to AR(1) here. nval (default zero) adds independent validation data from the same model, for any Gaussian-design setup, without changing the training sample. Multichannel and non-Gaussian outcomes extend the source's Gaussian model. For multinomial data the patterns multiply centered category contrasts.

Noise and responses. dispersion is a scalar or channel vector, default one: Gaussian variance, Gamma squared coefficient of variation, or inverse Gaussian dispersion. Gaussian error_sd, sd_noise or sigma specifies SD instead. Alternatively snr sets Gaussian variance to \(\beta^T\Sigma\beta/\mathrm{snr}\) for identity-link Gaussian models (defaults 10 for ABESS, 1 for bestsubset). SNR is not applied to non-Gaussian outcomes. theta is positive NB shape (default 10, or the supplied fixed family shape); conflicting fixed family and theta values are rejected. NB models have \(\mathrm{Var}(Y)=\mu+\mu^2/\theta\). trials (default one) sets binomial or multinomial totals; trialsval optionally supplies validation totals. Binomial responses are proportions with trial totals returned as weights; multinomial responses are count matrices, or factors with response="factor" and one trial. Multinomial coefficients use CLR. Quasipoisson uses Poisson/NB draws matching \(\phi\mu\); quasibinomial uses binomial/beta-binomial draws matching \(\phi\mu(1-\mu)/\mathrm{trials}\), requiring attainable overdispersion. Unrestricted quasi() is rejected because it does not specify a sampling distribution.

Cox data have \(H_0(t)=a t^b\): baseline.rate=a, weibull.shape=b, both default one. Censoring is independently uniform on [0, uniform.max], with default upper bound one. An intercept multiplies a by exp(intercept). y_true contains relative risks; event_time contains sampled uncensored times. No Cox intercept row is added to beta_true. Invalid means for a requested link are rejected rather than clipped; choose a valid link, intercept and coefficients.

References

ABESS source: https://github.com/abess-team/abess/blob/master/R-package/R/generate.data.R. Hastie, Tibshirani and Tibshirani, best-subset source: https://github.com/ryantibs/best-subset/blob/master/bestsubset/R/sim.R.

See also

simulate_spike_train(), L0glm.fit(), sparseglm()

Examples

# Blurred spike trains with a known PSF, sparse peaks and an optional baseline.
single <- simul("deconvolution", n=200, p=200, n_nonzero=20,
  family="poisson", intercept=1, seed=2)
multi <- simul("deconvolution", n=100, p=120, n_nonzero=10,
  family="gaussian", m=4, intercept=c(1, 2, 3, 4), seed=2,
  options=list(mean_beta=100, sd_noise=2))
dim(multi$beta_true)                  # 120 slopes plus an intercept, four channels
if (requireNamespace("mgcv", quietly=TRUE)) {
  nb_spikes <- simul("deconvolution", n=100, p=100, n_nonzero=10,
    family="nb", m=4, intercept=1, options=list(theta=5))
}

# Standard normal predictors, arbitrary n/p, and signed sparse coefficients.
regression <- simul("normal", n=150, p=30, n_nonzero=4, intercept=1.5,
  options=list(beta_range=c(-2, 2), error_sd=.5), seed=10)
fit <- L0glm.fit(regression$X, regression$y, family=regression$family,
  intercept=regression$intercept, nonnegative=FALSE, parallelize=FALSE)
L0glm.metrics(fit, beta_true=regression$beta_true)[c("gic", "F1_score", "MCC")]

# sparseglm reads family, intercept and coefficient truth from the simulation.
comparison <- sparseglm(regression, method=c("L0glm", "lasso"), nlambda=5)
L0glm.metrics(comparison)[c("method", "seconds", "gic", "F1_score", "MCC")]

# ABESS-style correlated predictors and coefficient magnitudes.
abess_data <- simul("abess", n=200, p=100, n_nonzero=5,
  family="binomial", intercept=-1, options=list(cortype=2, rho=.5))

# Best-subset pattern 3: decreasing effects, population SNR=2, validation data.
subset_data <- simul("bestsubset", n=100, p=50, n_nonzero=5, intercept=2,
  options=list(beta.type=3, rho=.5, snr=2, nval=50))
subset_data$dispersion
dim(subset_data$Xval)

# Other standard links and special families use the same normal-design setup.
probit_data <- simul("normal", n=100, p=20, n_nonzero=3,
  family=binomial("probit"), intercept=-.5)
gamma_data <- simul("normal", n=100, p=20, n_nonzero=3,
  family=Gamma("log"), intercept=1, options=list(dispersion=.5))
if (requireNamespace("mgcv", quietly=TRUE) && requireNamespace("survival", quietly=TRUE)) {
  nb_data <- simul("normal", n=150, p=20, n_nonzero=3,
    family=mgcv::nb(), intercept=1, options=list(theta=3))
  multinomial_data <- simul("normal", n=150, p=20, n_nonzero=3,
    family=mgcv::multinom(K=2), intercept=c(.3, 0, -.3))
  cox_data <- simul("normal", n=150, p=20, n_nonzero=3,
    family=mgcv::cox.ph(), intercept=log(2)) # doubles the baseline hazard
  cox_data$intercept                    # FALSE: no estimable Cox intercept
}