rb_dplr {rBahadur} | R Documentation |
Binary random variates with Diagonal Plus Low Rank (dplr) correlations
Description
Generate second Bahadur order multivariate Bernoulli random variates with Diagonal Plus Low Rank (dplr) correlation structures.
Usage
rb_dplr(n, mu, U)
Arguments
n |
number of observations |
mu |
vector of means |
U |
outer product component matrix |
Details
This generates multivariate Bernoulli (MVB) random vectors with mean
vector 'mu' and correlation matrix C = D + U U^T
where D
is a diagonal
matrix with values dictated by 'U'. 'mu' must take values in the open unit interval
and 'U' must induce a valid second Bahadur order probability distribution. That is,
there must exist an MVB probability distribution with first moments 'mu' and
standardized central second moments C
such that all higher order central
moments are zero.
Value
An n
-by-m
matrix of binary random variates, where m
is
the length of 'mu'.
Examples
set.seed(1)
h2_0 = .5; m = 200; n = 1000; r =.5; min_MAF=.1
## draw standardized diploid allele substitution effects
beta <- scale(rnorm(m))*sqrt(h2_0 / m)
## draw allele frequencies
AF <- runif(m, min_MAF, 1 - min_MAF)
## compute unstandardized effects
beta_unscaled <- beta/sqrt(2*AF*(1-AF))
## generate corresponding haploid quantities
beta_hap <- rep(beta, each=2)
AF_hap <- rep(AF, each=2)
## compute equilibrium outer product covariance component
U <- am_covariance_structure(beta, AF, r)
## draw multivariate Bernoulli haplotypes
H <- rb_dplr(n, AF_hap, U)
## convert to diploid genotypes
G <- H[,seq(1,ncol(H),2)] + H[,seq(2,ncol(H),2)]
## empirical allele frequencies vs target frequencies
emp_afs <- colMeans(G)/2
plot(AF, emp_afs)
## construct phenotype
heritable_y <- G%*%beta_unscaled
nonheritable_y <- rnorm(n, 0, sqrt(1-h2_0))
y <- heritable_y + nonheritable_y
## empirical h2 vs expected equilibrium h2
(emp_h2 <- var(heritable_y)/var(y))
h2_eq(r, h2_0)
[Package rBahadur version 1.0.0 Index]