il: il

In jamesliley/cfdr: Implements cFDR testing with overall FDR control

Description

Integrate over L (general). Assumes Zq| H^p=0 has a mixture distribution, being N(0,1) with probability pi0, and taking some other distribution distx with probability (1-pi0)

Usage

il(X, Y = NULL, pi0_null = NULL, sigma_null = rep(1,
  length(pi0_null)), rho_null = 0, distx = "norm", ...)

Arguments

X	either output from functions vl, vlx, or vly, or matrix nk x nc of x-co-ordinates (Z-plane) of regions to integrate over. X[k,] is the set of co-ordinates for the kth region.
Y	vector of length nc of y-co-ordinates to integrate over. Assumed to be constant for all columns of X. Leave as NULL if X is output from vl, vlx, or vly.
pi0_null	parameter for distribution of Q|H^p=0. Can be a vector of parameters of length np.
sigma_null	scale parameter for distribution of Q|H^p=0. Can be a vector of parameters
rho_null	optional parameter governing covariance between Z scores under the null; for instance, from shared controls
distx	distribution type for distribution of Q|H^p=0. Can be 1,2,or 3 (symbolising normal, t (3df) and Cauchy respectively) or a text string which can be appended to 'd' to get PDF and 'p' to get CDF
...	additional parameters passed to CDF and PDF functions, ie df=3

Details

We generally assume that distx is a Gaussian centred at 0.

Value

matrix of dimension nk x np; [k,p]th element is the integral for the kth region using the pth parameter values.

Author(s)

James Liley

Examples

# Generate standardised simulated dataset
set.seed(1); n=10000; n1p=100; n1pq=100; n1q=100
zp=c(rnorm(n1p,sd=2), rnorm(n1q,sd=1),rnorm(n1pq,sd=2), rnorm(n-n1p-n1q-n1pq,sd=1))
zq=c(rnorm(n1p,sd=1), rnorm(n1q,sd=2),rnorm(n1pq,sd=2), rnorm(n-n1p-n1q-n1pq,sd=1))
p=2*pnorm(-abs(zp)); q=2*pnorm(-abs(zq))

# Generate some L regions
example_indices=c(70,67,226)
v1=vl(p,q,indices=example_indices,mode=0); 


# The true distribution of Zq|H^p=0 is N(0,1) with weight n1q/(n-n1p-n1pq), and N(0,2^2) with weight 1- (n1q/(n-n1p-n1pq))
true_q0_pars=c(n1q/(n-n1p-n1pq),2)

# Estimate parameters:
est_q0_pars=fit.2g(q[which(p>0.5)])$pars

############### Integrals ################
int_true=il(v1,pi0_null=true_q0_pars[1],sigma_null=true_q0_pars[2],distx="norm")
int_est=il(v1,pi0_null=est_q0_pars[1],sigma_null=est_q0_pars[2],distx="norm")
##########################################

############# Check integral #############
# Sample values zp0,zq0 and p0,q0 following distribution of Zp,Zq|H^p=0
nsc=1000000 # generate nsc sample values
n0=round(nsc*true_q0_pars[1]); n1=nsc-n0
zp0=rnorm(nsc,sd=1)
zq0=c(rnorm(n0,sd=1), rnorm(n1,sd=true_q0_pars[2]))
p0=2*pnorm(-abs(zp0)); q0=2*pnorm(-abs(zq0))


# Proportion of values p0,q0 falling in region with co-ordinates v1$x[2,],v1$y
length(which(in.out(cbind(v1$x[2,],v1$y),cbind(p0,q0))))/nsc

# Value of integral over the region
int_true[2] 
##########################################

jamesliley/cfdr documentation built on July 31, 2020, 9:42 a.m.