关于BayesCpi
本文在进行基因型系数矩阵估计时采用的是 BayesCpi 方法
理论部分
第一部分：关于固定效应和随机效应系数矩阵的采样分布定义
根据混合线性模型的系数方程求解的方式，有：
公式1
其中∑代表协方差矩阵，X代表固定效应的设计矩阵，R代表随机效应的设计矩阵，β和r分别代表固定效应和随机效应的系数矩阵，A-1λ 代表遗传力矩阵（在该例子中 A-1λ = I(σe2/σr2)） ，如果协方差矩阵为单位矩阵，则可以得到如下简化：
公式2
其中∑代表协方差矩阵，X代表固定效应的设计矩阵，R代表随机效应的设计矩阵，β和r分别代表固定效应和随机效应的系数矩阵，A-1λ 代表遗传力矩阵（在该例子中 A-1λ = I(σe2/σr2)），此时根据 blue 和 blup 框架下可求解系数：
根据MCMC采样的原理，我们可以实现设置一个初始值：β0，r0，α0 和 μ0，并且定义 y* 如下：
公式3
，因此公式2可以改写为：
公式4
结合MCMC采样迭代的理论，第一轮采样就可以基于如下分布（ε~N(0,σe2)），上角标 0 和 1 代表mcmc迭代的次数，
公式5
更一般的：
公式6
X代表固定效应的设计矩阵，R代表随机效应的设计矩阵，β和r分别代表固定效应和随机效应的系数矩阵，上角标 i 代表mcmc迭代的次数，A-1λ 代表遗传力矩阵（在该例子中 A-1λ = I(σe2/σr2)）
如果把固定效应或者随机效应设计矩阵的每一个变量看作是一个随机变量
如果把固定效应或者随机效应设计矩阵的每一个变量（每一列）看作是一个随机变量，那么更一般地有：
公式7
X代表固定效应的设计矩阵，R代表随机效应的设计矩阵，β和r分别代表固定效应和随机效应的系数矩阵，上角标 i 代表mcmc迭代的次数，j
代表设计矩阵的第 j 列，A-1λ 代表遗传力矩阵（在该例子中 A-1λ = I(σe2/σr2)）
因此，每次采样完做最优化的目标函数为：
公式8
即如果系数 β 和 r 在mcmc采样的过程中达到稳定，那么 βji-1 - βji 和 rji-1 - rji 的差值将会为 0（或者很小）
第二部分：关于 genetic marker effect 系数矩阵的采样分布定义
笔者认为genetic marker effect相当于一项随机效应项，因此按照随机效应系数矩阵计算公式BLUP进行求解，根据第一部分公式1-6的计算，根据随机效应系数矩阵的计算方法，同理可推得genetic marker effect的系数矩阵 α 的采样分布为：
公式9
M代表genetic marker effect的设计矩阵，α 代表 genetic marker effect 的系数矩阵，K-1λ 代表遗传力矩阵（在该例子中 K-1λ = I(σe2/σα2)）有：
公式11
值得注意的是，genetic marker effect的系数矩阵的估计值采用的是拒绝-接收的采样方式，类似于M-H采样，因此需要构造接收函数或接收值，来方便判断该轮采样的值是否被接收
构造接收函数的方法如下：
极大似然的目标参数为Πk（Πk：the proportion of markers in zero effect size），依据极大似然估计的原理，我们构造参数的似然函数 LΠk 如下：（参考文章：https://www.nature.com/articles/s41467-019-12653-0#Sec18，补充材料的45页）
公式11
构造该似然函数的目的是为了构造拒绝-接收采样中的接收函数或接受率：
公式12
这个概率值就为接收函数或者接收值；一般情况下，nΠ = 2，第一个元素表征具有零效应的SNP比例，第二个元素表征具有非零效应的SNP比例。每次采样的时候从均匀分布中随机取一个数字 num 与上图的 acceptProb 作比较，若 num < acceptProb 接收采样 ，否则拒绝采样
在mcmc过程中，Πk（k > 0）服从 dirichlet 分布，η 代表当前轮次下拒绝-接收采样中，接收采样的次数（当前轮次，比方说 iter = 3，接收采样的次数为 2，则 η = 2；iter = 10，接收采样的次数为 6，则 η = 6）
公式13
代码部分
首先介绍下所用到的数据：
1.pheno： 一些表型信息和metadata的内容，
pheno 包括测序样品的id，性别，采样季节等基本信息，T1可能是表型信息（连续型）
2.geno：数字基因型矩阵（n × m, n为个体数，m为snp数，本例中为 600 × 1000），接收 0, 1, 2 或者 -1, 0, 1 这两种形式的矩阵，但矩阵中不能含有NAgeno矩阵
3.geno.id：代表的是geno矩阵的样品id名称，该例子中geno.id的长度为600geno.id
4.map：类似于VCF文件，列名依次表示为SNP名称，染色体编号，SNP位置和基因型map 该例子中基因型的编码规则为A1A1编码为2，A1A2编码为1，A2A2编码为0，其中A1代表SNP上第一个等位基因，即参考等位基因（通常出现频率较多）；A2代表SNP上第二个等位基因，即替代等位基因（通常出现频率较少）
第二部分，输入的介绍：
# 建立的混合线性模型
formula = T1 ~ season + bwt + (1 | loc) + (1 | dam)
# 一些表型信息和metadata的内容
data = pheno
# 数字基因型矩阵
M = geno
# geno矩阵的样品id名称
M.id = geno.id
# 方法选择
method = "BayesCpi"
# 类似于VCF文件，见 4
map = map
# 所设定的零效应SNP和有效应的SNP的比例（零效应的SNP比例为0.98，有效应的SNP的比例为0.02）
Pi = c(0.98, 0.02)
# MCMC 估计参数时的最大迭代次数
niter = 20000
# nburn 与 thin 和估计参数有关
nburn = 16000
thin = 5
lambda = 0.0
printfreq = 100
seed = 666666
threads = 4
verbose = TRUE
fold = NULL
windindx <- NULL
windsize = NULL
windnum = NULL
dfvr = NULL
s2vr = NULL
vg = NULL
dfvg = NULL
s2vg = NULL
ve = NULL
dfve = NULL
s2ve = NULL
正式的算法代码：
set.seed(seed)
# 整理格式
M.id <- as.character(as.matrix(M.id)[, 1, drop=TRUE])
data <- data[match(M.id, as.character(data[, 1, drop=TRUE])), ]
myformula <- paste(formula[2], formula[3], sep='~')
library(stringr)
# 将随机效应项 loc 和 dam 提取出来
rand_term <- unlist(str_extract_all(myformula, "(?<=(\\+ |~)\\(1 \\| )[:\\w\\d]+(?=\\))"))
R <- NULL
# 将随机效应项对应的内容提取出来
for(r in rand_term){
split_str = unlist(strsplit(r,":"))
if(length(split_str) != 1){
Ri <- as.matrix(apply(data[, split_str], 1, function(x){if(sum(is.na(x)) != 0){NA}else{paste(x, collapse = ":")}}))
}else{
Ri <- as.matrix(as.character(data[, split_str]))
}
R <- cbind(R, Ri)
}
# 将固定效应项对应的内容提取出来
fixed_formula <- str_replace_all(myformula, pattern = "(\\+ |(?<=~))\\(1 \\| [:\\w\\d]+\\)", replacement = "" )
fixed_formula <- str_replace_all(fixed_formula, pattern = "~ *\\+ ", replacement = "~" )
fixed_formula <- str_replace_all(fixed_formula, pattern = "~ *\\- ", replacement = "~" )
fixed_formula <- str_replace_all(fixed_formula, pattern = "~ *$", replacement = "~1" )
# warning
warn_pattern = "(. |~)\\(.*? \\| .*?\\)"
fixed_formula <- formula(fixed_formula)
其中：
myformula在该例子中代表T1 ~ season + bwt + (1 | loc) + (1 | dam)，这个式子表示 表型T1作为响应变量，固定效应项为season和bwt，随机效应则不考虑随机效应项的斜率，只考虑随机效应项的截距，分组的随机因子分别是loc和dam
2.上述循环中，将随机效应项对应的内容提取出来在该例子中指的是!pheno 将红框中的内容提取出来
3.最终的 fixed_formula 提取出来为字符串T1 ~ season + bwt
# 对应多张表的id
yNA <- union(attr(model.frame(fixed_formula, data = data), "na.action"), attr(na.omit(R), "na.action"))
yNA <- union(yNA, which(is.na(data[, as.character(formula[2]), drop = TRUE])))
yNA <- c(1:length(M.id)) %in% yNA
if(all(yNA)) stop("no effective data left.")
# 将固定效应按条件转换为设计矩阵
X <- model.matrix(fixed_formula, data = data[!yNA, , drop = FALSE])
X <- X[, !apply(X, 2, function(x){all(x == 1)}), drop = FALSE]
if(!ncol(X)) X <- NULL
# 获得随机效应项的内容
R <- R[!yNA, , drop = FALSE]
# 提取对应的响应变量值
y <- data[!yNA, as.character(formula[2]), drop = TRUE]
if(!is.matrix(M)){M <- as.matrix(M); gc()}
# 将id对不上的 geno：数字基因型矩阵赋给 Mp
Mp <- M[yNA, , drop=FALSE]
# 将id对上的 geno：数字基因型矩阵赋给 M
M <- M[!yNA, , drop=FALSE]
library(Rcpp)
Rcpp::sourceCpp('C:/Users/lenovo/Downloads/hibayes-master/src/Bayes.cpp')
windindx <- NULL
res <- Bayes(y=y, X=M, model=method, Pi=Pi,
fold=fold, C=X, R=R, niter=niter,
nburn=nburn, thin=thin, windindx=windindx,
dfvr=dfvr, s2vr=s2vr, vg=vg, dfvg=dfvg,
s2vg=s2vg, ve=ve, dfve=dfve, s2ve=s2ve,
outfreq=printfreq, threads=threads,
verbose=verbose)
其中：
- 固定效应的设计矩阵 X 为
X
其中season为因子型变量，且将Autumn作为参考物，Spring，Summer和Winter都与它相比；而bwt为连续型变量，直接用一列向量表示
- y 为表型值，在该例子中为一列长度是300的向量
- 此时经过id对应后的M矩阵（经过id对应筛选后的geno：数字基因型矩阵；300 × 1000）和表型值y（长度为300），就可以做后续贝叶斯分析了
接下来重点解析函数 Bayes()，这一部分是用cpp写的
// 首先是定义变量
# y 值得是表型数据，数据类型为一列向量
arma::vec &y,
# X 是 geno：数字基因型矩阵，数据类型为矩阵
arma::mat &X,
# 定义模型类型，数据类型为字符串，该例子为 BayesCpi
std::string model,
# 所设定的零效应SNP和有效应的SNP的比例（零效应的SNP比例为0.98，有效应的SNP的比例为0.02），数据类型为一列向量
arma::vec Pi,
# 以下这些是计算时所需要的数值型变量
const Nullable<arma::vec> Kival = R_NilValue,
const Nullable<arma::mat> Ki = R_NilValue,
const Nullable<arma::mat> C = R_NilValue,
const Nullable<CharacterMatrix> R = R_NilValue,
# fold proportion of variance explained for groups of SNPs，即SNPs的方差解释性，本例未涉及
const Nullable<arma::vec> fold = R_NilValue,
const int niter = 50000,
const int nburn = 20000,
const int thin = 5,
const Nullable<arma::vec> epsl_y_J = R_NilValue,
const Nullable<arma::sp_mat> epsl_Gi = R_NilValue,
const Nullable<arma::uvec> epsl_index = R_NilValue,
const Nullable<double> dfvr = R_NilValue,
const Nullable<double> s2vr = R_NilValue,
const Nullable<double> vg = R_NilValue,
const Nullable<double> dfvg = R_NilValue,
const Nullable<double> s2vg = R_NilValue,
const Nullable<double> ve = R_NilValue,
const Nullable<double> dfve = R_NilValue,
const Nullable<double> s2ve = R_NilValue,
const Nullable<arma::uvec> windindx = R_NilValue,
const int outfreq = 100,
const int threads = 0,
const bool verbose = true
接下来，cpp代码中对一些基本参数进行了判断
omp_setup(threads);
# 判断表型数据 y 中是否含有 NA值
if(y.has_nan()) throw Rcpp::exception("NAs are not allowed in y.");
# 判断表型数据 y 的长度是否与geno：数字基因型矩阵的行数相同
if(y.n_elem != X.n_rows) throw Rcpp::exception("Number of individuals not equals.");
# 判断模型类型，并赋予数字index
int model_index = (model == "BayesRR" ? 1 : (model == "BayesA" ? 2 : (model == "BayesB" || model == "BayesBpi" ? 3 : (model == "BayesC" || model == "BayesCpi" || model == "BSLMM" ? 4 : (model == "BayesL" ? 5 : 6)))));
bool fixpi = false;
# 若模型类型为 "BayesB" 或 "BayesC" 则将 fixpi = true
if(model == "BayesB" || model == "BayesC") fixpi = true;
# 判断所设定的零效应SNP和有效应的SNP的比例中的元素是否小于2个
if(Pi.n_elem < 2) throw Rcpp::exception("Pi should be a vector.");
# 判断所设定的零效应SNP和有效应的SNP的比例中的元素之和是否等于1
if(sum(Pi) != 1) throw Rcpp::exception("sum of Pi should be 1.");
if(Pi[0] == 1) throw Rcpp::exception("all markers have no effect size.");
# 判断所设定的零效应SNP和有效应的SNP的比例，其数值是否介于0-1之间
for(int i = 0; i < Pi.n_elem; i++){
if(Pi[i] < 0 || Pi[i] > 1){
throw Rcpp::exception("elements of Pi should be at the range of [0, 1]");
}
}
# fold proportion of variance explained for groups of SNPs
# 即SNPs的方差解释性，本例未涉及，因此更改向量 fold 的容器大小，为两个长度的向量
vec fold_;
if(fold.isNotNull()){
fold_ = as<arma::vec>(fold);
}else{
if(model == "BayesR") throw Rcpp::exception("'fold' should be provided for BayesR model.");
# 更改向量 fold 的容器大小，为两个长度的向量
fold_.resize(2);
}
# 判断SNPs方差解释的向量是否与零效应SNP和有效应的SNP比例向量等长
if(fold_.n_elem != Pi.n_elem){
throw Rcpp::exception("length of Pi and fold not equals.");
}
# 获得表型向量 y 的元素个数
int n = y.n_elem;
# 获得 geno：数字基因型矩阵的列数
int m = X.n_cols;
# 计算表型向量 y 的方差
double vary = var(y);
double h2 = 0.5;
int n_records = (niter - nburn) / thin;
# 定义参数
double* dci;
int nc = 0;
# C_ 为 X 矩阵，即固定效应的设计矩阵
mat C_;
vec beta;
vec betasd;
mat beta_store;
# 此时的C_即为X矩阵，即固定效应的设计矩阵，以下代码的目的是对矩阵进行基本判断，
## 1.判断C_与X矩阵的行数是否相等；
## 2.C_矩阵是否含有NA；
## 3. nc值为C_矩阵（即X矩阵的列数）的列数，更改向量 beta的大小为 nc 个长度的向量
if(C.isNotNull()){
C_ = as<arma::mat>(C);
if(C_.n_rows != X.n_rows) throw Rcpp::exception("Number of individuals does not match for covariates.");
if(C_.has_nan()) throw Rcpp::exception("Individuals with phenotypic value should not have missing covariates.");
nc = C_.n_cols;
beta.resize(nc);
beta_store.resize(nc, n_records);
}
# 定义一个长度为nc的向量cpc
vec cpc;
if(nc){
cpc.resize(nc);
#pragma omp parallel for
# cpc 的元素为C_矩阵（即X矩阵）第 i 列的列向量的点乘
for(int i = 0; i < nc; i++){
# C_.col(i) 相当于第 i 列的列向量
cpc[i] = dot(C_.col(i), C_.col(i));
}
}
int nr = 0;
vec vr;
vec vrsd;
vec vrtmp;
vec estR;
vec estR_tmp;
vec R_idx;
vec r_RHS;
vec r_LHS;
int qr;
double dfr;
double s2r;
if(dfvr.isNotNull()){
dfr = as<double>(dfvr);
}else{
dfr = -1;
}
if(s2vr.isNotNull()){
s2r = as<double>(s2vr);
}else{
s2r = 0;
}
mat vr_store;
mat estR_store;
vector<sp_mat> Z;
vector<sp_mat> ZZ;
vector< vector<string> > Z_levels;
# 定义dfvr值，dfvf 是环境变量的自由度，默认为 -1
if(dfvr.isNotNull()){
dfr = as<double>(dfvr);
}else{
dfr = -1;
}
# 定义s2vr值，s2vr 是环境变量的标准化参数，默认为 0
if(s2vr.isNotNull()){
s2r = as<double>(s2vr);
}else{
s2r = 0;
}
mat vr_store;
mat estR_store;
# 定义三个空向量，下面有用
## sp_mat 矩阵型向量；vector<string> 字符型向量
vector<sp_mat> Z;
vector<sp_mat> ZZ;
vector< vector<string> > Z_levels;
# 定义R矩阵，R矩阵是随机效应项的设计矩阵
if(R.isNotNull()){
# 定义一个字符矩阵 R_
CharacterMatrix R_ = as<CharacterMatrix>(R);
# 判断 R_ 矩阵的行数（随机效应项内容）是否等于 X 矩阵的行数（固定效应项设计矩阵），即样本量是否相等
if(R_.nrow() != X.n_rows) throw Rcpp::exception("Number of individuals does not match for environmental random effects.");
# 判断 X 矩阵是否有 NA 值
if(xhasNA(R_)) throw Rcpp::exception("Individuals with phenotypic value should not have missing environmental random effects.");
# nr 定义为 R_ 矩阵的列数
nr = R_.ncol();
# vr 向量的长度为 nr，与R_ 矩阵的列数相等
vr.resize(nr);
# vrtmp 向量的长度为 nr，与R_ 矩阵的列数相等
vrtmp.resize(nr);
# R_idx 向量的长度为 nr + 1，且R_idx 向量的第一个元素为 -1
R_idx.resize(nr + 1); R_idx[0] = -1;
# vr_store 矩阵的长为 nr，宽为 n_records
vr_store.resize(nr, n_records);
# vary 为表型向量 y 的方差，h2 = 0.5
# fill() 函数的意义是对容器 vrtmp 所有位置都填充元素 (vary * (1 - h2) / (nr + 1))
vrtmp.fill(vary * (1 - h2) / (nr + 1));
int n_levels = 0;
for(int i = 0; i < nr; i++){
# 将R_矩阵（随机效应项的内容）的每一列向量提取出来，赋予Ri
CharacterVector Ri = R_(_, i);
# makeZ() 的作用是将随机效应项的内容转换为稀疏的设计矩阵
# 本例中为两列，将这两列随机效应项的内容转换为稀疏的设计矩阵全部赋予 Z_info
List Z_info = makeZ(Ri);
# Z_info[0] 代表稀疏设计矩阵，sp_mat 即为 sparse matrix
sp_mat z = Z_info[0];
# Z_info[1] 代表因子变量的水平
vector<string> z_levels = Z_info[1];
Z.push_back(z);
ZZ.push_back(z.t() * z);
Z_levels.push_back(z_levels);
n_levels += z_levels.size();
R_idx[i + 1] = n_levels - 1;
}
# 建立两个零矩阵estR，estR_store为矩阵
estR.zeros(n_levels);
estR_store.zeros(n_levels, n_records);
}
其中
- 函数
resize(size)作用是更改容器大小，vec a; a.resize(size)相当于将向量a的容积改为2；fill(value)的作用是对容器中的所有位置都填充元素 value；函数push_back(value)作用是将一个新的元素加到vector的最后面- 对于矩阵 X，
X.col(i)代表将 X 矩阵的第 i 个列向量取出来X矩阵 固定效应设计矩阵
若X.col(2)，相当于图中的第二列列向量（红框），例子中的X矩阵为固定效应的设计矩阵 X，该例子中固定效应矩阵 X 点乘后的结果为
- nc 为 C_ 矩阵或 X 矩阵的总列数，而C_或X矩阵是固定效应项的设计矩阵；dfvf 是环境变量的自由度，默认为 -1；s2vr 是环境变量的标准化参数，默认为 0；R矩阵是随机效应项的内容：
R矩阵 本例中，R矩阵第一列表示 loc，第二列表示 dam
makeZ()（见补充函数 1） 的作用是将随机效应项的内容（R矩阵或R_矩阵）随机效应项的内容 R 矩阵
转换为稀疏的设计矩阵随机效应项稀疏设计矩阵 上图为R矩阵或R_矩阵的第一列（也就是第一个随机效应变量loc，该例子有两个loc和dam）作为例子（上图点代表 0）；Z_info 这个 list 有两个元素：Z_info[0] 为随机效应loc的稀疏设计矩阵，Z_info[1] 为因子变量的 level
- 由于
sp_mat z = Z_info[0]; Z.push_back(z); ZZ.push_back(z.t() * z);Z 代表随机效应稀疏的设计矩阵（loc和dam）；ZZ代表每一个随机效应变量的稀疏设计矩阵（loc和dam）的点乘：
z.t() * z
接下来需要对方差-协方差矩阵进行定义
int nk = 0;
double va, vb;
double vbtmp;
double vasd, vbsd;
vec va_store;
vec vb_store;
// sp_mat k_Z, k_ZZ;
mat K;
vec Kval;
vec k_estR;
vec k_estR_store;
mat k_LHS;
vec k_estR_tmp;
vec k_RHS;
// mat ghat;
# 判断方差-协方差矩阵
if(Ki.isNotNull()){
K = as<arma::mat>(Ki);
Kval = as<arma::vec>(Kival);
if(K.n_cols != K.n_rows) throw Rcpp::exception("variance-covariance matrix should be in square.");
nk = K.n_cols;
// K_index -= 1;
va_store.resize(n_records);
vb_store.resize(n_records);
// ghat.resize(m, n_records);
// k_Z.resize(K_index.n_elem, nk);
k_estR.zeros(nk);
k_estR_store.zeros(nk);
k_estR_tmp.zeros(nk);
// for(int i = 0; i < K_index.n_elem; i++){k_Z(i, K_index[i]) = 1.0;}
// k_ZZ = k_Z.t() * k_Z;
}
int ne = 0;
double veps;
double vepstmp;
double vepssd;
double JtJ;
vec veps_store;
sp_mat epsl_Gi_, epsl_Z, epsl_ZZ;
uvec epsl_index_;
vec epsl_estR;
mat epsl_estR_store;
vec epsl_y_J_;
sp_mat epsl_LHS;
vec epsl_yadj;
vec epsl_estR_tmp;
vec epsl_RHS;
int qe = 0;
double epsl_J_beta = 0;
double epsl_J_betasd = 0;
vec epsl_J_beta_store;
if(epsl_index.isNotNull()){
epsl_index_ = as<uvec>(epsl_index);
epsl_index_ -= 1;
ne = epsl_index_.n_elem;
}
if(ne){
if(!epsl_Gi.isNotNull()) throw Rcpp::exception("variance-covariance matrix should be provided for epsilon term.");
epsl_Gi_ = as<sp_mat>(epsl_Gi);
epsl_y_J_ = as<vec>(epsl_y_J);
if(epsl_Gi_.n_cols != epsl_Gi_.n_rows) throw Rcpp::exception("variance-covariance matrix should be in square.");
JtJ = dot(epsl_y_J_, epsl_y_J_);
veps_store.resize(n_records);
epsl_J_beta_store.resize(n_records);
epsl_Z.resize(ne, epsl_Gi_.n_cols);
epsl_estR.zeros(epsl_Gi_.n_cols);
epsl_estR_tmp.zeros(epsl_Gi_.n_cols);
epsl_estR_store.resize(epsl_Gi_.n_cols, n_records);
epsl_yadj.resize(ne);
qe = epsl_Gi_.n_cols;
for(int i = 0; i < ne; i++){epsl_Z(i, epsl_index_[i]) = 1.0;}
epsl_ZZ = epsl_Z.t() * epsl_Z;
}
从这里开始，将介绍算法的核心部分：
int count = 0;
int nzct = 0;
int inc = 1;
int n_fold = fold_.n_elem;
int NnzSnp, indistflag;
double doc = 1.0;
double xx, oldgi, gi, gi_, rhs, lhs, logdetV, acceptProb, uhat, v;
double vara_, dfvara_, s2vara_, vare_, dfvare_, s2vare_, vargi, hsq, s2varg_;
// double sumvg;
vec snptracker;
vec nzrate;
# 对模型进行判断
if(model == "BayesRR" || model == "BayesA" || model == "BayesL"){
NnzSnp = m;
Pi[0] = 0; Pi[1] = 1;
fixpi = true;
}else{
if(model != "BayesR" && Pi.n_elem != 2) throw Rcpp::exception("length of Pi should be 2, the first value is the proportion of non-effect markers.");
nzrate.zeros(m);
snptracker.zeros(m);
}
vec g = zeros(m);
mat g_store = zeros(m, n_records);
vec u = zeros(n);
vec xpx = zeros(m);
vec vx = zeros(m);
vec mu_store = zeros(n_records);
// vec sumvg_store = zeros(niter - nburn + 1);
vec vara_store = zeros(n_records);
vec vare_store = zeros(n_records);
vec hsq_store = zeros(n_records);
mat pi_store;
pi_store.resize(n_fold, n_records);
#pragma omp parallel for
for(int i = 0; i < m; i++){
# 计算 X 矩阵每一个列向量的平方和与方差，
# X 矩阵为固定效应的设计矩阵
vec Xi = X.col(i);
xpx[i] = sum(square(Xi));
vx[i] = var(Xi);
}
# 计算 X 矩阵列向量方差和
double sumvx = sum(vx);
int nvar0 = sum(vx == 0);
# 判断其他参数，并赋值
if(dfvg.isNotNull()){
dfvara_ = as<double>(dfvg);
}else{
dfvara_ = 4;
}
if(dfvara_ <= 2){
throw Rcpp::exception("dfvg should not be less than 2.");
}
if(vg.isNotNull()){
vara_ = as<double>(vg);
}else{
vara_ = ((dfvara_ - 2) / dfvara_) * vary * h2;
}
vepstmp = vara_;
vbtmp = vara_;
if(ve.isNotNull()){
vare_ = as<double>(ve);
}else{
vare_ = vary * (1 - h2) / (nr + 1);
}
if(dfve.isNotNull()){
dfvare_ = as<double>(dfve);
}else{
dfvare_ = -2;
}
if(s2vg.isNotNull()){
s2vara_ = as<double>(s2vg);
}else{
s2vara_ = vara_ * (dfvara_ - 2) / dfvara_;
}
double varg = vara_ / ((1 - Pi[0]) * sumvx);
s2varg_ = s2vara_ / ((1 - Pi[0]) * sumvx);
if(s2ve.isNotNull()){
s2vare_ = as<double>(s2ve);
}else{
s2vare_ = 0;
}
if(niter < nburn){
throw Rcpp::exception("Number of total iteration ('niter') shold be larger than burn-in ('nburn').");
}
double R2 = (dfvara_ - 2) / dfvara_;
double lambda2 = 2 * (1 - R2) / (R2) * sumvx;
double lambda = sqrt(lambda2);
double shape, shape0 = 1.1;
double rate, rate0 = (shape0 - 1) / lambda2;
vec vargL;
if(model == "BayesL"){
vargL.resize(m);
vargL.fill(varg);
}
vec stemp = zeros(n_fold);
vec fold_snp_num = zeros(n_fold);
vec logpi = zeros(n_fold);
vec s = zeros(n_fold);
vec vara_fold = (vara_ / ((1 - Pi[0]) * sumvx)) * fold_;
vec vare_vara_fold = zeros(n_fold);
# 输出的日志文件，分别解释每个参数的意义
if(verbose){
Rcpp::Rcout.precision(4);
Rcpp::Rcout << "Prior parameters:" << std::endl;
Rcpp::Rcout << " Model fitted at [" << (model == "BayesRR" ? "Bayes Ridge Regression" : model) << "]" << std::endl;
Rcpp::Rcout << " Number of observations " << n << std::endl;
if(epsl_index.isNotNull()){
Rcpp::Rcout << " Observations with genotype " << (n - ne) << std::endl;
Rcpp::Rcout << " Observations with imputed genotype " << ne << std::endl;
}
Rcpp::Rcout << " Number of covariates " << (nc + 1) << std::endl;
Rcpp::Rcout << " Number of envir-random effects " << nr << std::endl;
Rcpp::Rcout << " Number of markers " << m << std::endl;
for(int i = 0; i < Pi.n_elem; i++){
if(i == 0){
Rcpp::Rcout << " π for markers in zero effect size " << Pi[i] << endl;
}else{
if(i == 1){
Rcpp::Rcout << " π for markers in non-zero effect size " << Pi[i] << " ";
}else{
Rcpp::Rcout << Pi[i] << " ";
}
}
}
Rcpp::Rcout << std::endl;
if(model == "BayesR"){
Rcpp::Rcout << " Group fold ";
for(int i = 0; i < fold_.n_elem; i++){
Rcpp::Rcout << fold_[i] << " ";
}
Rcpp::Rcout << std::endl;
}
Rcpp::Rcout << " Total number of iteration " << niter << std::endl;
Rcpp::Rcout << " Total number of burn-in " << nburn << std::endl;
Rcpp::Rcout << " Frequency of collecting " << thin << std::endl;
Rcpp::Rcout << " Phenotypic var " << std::fixed << vary << std::endl;
if(nr){
Rcpp::Rcout << " Environmental var ";
for(int i = 0; i < nr; i++){
Rcpp::Rcout << std::fixed << vrtmp[i] << " ";
}
Rcpp::Rcout << std::endl;
}
Rcpp::Rcout << " Genetic var " << std::fixed << vara_ << std::endl;
Rcpp::Rcout << " Inv-Chisq gpar " << std::fixed << dfvara_ << " " << std::fixed << s2vara_ << std::endl;
Rcpp::Rcout << " Residual var " << std::fixed << vare_ << std::endl;
Rcpp::Rcout << " Inv-Chisq epar " << std::fixed << dfvare_ << " " << std::fixed << s2vare_ << std::endl;
Rcpp::Rcout << " Marker var " << std::fixed << varg << std::endl;
Rcpp::Rcout << " Inv-Chisq alpar " << std::fixed << dfvara_ << " " << std::fixed << s2varg_ << std::endl;
if(WPPA){
Rcpp::Rcout << " Number of windows for GWAS analysis " << nw << std::endl;
}
Rcpp::Rcout << "MCMC started: " << std::endl;
Rcpp::Rcout << " Iter" << " ";
Rcpp::Rcout << "NumNZSnp" << " ";
for(int i = 0; i < n_fold; i++){
Rcpp::Rcout << "π" << i + 1 << " ";
}
if(model == "BayesL") Rcpp::Rcout << "Lambda" << " ";
if(model == "BSLMM") Rcpp::Rcout << "Va" << " " << "Vb" << " ";
for(int i = 0; i < nr; i++){
Rcpp::Rcout << "Vr" << (i + 1) << " ";
}
if(ne) Rcpp::Rcout << "Vε" << " ";
Rcpp::Rcout << "Vg" << " ";
Rcpp::Rcout << "Ve" << " ";
Rcpp::Rcout << "h2" << " ";
Rcpp::Rcout << "Timeleft" << std::endl;
}
其中在
if(verbose)的参数：
- Number of observations的参数 n 代表所观察到有表型值的个体数
int n = y.n_elem;- Observations with genotype的参数 (n - ne) 代表所观察到的基因型数量，默认下的 ne = 0，因此基因型数目与所观察到有表型值的个体数相等；
int ne = 0;- Observations with imputed genotype的参数 ne 代表的是由于缺少基因型而做的impute，默认下的 ne = 0，
int ne = 0;- Number of covariates 的参数 (nc + 1) 代表的是协方差（协变量）的数量，默认下的 nc = 0，即
int nc = 0;，相当于协方差的数量为 1- Number of envir-random effects 的参数 nr 代表随机效应（或者环境效应）的数量，默认下 nr 等于随机效应内容矩阵 R 矩阵的列数（R 矩阵见下）
- Number of markers 的参数 m 代表的是固定效应矩阵的设计矩阵变量个数；
int m = X.n_cols;- π for markers in zero effect size 代表零SNP效应的比例
Pi[0]- π for markers in non-zero effect size 代表非零SNP效应的比例
Pi[1]- Total number of iteration 代表的是MCMC的迭代次数，默认 niter = 50000
- Total number of burn-in 代表的是
- Frequency of collecting 代表的是采样频率，每相隔thin轮采一次样
- Phenotypic var 代表的是表型值的方差，
double vary = var(y);- Environmental var 代表的是环境变量的方差：
vec vrtmp; vrtmp.resize(nr); vrtmp.fill(vary * (1 - h2) / (nr + 1)); ## h2默认等于 0.5环境变量的方差 vrtmp 为长度为 nr （见第5条）的向量，默认的环境变量所有的值相等（vrtmp与表型值方差vary有关），等于
vary * (1 - h2) / (nr + 1)
- Genetic var 的参数 vara_ 代表遗传方差，是一个双精型数值：
double vara_ vara_ = ((dfvara_ - 2) / dfvara_) * vary * h2; # 默认下h2为0.5，dfvara_ 为 4显然遗传方差与表型值方差vary有关
- Inv-Chisq gpar 的参数 dfvara_ 代表遗传方差的自由度，默认为 4，
dfvara_ = 4- Residual var 的参数 vare_ 代表残差方差的先验值，
vare_ = vary * (1 - h2) / (nr + 1);，其中 vary 代表表型值的方差，h2 = 0.5，nr 为随机效应内容矩阵 R 矩阵的列数（见第5条），该例子中 nr = 2- Inv-Chisq epar 代表残差方差的自由度，默认为 -2，
dfvare_ = -2;- Marker var 的参数 varg 代表固定效应的方差，默认
double varg = vara_ / ((1 - Pi[0]) * sumvx);，其中 vara_ 代表遗传方差（参见第14条），Pi[0] 代表零SNP效应的比例（见第7条），sumvx 代表固定效应矩阵列向量的方差之和for(int i = 0; i < m; i++){ vec Xi = X.col(i); xpx[i] = sum(square(Xi)); vx[i] = var(Xi); } double sumvx = sum(vx);因此固定效应的方差与遗传方差和固定效应的设计矩阵列向量方差有关
MCMC 估计参数部分（根据示例数据的内容，代码做了部分删减）
# 进行时间的计算
MyTimer timer;
timer.step("start");
MyTimer timer_loop;
timer_loop.step("start");
double tt0;
int tt, hor, min, sec;
# 定义表型观测值变量 y 的平均值 mu_, mu
double mu_, mu = mean(y);
# 定义一个与表型值的个体数长度相等的向量，向量元素值全为 1
vec one(n); one.fill(1);
# 定义一个减去表型均值后的表型观察值向量，即每个表型观测者减去其均值
vec yadj = y - mu;
double* dyadj = yadj.memptr();
double* dxi;
# 这里的niter代表总迭代轮数，iter代表迭代到第几轮
for(int iter = 0; iter < niter; iter++){
// sample intercept
# 对截距 μ 进行采样
## 每迭代一轮需要对表型值的均值重新进行正态分布采样
## 函数 norm_sample() 见补充函数 2
## 这一步的作用是想通过多次采样获得稳定的表型值均值
## mu_ 满足均值为 sum(yadj) / n，方差为 sqrt(vare_ / n) 的正态分布
## 代码 sum(yadj) / n 为 yadj（减去表型均值后的表型观察值向量）的均值
# 以下的 mu_ 可以理解为表型向量均值的估计值
mu_ = - norm_sample(sum(yadj) / n, sqrt(vare_ / n));
mu -= mu_;
# 初始的截距值为 dyadj，初始截距 dyadj = mu_ + y - mu (dyadj, 减去表型均值后的表型观察值向量)，即 y*
# 本质上此时的 dyadj 代表表型向量的估计值，详见下面第 7 条
daxpy_(&n, &mu_, one.memptr(), &inc, dyadj, &inc);
# 对固定效应的系数进行MCMC采样，参见理论的第一部分和公式 7
for(int i = 0; i < nc; i++){
# 这里的C_矩阵（或C矩阵）为固定效应设计矩阵 X
# 将 C_矩阵（或C矩阵，即固定效应设计矩阵 X）的每一列转换为二进制
dci = C_.colptr(i);
# 向量 beta 的中所有元素的初始值为 0，这里的 oldgi 代表上一轮采样的系数 β
oldgi = beta[i];
# 即cpc 的元素为C_矩阵（即X矩阵）第 i 列的列向量的点乘
# 即cpc[i] 代表C_矩阵（即X矩阵）第 i 列向量的点乘
v = cpc[i];
# 计算固定效应设计矩阵 X （C 矩阵）的每一列向量与初始截距（即为表型值向量）向量的点积
rhs = ddot_(&n, dci, &inc, dyadj, &inc);
rhs += v * oldgi;
# 这里的变量 gi 代表新采样的 β 值
gi = norm_sample(rhs / v, sqrt(vare_ / v));
# 计算前后两次采样的 β 值的误差，gi 代表新采样的 β 值；oldgi 代表上一次采样的 β 值
# gi_ 代表前后两次采样的 β 值的误差
gi_ = oldgi - gi;
# 计算最优化的目标函数，见理论的第一部分和公式 8
# 和下面补充说明第 8 条
daxpy_(&n, &gi_, dci, &inc, dyadj, &inc);
# 新获得的 gi 值赋予给向量 beta 的对应元素
beta[i] = gi;
}
# 对随机效应的系数进行MCMC采样，参见理论的第一部分
for(int i = 0; i < nr; i++){
estR_tmp = estR.subvec(R_idx[i] + 1, R_idx[i + 1]);
qr = estR_tmp.n_elem;
# 计算公式6中，随机效应正态分布均值的分子部分，Z 和 ZZ 详细见下面第 9 条
r_RHS = Z[i].t() * yadj;
# 其中这里的 estR_tmp 相当于前一轮采样的随机效应系数 rj 的估计值
r_RHS += ZZ[i] * estR_tmp;
# 计算公式6中，随机效应正态分布均值的分母部分
# vare_ / vrtmp[i] 可以解释为遗传力
r_LHS = diagvec(ZZ[i]) + vare_ / vrtmp[i];
# 进行随机效应系数的采样，estR_tmp 代表待采样的随机效应系数 rj 的估计值
# 参数 qr 代表随机效应内容 R 矩阵的列数（即一共有 qr 个随机效应项）
# estR_tmp[qi] 代表第 qi 个随机效应的系数 ( r )
for(int qi = 0; qi < qr; qi++){
# r_RHS[qi] / r_LHS[qi] 为正态分布的方差
estR_tmp[qi] = norm_sample(r_RHS[qi] / r_LHS[qi], sqrt(vare_ / r_LHS[qi]));
}
# 计算最优化的目标函数，见理论的第一部分和公式 8
vec diff = Z[i] * (estR.subvec(R_idx[i] + 1, R_idx[i + 1]) - estR_tmp);
# 计算最优化的目标函数，详细说明见下面第10条
daxpy_(&n, &doc, diff.memptr(), &inc, dyadj, &inc);
vrtmp[i] = (ddot_(&qr, estR_tmp.memptr(), &inc, estR_tmp.memptr(), &inc) + s2r * dfr) / chisq_sample(qr + dfr);
# 计算随机效应系数 r 的方差
vr[i] = var(estR_tmp);
estR.subvec(R_idx[i] + 1, R_idx[i + 1]) = estR_tmp;
}
# 对Genetic marker effects的系数 α 进行采样，即基因型矩阵 M
# 同时计算遗传方差 Vg
# 这里先考虑 BayesCpi 方法的参数估计
switch(model_index){
case 4:
logpi = log(Pi);
# 这里的初始 s 为：vec s = zeros(n_fold)；默认的 n_fold = 2
# 其中 logpi[0] 代表 logpi 的第一个元素零效应SNP的比例
s[0] = logpi[0];
// sumvg = 0;
vargi = 0;
for(int i = 0; i < m; i++){
if(!vx[i]) continue;
# 这一块的变量 X 为基因型矩阵 M，只不过变量定义为X
# 切勿与固定效应的设计矩阵 X 混淆
# 将基因型矩阵 M 的第 i 列取出
dxi = X.colptr(i);
# 变量 xpx 为基因型矩阵 M 的每一个列向量的平方和所构成的集合
xx = xpx[i];
# 这里的 g[i] 初始值为 vec g = zeros(m); 长度为 m (m 为 M 矩阵的列长，元素全为 0)
oldgi = g[i];
# 这一块详细见下面的第 12 条
rhs = ddot_(&n, dxi, &inc, dyadj, &inc);
# 这里的 oldgi 相当于上一轮采样的 α，初始的 α = 0
if(oldgi) rhs += xx * oldgi; // 首次采样时，oldgi = 0，拒绝采样时 oldgi = 0
lhs = xx / vare_;
logdetV = log(varg * lhs + 1);
uhat = rhs / (xx + vare_ / varg);
# 这一块计算的是公式 11 的内容，目的是构造接收概率 acceptProb
s[1] = -0.5 * (logdetV - (rhs * uhat / vare_)) + logpi[1];
# 拒绝-接收采样，定义接收率
acceptProb = 1 / sum(exp(s - s[0]));
# 参考下面第 12 条
# indistflag 作用是决定本次采样是否被接收
# 从均匀分布随机取一个值与 acceptProb 比大小，acceptProb 的计算参考公式 12
# 若从均匀分布中随机取的一个值小于 acceptProb，indistflag 值被则赋予 0，否则被赋予 1
indistflag = (uniform_sample()) < acceptProb ? 0 : 1;
snptracker[i] = indistflag;
# 接收采样
if(indistflag){
# 这里对 genetic marker effect 的系数进行采样
v = xx + vare_ / varg;
# 采样分布参考公式 9 采样分布部分，详细介绍参考下面第 13 条
# 这里的 gi 相当于本轮采样的 α
gi = norm_sample(rhs / v, sqrt(vare_ / v));
# 定义最优化函数参考公式 9 的最优化函数部分
gi_ = oldgi - gi;
daxpy_(&n, &gi_, dxi, &inc, dyadj, &inc);
// sumvg += (gi * gi) * vx[i];
vargi += (gi * gi);
# 拒绝采样
}else{
# 拒绝采样的最优化函数
gi = 0;
if(oldgi){
gi_ = oldgi;
daxpy_(&n, &gi_, dxi, &inc, dyadj, &inc);
}
}
g[i] = gi;
}
# 详细解释查看下面第 14 条
fold_snp_num[1] = sum(snptracker);
fold_snp_num[0] = m - nvar0 - fold_snp_num[1];
NnzSnp = fold_snp_num[1];
varg = (vargi + s2varg_ * dfvara_) / chisq_sample(dfvara_ + NnzSnp);
// varg = invchisq_sample(NnzSnp + dfvara_, varg);
if(nk) va = varg;
# BayesCpi 的方法执行下面语句，!fixpi = TRUE，n_fold = 2
# fold_snp_num 代表当前轮次下拒绝-接收采样中，接收采样的次数
# 设 Πk 满足 dirichlet 分布
if(!fixpi) Pi = rdirichlet_sample(n_fold, (fold_snp_num + 1));
break;
}
// collect the parameters to obtain posterior estimation
# 获得采样的数据，并不是每一个 iter 都采样，而是隔 nburn 采样一次
if(iter >= nburn && (iter + 1 - nburn) % thin == 0){
# 储存采样的均值（截距） μ，count 为记录采样的次数
mu_store[count] = mu;
if(!fixpi) pi_store.col(count) = Pi;
# 储存采样的基因型系数矩阵 α 的方差
vara_store[count] = vara_;
# 储存采样的基因型系数矩阵 α 的值
g_store.col(count) = g;
# 储存采样的固定效应系数矩阵 β 的值
if(nc){
beta_store.col(count) = beta;
}
# 储存采样的随机效应系数矩阵 r 的值
if(nr){
# 储存采样的随机效应系数矩阵 r 值的方差
vr_store.col(count) = vr;
# 储存采样的随机效应系数矩阵 r 的值
estR_store.col(count) = estR;
}
count++;
}
}
上述代码中出现的 oldgi 代表上一轮采样的参数（β，α等；固定效应，genetic marker effect 的系数）；oldgi 代表本轮采样的参数（β，α等；固定效应，genetic marker effect 的系数）； estR_tmp[qi] 代表第 qi 个随机效应的系数 ( r )
其中：
- 上述代码中出现的 C_或C 代表固定效应设计矩阵 X
- BLAS 库是cpp中一个进行线性代数计算的库，用这个库进行计算时首先要将矩阵或者向量转换为二进制：
# 将向量转换为二进制 double* dyadj = yadj.memptr(); dci = C_.colptr(i); # 计算点积 rhs = ddot_(&n, dci, &inc, dyadj, &inc);而函数
ddot_(N, X, INCX, Y, INCY)的作用是用于计算两个向量的点积，参数 N 代表向量的长度；参数 X 代表待计算的向量 A；参数 INCX 代表待计算的向量 A 元素之间的间距（默认为1）；参数 Y 代表待计算的向量 B；参数 INCY 代表待计算的向量 B 元素之间的间距（默认为1）；
- 函数
daxpy_()的作用是线性化相加，daxpy_(const MKL_INT n, const double a, const double *x, const MKL_INT incx, double *y, const MKL_INT incy);，其中 const MKL_INT n 代表待输入向量的长度，const double a 代表系数 a，const double * x 代表待输入的向量 x，const MKL_INT incx 代表参数 INCX 代表待计算的向量 x 元素之间的间距（默认为1），double * y 代表待计算的向量 y，const MKL_INT incy 代表参数 INCY 代表待计算的向量 y 元素之间的间距（默认为1）；因此函数所作运算为：ynew = a*x + y，然后将新的 y（ynew）赋予 y- 在对计算随机效应系数的mcmc采样过程中
相当于计算公式6中正态分布均值的分子部分，即r_RHS = Z[i].t() * yadj; r_RHS += ZZ[i] * estR_tmp;
- 在对计算随机效应系数的mcmc采样过程中
相当于计算公式6中正态分布均值的分母部分，r_LHS = diagvec(ZZ[i]) + vare_ / vrtmp[i];
- 在对计算随机效应系数的mcmc采样过程中，
vare_ / r_LHS[qi]相当于
- 为什么进行以下操作：
double mu_, mu = mean(y); vec yadj = y - mu; double* dyadj = yadj.memptr(); mu_ = - norm_sample(sum(yadj) / n, sqrt(vare_ / n)); mu -= mu_; daxpy_(&n, &mu_, one.memptr(), &inc, dyadj, &inc);笔者认为经过上述代码运行后，初始截距
dyadj = mu_ + y - mu（而 mu_new = mu - mu_采样后，则mu_new表示表型向量均值的真实值与估计值之间的误差）的目的是创造一个表型向量的估计值，即 y*
- 对于固定效应系数的最优化目标函数：
daxpy_(&n, &gi_, dci, &inc, dyadj, &inc);的目的是计算最优化的目标函数，见理论的第一部分和公式 8，相当于执行了公式 8，即dyadj = dyadj + dci(oldgi - gi)，其中 oldgi - gi 代表 βji-1 - βji，dci 代表 Xj，dyadj 代表 y*- Z 和 ZZ 分别代表什么？Z 代表随机效应稀疏的设计矩阵的list（loc和dam），例如 Zloc为
是一个 300 X 50 的稀疏矩阵（300 为随机效应的R矩阵的行数，50为随机效应项 loc 的因子水平为 50），点代表 0 ，loc 的因子水平为 50：
loc 的因子水平
而 Zdam 为是一个 300 X 150 的稀疏矩阵（300 为随机效应的R矩阵的行数，150为随机效应项 dam 的因子水平为 150），点代表 0 ，dam 的因子水平为 150：
dam 的因子水平 而
ZZ_loc = t(Z_loc) %*% Z_loc；ZZ_dam = t(Z_dam) %*% Z_dam- 对于随机效应系数的最优化目标函数：
vec diff = Z[i] * (estR.subvec(R_idx[i] + 1, R_idx[i + 1]) - estR_tmp); daxpy_(&n, &doc, diff.memptr(), &inc, dyadj, &inc); # doc = 1.0其中
vec diff代表而 dyadj 代表 y*
- 基因型矩阵M（为了区别于固定效应的设计矩阵 X，基因型矩阵在全文中仅称为 M 矩阵），M 矩阵，即geno：数字基因型矩阵（n × m, n为个体数，m为snp数，本例中为 600 × 1000），接收 0, 1, 2 或者 -1, 0, 1 这两种形式的矩阵，但矩阵中不能含有NA，
- 有关于公式11的代码部分解释：
rhs = ddot_(&n, dxi, &inc, dyadj, &inc); if(oldgi) rhs += xx * oldgi; # xx 为 M 矩阵其中一列的平方和上述代码计算的 rhs
这里的
oldgi = g[i]表示：g[i] 初始值为vec g = zeros(m);长度为 m（m 为 M 矩阵的列长，元素全为 0）；uhat = rhs / (xx + vare_ / varg);代表的是
logdetV = log(varg * lhs + 1);计算的是而
s[1] = -0.5 * (logdetV - (rhs * uhat / vare_)) + logpi[1];计算的是公式11
这个概率值就为接收函数或者接收值；一般情况下，nΠ = 2，第一个元素表征具有零效应的SNP比例，第二个元素表征具有非零效应的SNP比例。每次采样的时候从均匀分布中随机取一个数字 num 与上图的 acceptProb 作比较，若 num < acceptProb 接收采样 ，否则拒绝采样；该例子由于Pi = c(0.98, 0.02)，因此nk = 2;Π0 = 0.98;Π1 = 0.02，Π0 = 0.98代表具有零效应的SNP比例；Π1 = 0.02代表具有非零效应的SNP比例
- 对系数 α 的采样介绍：
gi = norm_sample(rhs / v, sqrt(vare_ / v));计算的采样分布为接下来计算最优化的函数
上述代码则表示最优化目标的函数：gi_ = oldgi - gi; daxpy_(&n, &gi_, dxi, &inc, dyadj, &inc);
- 其中 Πk（ k > 0，Π0 在该例子中为 0.98），n_fold = 2，fold_snp_num 代表当前轮次下拒绝-接收采样中，接收采样的次数（当前轮次，比方说 iter = 3，接收采样的次数为 2，则 fold_snp_num = 2；iter = 10，接收采样的次数为 6，则 fold_snp_num = 6），且 Πk 满足 dirichlet 分布
η 即为fold_snp_num，代表当前轮次下拒绝-接收采样中，接收采样的次数（当前轮次，比方说 iter = 3，接收采样的次数为 2，则 η = 2；iter = 10，接收采样的次数为 6，则 η = 6）
收集数据，返回结果（根据示例数据做了一定删减）
# 定义结果的 list
List results;
List MCMCsample;
# 储存随机效应系数矩阵的方差
if(nr){
# 取多次迭代的均值
vr = conv_to<vec>::from(mean(vr_store, 1));
vrsd = conv_to<vec>::from(stddev(vr_store, 0, 1));
results["Vr"] = vr;
MCMCsample["Vr"] = vr_store;
}
# 储存基因型系数矩阵的方差，取多次迭代的均值
vara_ = mean(vara_store);
double varasd = stddev(vara_store);
results["Vg"] = vara_;
MCMCsample["Vg"] = vara_store.t();
# 储存均值（截距）的方差，取多次迭代的均值
double Mu = mean(mu_store);
vec e = y - (Mu * one);
double musd = stddev(mu_store);
results["mu"] = Mu;
MCMCsample["mu"] = mu_store.t();
if(nc){
# 储存随机效应系数矩阵的方差，取多次迭代的均值
beta = conv_to<vec>::from(mean(beta_store, 1));
betasd = conv_to<vec>::from(stddev(beta_store, 0, 1));
# 储存随机效应系数矩阵 β 值
results["beta"] = beta;
MCMCsample["beta"] = beta_store;
}
# 储存基因型系数矩阵 α 值，取多次迭代的均值
g = conv_to<vec>::from(mean(g_store, 1));
results["alpha"] = g;
MCMCsample["alpha"] = g_store;
# # 储存非零效应SNP的比例 Π 值，取多次迭代的均值
vec pisd;
if(!fixpi){
Pi = conv_to<vec>::from(mean(pi_store, 1));
pisd = conv_to<vec>::from(stddev(pi_store, 0, 1));
}else{
pisd = zeros(Pi.n_elem);
pi_store.row(0).fill(Pi[0]);
pi_store.row(1).fill(Pi[1]);
}
results["pi"] = Pi;
MCMCsample["pi"] = pi_store;
# 储存随机效应系数矩阵 r 值，取多次迭代的均值
if(nr){
List listr(2);
estR = conv_to<vec>::from(mean(estR_store, 1));
vector<string> r_levers;
for(int i = 0; i < nr; i++){
for(int j = 0; j < (Z_levels[i]).size(); j++){
r_levers.push_back(Z_levels[i][j]);
}
}
listr[0] = wrap(r_levers.begin(), r_levers.end());
listr[1] = wrap(estR);
DataFrame r = listr;
Rcpp::CharacterVector names(2);
names[0] = "Levels";
names[1] = "Estimation";
r.attr("names") = names;
results["r"] = r;
MCMCsample["r"] = estR_store;
}
补充函数
- makeZ()：makeZ() 的作用是将随机效应项的内容（R矩阵或R_矩阵），转换为稀疏的设计矩阵
Rcpp::List makeZ(
const CharacterVector &R
)
{
int n = R.length();
std::vector<std::string> levels;
std::vector<std::string> value = Rcpp::as<std::vector<std::string> >(R);
stable_sort(value.begin(), value.end());
value.erase(unique(value.begin(), value.end()), value.end());
if (value.size() == n){
throw Rcpp::exception("number of class of environmental random effects should be less than population size.");
}
if(value.size() == 1){
throw Rcpp::exception("number of class of environmental random effects should be bigger than 1.");
}
std::map<string, int> val_map;
for (int j = 0; j < value.size(); j++){
levels.push_back(value[j]);
val_map.insert(pair<string, int>(value[j], j));
}
map<string, int>::iterator iter;
arma::sp_mat z(n, value.size());
for (int j = 0; j < n; j++) {
std::string v = Rcpp::as<string>(R[j]);
iter = val_map.find(v);
z(j, iter->second) = 1.0;
}
return List::create(Named("z") = z, Named("levels") = levels);
}
- norm_sample() 的作用是进行正态分布的采样：
double norm_sample(double mean, double sd){
// return R::rnorm(mean, sd);
return mean + sd * norm_rand();
}
补充的代数知识：
如果将每一个 βj，rj 或 αj 看作是随机变量（假设 j 代表第 j 列列向量）：
一般地推导得：
网友评论