Linear Mixed Models using Gibbs Sampling
Description
This function runs linear mixed models of the following form:
\mathbf{y} = \mathbf{Xb} + \mathbf{Z}_{1}\mathbf{u}_1 + \mathbf{Z}_{2}\mathbf{u}_2 + \mathbf{Z}_{3}\mathbf{u}_3 + ... + \mathbf{Z}_{k}\mathbf{u}_k + \mathbf{e}
The function allows to include an arbitrary number of independent random effects with each of them being assumed to follow: MVN(\mathbf{0},\mathbf{I}σ^2_{u_k}). If the covariance structure of one random effect is assumed to follow some \mathbf{G} then the design matrix for
that random effect can be constructed as described in Waldmann et al. (2008): \mathbf{F} = \mathbf{ZG}^{1/2}.
Alternatively, the argument ginverse
can be passed.
Usage
1 2 3 
Arguments
y 
vector or list of phenotypes 
X 
fixed effects design matrix of type: 
Z 
list of design matrices for random effects  every element of the list represents one random effect and may be of type: 
ginverse 
list of inverse covariance matrices for random effects in 
par_random 
list of options for random effects. If passed, the list must have as many elements as

niter 
number of iterations 
burnin 
number of iterations to be discarded as burnin 
verbose 
prints progress to the screen 
beta_posterior 
save all posterior samples of regression coefficients 
timings 
prints time per iteration to the screen  sets 
scale_e 
scale parameter for the inverse chisquare prior for the residuals 
df_e 
degrees of freedom for the inverse chisquare prior for the residuals 
seed 
seed for the random number generator. If omitted, a seed will be generated based on machine and time 
use_BLAS 
use BLAS library instead of Eigen 
Details
Single Model run
At this point the function allows to specify the method for any random term as: 'ridge' or 'BayesA'. In Ridge Regression the prior distribution of the random effects is assumed to be normal with a constant variance component, while in BayesA the marginal prior is a tdistribution, resulting from locus specific variances with inverse chisquare priors (Gianola et al., 2009). A wider range of methods is available in the excellent BGLRpackage, which also allows phenotypes to be discrete (de los Campos et al. 2013).
The focus of this function is to allow solving highdimensional problems that are mixtures of sparse and dense features in the design matrices.
The computational expensive parts of the Gibbs Sampler are parallelized as described in Fernando et al. (2014).
Note that the parallel performance highly depends on the number of observations and features present in the design matrices.
It is highly recommended to set the number of threads for less than 10000 observations (length of phenotype vector) to 1 using: set_num_threads(1)
before
running a model. Even for larger sample sizes the parallel performance still depends on the dimension of the feature matrices.
Good results in terms of parallel scaling were observed starting from 50000 observations and more than 10000 features (i.e. number of markers).
Single threaded performance is very good thanks to smart computations during gibbs sampling (Fernando, 2013 (personal communication), de los Campos et al., 2009)
and the use of efficient Eigenmethods for dense and sparse algebra.
Parallel Model runs
In the case of multiple phenotypes passed to the function as a list, the main advantage of the function is that several threads can access the very same data once assigned, which means that the design matrices only have to be allocated once. The parallel scaling of this function using multiple phenotypes is almost linear.
In C++:
For every element of the phenotype list a new instance of an MCMCobject will be created. All the memory allocation needed for running the model is done by the major thread. The function then iterates over all objects and runs the gibbs sampler. This step is parallelized, which means that as many models are being run at the same time as threads available. All MCMCobjects are totally independent from each other, they only share the same designmatrices. Every object has its own randomnumber generator with its own seed which allows perfectly reproducible results.
GWAS using genomic windows
The function allows to specify options to any random effect for conducting genomewide association studies using prediction vector variances of marker windows as described in Fernando et al. (2013). In every effective sample of the Gibbs Sampler the sampling variance of the genotypic value vector \mathbf{g} = \mathbf{Zu} of the particular random effect is computed as: \tilde{σ}_{g}^2 = ≤ft(∑_{j=1}^{n} (g_j  μ_g)^2\right) (n 1)^{1} , with μ_g being the mean of \mathbf{g} and n the number of observations. Then for any window w the sampling variance of \mathbf{g_w} = \mathbf{Z}_w\mathbf{u}_w is obtained as: \tilde{σ}_{g_{w}}^2 = ≤ft(∑_{j=1}^{n} (g_{w_j}  μ_{g_w})^2\right) (n 1)^{1}, where w indicates the range over the columns of \mathbf{Z} that forms the window w. The posterior probability that a window exceeds a specified proportion δ of the total variance is estimated by the number of samples in which \frac{\tilde{σ}_{g_{w}}^2}{\tilde{σ}_{g}^2} > δ divided by the total number of samples. It can be shown that among marker windows that have a posterior probability p or greater for having a variance greater than δ of the total variance, the proportion of false positive signals (PFP) are expected to be lower than 1p (Fernando et al. 2004, Fernando et al., 2013).
Value
List of 4 + number of random effects:
Residual_Variance 
List of 4:

