createData {ptycho} | R Documentation |
Simulate Data
Description
Create a data object suitable for ptycho.all
.
Usage
createData(X, y, omega = NULL, beta = NULL)
createDataBayesModel(mode = c("exchange","pleiotropy","gene"), n, p, q,
nreps, tau.min, tau.max, G)
createPubData(mode = c("tinysim","ptychoIn",
"exchange","pleiotropy","gene",
"actualGeno","actualPheno","corTest",
"fixedOmega","uniformEffects"),
X=NULL, y=NULL, var.detail=NULL, variants=NULL)
Arguments
X |
Design matrix or alist specifying how to generate such a matrix. If
a list, the first entry is a function name and the second is a list of
arguments to the function. In |
y |
Numeric vector or matrix or list with the following components:
In |
omega |
Numeric vector or matrix or list specifying how to generate a
list with the component |
beta |
List specifying how to generate a matrix of effect sizes. The
first entry of the list is a function name and the second is a list of
arguments to the function, which will be prepended by a matrix specifying
the variables selected and |
n |
Number of observations to simulate |
p |
Number of covariates to simulate |
q |
Number of responses to simulate for each replicate |
nreps |
Number of replicates to simulate |
mode |
String specifying type of dataset to create:
For |
tau.min , tau.max |
Endpoints of uniform distribution from which to draw
|
G |
Number of groups of covariates; unused if |
var.detail |
Data frame with row names same as column names of |
variants |
Character vector containing names of two columns of |
Details
We describe createData
and then describe its wrappers
createDataBayesModel
and createPubData
.
Although createData
can form the data object required by
ptycho.all
when X
and y
are input, it primarily exists to
simplify simulating data from Y=X\beta+\epsilon
, where
\epsilon
is normal with mean zero and specified standard deviation and
\beta
is sparse with entries simulated as specified.
The function generates a specified number of replicates, all of which use the
same design matrix X
. If this matrix is not input, then its argument
must specify a function call to generate it. In either case, suppose X
has n
rows and p
columns.
If the input y
is numeric, then it will be used for the lone replicate.
If it is a matrix, it must have n
rows; let q
be its number of
columns. If input y
is a numeric vector, it must have n
entries
and will be cast as a matrix with q=1
column. Otherwise, input y
is a list specifying, along with the arguments omega
and beta
,
how to simulate the response(s). Because it is useful in analysis of the
estimation of the marginal posterior distribution, the returned object always
contains, regardless of how X
and y
are specified, a matrix
eta2
with (j,k)
entry equal to
\mathbf{x}_j^T \mathbf{y}_k / (n \mathbf{y}_k^T \mathbf{y}_k)
If y
is to be simulated, the first step is to choose the probability
that each covariate is associated with each reponse as specified by the input
argument omega
. If this argument is a matrix, it must have size
p
-by-q
. If it is not a matrix but is numeric, it will be passed
to matrix
to create a matrix of the correct size. Otherwise,
the matrix for each replicate will be generated by calling the function whose
name is given by omega[[1]]
with argument list
(p, q, omega[[2]])
. This function must return a list with component
omega
set to a p
-by-q
matrix; the list may also contain
additional components. The package contains several functions whose names
start with “createOmega” that might guide users in writing their own
functions.
The next step is to draw a p
-by-q
matrix indic.var
whose
(j,k)
entry is equal to one with probability omega[j,k]
and zero
otherwise. This matrix will be drawn until all column sums are positive.
For each entry in indic.var
that is equal to one, the effect size must
be drawn. This is done by calling the function whose name is given by
beta[[1]]
with argument list (indic.var, y$sd, beta[[2]])
. This
function must return a list with component beta
set to a
p
-by-q
matrix; the list may also contain additional components.
If indic.var[j,k]
is zero, then beta[j,k]
should be zero. The
package contains functions whose names start with “createBeta” that
might guide users in writing their own functions.
Finally, an n
-by-q
matrix of noise is drawn from
N(0,\sigma^2)
, where \sigma
is the input noise.sd
, and
added to X\beta
to obtain y
. The column names of each
response matrix generated will be y1
, y2
, and so forth.
The function createPubData
generates the data sets used in Stell and
Sabatti (2015). For mode
equal to “exchange”,
“pleiotropy”, or “geno”, it calls createData
via
createDataBayesModel
; otherwise, it calls createData
directly.
These functions also serve as additional examples of the use of
createData
. For reproducibility, createPubData
first sets the
random seed to 1234, except that it is set to 4 when mode
equals
“ptychoIn” and it does not set it when mode
equals
“corTest”.
In createDataBayesModel
, if mode
is “exchange”, then
one \omega \sim \mbox{Beta}(12,48)
is drawn
independently for each trait. If mode
is “pleiotropy”, then one
probability of association for a trait is drawn from Beta(16,55) for each data
set, that probability is used to draw indic.grp
for each variant, and
then the probability of nonzero indic.var[j,k]
is drawn from
Beta(48,12) for each nonzero indic.grp[j]
. Finally, if mode
is
“gene”, the process is analogous to pleiotropy except that each trait
is simulated independently.
Value
List containing:
X |
Design matrix |
q |
Number of columns in each response |
noise.sd |
Standard deviation of the simulated noise; |
omega |
Input |
beta |
Input |
replicates |
List of length
For |
Author(s)
Laurel Stell and Chiara Sabatti
Maintainer: Laurel Stell <lstell@stanford.edu>
References
Stell, L. and Sabatti, C. (2015) Genetic variant selection: learning across traits and sites, arXiv:1504.00946.
See Also
createOrthogonalX
, createGroupsSim
;
also Data describes tinysim
in example below as well as another
object output by createData
Examples
### EXAMPLE 1
data(tinysim)
# Data generated with mode equal to pleiotropy, so indic.grp exists and
# has an entry for each column in X.
colnames(tinysim$X)
tinysim$replicates[[5]]$indic.grp
# X4, X6, and X9 are associated with some responses.
tinysim$replicates[[5]]$indic.var
### EXAMPLE 2
# Generate miniature data set with information shared across covariates.
set.seed(1234)
tiny1 <- createDataBayesModel(mode="gene", n=100, p=10, q=5, nreps=10,
tau.min=0.045, tau.max=0.063, G=2)
# A covariate can only have indic.var=1 if the group it belongs to has
# indic.grp=1. For example,indic.grp[1,4]=0 implies
# indic.var[groups$group2var[1],4]=0.
tiny1$replicates[[1]]$indic.grp
tiny1$omega[[2]]$groups$group2var[1]
tiny1$replicates[[1]]$indic.var
### EXAMPLE 3
# Alternatively, call createData directly
groups <- createGroupsSim(G=2, p=10)
omegaargs <- list(indic.grp.shape1=16, indic.grp.shape2=55,
shape1=48, shape2=12, groups=groups)
betaargs <- list(tau.min=0.045, tau.max=0.063)
set.seed(1234)
tiny2 <- createData(X=list("createOrthogonalX", list(n=100, p=10)),
y=list(nreps=10, q=5, sd=1),
omega=list("createOmegaCrossVars", omegaargs),
beta=list("createBetaNormal", betaargs))
identical(tiny1, tiny2)
### SEE THE CODE FOR createPubData FOR MORE EXAMPLES.