qpAvgNrr: Average non-rejection rate estimation

In rcastelo/qpgraph: Estimation of genetic and molecular regulatory networks from high-throughput genomics data

Average non-rejection rate estimation

Description

Estimates average non-rejection rates for every pair of variables.

Usage

## S4 method for signature 'ExpressionSet'
qpAvgNrr(X, qOrders=4, I=NULL, restrict.Q=NULL,
                                   fix.Q=NULL, nTests=100, alpha=0.05,
                                   pairup.i=NULL, pairup.j=NULL, type=c("arith.mean"),
                                   verbose=TRUE, identicalQs=TRUE,
                                   exact.test=TRUE, use=c("complete.obs", "em"),
                                   tol=0.01, R.code.only=FALSE, clusterSize=1,
                                   estimateTime=FALSE, nAdj2estimateTime=10)
## S4 method for signature 'data.frame'
qpAvgNrr(X, qOrders=4, I=NULL, restrict.Q=NULL,
                                fix.Q=NULL, nTests=100, alpha=0.05, pairup.i=NULL,
                                pairup.j=NULL, long.dim.are.variables=TRUE,
                                type=c("arith.mean"), verbose=TRUE,
                                identicalQs=TRUE, exact.test=TRUE,
                                use=c("complete.obs", "em"), tol=0.01, R.code.only=FALSE,
                                clusterSize=1, estimateTime=FALSE, nAdj2estimateTime=10)
## S4 method for signature 'matrix'
qpAvgNrr(X, qOrders=4, I=NULL, restrict.Q=NULL, fix.Q=NULL,
                            nTests=100, alpha=0.05, pairup.i=NULL,
                            pairup.j=NULL, long.dim.are.variables=TRUE,
                            type=c("arith.mean"), verbose=TRUE,
                            identicalQs=TRUE, exact.test=TRUE,
                            use=c("complete.obs", "em"), tol=0.01, R.code.only=FALSE,
                            clusterSize=1, estimateTime=FALSE, nAdj2estimateTime=10)

Arguments

X	data set from where to estimate the average non-rejection rates. It can be an ExpressionSet object, a data frame or a matrix.
qOrders	either a number of partial-correlation orders or a vector of vector of particular orders to be employed in the calculation.
I	indexes or names of the variables in X that are discrete. When X is an ExpressionSet then I may contain only names of the phenotypic variables in X. See details below regarding this argument.
restrict.Q	indexes or names of the variables in X that restrict the sample space of conditioning subsets Q.
fix.Q	indexes or names of the variables in X that should be fixed within every conditioning conditioning subsets Q.
nTests	number of tests to perform for each pair for variables.
alpha	significance level of each test.
pairup.i	subset of vertices to pair up with subset pairup.j
pairup.j	subset of vertices to pair up with subset pairup.i
long.dim.are.variables	logical; if TRUE it is assumed that when the data is a data frame or a matrix, the longer dimension is the one defining the random variables; if FALSE, then random variables are assumed to be at the columns of the data frame or matrix.
type	type of average. By now only the arithmetic mean is available.
verbose	show progress on the calculations.
identicalQs	use identical conditioning subsets for every pair of vertices (default), otherwise sample a new collection of nTests subsets for each pair of vertices.
exact.test	logical; if FALSE an asymptotic conditional independence test is employed with mixed (i.e., continuous and discrete) data; if TRUE (default) then an exact conditional independence test with mixed data is employed.
use	a character string defining the way in which calculations are done in the presence of missing values. It can be either "complete.obs" (default) or "em".
tol	maximum tolerance controlling the convergence of the EM algorithm employed when the argument use="em".
R.code.only	logical; if FALSE then the faster C implementation is used (default); if TRUE then only R code is executed.
clusterSize	size of the cluster of processors to employ if we wish to speed-up the calculations by performing them in parallel. A value of 1 (default) implies a single-processor execution. The use of a cluster of processors requires having previously loaded the packages snow and rlecuyer.
estimateTime	logical; if TRUE then the time for carrying out the calculations with the given parameters is estimated by calculating for a limited number of adjacencies, specified by nAdj2estimateTime, and extrapolating the elapsed time; if FALSE (default) calculations are performed normally till they finish.
nAdj2estimateTime	number of adjacencies to employ when estimating the time of calculations (estimateTime=TRUE). By default this has a default value of 10 adjacencies and larger values should provide more accurate estimates. This might be relevant when using a cluster facility.

Details

Note that when specifying a vector of particular orders q, these values should be in the range 1 to min(p, n-3), where p is the number of variables and n the number of observations. The computational cost increases linearly within each q value and quadratically in p. When setting identicalQs to FALSE the computational cost may increase between 2 times and one order of magnitude (depending on p and q) while asymptotically the estimation of the non-rejection rate converges to the same value.

When I is set different to NULL then mixed graphical model theory is employed and, concretely, it is assumed that the data comes from an homogeneous conditional Gaussian distribution. In this setting further restrictions to the maximum value of q apply, concretely, it cannot be smaller than p plus the number of levels of the discrete variables involved in the marginal distributions employed by the algorithm. By default, with exact.test=TRUE, an exact test for conditional independence is employed, otherwise an asymptotic one will be used. Full details on these features can be found in Tur, Roverato and Castelo (2014).

Value

A dspMatrix-class symmetric matrix of estimated average non-rejection rates with the diagonal set to NA values. When using the arguments pairup.i and pairup.j, those cells outside the constraint pairs will get also a NA value.

Note, however, that when estimateTime=TRUE, then instead of the matrix of estimated average non-rejection rates, a vector specifying the estimated number of days, hours, minutes and seconds for completion of the calculations is returned.

Author(s)

R. Castelo and A. Roverato

References

Castelo, R. and Roverato, A. Reverse engineering molecular regulatory networks from microarray data with qp-graphs. J. Comp. Biol., 16(2):213-227, 2009.

Tur, I., Roverato, A. and Castelo, R. Mapping eQTL networks with mixed graphical Markov models. Genetics, 198:1377-1393, 2014.

See Also

qpNrr qpEdgeNrr qpHist qpGraphDensity qpClique

Examples

require(mvtnorm)

nVar <- 50  ## number of variables
maxCon <- 3 ## maximum connectivity per variable
nObs <- 30  ## number of observations to simulate

set.seed(123)

A <- qpRndGraph(p=nVar, d=maxCon)
Sigma <- qpG2Sigma(A, rho=0.5)
X <- rmvnorm(nObs, sigma=as.matrix(Sigma))

avgnrr.estimates <- qpAvgNrr(X, verbose=FALSE)

## distribution of average non-rejection rates for the present edges
summary(avgnrr.estimates[upper.tri(avgnrr.estimates) & A])

## distribution of average non-rejection rates for the missing edges
summary(avgnrr.estimates[upper.tri(avgnrr.estimates) & !A])

## Not run: 
library(snow)
library(rlecuyer)

## only for moderate and large numbers of variables the
## use of a cluster of processors speeds up the calculations

nVar <- 500
maxCon <- 3
A <- qpRndGraph(p=nVar, d=maxCon)
Sigma <- qpG2Sigma(A, rho=0.5)
X <- rmvnorm(nObs, sigma=as.matrix(Sigma))

system.time(avgnrr.estimates <- qpAvgNrr(X, q=10, verbose=TRUE))
system.time(avgnrr.estimates <- qpAvgNrr(X, q=10, verbose=TRUE, clusterSize=4))

## End(Not run)

rcastelo/qpgraph documentation built on Oct. 28, 2024, 5:15 a.m.