ashr (version 0.9)

ash.workhorse: Detailed Adaptive Shrinkage function

Description

Takes vectors of estimates (betahat) and their standard errors (sebetahat), and applies shrinkage to them, using Empirical Bayes methods, to compute shrunk estimates for beta. This is the more detailed version of ash for "research" use. Most users will be happy with the ash function, which provides the same usage, but documents only the main options for simplicity.

Usage

ash.workhorse(betahat, sebetahat, method = c("fdr", "shrink"), mixcompdist = c("uniform", "halfuniform", "normal"), lambda1 = 1, lambda2 = 0, nullcheck = TRUE, df = NULL, randomstart = FALSE, nullweight = 10, nonzeromode = FALSE, pointmass = TRUE, onlylogLR = FALSE, prior = c("nullbiased", "uniform"), mixsd = NULL, VB = FALSE, gridmult = sqrt(2), minimaloutput = FALSE, multiseqoutput = FALSE, g = NULL, cxx = FALSE, model = c("EE", "ES"), control = list())

Arguments

betahat

a p vector of estimates

sebetahat

a p vector of corresponding standard errors

method

specifies how ash is to be run. Can be "shrinkage" (if main aim is shrinkage) or "fdr" (if main aim is to assess fdr or fsr) This is simply a convenient way to specify certain combinations of parameters: "shrinkage" sets pointmass=FALSE and prior="uniform"; "fdr" sets pointmass=TRUE and prior="nullbiased".

mixcompdist

distribution of components in mixture ( "uniform","halfuniform" or "normal"), the default value is "uniform"

lambda1

multiplicative "inflation factor" for standard errors (like Genomic Control)

lambda2

additive "inflation factor" for standard errors (like Genomic Control)

nullcheck

whether to check that any fitted model exceeds the "null" likelihood in which all weight is on the first component

df

appropriate degrees of freedom for (t) distribution of betahat/sebetahat, default is NULL(Gaussian)

randomstart

logical, indicating whether to initialize EM randomly. If FALSE, then initializes to prior mean (for EM algorithm) or prior (for VBEM)

nullweight

scalar, the weight put on the prior under "nullbiased" specification, see prior

nonzeromode

logical, indicating whether to use a non-zero unimodal mixture(default is "FALSE")

pointmass

logical, indicating whether to use a point mass at zero as one of components for a mixture distribution

onlylogLR

logical, indicating whether to use this function to get logLR. Skip posterior prob, posterior mean, lfdr...

prior

string, or numeric vector indicating Dirichlet prior on mixture proportions (defaults to "uniform", or (1,1...,1); also can be "nullbiased" (nullweight,1,...,1) to put more weight on first component)

mixsd

vector of sds for underlying mixture components

VB

whether to use Variational Bayes to estimate mixture proportions (instead of EM to find MAP estimate), see mixVBEM and mixEM

gridmult

the multiplier by which the default grid values for mixsd differ by one another. (Smaller values produce finer grids)

minimaloutput

if TRUE, just outputs the fitted g and the lfsr (useful for very big data sets where memory is an issue)

multiseqoutput

if TRUE, just outputs the fitted g, logLR, PosteriorMean, PosteriorSD, function call and df

g

the prior distribution for beta (usually estimated from the data; this is used primarily in simulated data to do computations with the "true" g)

cxx

flag to indicate whether to use the c++ (Rcpp) version. After application of Squared extrapolation methods for accelerating fixed-point iterations (R Package "SQUAREM"), the c++ version is no longer faster than non-c++ version, thus we do not recommend using this one, and might be removed at any point.

model

c("EE","ES") specifies whether to assume exchangeable effects (EE) or exchangeable standardized effects (ES).

control

A list of control parameters for the SQUAREM algorithm, default value is set to be control.default=list(K = 1, method=3, square=TRUE, step.min0=1, step.max0=1, mstep=4, kr=1, objfn.inc=1,tol=1.e-07, maxiter=5000, trace=FALSE). User may supply changes to this list of parameter, say, control=list(maxiter=10000,trace=TRUE)

Value

ash returns an object of class "ash", a list with the following elements(or a simplified list, if $onlylogLR=TRUE$, $minimaloutput=TRUE$ or $multiseqoutput=TRUE$)

Details

See readme for more details.

See Also

ash for simplified specification of ash function

ashci for computation of credible intervals after getting the ash object return by ash()

ashm for Multi-model Adaptive Shrinkage function

Examples

beta = c(rep(0,100),rnorm(100))
sebetahat = abs(rnorm(200,0,1))
betahat = rnorm(200,beta,sebetahat)
beta.ash = ash(betahat, sebetahat)
summary(beta.ash)
plot(betahat,beta.ash$PosteriorMean,xlim=c(-4,4),ylim=c(-4,4))

CIMatrix=ashci(beta.ash,level=0.95)
print(CIMatrix)

#Testing the non-zero mode feature
betahat=betahat+5
beta.ash = ash(betahat, sebetahat)
plot(betahat,beta.ash$PosteriorMean)
summary(beta.ash)
betan.ash=ash(betahat, sebetahat,nonzeromode=TRUE)
plot(betahat, betan.ash$PosteriorMean)
summary(betan.ash)