rBahadur (version 1.0.0)

rb_unstr: Binary random variates with unstructured correlations

Description

Generate Bahadur order-2 multivariate Bernoulli random variates with unstructured correlations.

Usage

rb_unstr(n, mu, C)

Value

An \(n\)-by-\(m\) matrix of binary random variates, where \(m\) is the length of 'mu'.

Arguments

n

number of observations

mu

vector of means

C

correlation matrix

Details

This generates multivariate Bernoulli (MVB) random vectors with mean vector 'mu' and correlation matrix 'C'. 'mu' must take values in the open unit interval and 'C' 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.

Examples

set.seed(1)
h2_0 = .5; m = 200; n = 500; 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)

## construct Correlation matrix
S <- diag(1/sqrt(AF_hap*(1-AF_hap)))
DPLR <- U%o%U
diag(DPLR) <- 1
C <- cov2cor(S%*%DPLR%*%S)

## draw multivariate Bernoulli haplotypes
H <- rb_unstr(n, AF_hap, C)

## 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)