rb_dplr {rBahadur}

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)