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", "+uniform", "-uniform"), optmethod = c("mixIP", "cxxMixSquarem", "mixEM", "mixVBEM"), df = NULL, randomstart = FALSE, nullweight = 10, nonzeromode = FALSE, pointmass = TRUE, prior = c("nullbiased", "uniform", "unit"), mixsd = NULL, gridmult = sqrt(2), outputlevel = 2, g = NULL, fixg = FALSE, cxx = FALSE, VB = FALSE, model = c("EE", "ET"), 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","normal" or "+uniform"), the default value is "uniform" use "halfuniform" to allow for assymetric g, and "+uniform"/"-uniform" to constrain g to be positive/negative.

optmethod

specifies optimization method used. Default is "mixIP", an interior point method, if REBayes is installed; otherwise an EM algorithm is used. The interior point method is faster for large problems (n>2000).

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

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), or "unit" (1/K,...,1/K) [for optmethod=mixVBEM version only]

mixsd

vector of sds for underlying mixture components

gridmult

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

outputlevel

determines amount of output [0=just fitted g; 1=also PosteriorMean and PosteriorSD; 2= all but data; 3=also include copy of input data; 4= output additional things required by flash (flash.data)]

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

fixg

if TRUE, don't estimate g but use the specified g - useful for computations under the "true" g in simulations

cxx

flag (deprecated, use optmethod) 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.

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

model

c("EE","ET") specifies whether to assume exchangeable effects (EE) or exchangeable T stats (ET).

control

A list of control parameters for the optmization 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 some or all of the following elements (determined by outputlevel)

Details

See readme for more details.

Examples

Run this code

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,betahat,sebetahat,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)

#Running ash with a pre-specified g, rather than estimating it
beta = c(rep(0,100),rnorm(100))
sebetahat = abs(rnorm(200,0,1))
betahat = rnorm(200,beta,sebetahat)
true_g = normalmix(c(0.5,0.5),c(0,0),c(0,1)) # define true g 
#Passing this g into ash causes it to i) take the sd and the means for each component from this g, and ii) initialize pi to the value from this g.
beta.ash = ash(betahat, sebetahat,g=true_g,fixg=TRUE)

Run the code above in your browser using DataLab