simulBslmm: BSLMM

In timflutre/rutilstimflutre: Timothee Flutre's personal R code

BSLMM

Description

Simulate phenotypes according to the Bayesian sparse linear mixed model (Zhou, Carbonetto & Stephens, 2013)): y = W alpha + Z X_c beta-tilde + Z u + epsilon where y is N x 1; W is N x Q; alpha is Q x 1; Z is N x I; X_c is I x P; u is I x 1 and epsilon is N x 1. For all p, beta-tilde_p ~ pi Norm_1(0, sigma_betat^2/sigma^2) + (1 - pi) delta_0; u ~ Norm_I(0, (sigma_u^2/sigma^2) K) and epsilon ~ Norm_N(0, sigma^2 I).

Usage

simulBslmm(
  Q = 3,
  mu = 50,
  mean.a = 5,
  sd.a = 2,
  X,
  pi = NULL,
  h = NULL,
  rho = NULL,
  tau = 1,
  enforce.zhou = TRUE,
  perc.NA = 0,
  err.df = Inf,
  seed = NULL
)

Arguments

Q	number of fixed effects, i.e. the intercept plus the number of years during which genotypes are phenotyped (starting in 2010)
mu	overall mean
mean.a	mean of the prior on alpha[2:Q]
sd.a	std dev of the prior on alpha[2:Q]
X	matrix of bi-allelic SNP genotypes encoded in allele doses in [0,2], with genotypes in rows and SNPs in columns (SNPs with missing values or low MAF should be discarded beforehand).
pi	proportion of beta-tilde values that are non-zero; if NULL, log(pi) is drawn from Unif(log(1/P),log(1))
h	approximation to E[PVE]; if NULL, drawn from Unif(0,1)
rho	approximation to E[PGE]; if NULL and pi=0, then rho=0, else rho is drawn from Unif(0,1)
tau	precision of the errors (1 / sigma^2)
enforce.zhou	if TRUE, X will be transformed into X_c by centering columns, and K will be calculated as X_c X_c^T / P; otherwise, X will be transformed into (-1,0,1) and K will be calculated as estimGenRel with method="vanraden1"
perc.NA	percentage of missing phenotypes, at random
err.df	degrees of freedom of errors' Student's t-distribution
seed	seed for the pseudo-random number generator

Value

list

Author(s)

Timothee Flutre

Examples

## Not run: ## simulate genotypes
set.seed(1859)
X <- simulGenosDose(nb.genos=200, nb.snps=2000)
afs <- estimSnpAf(X)

## particular case: LMM (only u contributes)
mod.lmm <- simulBslmm(X=X, pi=0, h=0.7, seed=3591)

## particular case: BVSR (only beta-tilde contributes)
mod.bvsr <- simulBslmm(X=X, pi=0.1, h=0.7, rho=1, seed=3591)

## general case: BSLMM (both u and beta-tilde contribute)
mod.bslmm <- simulBslmm(X=X, pi=0.1, h=0.7, rho=0.7, seed=3591)

library(rrBLUP)
mod.lmmC <- simulBslmm(X=X, pi=0, h=0.7, enforce.zhou=FALSE, seed=3591)
fit.lmmC1 <- mixed.solve(y=mod.lmmC$y, Z=mod.lmmC$Z,
                         K=mod.lmmC$X.c %*% t(mod.lmmC$X.c),
                         X=mod.lmmC$W)
cbind(truth=c(mod.lmmC$alpha), estim=fit.lmmC1$beta)
cbind(truth=c(mod.lmmC$sigma2, mod.lmmC$sigma.u2 / sum(2 * afs * (1-afs))),
      estim=c(fit.lmmC1$Ve, fit.lmmC1$Vu))
fit.lmmC2 <- mixed.solve(y=mod.lmmC$y, Z=mod.lmmC$Z %*% mod.lmmC$X.c,
                         X=mod.lmmC$W)
cbind(truth=c(mod.lmmC$alpha), estim=fit.lmmC2$beta)
cbind(truth=c(mod.lmmC$sigma2, mod.lmmC$sigma.u2 / sum(2 * afs * (1-afs))),
      estim=c(fit.lmmC2$Ve, fit.lmmC2$Vu))

## End(Not run)

timflutre/rutilstimflutre documentation built on Feb. 7, 2024, 8:17 a.m.