Relevant One-step Selective Inference for the LASSO
Description
Compute p-values and confidence intervals for the lasso estimate, at a fixed value of the tuning parameter lambda using the "relevant" conditioning event of arxiv.org/1801.09037.
Usage
ROSI(X,
y,
soln,
lambda,
penalty_factor=NULL,
dispersion=1,
family=c('gaussian', 'binomial'),
solver=c('QP', 'glmnet'),
construct_ci=TRUE,
debiasing_method=c("JM", "BN"),
verbose=FALSE,
level=0.9,
use_debiased=TRUE)
Arguments
X
Matrix of predictors (n by p);
y
Vector of outcomes (length n)
soln
Estimated lasso coefficients (e.g., from glmnet). This is of length p (so the intercept is not included as the first component).
Be careful! This function uses the "standard" lasso objective
1/2 \|y - X \beta\|_2^2 + \lambda \|\beta\|_1.
In contrast, glmnet multiplies the first term by a factor of 1/n.
So after running glmnet, to extract the beta corresponding to a value lambda,
you need to use beta = coef(obj, s=lambda/n)[-1],
where obj is the object returned by glmnet (and [-1] removes the intercept,
which glmnet always puts in the first component)
lambda
Value of lambda used to compute beta. See the above warning
penalty_factor
Penalty factor as used by glmnet. Actual penalty used in solving the problem is
\lambda \cdot \sum_{i=1}^p f_i |\beta_i|
with f being the penalty_factor. Defaults to vector of 1s.
dispersion
Estimate of dispersion in the GLM. Can be taken to be 1 for logisitic and should be an estimate of the error variance for the Gaussian.
family
Family used for likelihood.
solver
Solver used to solve restricted problems needed to find truncation set. Each active variable requires solving a new LASSO problem obtained by zeroing out one coordinate of original problem. The "QP" choice uses coordinate descent for a specific value of lambda, rather than glmnet which would solve for a new path each time.
construct_ci
Report confidence intervals or just p-values?
debiasing_method
Which method should be used for debiasing? Choices are "JM" (Javanmard, Montanari) or "BN" (method described in arxiv.org/1703.03282).
verbose
Print out progress along the way? Default is FALSE.
level
Confidence level for intervals.
use_debiased
Use the debiased estimate of the parameter or not. When FALSE, this is the method desribed in arxiv.org/1801.09037. The default TRUE often produces noticably shorter intervals and more powerful tests when p is comparable to n. Ignored when n<p and set to TRUE. Also note that with "BN" as debiasing method and n > p, this agrees with method in arxiv.org/1801.09037.
Details
???
Value
active_set
Active set of LASSO.
pvalues
Two-sided P-values for active variables.
intervals
Confidence intervals
estimate
Relaxed (i.e. unshrunk) selected estimates.
std_err
Standard error of relaxed estimates (pre-selection).
dispersion
Dispersion parameter.
lower_trunc
Lower truncation point. The estimates should be outside the interval formed by the lower and upper truncation poitns.
upper_trunc
Lower truncation point. The estimates should be outside the interval formed by the lower and upper truncation poitns.
lambda
Value of tuning parameter lambda used.
penalty_factor
Penalty factor used for solving problem.
level
Confidence level.
call
The call to fixedLassoInf.
Author(s)
Jelena Markovic, Jonathan Taylor
References
Keli Liu, Jelena Markovic, Robert Tibshirani. More powerful post-selection inference, with application to the Lasso. arXiv:1801.09037
Tom Boot, Didier Nibbering. Inference in high-dimensional linear regression models. arXiv:1703.03282
Examples
library(selectiveInference)
library(glmnet)
set.seed(43)
n = 100
p = 200
s = 2
sigma = 1
x = matrix(rnorm(n*p),n,p)
x = scale(x,TRUE,TRUE)
beta = c(rep(10, s), rep(0,p-s)) / sqrt(n)
y = x %*% beta + sigma*rnorm(n)
# first run glmnet
gfit = glmnet(x,y,standardize=FALSE)
# extract coef for a given lambda; note the 1/n factor!
# (and we don't save the intercept term)
lambda = 4 * sqrt(n)
lambda_glmnet = 4 / sqrt(n)
beta = selectiveInference:::solve_problem_glmnet(x,
y,
lambda_glmnet,
penalty_factor=rep(1, p),
family="gaussian")
# compute fixed lambda p-values and selection intervals
out = ROSI(x,
y,
beta,
lambda,
dispersion=sigma^2)
out
# an alternate approximate inverse from Boot and Nibbering
out = ROSI(x,
y,
beta,
lambda,
dispersion=sigma^2,
debiasing_method="BN")
out
Truncated Gaussian confidence interval.
Description
Compute truncated Gaussian interval of Lee et al. (2016) with arbitrary affine selection and covariance. Z should satisfy A
Usage
TG.interval(Z, A, b, eta, Sigma=NULL, alpha=0.1,
gridrange=c(-100,100),
gridpts=100,
griddepth=2,
flip=FALSE,
bits=NULL)
Arguments
Z
Observed data (assumed to follow N(mu, Sigma) with sum(eta*mu)=null_value)
A
Matrix specifiying affine inequalities AZ <= b
b
Offsets in the affine inequalities AZ <= b.
eta
Determines the target sum(eta*mu) and estimate sum(eta*Z).
Sigma
Covariance matrix of Z. Defaults to identity.
alpha
Significance level for confidence intervals (target is miscoverage alpha/2 in each tail)
gridrange
Grid range for constructing confidence intervals, on the standardized scale.
gridpts
???????
griddepth
???????
flip
???????
bits
Number of bits to be used for p-value and confidence interval calculations. Default is
NULL, in which case standard floating point calculations are performed. When not NULL,
multiple precision floating point calculations are performed with the specified number
of bits, using the R package Rmpfr (if this package is not installed, then a
warning is thrown, and standard floating point calculations are pursued).
Note: standard double precision uses 53 bits
so, e.g., a choice of 200 bits uses about 4 times double precision. The confidence
interval computation is sometimes numerically challenging, and the extra precision can be
helpful (though computationally more costly). In particular, extra precision might be tried
if the values in the output columns of tailarea differ noticeably from alpha/2.
Details
This function computes selective confidence intervals based on the polyhedral lemma of Lee et al. (2016).
Value
int
Selective confidence interval.
tailarea
Realized tail areas (lower and upper) for each confidence interval.
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid
References
Jason Lee, Dennis Sun, Yuekai Sun, and Jonathan Taylor (2016). Exact post-selection inference, with application to the lasso. Annals of Statistics, 44(3), 907-927.
Jonathan Taylor and Robert Tibshirani (2017) Post-selection inference for math L1-penalized likelihood models. Canadian Journal of Statistics, xx, 1-21. (Volume still not posted)
Examples
A = diag(5)
b = rep(1, 5)
Z = rep(0, 5)
Sigma = diag(5)
eta = as.numeric(c(1, 1, 0, 0, 0))
TG.interval(Z, A, b, eta, Sigma)
Truncation limits and standard deviation.
Description
Compute truncated limits and SD for use in computing p-values or confidence intervals of Lee et al. (2016). Z should satisfy A
Usage
TG.limits(Z, A, b, eta, Sigma)
Arguments
Z
Observed data (assumed to follow N(mu, Sigma) with sum(eta*mu)=null_value)
A
Matrix specifiying affine inequalities AZ <= b
b
Offsets in the affine inequalities AZ <= b.
eta
Determines the target sum(eta*mu) and estimate sum(eta*Z).
Sigma
Covariance matrix of Z. Defaults to identity.
Details
This function computes the limits of truncation and the implied standard deviation in the polyhedral lemma of Lee et al. (2016).
Value
vlo
Lower truncation limits for statistic
vup
Upper truncation limits for statistic
sd
Standard error of sum(eta*Z)
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid
References
Jason Lee, Dennis Sun, Yuekai Sun, and Jonathan Taylor (2016). Exact post-selection inference, with application to the lasso. Annals of Statistics, 44(3), 907-927.
Jonathan Taylor and Robert Tibshirani (2017) Post-selection inference for math L1-penalized likelihood models. Canadian Journal of Statistics, xx, 1-21. (Volume still not posted)
Examples
A = diag(5)
b = rep(1, 5)
Z = rep(0, 5)
Sigma = diag(5)
eta = as.numeric(c(1, 1, 0, 0, 0))
TG.limits(Z, A, b, eta, Sigma)
Truncated Gaussian p-value.
Description
Compute truncated Gaussian p-value of Lee et al. (2016) with arbitrary affine selection and covariance. Z should satisfy A
Usage
TG.pvalue(Z, A, b, eta, Sigma, null_value=0, bits=NULL)
Arguments
Z
Observed data (assumed to follow N(mu, Sigma) with sum(eta*mu)=null_value)
A
Matrix specifiying affine inequalities AZ <= b
b
Offsets in the affine inequalities AZ <= b.
eta
Determines the target sum(eta*mu) and estimate sum(eta*Z).
Sigma
Covariance matrix of Z. Defaults to identity.
null_value
Hypothesized value of sum(eta*mu) under the null.
bits
Number of bits to be used for p-value and confidence interval calculations. Default is
NULL, in which case standard floating point calculations are performed. When not NULL,
multiple precision floating point calculations are performed with the specified number
of bits, using the R package Rmpfr (if this package is not installed, then a
warning is thrown, and standard floating point calculations are pursued).
Note: standard double precision uses 53 bits
so, e.g., a choice of 200 bits uses about 4 times double precision. The confidence
interval computation is sometimes numerically challenging, and the extra precision can be
helpful (though computationally more costly). In particular, extra precision might be tried
if the values in the output columns of tailarea differ noticeably from alpha/2.
Details
This function computes selective p-values based on the polyhedral lemma of Lee et al. (2016).
Value
pv
One-sided P-values for active variables, uses the fact we have conditioned on the sign.
vlo
Lower truncation limits for statistic
vup
Upper truncation limits for statistic
sd
Standard error of sum(eta*Z)
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid
References
Jason Lee, Dennis Sun, Yuekai Sun, and Jonathan Taylor (2016). Exact post-selection inference, with application to the lasso. Annals of Statistics, 44(3), 907-927.
Jonathan Taylor and Robert Tibshirani (2017) Post-selection inference for math L1-penalized likelihood models. Canadian Journal of Statistics, xx, 1-21. (Volume still not posted)
Examples
A = diag(5)
b = rep(1, 5)
Z = rep(0, 5)
Sigma = diag(5)
eta = as.numeric(c(1, 1, 0, 0, 0))
TG.pvalue(Z, A, b, eta, Sigma)
TG.pvalue(Z, A, b, eta, Sigma, null_value=1)
Find an approximate inverse of a non-negative definite matrix.
Description
Find some rows of an approximate inverse of a non-negative definite symmetric matrix by solving optimization problem described in Javanmard and Montanari (2013). Can be used for approximate Newton step from some consistent estimator (such as the LASSO) to find a debiased solution.
Usage
debiasingMatrix(Xinfo,
is_wide,
nsample,
rows,
verbose=FALSE,
bound=NULL,
linesearch=TRUE,
scaling_factor=1.5,
max_active=NULL,
max_try=10,
warn_kkt=FALSE,
max_iter=50,
kkt_stop=TRUE,
parameter_stop=TRUE,
objective_stop=TRUE,
kkt_tol=1.e-4,
parameter_tol=1.e-4,
objective_tol=1.e-4)
Arguments
Xinfo
Either a non-negative definite matrix S=t(X) is_wide is TRUE, then Xinfo should be X, otherwise it should be S.
is_wide
Are we solving for rows of the debiasing matrix assuming it is a wide matrix so that Xinfo=X and the non-negative definite matrix of interest is t(X)
nsample
Number of samples used in forming the cross-covariance matrix. Used for default value of the bound parameter.
rows
Which rows of the approximate inverse to compute.
verbose
Print out progress as rows are being computed.
bound
Initial bound parameter for each row. Will be changed if linesearch is TRUE.
linesearch
Run a line search to find as small as possible a bound parameter for each row?
scaling_factor
In the linesearch, the bound parameter is either multiplied or divided by this factor at each step.
max_active
How large an active set to consider in solving the problem with coordinate descent. Defaults to max(50, 0.3*nsample).
max_try
How many tries in the linesearch.
warn_kkt
Warn if the problem does not seem to be feasible after running the coordinate descent algorithm.
max_iter
How many full iterations to run of the coordinate descent for each value of the bound parameter.
kkt_stop
If TRUE, check to stop coordinate descent when KKT conditions are approximately satisfied.
parameter_stop
If TRUE, check to stop coordinate descent based on relative convergence of parameter vector, checked at geometrically spaced iterations 2^k.
objective_stop
If TRUE, check to stop coordinate descent based on relative decrease of objective value, checked at geometrically spaced iterations 2^k.
kkt_tol
Tolerance value for assessing whether KKT conditions for solving the dual problem and feasibility of the original problem.
parameter_tol
Tolerance value for assessing convergence of the problem using relative convergence of the parameter.
objective_tol
Tolerance value for assessing convergence of the problem using relative decrease of the objective.
Details
This function computes an approximate inverse as described in Javanmard and Montanari (2013), specifically display (4). The problem is solved by considering a dual problem which has an objective similar to a LASSO problem and is solvable by coordinate descent. For some values of bound the original problem may not be feasible, in which case the dual problem has no solution. An attempt to detect this is made by stopping when the active set grows quite large, determined by max_active.
Value
M
Rows of approximate inverse of Sigma.
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid
References
Adel Javanmard and Andrea Montanari (2013). Confidence Intervals and Hypothesis Testing for High-Dimensional Regression. Arxiv: 1306.3171
Examples
set.seed(10)
n = 50
p = 100
X = matrix(rnorm(n * p), n, p)
S = t(X) %*% X / n
M = debiasingMatrix(S, FALSE, n, c(1,3,5))
M2 = debiasingMatrix(X, TRUE, n, c(1,3,5))
max(M - M2)
Estimate the noise standard deviation in regression
Description
Estimates the standard deviation of the noise, for use in the selectiveInference package
Usage
estimateSigma(x, y, intercept=TRUE, standardize=TRUE)
Arguments
x
Matrix of predictors (n by p)
y
Vector of outcomes (length n)
intercept
Should glmnet be run with an intercept? Default is TRUE
standardize
Should glmnet be run with standardized predictors? Default is TRUE
Details
This function estimates the standard deviation of the noise, in a linear regresion setting.
A lasso regression is fit, using cross-validation to estimate the tuning parameter lambda.
With sample size n, yhat equal to the predicted values and df being the number of nonzero
coefficients from the lasso fit, the estimate of sigma is sqrt(sum((y-yhat)^2) / (n-df-1)).
Important: if you are using glmnet to compute the lasso estimate, be sure to use the settings
for the "intercept" and "standardize" arguments in glmnet and estimateSigma. Same applies to fs
or lar, where the argument for standardization is called "normalize".
Value
sigmahat
The estimate of sigma
df
The degrees of freedom of lasso fit used
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid
References
Stephen Reid, Jerome Friedman, and Rob Tibshirani (2014). A study of error variance estimation in lasso regression. arXiv:1311.5274.
Examples
set.seed(33)
n = 50
p = 10
sigma = 1
x = matrix(rnorm(n*p),n,p)
beta = c(3,2,rep(0,p-2))
y = x%*%beta + sigma*rnorm(n)
# run forward stepwise
fsfit = fs(x,y)
# estimate sigma
sigmahat = estimateSigma(x,y)$sigmahat
# run sequential inference with estimated sigma
out = fsInf(fsfit,sigma=sigmahat)
out
Expand a data frame with factors to form a design matrix with the full binary encoding of each factor.
Description
When using groupfs with factor variables call this function first to create a design matrix.
Usage
factorDesign(df)
Arguments
df
Data frame containing some columns which are factors.
Value
List containing
- x
Design matrix, the first columns contain any numeric variables from the original date frame.
- index
Group membership indicator for expanded matrix.
Examples
## Not run:
fd = factorDesign(warpbreaks)
y = rnorm(nrow(fd$x))
fit = groupfs(fd$x, y, fd$index, maxsteps=2, intercept=F)
pvals = groupfsInf(fit)
## End(Not run)
Inference for the lasso, with a fixed lambda
Description
Compute p-values and confidence intervals for the lasso estimate, at a fixed value of the tuning parameter lambda
Usage
fixedLassoInf(x,
y,
beta,
lambda,
family = c("gaussian", "binomial", "cox"),
intercept=TRUE,
add.targets=NULL,
status=NULL,
sigma=NULL,
alpha=0.1,
type=c("partial","full"),
tol.beta=1e-5,
tol.kkt=0.1,
gridrange=c(-100,100),
bits=NULL,
verbose=FALSE,
linesearch.try=10)
Arguments
x
Matrix of predictors (n by p);
y
Vector of outcomes (length n)
beta
Estimated lasso coefficients (e.g., from glmnet). This is of length p (so the intercept is not included as the first component).
Be careful! This function uses the "standard" lasso objective
1/2 \|y - x \beta\|_2^2 + \lambda \|\beta\|_1.
In contrast, glmnet multiplies the first term by a factor of 1/n.
So after running glmnet, to extract the beta corresponding to a value lambda,
you need to use beta = coef(obj, s=lambda/n)[-1],
where obj is the object returned by glmnet (and [-1] removes the intercept,
which glmnet always puts in the first component)
lambda
Value of lambda used to compute beta. See the above warning
family
Response type: "gaussian" (default), "binomial", or "cox" (for censored survival data)
sigma
Estimate of error standard deviation. If NULL (default), this is estimated
using the mean squared residual of the full least squares fit when n >= 2p, and
using the standard deviation of y when n < 2p. In the latter case, the user
should use estimateSigma function for a more accurate estimate.
Not used for family= "binomial", or "cox"
alpha
Significance level for confidence intervals (target is miscoverage alpha/2 in each tail)
intercept
Was the lasso problem solved (e.g., by glmnet) with an intercept in the model? Default is TRUE. Must be TRUE for "binomial" family. Not used for 'cox" family, where no intercept is assumed.
add.targets
Optional vector of predictors to be included as targets of inference, regardless of whether or not they are selected by the lasso. Default is NULL.
status
Censoring status for Cox model; 1=failurem 0=censored
type
Contrast type for p-values and confidence intervals: default is "partial"—meaning that the contrasts tested are the partial population regression coefficients, within the active set of predictors; the alternative is "full"—meaning that the full population regression coefficients are tested. The latter does not make sense when p > n.
tol.beta
Tolerance for determining if a coefficient is zero
tol.kkt
Tolerance for determining if an entry of the subgradient is zero
gridrange
Grid range for constructing confidence intervals, on the standardized scale
bits
Number of bits to be used for p-value and confidence interval calculations. Default is
NULL, in which case standard floating point calculations are performed. When not NULL,
multiple precision floating point calculations are performed with the specified number
of bits, using the R package Rmpfr (if this package is not installed, then a
warning is thrown, and standard floating point calculations are pursued).
Note: standard double precision uses 53 bits
so, e.g., a choice of 200 bits uses about 4 times double precision. The confidence
interval computation is sometimes numerically challenging, and the extra precision can be
helpful (though computationally more costly). In particular, extra precision might be tried
if the values in the output columns of tailarea differ noticeably from alpha/2.
verbose
Print out progress along the way? Default is FALSE
linesearch.try
When running type="full" (i.e. debiased LASSO) how many attempts in the line search?
Details
This function computes selective p-values and confidence intervals for the lasso,
given a fixed value of the tuning parameter lambda.
Three different response types are supported: gaussian, binomial and Cox.
The confidence interval construction involves numerical search and can be fragile:
if the observed statistic is too close to either end of the truncation interval
(vlo and vup, see references), then one or possibly both endpoints of the interval of
desired coverage cannot be computed, and default to +/- Inf. The output tailarea
gives the achieved Gaussian tail areas for the reported intervals—these should be close
to alpha/2, and can be used for error-checking purposes.
Important!: Before running glmnet (or some other lasso-solver) x should be centered, that is x <- scale(X,TRUE,FALSE). In addition, if standardization of the predictors is desired, x should be scaled as well: x <- scale(x,TRUE,TRUE). Then when running glmnet, set standardize=F. See example below.
The penalty.factor facility in glmmet– allowing different penalties lambda for each predictor, is not yet implemented in fixedLassoInf. However you can finesse this— see the example below. One caveat- using this approach, a penalty factor of zero (forcing a predictor in) is not allowed.
Note that the coefficients and standard errors reported are unregularized. Eg for the Gaussian, they are the usual least squares estimates and standard errors for the model fit to the active set from the lasso.
Value
type
Type of coefficients tested (partial or full)
lambda
Value of tuning parameter lambda used
pv
One-sided P-values for active variables, uses the fact we have conditioned on the sign.
ci
Confidence intervals
tailarea
Realized tail areas (lower and upper) for each confidence interval
vlo
Lower truncation limits for statistics
vup
Upper truncation limits for statistics
vmat
Linear contrasts that define the observed statistics
y
Vector of outcomes
vars
Variables in active set
sign
Signs of active coefficients
alpha
Desired coverage (alpha/2 in each tail)
sigma
Value of error standard deviation (sigma) used
call
The call to fixedLassoInf
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid
References
Jason Lee, Dennis Sun, Yuekai Sun, and Jonathan Taylor (2013). Exact post-selection inference, with application to the lasso. arXiv:1311.6238.
Jonathan Taylor and Robert Tibshirani (2016) Post-selection inference for L1-penalized likelihood models. arXiv:1602.07358
Examples
set.seed(43)
n = 50
p = 10
sigma = 1
x = matrix(rnorm(n*p),n,p)
x = scale(x,TRUE,TRUE)
beta = c(3,2,rep(0,p-2))
y = x%*%beta + sigma*rnorm(n)
# first run glmnet
gfit = glmnet(x,y,standardize=FALSE)
# extract coef for a given lambda; note the 1/n factor!
# (and we don't save the intercept term)
lambda = .8
beta = coef(gfit, x=x, y=y, s=lambda/n, exact=TRUE)[-1]
# compute fixed lambda p-values and selection intervals
out = fixedLassoInf(x,y,beta,lambda,sigma=sigma)
out
## as above, but use lar function instead to get initial
## lasso fit (should get same results)
lfit = lar(x,y,normalize=FALSE)
beta = coef(lfit, s=lambda, mode="lambda")
out2 = fixedLassoInf(x, y, beta, lambda, sigma=sigma)
out2
## mimic different penalty factors by first scaling x
set.seed(43)
n = 50
p = 10
sigma = 1
x = matrix(rnorm(n*p),n,p)
x=scale(x,TRUE,TRUE)
beta = c(3,2,rep(0,p-2))
y = x%*%beta + sigma*rnorm(n)
pf=c(rep(1,7),rep(.1,3)) #define penalty factors
pf=p*pf/sum(pf) # penalty factors should be rescaled so they sum to p
xs=scale(x,FALSE,pf) #scale cols of x by penalty factors
# first run glmnet
gfit = glmnet(xs, y, standardize=FALSE)
# extract coef for a given lambda; note the 1/n factor!
# (and we don't save the intercept term)
lambda = .8
beta_hat = coef(gfit, x=xs, y=y, s=lambda/n, exact=TRUE)[-1]
# compute fixed lambda p-values and selection intervals
out = fixedLassoInf(xs,y,beta_hat,lambda,sigma=sigma)
#rescale conf points to undo the penalty factor
out$ci=t(scale(t(out$ci),FALSE,pf[out$vars]))
out
#logistic model
set.seed(43)
n = 50
p = 10
sigma = 1
x = matrix(rnorm(n*p),n,p)
x=scale(x,TRUE,TRUE)
beta = c(3,2,rep(0,p-2))
y = x%*%beta + sigma*rnorm(n)
y=1*(y>mean(y))
# first run glmnet
gfit = glmnet(x,y,standardize=FALSE,family="binomial")
# extract coef for a given lambda; note the 1/n factor!
# (and here we DO include the intercept term)
lambda = .8
beta_hat = coef(gfit, x=x, y=y, s=lambda/n, exact=TRUE)
# compute fixed lambda p-values and selection intervals
out = fixedLassoInf(x,y,beta_hat,lambda,family="binomial")
out
# Cox model
set.seed(43)
n = 50
p = 10
sigma = 1
x = matrix(rnorm(n*p), n, p)
x = scale(x, TRUE, TRUE)
beta = c(3,2,rep(0,p-2))
tim = as.vector(x%*%beta + sigma*rnorm(n))
tim= tim-min(tim)+1
status=sample(c(0,1),size=n,replace=TRUE)
# first run glmnet
y = Surv(tim,status)
gfit = glmnet(x, y, standardize=FALSE, family="cox")
# extract coef for a given lambda; note the 1/n factor!
lambda = 1.5
beta_hat = as.numeric(coef(gfit, x=x, y=y, s=lambda/n, exact=TRUE))
# compute fixed lambda p-values and selection intervals
out = fixedLassoInf(x, tim, beta_hat, lambda, status=status, family="cox")
out
# Debiased lasso or "full"
n = 50
p = 100
sigma = 1
x = matrix(rnorm(n*p),n,p)
x = scale(x,TRUE,TRUE)
beta = c(3,2,rep(0,p-2))
y = x%*%beta + sigma*rnorm(n)
# first run glmnet
gfit = glmnet(x, y, standardize=FALSE, intercept=FALSE)
# extract coef for a given lambda; note the 1/n factor!
# (and we don't save the intercept term)
lambda = 2.8
beta = coef(gfit, x=x, y=y, s=lambda/n, exact=TRUE)[-1]
# compute fixed lambda p-values and selection intervals
out = fixedLassoInf(x, y, beta, lambda, sigma=sigma, type='full', intercept=FALSE)
out
# When n > p and "full" we use the full inverse
# instead of Javanmard and Montanari's approximate inverse
n = 200
p = 50
sigma = 1
x = matrix(rnorm(n*p),n,p)
x = scale(x,TRUE,TRUE)
beta = c(3,2,rep(0,p-2))
y = x%*%beta + sigma*rnorm(n)
# first run glmnet
gfit = glmnet(x, y, standardize=FALSE, intercept=FALSE)
# extract coef for a given lambda; note the 1/n factor!
# (and we don't save the intercept term)
lambda = 2.8
beta = coef(gfit, x=x, y=y, s=lambda/n, exact=TRUE)[-1]
# compute fixed lambda p-values and selection intervals
out = fixedLassoInf(x, y, beta, lambda, sigma=sigma, type='full', intercept=FALSE)
out
ForwardStop rule for sequential p-values
Description
Computes the ForwardStop sequential stopping rule of G'Sell et al (2014)
Usage
forwardStop(pv, alpha=0.1)
Arguments
pv
Vector of **sequential** p-values, for example from fsInf or larInf
alpha
Desired type FDR level (between 0 and 1)
Details
Computes the ForwardStop sequential stopping rule of G'Sell et al (2014). Guarantees FDR control at the level alpha, for independent p-values.
Value
Step number for sequential stop.
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid
References
Max Grazier G'Sell, Stefan Wager, Alexandra Chouldechova, and Rob Tibshirani (2014). Sequential selection procedures and Fflse Discovery Rate Control. arXiv:1309.5352. To appear in Journal of the Royal Statistical Society: Series B.
Examples
set.seed(33)
n = 50
p = 10
sigma = 1
x = matrix(rnorm(n*p),n,p)
beta = c(3,2,rep(0,p-2))
y = x%*%beta + sigma*rnorm(n)
# run forward stepwise
fsfit = fs(x,y)
# compute sequential p-values and confidence intervals
# (sigma estimated from full model)
out = fsInf(fsfit)
out
# estimate optimal stopping point
forwardStop(out$pv, alpha=.10)
Forward stepwise regression
Description
This function implements forward stepwise regression, for use in the selectiveInference package
Usage
fs(x, y, maxsteps=2000, intercept=TRUE, normalize=TRUE, verbose=FALSE)
Arguments
x
Matrix of predictors (n by p)
y
Vector of outcomes (length n)
maxsteps
Maximum number of steps to take
intercept
Should an intercept be included on the model? Default is TRUE
normalize
Should the predictors be normalized? Default is TRUE. (Note:
this argument has no real effect on model selection since forward stepwise is
scale invariant already; however, it is included for completeness, and to match
the interface for the lar function)
verbose
Print out progress along the way? Default is FALSE
Details
This function implements forward stepwise regression, adding the predictor at each
step that maximizes the absolute correlation between the predictors—once
orthogonalized with respect to the current model—and the residual. This entry
criterion is standard, and is equivalent to choosing the variable that achieves
the biggest drop in RSS at each step; it is used, e.g., by the step function
in R. Note that, for example, the lars package implements a stepwise option
(with type="step"), but uses a (mildly) different entry criterion, based on maximal
absolute correlation between the original (non-orthogonalized) predictors and the
residual.
Value
action
Vector of predictors in order of entry
sign
Signs of coefficients of predictors, upon entry
df
Degrees of freedom of each active model
beta
Matrix of regression coefficients for each model along the path, one column per model
completepath
Was the complete stepwise path computed?
bls
If completepath is TRUE, the full least squares coefficients
Gamma
Matrix that captures the polyhedral selection at each step
nk
Number of polyhedral constraints at each step in path
vreg
Matrix of linear contrasts that gives coefficients of variables to enter along the path
x
Matrix of predictors used
y
Vector of outcomes used
bx
Vector of column means of original x
by
Mean of original y
sx
Norm of each column of original x
intercept
Was an intercept included?
normalize
Were the predictors normalized?
call
The call to fs
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid
See Also
fsInf , predict.fs ,coef.fs , plot.fs
Examples
set.seed(33)
n = 50
p = 10
sigma = 1
x = matrix(rnorm(n*p),n,p)
beta = c(3,2,rep(0,p-2))
y = x%*%beta + sigma*rnorm(n)
# run forward stepwise, plot results
fsfit = fs(x,y)
plot(fsfit)
# compute sequential p-values and confidence intervals
# (sigma estimated from full model)
out = fsInf(fsfit)
out
Selective inference for forward stepwise regression
Description
Computes p-values and confidence intervals for forward stepwise regression
Usage
fsInf(obj, sigma=NULL, alpha=0.1, k=NULL, type=c("active","all","aic"),
gridrange=c(-100,100), bits=NULL, mult=2, ntimes=2, verbose=FALSE)
Arguments
obj
Object returned by fs function
sigma
Estimate of error standard deviation. If NULL (default), this is estimated
using the mean squared residual of the full least squares fit when n >= 2p, and
using the standard deviation of y when n < 2p. In the latter case, the user
should use estimateSigma function for a more accurate estimate
alpha
Significance level for confidence intervals (target is miscoverage alpha/2 in each tail)
k
See "type" argument below. Default is NULL, in which case k is taken to be the the number of steps computed in the forward stepwise path
type
Type of analysis desired: with "active" (default), p-values and confidence intervals are computed for each predictor as it is entered into the active step, all the way through k steps; with "all", p-values and confidence intervals are computed for all variables in the active model after k steps; with "aic", the number of steps k is first estimated using a modified AIC criterion, and then the same type of analysis as in "all" is carried out for this particular value of k.
Note that the AIC scheme is defined to choose a number of steps k after which the AIC criterion
increases ntimes in a row, where ntimes can be specified by the user (see below).
Under this definition, the AIC selection event is characterizable as a polyhedral set, and hence
the extra conditioning can be taken into account exactly. Also note that an analogous BIC scheme
can be specified through the mult argument (see below)
gridrange
Grid range for constructing confidence intervals, on the standardized scale
bits
Number of bits to be used for p-value and confidence interval calculations. Default is
NULL, in which case standard floating point calculations are performed. When not NULL,
multiple precision floating point calculations are performed with the specified number
of bits, using the R package Rmpfr (if this package is not installed, then a
warning is thrown, and standard floating point calculations are pursued).
Note: standard double precision uses 53 bits
so, e.g., a choice of 200 bits uses about 4 times double precision. The confidence
interval computation is sometimes numerically challenging, and the extra precision can be
helpful (though computationally more costly). In particular, extra precision might be tried
if the values in the output columns of tailarea differ noticeably from alpha/2.
mult
Multiplier for the AIC-style penalty. Hence a value of 2 (default) gives AIC, whereas a value of log(n) would give BIC
ntimes
Number of steps for which AIC-style criterion has to increase before minimizing point is declared
verbose
Print out progress along the way? Default is FALSE
Details
This function computes selective p-values and confidence intervals (selection intervals)
for forward stepwise regression. The default is to report the results for
each predictor after its entry into the model. See the "type" argument for other options.
The confidence interval construction involves numerical search and can be fragile:
if the observed statistic is too close to either end of the truncation interval
(vlo and vup, see references), then one or possibly both endpoints of the interval of
desired coverage cannot be computed, and default to +/- Inf. The output tailarea
gives the achieved Gaussian tail areas for the reported intervals—these should be close
to alpha/2, and can be used for error-checking purposes.
Value
type
Type of analysis (active, all, or aic)
k
Value of k specified in call
khat
When type is "active", this is an estimated stopping point
declared by forwardStop ; when type is "aic", this is the
value chosen by the modified AIC scheme
pv
One sided P-values for active variables, uses the sign that a variable entered the model with.
ci
Confidence intervals
tailarea
Realized tail areas (lower and upper) for each confidence interval
vlo
Lower truncation limits for statistics
vup
Upper truncation limits for statistics
vmat
Linear contrasts that define the observed statistics
y
Vector of outcomes
vars
Variables in active set
sign
Signs of active coefficients
alpha
Desired coverage (alpha/2 in each tail)
sigma
Value of error standard deviation (sigma) used
call
The call to fsInf
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid
References
Ryan Tibshirani, Jonathan Taylor, Richard Lockhart, and Rob Tibshirani (2014). Exact post-selection inference for sequential regression procedures. arXiv:1401.3889.
Joshua Loftus and Jonathan Taylor (2014). A significance test for forward stepwise model selection. arXiv:1405.3920.
See Also
Examples
set.seed(33)
n = 50
p = 10
sigma = 1
x = matrix(rnorm(n*p),n,p)
beta = c(3,2,rep(0,p-2))
y = x%*%beta + sigma*rnorm(n)
# run forward stepwise
fsfit = fs(x,y)
# compute sequential p-values and confidence intervals
# (sigma estimated from full model)
out.seq = fsInf(fsfit)
out.seq
# compute p-values and confidence intervals after AIC stopping
out.aic = fsInf(fsfit,type="aic")
out.aic
# compute p-values and confidence intervals after 5 fixed steps
out.fix = fsInf(fsfit,type="all",k=5)
out.fix
Select a model with forward stepwise.
Description
This function implements forward selection of linear models almost identically to step with direction = "forward". The reason this is a separate function from fs is that groups of variables (e.g. dummies encoding levels of a categorical variable) must be handled differently in the selective inference framework.
Usage
groupfs(x, y, index, maxsteps, sigma = NULL, k = 2, intercept = TRUE,
center = TRUE, normalize = TRUE, aicstop = 0, verbose = FALSE)
Arguments
x
Matrix of predictors (n by p).
y
Vector of outcomes (length n).
index
Group membership indicator of length p. Check that sort(unique(index)) = 1:G where G is the number of distinct groups.
maxsteps
Maximum number of steps for forward stepwise.
sigma
Estimate of error standard deviation for use in AIC criterion. This determines the relative scale between RSS and the degrees of freedom penalty. Default is NULL corresponding to unknown sigma. When NULL, link{groupfsInf} performs truncated F inference instead of truncated \chi. See extractAIC for details on the AIC criterion.
k
Multiplier of model size penalty, the default is k = 2 for AIC. Use k = log(n) for BIC, or k = 2log(p) for RIC (best for high dimensions, when p > n). If G < p then RIC may be too restrictive and it would be better to use log(G) < k < 2log(p).
intercept
Should an intercept be included in the model? Default is TRUE. Does not count as a step.
center
Should the columns of the design matrix be centered? Default is TRUE.
normalize
Should the design matrix be normalized? Default is TRUE.
aicstop
Early stopping if AIC increases. Default is 0 corresponding to no early stopping. Positive integer values specify the number of times the AIC is allowed to increase in a row, e.g. with aicstop = 2 the algorithm will stop if the AIC criterion increases for 2 steps in a row. The default of step corresponds to aicstop = 1.
verbose
Print out progress along the way? Default is FALSE.
Value
An object of class "groupfs" containing information about the sequence of models in the forward stepwise algorithm. Call the function groupfsInf on this object to compute selective p-values.
See Also
Examples
x = matrix(rnorm(20*40), nrow=20)
index = sort(rep(1:20, 2))
y = rnorm(20) + 2 * x[,1] - x[,4]
fit = groupfs(x, y, index, maxsteps = 5)
out = groupfsInf(fit)
out
Compute selective p-values for a model fitted by groupfs.
Description
Computes p-values for each group of variables in a model fitted by groupfs . These p-values adjust for selection by truncating the usual \chi^2 statistics to the regions implied by the model selection event. If the sigma to groupfs was NULL then groupfsInf uses truncated F statistics instead of truncated \chi. The sigma argument to groupfsInf allows users to override and use \chi, but this is not recommended unless \sigma can be estimated well (i.e. n > p).
Usage
groupfsInf(obj, sigma = NULL, verbose = TRUE)
Arguments
obj
Object returned by groupfs function
sigma
Estimate of error standard deviation. Default is NULL and in this case groupfsInf uses the value of sigma specified to groupfs .
verbose
Print out progress along the way? Default is TRUE.
Value
An object of class "groupfsInf" containing selective p-values for the fitted model obj. For comparison with fsInf , note that the option type = "active" is not available.
- vars
Labels of the active groups in the order they were included.
- pv
Selective p-values computed from appropriate truncated distributions.
- sigma
Estimate of error variance used in computing p-values.
- TC or TF
Observed value of truncated
\chiorF.- df
Rank of group of variables when it was added to the model.
- support
List of intervals defining the truncation region of the corresponding statistic.
Least angle regression
Description
This function implements least angle regression, for use in the selectiveInference package
Usage
lar(x, y, maxsteps=2000, minlam=0, intercept=TRUE, normalize=TRUE,
verbose=FALSE)
Arguments
x
Matrix of predictors (n by p)
y
Vector of outcomes (length n)
maxsteps
Maximum number of steps to take
minlam
Minimum value of lambda to consider
intercept
Should an intercept be included on the model? Default is TRUE
normalize
Should the predictors be normalized? Default is TRUE
verbose
Print out progress along the way? Default is FALSE
Details
The least angle regression algorithm is described in detail by Efron et al. (2002).
This function should match (in terms of its output) that from the lars package,
but returns additional information (namely, the polyhedral constraints) needed for the
selective inference calculations.
Value
lambda
Values of lambda (knots) visited along the path
action
Vector of predictors in order of entry
sign
Signs of coefficients of predictors, upon entry
df
Degrees of freedom of each active model
beta
Matrix of regression coefficients for each model along the path, one model per column
completepath
Was the complete stepwise path computed?
bls
If completepath is TRUE, the full least squares coefficients
Gamma
Matrix that captures the polyhedral selection at each step
nk
Number of polyhedral constraints at each step in path
vreg
Matrix of linear contrasts that gives coefficients of variables to enter along the path
mp
Value of M+ (for internal use with the spacing test)
x
Matrix of predictors used
y
Vector of outcomes used
bx
Vector of column means of original x
by
Mean of original y
sx
Norm of each column of original x
intercept
Was an intercept included?
normalize
Were the predictors normalized?
call
The call to lar
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Max G'Sell, Joshua Loftus, Stephen Reid
References
Brad Efron, Trevor Hastie, Iain Johnstone, and Rob Tibshirani (2002). Least angle regression. Annals of Statistics (with discussion).
See also the descriptions in Trevor Hastie, Rob Tibshirani, and Jerome Friedman (2002, 2009). Elements of Statistical Learning.
See Also
larInf , predict.lar , coef.lar , plot.lar
Examples
set.seed(43)
n = 50
p = 10
sigma = 1
x = matrix(rnorm(n*p),n,p)
beta = c(3,2,rep(0,p-2))
y = x%*%beta + sigma*rnorm(n)
# run LAR, plot results
larfit = lar(x,y)
plot(larfit)
# compute sequential p-values and confidence intervals
# (sigma estimated from full model)
out = larInf(larfit)
out
Selective inference for least angle regression
Description
Computes p-values and confidence intervals for least angle regression
Usage
larInf(obj, sigma=NULL, alpha=0.1, k=NULL, type=c("active","all","aic"),
gridrange=c(-100,100), bits=NULL, mult=2, ntimes=2, verbose=FALSE)
Arguments
obj
Object returned by lar function (not the lars function!)
sigma
Estimate of error standard deviation. If NULL (default), this is estimated
using the mean squared residual of the full least squares fit when n >= 2p, and
using the standard deviation of y when n < 2p. In the latter case, the user
should use estimateSigma function for a more accurate estimate
alpha
Significance level for confidence intervals (target is miscoverage alpha/2 in each tail)
k
See "type" argument below. Default is NULL, in which case k is taken to be the the number of steps computed in the least angle regression path
type
Type of analysis desired: with "active" (default), p-values and confidence intervals are computed for each predictor as it is entered into the active step, all the way through k steps; with "all", p-values and confidence intervals are computed for all variables in the active model after k steps; with "aic", the number of steps k is first estimated using a modified AIC criterion, and then the same type of analysis as in "all" is carried out for this particular value of k.
Note that the AIC scheme is defined to choose a number of steps k after which the AIC criterion
increases ntimes in a row, where ntimes can be specified by the user (see below).
Under this definition, the AIC selection event is characterizable as a polyhedral set, and hence
the extra conditioning can be taken into account exactly. Also note that an analogous BIC scheme
can be specified through the mult argument (see below)
gridrange
Grid range for constructing confidence intervals, on the standardized scale
bits
Number of bits to be used for p-value and confidence interval calculations. Default is
NULL, in which case standard floating point calculations are performed. When not NULL,
multiple precision floating point calculations are performed with the specified number
of bits, using the R package Rmpfr (if this package is not installed, then a
warning is thrown, and standard floating point calculations are pursued).
Note: standard double precision uses 53 bits
so, e.g., a choice of 200 bits uses about 4 times double precision. The confidence
interval computation is sometimes numerically challenging, and the extra precision can be
helpful (though computationally more costly). In particular, extra precision might be tried
if the values in the output columns of tailarea differ noticeably from alpha/2.
mult
Multiplier for the AIC-style penalty. Hence a value of 2 (default) gives AIC, whereas a value of log(n) would give BIC
ntimes
Number of steps for which AIC-style criterion has to increase before minimizing point is declared
verbose
Print out progress along the way? Default is FALSE
Details
This function computes selective p-values and confidence intervals (selection intervals)
for least angle regression. The default is to report the results for
each predictor after its entry into the model. See the "type" argument for other options.
The confidence interval construction involves numerical search and can be fragile:
if the observed statistic is too close to either end of the truncation interval
(vlo and vup, see references), then one or possibly both endpoints of the interval of
desired coverage cannot be computed, and default to +/- Inf. The output tailarea
gives the achieved Gaussian tail areas for the reported intervals—these should be close
to alpha/2, and can be used for error-checking purposes.
Value
type
Type of analysis (active, all, or aic)
k
Value of k specified in call
khat
When type is "active", this is an estimated stopping point
declared by forwardStop ; when type is "aic", this is the
value chosen by the modified AIC scheme
pv
P-values for active variables
ci
Confidence intervals
tailarea
Realized tail areas (lower and upper) for each confidence interval
vlo
Lower truncation limits for statistics
vup
Upper truncation limits for statistics
vmat
Linear contrasts that define the observed statistics
y
Vector of outcomes
pv.spacing
P-values from the spacing test (here M+ is used)
pv.modspac
P-values from the modified form of the spacing test (here M+ is replaced by the next knot)
pv.covtest
P-values from covariance test
vars
Variables in active set
sign
Signs of active coefficients
alpha
Desired coverage (alpha/2 in each tail)
sigma
Value of error standard deviation (sigma) used
call
The call to larInf
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid
References
Ryan Tibshirani, Jonathan Taylor, Richard Lockhart, and Rob Tibshirani (2014). Exact post-selection inference for sequential regression procedures. arXiv:1401.3889.
See Also
Examples
set.seed(43)
n = 50
p = 10
sigma = 1
x = matrix(rnorm(n*p),n,p)
beta = c(3,2,rep(0,p-2))
y = x%*%beta + sigma*rnorm(n)
# run LAR
larfit = lar(x,y)
# compute sequential p-values and confidence intervals
# (sigma estimated from full model)
out.seq = larInf(larfit)
out.seq
# compute p-values and confidence intervals after AIC stopping
out.aic = larInf(larfit,type="aic")
out.aic
# compute p-values and confidence intervals after 5 fixed steps
out.fix = larInf(larfit,type="all",k=5)
out.fix
Selective inference for many normal means
Description
Computes p-values and confidence intervals for the largest k among many normal means
Usage
manyMeans(y, alpha=0.1, bh.q=NULL, k=NULL, sigma=1, verbose=FALSE)
Arguments
y
Vector of outcomes (length n)
alpha
Significance level for confidence intervals (target is miscoverage alpha/2 in each tail)
bh.q
q parameter for BH(q) procedure
k
Number of means to consider
sigma
Estimate of error standard deviation
verbose
Print out progress along the way? Default is FALSE
Details
This function compute p-values and confidence intervals for the largest k among many normal means. One can specify a fixed number of means k to consider, or choose the number to consider via the BH rule.
Value
mu.hat
Vector of length n containing the estimated signal sizes. If a sample element is not selected, then its signal size estimate is 0
selected.set
Indices of the vector y of the sample elements that were selected by the procedure (either BH(q) or top-K). Labelled "Selind" in output table.
pv
P-values for selected signals
ci
Confidence intervals
method
Method used to choose number of means
sigma
Value of error standard deviation (sigma) used
bh.q
BH q-value used
k
Desired number of means
threshold
Computed cutoff
call
The call to manyMeans
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid
References
Stephen Reid, Jonathan Taylor, and Rob Tibshirani (2014). Post-selection point and interval estimation of signal sizes in Gaussian samples. arXiv:1405.3340.
Examples
set.seed(12345)
n = 100
mu = c(rep(3,floor(n/5)), rep(0,n-floor(n/5)))
y = mu + rnorm(n)
out = manyMeans(y, bh.q=0.1)
out
Plot function for forward stepwise regression
Description
Plot coefficient profiles along the forward stepwise path
Usage
## S3 method for class 'fs'
plot(x, breaks=TRUE, omit.zeros=TRUE, var.labels=TRUE, ...)
Arguments
x
Object returned by a call to fs function
breaks
Should vertical lines be drawn at each break point in the piecewise linear coefficient paths? Default is TRUE
omit.zeros
Should segments of the coefficients paths that are equal to zero be omitted (to avoid clutter in the figure)? Default is TRUE
var.labels
Should paths be labelled with corresponding variable numbers? Default is TRUE
...
Additional arguments for plotting
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid
Examples
set.seed(33)
n = 50
p = 10
sigma = 1
x = matrix(rnorm(n*p),n,p)
beta = c(3,2,rep(0,p-2))
y = x%*%beta + sigma*rnorm(n)
# run forward stepwise, plot results
fsfit = fs(x,y)
plot(fsfit)
Plot function for least angle regression
Description
Plot coefficient profiles along the LAR path
Usage
## S3 method for class 'lar'
plot(x, xvar=c("norm","step","lambda"), breaks=TRUE,
omit.zeros=TRUE, var.labels=TRUE, ...)
Arguments
x
Object returned by a call to lar function
(not the lars function!)
xvar
Either "norm" or "step" or "lambda", determining what is plotted on the x-axis
breaks
Should vertical lines be drawn at each break point in the piecewise linear coefficient paths? Default is TRUE
omit.zeros
Should segments of the coefficients paths that are equal to zero be omitted (to avoid clutter in the figure)? Default is TRUE
var.labels
Should paths be labelled with corresponding variable numbers? Default is TRUE
...
Additional arguments for plotting
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid
Examples
set.seed(43)
n = 50
p = 10
sigma = 1
x = matrix(rnorm(n*p),n,p)
beta = c(3,2,rep(0,p-2))
y = x%*%beta + sigma*rnorm(n)
# run LAR, plot results
larfit = lar(x,y)
plot(larfit)
Prediction and coefficient functions for forward stepwise regression
Description
Make predictions or extract coefficients from a forward stepwise object
Usage
## S3 method for class 'fs'
predict(object, newx, s, ...)
## S3 method for class 'fs'
coef(object, s, ...)
Arguments
object
Object returned by a call to fs function
newx
Matrix of x values at which the predictions are desired. If NULL, the x values from forward stepwise fitting are used
s
Step number(s) at which predictions or coefficients are desired
...
Additional arguments
Value
Either a vector/matrix of predictions, or a vector/matrix of coefficients.
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid
Examples
set.seed(33)
n = 200
p = 20
sigma = 1
x = matrix(rnorm(n*p),n,p)
beta = c(rep(3,10),rep(0,p-10))
y = x%*%beta + sigma*rnorm(n)
# run forward stepwise and predict functions
obj = fs(x,y)
fit = predict(obj,x,s=3)
Prediction and coefficient functions for groupfs .
Description
Make predictions or extract coefficients from a groupfs forward stepwise object.
Usage
## S3 method for class 'groupfs'
predict(object, newx)
Arguments
object
Object returned by a call to groupfs .
newx
Matrix of x values at which the predictions are desired. If NULL, the x values from groupfs fitting are used.
Value
A vector of predictions or a vector of coefficients.
Prediction and coefficient functions for least angle regression
Description
Make predictions or extract coefficients from a least angle regression object
Usage
## S3 method for class 'lar'
predict(object, newx, s, mode=c("step","lambda"), ...)
## S3 method for class 'lar'
coef(object, s, mode=c("step","lambda"), ...)
Arguments
object
Object returned by a call to lar function
(not the lars function!)
newx
Matrix of x values at which the predictions are desired. If NULL, the x values from least angle regression fitting are used
s
Step number(s) or lambda value(s) at which predictions or coefficients are desired
mode
Either "step" or "lambda", determining the role of s (above)
...
Additional arguments
Value
Either a vector/matrix of predictions, or a vector/matrix of coefficients.
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid
Examples
set.seed(33)
n = 200
p = 20
sigma = 1
x = matrix(rnorm(n*p),n,p)
beta = c(rep(3,10),rep(0,p-10))
y = x%*%beta + sigma*rnorm(n)
# run lar and predict functions
obj = lar(x,y)
fit = predict(obj,x,s=3)
Inference for the randomized lasso, with a fixed lambda
Description
Solve a randomly perturbed LASSO problem.
Usage
randomizedLasso(X,
y,
lam,
family=c("gaussian", "binomial"),
noise_scale=NULL,
ridge_term=NULL,
max_iter=100,
kkt_tol=1.e-4,
parameter_tol=1.e-8,
objective_tol=1.e-8,
objective_stop=FALSE,
kkt_stop=TRUE,
parameter_stop=TRUE)
Arguments
X
Matrix of predictors (n by p);
y
Vector of outcomes (length n)
lam
Value of lambda used to compute beta. See the above warning Be careful! This function uses the "standard" lasso objective
1/2 \|y - x \beta\|_2^2 + \lambda \|\beta\|_1.
In contrast, glmnet multiplies the first term by a factor of 1/n.
So after running glmnet, to extract the beta corresponding to a value lambda,
you need to use beta = coef(obj, s=lambda/n)[-1],
where obj is the object returned by glmnet (and [-1] removes the intercept,
which glmnet always puts in the first component)
family
Response type: "gaussian" (default), "binomial".
noise_scale
Scale of Gaussian noise added to objective. Default is 0.5 * sd(y) times the sqrt of the mean of the trace of X^TX.
ridge_term
A small "elastic net" or ridge penalty is added to ensure the randomized problem has a solution. 0.5 * sd(y) times the sqrt of the mean of the trace of X^TX divided by sqrt(n).
max_iter
How many rounds of updates used of coordinate descent in solving randomized LASSO.
kkt_tol
Tolerance for checking convergence based on KKT conditions.
parameter_tol
Tolerance for checking convergence based on convergence of parameters.
objective_tol
Tolerance for checking convergence based on convergence of objective value.
kkt_stop
Should we use KKT check to determine when to stop?
parameter_stop
Should we use convergence of parameters to determine when to stop?
objective_stop
Should we use convergence of objective value to determine when to stop?
Details
For family="gaussian" this function uses the "standard" lasso objective
1/2 \|y - x \beta\|_2^2 + \lambda \|\beta\|_1
and adds a term
- \omega^T\beta + \frac{\epsilon}{2} \|\beta\|^2_2
where omega is drawn from IID normals with standard deviation
noise_scale and epsilon given by ridge_term.
See below for default values of noise_scale and ridge_term.
For family="binomial", the squared error loss is replaced by the
negative of the logistic log-likelihood.
Value
X
Design matrix.
y
Response vector.
lam
Vector of penalty parameters.
family
Family: "gaussian" or "binomial".
active_set
Set of non-zero coefficients in randomized solution that were penalized. Integers from 1:p.
inactive_set
Set of zero coefficients in randomized solution. Integers from 1:p.
unpenalized_set
Set of non-zero coefficients in randomized solution that were not penalized. Integers from 1:p.
sign_soln
The sign pattern of the randomized solution.
full_law
List describing sampling parameters for conditional law of all optimization variables given the data in the LASSO problem.
conditional_law
List describing sampling parameters for conditional law of only the scaling variables given the data and the observed subgradient in the LASSO problem.
internal_transform
Affine transformation describing relationship between internal representation of the data and the data compontent of score of the likelihood at the unregularized MLE based on the sign_vector (a.k.a. relaxed LASSO).
observed_raw
Data component of the score at the unregularized MLE.
noise_scale
SD of Gaussian noise used to draw the perturbed objective.
soln
The randomized solution. Inference is made conditional on its sign vector (so no more snooping of this value is formally permitted.)
If condition_subgrad == TRUE when sampling, then we may snoop on the observed subgradient.
perturb
The random vector in the linear term added to the objective.
Author(s)
Jelena Markovic, Jonathan Taylor
References
Xiaoying Tian, and Jonathan Taylor (2015). Selective inference with a randomized response. arxiv.org:1507.06739
Xiaoying Tian, Snigdha Panigrahi, Jelena Markovic, Nan Bi and Jonathan Taylor (2016). Selective inference after solving a convex problem. arxiv:1609.05609
Examples
set.seed(43)
n = 50
p = 10
sigma = 0.2
lam = 0.5
X = matrix(rnorm(n*p), n, p)
X = scale(X, TRUE, TRUE) / sqrt(n-1)
beta = c(3,2,rep(0,p-2))
y = X%*%beta + sigma*rnorm(n)
result = randomizedLasso(X, y, lam)
Inference for the randomized lasso, with a fixed lambda
Description
Compute p-values and confidence intervals based on selecting an active set with the randomized lasso, at a fixed value of the tuning parameter lambda and using Gaussian randomization.
Usage
randomizedLassoInf(rand_lasso_soln,
targets=NULL,
level=0.9,
sampler=c("norejection", "adaptMCMC"),
nsample=10000,
burnin=2000,
opt_samples=NULL)
Arguments
rand_lasso_soln
A randomized lasso solution as returned by randomizedLasso.
targets
If not NULL, should be a list with entries observed_target, cov_target, crosscov_target_internal.
The observed_target should be (pre-selection) asymptotically Gaussian around targeted
parameters. The quantity cov_target should be an estimate of the (pre-selection) covariance
of observed_target. Finally, crosscov_target_internal should be an estimate of
the (pre-selection) covariance of observed_target and the internal representation of the
data of the LASSO. For both "gaussian" and "binomial", this is the vector
\hat{\beta}_{E,MLE}, X_{-E}^T(y - \mu(X_E\hat{\beta}_{E,MLE}))
For example, this cross-covariance could be estimated by jointly bootstrapping the target of interest and the above vector.
level
Level for confidence intervals.
sampler
Which sampler to use – default is a no-rejection sampler. Otherwise use MCMC from the adaptMCMC package.
nsample
Number of samples of optimization variables to sample.
burnin
How many samples of optimization variable to discard (should be less than nsample).
opt_samples
Optional sample of optimization variables. If not NULL then no MCMC will be run.
Details
This function computes selective p-values and confidence intervals for a randomized version of the lasso, given a fixed value of the tuning parameter lambda.
Value
targets
A list with entries observed_target, cov_target, crosscov_target_internal. See argument description above.
pvalues
P-values testing hypotheses that each specific target is 0.
ci
Confidence interval for parameters determined by targets.
Author(s)
Jelena Markovic, Jonathan Taylor
References
Jelena Markovic and Jonathan Taylor (2016). Bootstrap inference after using multiple queries for model selection. arxiv.org:1612.07811
Xiaoying Tian and Jonathan Taylor (2015). Selective inference with a randomized response. arxiv.org:1507.06739
Xiaoying Tian, Snigdha Panigrahi, Jelena Markovic, Nan Bi and Jonathan Taylor (2016). Selective inference after solving a convex problem. arxiv.org:1609.05609
Examples
set.seed(43)
n = 50
p = 10
sigma = 0.2
lam = 0.5
X = matrix(rnorm(n*p), n, p)
X = scale(X, TRUE, TRUE) / sqrt(n-1)
beta = c(3,2,rep(0,p-2))
y = X%*%beta + sigma*rnorm(n)
result = randomizedLasso(X, y, lam)
inf_result = randomizedLassoInf(result)
Center and scale design matrix by groups
Description
For internal use by groupfs .
Usage
scaleGroups(x, index, center = TRUE, normalize = TRUE)
Arguments
x
Design matrix.
index
Group membership indicator of length p.
center
Center groups, default is TRUE.
normalize
Scale groups by Frobenius norm, default is TRUE.
Value
- x
Optionally centered/scaled design matrix.
- xm
Means of groups in original design matrix.
- xs
Frobenius norms of groups in original design matrix.
Tools for selective inference
Description
Functions to perform post-selection inference for forward stepwise regression, least angle regression, the lasso and the many normal means problem. The lasso function also supports logistic regression and the Cox model.
Details
This package provides tools for inference after selection, in forward stepwise regression, least angle regression, the lasso, and the many normal means problem. The functions compute p-values and selection intervals that properly account for the inherent selection carried out by the procedure. These have exact finite sample type I error and coverage under Gaussian errors. For the logistic and Cox familes (fixedLassoInf), the coverage is asymptotically valid
This R package was developed as part of the selective inference software project in Python and R:
https://github.com/selective-inference
Some of the R code in this work is a modification of Python code from this repository. Here is the current selective inference software team:
Yuval Benjamini, Leonard Blier, Will Fithian, Jason Lee, Joshua Loftus, Joshua Loftus, Stephen Reid, Dennis Sun, Yuekai Sun, Jonathan Taylor, Xiaoying Tian, Ryan Tibshirani, Rob Tibshirani
The main functions included in the package are:
fs ,
fsInf ,
lar ,
larInf ,
fixedLassoInf ,
manyMeans
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid
Maintainer: Rob Tibshirani <tibs@stanford.edu>
References
Ryan Tibshirani, Jonathan Taylor, Richard Lockhart, and Rob Tibshirani (2014). Exact post-selection inference for sequential regression procedures. arXiv:1401.3889.
Jason Lee, Dennis Sun, Yuekai Sun, and Jonathan Taylor (2013). Exact post-selection inference, with application to the lasso. arXiv:1311.6238.
Stephen Reid, Jonathan Taylor, and Rob Tibshirani (2014). Post-selection point and interval estimation of signal sizes in Gaussian samples. arXiv:1405.3340.
Jonathan Taylor and Robert Tibshirani (2016) Post-selection inference for L1-penalized likelihood models. arXiv:1602.07358
Examples
set.seed(33)
n = 50
p = 10
sigma = 1
x = matrix(rnorm(n*p),n,p)
beta = c(3,2,rep(0,p-2))
y = x%*%beta + sigma*rnorm(n)
# run forward stepwise
fsfit = fs(x,y)
# compute sequential p-values and confidence intervals
# (sigma estimated from full model)
out.seq = fsInf(fsfit)
out.seq
# compute p-values and confidence intervals after AIC stopping
out.aic = fsInf(fsfit,type="aic")
out.aic
# compute p-values and confidence intervals after 5 fixed steps
out.fix = fsInf(fsfit,type="all",k=5)
out.fix
## NOT RUN---lasso at fixed lambda- Gaussian family
## first run glmnet
# gfit = glmnet(x,y)
## extract coef for a given lambda; note the 1/n factor!
## (and we don't save the intercept term)
# lambda = .1
# beta = coef(gfit, s=lambda/n, exact=TRUE)[-1]
## compute fixed lambda p-values and selection intervals
# out = fixedLassoInf(x,y,beta,lambda,sigma=sigma)
# out
#lasso at fixed lambda- logistic family
#set.seed(43)
# n = 50
# p = 10
# sigma = 1
# x = matrix(rnorm(n*p),n,p)
x=scale(x,TRUE,TRUE)
#
# beta = c(3,2,rep(0,p-2))
# y = x%*%beta + sigma*rnorm(n)
# y=1*(y>mean(y))
# first run glmnet
# gfit = glmnet(x,y,standardize=FALSE,family="binomial")
# extract coef for a given lambda; note the 1/n factor!
# (and here we DO include the intercept term)
# lambda = .8
# beta = coef(gfit, s=lambda/n, exact=TRUE)
# # compute fixed lambda p-values and selection intervals
# out = fixedLassoInf(x,y,beta,lambda,family="binomial")
# out
##lasso at fixed lambda- Cox family
#set.seed(43)
# n = 50
# p = 10
# sigma = 1
# x = matrix(rnorm(n*p),n,p)
# x=scale(x,TRUE,TRUE)
# beta = c(3,2,rep(0,p-2))
# tim = as.vector(x%*%beta + sigma*rnorm(n))
# tim= tim-min(tim)+1
#status=sample(c(0,1),size=n,replace=T)
# first run glmnet
# gfit = glmnet(x,Surv(tim,status),standardize=FALSE,family="cox")
# extract coef for a given lambda; note the 1/n factor!
# lambda = 1.5
# beta = as.numeric(coef(gfit, s=lambda/n, exact=TRUE))
# compute fixed lambda p-values and selection intervals
# out = fixedLassoInf(x,tim,beta,lambda,status=status,family="cox")
# out
## NOT RUN---many normal means
# set.seed(12345)
# n = 100
# mu = c(rep(3,floor(n/5)), rep(0,n-floor(n/5)))
# y = mu + rnorm(n)
# out = manyMeans(y, bh.q=0.1)
# out
## NOT RUN---forward stepwise with groups
# set.seed(1)
# n = 20
# p = 40
# x = matrix(rnorm(n*p), nrow=n)
# index = sort(rep(1:(p/2), 2))
# y = rnorm(n) + 2 * x[,1] - x[,4]
# fit = groupfs(x, y, index, maxsteps = 5)
# out = groupfsInf(fit)
# out
## NOT RUN---estimation of sigma for use in fsInf
## (or larInf or fixedLassoInf)
# set.seed(33)
# n = 50
# p = 10
# sigma = 1
# x = matrix(rnorm(n*p),n,p)
# beta = c(3,2,rep(0,p-2))
# y = x%*%beta + sigma*rnorm(n)
## run forward stepwise
# fsfit = fs(x,y)
## estimate sigma
# sigmahat = estimateSigma(x,y)$sigmahat
## run sequential inference with estimated sigma
# out = fsInf(fit,sigma=sigmahat)
# out
Internal PMA functions
Description
Internal selectiveInference functions
Usage
## S3 method for class 'fs'
print(x, ...)
## S3 method for class 'fsInf'
print(x, tailarea = TRUE, ...)
## S3 method for class 'lar'
print(x,...)
## S3 method for class 'larInf'
print(x, tailarea = TRUE, ...)
## S3 method for class 'fixedLassoInf'
print(x, tailarea = TRUE, ...)
## S3 method for class 'manyMeans'
print(x, ...)
## S3 method for class 'ROSI'
print(x, ...)
Author(s)
Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor, Joshua Loftus, Stephen Reid