Predicted 
numeric vector of predicted values 
fixed_effects 
List of 4:

Susequently as many additional items as random effects of the following form
Effect_k 
List of 4 + 1 (if

mcmc 
List of 4 + 1 (if

Author(s)
Claas Heuer
Credits: Xiaochen Sun (Iowa State University, Ames) gave strong assistance in the theoretical parts and contributed in the very first implementation of the Gibbs Sampler. Essential parts were adopted from the BayesCimplementation of Rohan Fernando and the BLRpackage of Gustavo de los Campos. The idea of how to parallelize the single site Gibbs Sampler came from Rohan Fernando (2013).
References
Gianola, D., de Los Campos, G., Hill, W.G., Manfredi, E., Fernando, R.: "Additive genetic variability and the bayesian alphabet." Genetics 183(1), 347363 (2009)
de los Campos, G., H. Naya, D. Gianola, J. Crossa, A. Legarra, E. Manfredi, K. Weigel, and J. M. Cotes. "Predicting Quantitative Traits With Regression Models for Dense Molecular Markers and Pedigree." Genetics 182, no. 1 (May 1, 2009): 37585. doi:10.1534/genetics.109.101501.
Waldmann, Patrik, Jon Hallander, Fabian Hoti, and Mikko J. Sillanpaa. "Efficient Markov Chain Monte Carlo Implementation of Bayesian Analysis of Additive and Dominance Genetic Variances in Noninbred Pedigrees." Genetics 179, no. 2 (June 1, 2008): 110112. doi:10.1534/genetics.107.084160.
Meuwissen, T., B. J. Hayes, and M. E. Goddard. "Prediction of Total Genetic Value Using GenomeWide Dense Marker Maps." Genetics 157, no. 4 (2001): 181929.
de los Campos, Gustavo, Paulino Perez Rodriguez, and Maintainer Paulino Perez Rodriguez. "Package 'BGLR,'" 2013. ftp://128.31.0.28/pub/CRAN/web/packages/BGLR/BGLR.pdf.
Fernando, R.L., Dekkers, J.C., Garrick, D.J. "A class of bayesian methods to combine large numbers of genotyped and nongenotyped animals for wholegenome analyses." Genetics Selection Evolution 46(1), 50 (2014)
Fernando, R., Nettleton, D., Southey, B., Dekkers, J., Rothschild, M., Soller, M. "Controlling the proportion of false positives in multiple dependent tests." Genetics 166(1), 611619 (2004)
Fernando, Rohan L., and Dorian Garrick. "Bayesian methods applied to GWAS." GenomeWide Association Studies and Genomic Prediction. Humana Press, 2013. 237274.
See Also
cGBLUP, cSSBR, cGWAS.emmax
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219  ## Not run:
#############################################################
### Running a model with an additive and dominance effect ###
#############################################################
# generate random data
rand_data(500,5000)
### compute the relationship matrices
G.A < cgrm.A(M,lambda=0.01)
G.D < cgrm.D(M,lambda=0.01)
### generate the list of design matrices for clmm
Z_list = list(t(chol(G.A)),t(chol(G.D)))
### specify options
par_random = list(list(method="ridge",scale=var(y)/2 ,df=5, name="add"),
list(method="ridge",scale=var(y)/10,df=5, name="dom"))
### run
set_num_threads(1)
fit < clmm(y = y, Z = Z_list, par_random=par_random, niter=5000, burnin=2500)
### inspect results
str(fit)
########################
### Cross Validation ###
########################
### 4fold crossvalidation with one repetition:
# generate random data
rand_data(500,5000)
### compute the list of masked phenotypevectors for CV
y_CV < cCV(y, fold=4, reps=1)
### Cross Validation using GBLUP
G.A < cgrm.A(M,lambda=0.01)
### generate the list of design matrices for clmm
Z_list = list(t(chol(G.A)))
### specify options
h2 = 0.3
scale = unlist(lapply(y_CV,function(x)var(x,na.rm=T))) * h2
df = rep(5,length(y_CV))
par_random = list(list(method="ridge",scale=scale,df=df, name="animal"))
### run model with 4 threads
set_num_threads(4)
fit < clmm(y = y_CV, Z = Z_list, par_random=par_random, niter=5000, burnin=2500)
### inspect results
str(fit)
### obtain predictions
pred < get_pred(fit)
### prediction accuracy
get_cor(pred,y_CV,y)
########################################################
### GWAS using Bayesian Regression on marker windows ###
########################################################
# generate random data
rand_data(500,5000)
### generate the list of design matrices for clmm
Z_list = list(M)
### specify options
h2 = 0.3
scale = var(y) * h2
df = 5
# specifying the model
# Here we use ridge regression on the marker covariates
# and define a window size of 100 and a threshold of 0.01
# which defines the proportion of genetic variance accounted
# for by a single window
par_random = list(list(method="ridge",scale=scale,df=df,
GWAS=list(window_size=100, threshold=0.01), name="markers"))
### run
set_num_threads(1)
fit < clmm(y =y, Z = Z_list, par_random=par_random, niter=2000, burnin=1000)
### inspect results
str(fit)
### extract GWAS part
gwas < fit$markers$GWAS
# plot window variance proportions
plot(gwas$window_variance_proportion)
##########################################################
### Sparse Animal Model using the pedigreemm milk data ###
##########################################################
# cpgen offers two ways of running models with random effects that
# are assumed to follow some covariance structure.
# 1) Construct the Covariance matrix and pass the cholesky of that
# as design matrix for that random effect
# 2) Construct the inverse of the covariance matrix (ginverse) and pass the design
# matrix 'Z' that relates observations to factors in ginverse in conjunction
# with the symmetric ginverse.
# This is approach 2) which is more convenient for pedigree based
# animal models, as ginverse (Inverse of numerator relationship matrix) is
# very sparse and can be obtained very efficiently
set_num_threads(1)
# load the data
data(milk)
# get Ainverse
# Ainv < as(getAInv(pedCows),"dgCMatrix")
T_Inv < as(pedCows, "sparseMatrix")
D_Inv < Diagonal(x=1/Dmat(pedCows))
Ainv<t(T_Inv)
dimnames(Ainv)[[1]]<dimnames(Ainv)[[2]] <pedCows@label
Ainv < as(Ainv, "dgCMatrix")
# We need to construct the design matrix.
# Therefore we create a second id column with factor levels
# equal to the animals in the pedigree
milk$id2 < factor(as.character(milk$id), levels = pedCows@label)
# set up the design matrix
Z < sparse.model.matrix(~ 1 + id2, data = milk, drop.unused.levels=FALSE)
# run the model
niter = 5000
burnin = 2500
modAinv < clmm(y = as.numeric(milk$milk), Z = list(Z), ginverse = list(Ainv),
niter = niter, burnin = burnin)
# This is approach 1) run an equivalent model using the cholesky of A
# get L from A = LL'
L < as(t(relfactor(pedCows)),"dgCMatrix")
# match with ids
ZL < L[match(milk$id, pedCows@label),]
# run the model
modL < clmm(as.numeric(milk$milk), Z= list(ZL),
niter = niter, burnin = burnin)
### a more advanced model
# y = Xb + Zu + a + e
#
# u = permanent environment of animal
# a = additive genetic effect of animal
Zpe < sparse.model.matrix(~ 1 + id2, drop.unused.levels = TRUE, data = milk)
# make X and account for lactation and herd
X < sparse.model.matrix(~ 1 + lact + herd, data = milk)
niter = 10000
burnin = 2500
mod2 < clmm(as.numeric(milk$milk), X = X, Z = list(Zpe,Z), ginverse = list(NULL, Ainv),
niter = niter, burnin = burnin)
# run all phenotypes in the milk dataset at once in parallel
Y < list(as.numeric(milk$milk),as.numeric(milk$fat),as.numeric(milk$prot),as.numeric(milk$scs))
set_num_threads(4)
# ginverse version
model < clmm(Y, X = X, Z = list(Zpe,Z), ginverse = list(NULL, Ainv),
niter = niter, burnin = burnin)
# get heritabilities and repeatabilities with their standard deviations
heritabilities < array(0, dim=c(length(Y),2))
colnames(heritabilities) < c("h2","sd")
# only use postburnin samples
range < (burnin+1):niter
# h2
heritabilities[,1] < unlist(lapply(model, function(x)
mean(
x$Effect_2$posterior$variance[range] /
(x$Effect_1$posterior$variance[range] +
x$Effect_2$posterior$variance[range] +
x$Residual_Variance$Posterior[range]))
)
)
# standard deviation of h2
heritabilities[,2] < unlist(lapply(model, function(x)
sd(
x$Effect_2$posterior$variance[range] /
(x$Effect_1$posterior$variance[range] +
x$Effect_2$posterior$variance[range] +
x$Residual_Variance$Posterior[range]))
)
)
## End(Not run)
