關(guān)于BayesCpi

本文在進(jìn)行基因型系數(shù)矩陣估計(jì)時(shí)采用的是 BayesCpi 方法

理論部分

第一部分：關(guān)于固定效應(yīng)和隨機(jī)效應(yīng)系數(shù)矩陣的采樣分布定義
根據(jù)混合線性模型的系數(shù)方程求解的方式，有：

公式1

公式2

公式3

公式4

公式5

更一般的：

公式6

如果把固定效應(yīng)或者隨機(jī)效應(yīng)設(shè)計(jì)矩陣的每一個(gè)變量看作是一個(gè)隨機(jī)變量
如果把固定效應(yīng)或者隨機(jī)效應(yīng)設(shè)計(jì)矩陣的每一個(gè)變量（每一列）看作是一個(gè)隨機(jī)變量，那么更一般地有：

公式7

因此，每次采樣完做最優(yōu)化的目標(biāo)函數(shù)為：

公式8

第二部分：關(guān)于 genetic marker effect 系數(shù)矩陣的采樣分布定義
筆者認(rèn)為genetic marker effect相當(dāng)于一項(xiàng)隨機(jī)效應(yīng)項(xiàng)，因此按照隨機(jī)效應(yīng)系數(shù)矩陣計(jì)算公式BLUP進(jìn)行求解士葫，根據(jù)第一部分公式1-6的計(jì)算菩颖，根據(jù)隨機(jī)效應(yīng)系數(shù)矩陣的計(jì)算方法，同理可推得genetic marker effect的系數(shù)矩陣 α 的采樣分布為：

公式9

M代表genetic marker effect的設(shè)計(jì)矩陣为障，α 代表 genetic marker effect 的系數(shù)矩陣，K^-1λ 代表遺傳力矩陣（在該例子中 K^-1λ = I(σ_e²/σ_α²)）有：

公式11

構(gòu)造接收函數(shù)的方法如下：
極大似然的目標(biāo)參數(shù)為Π_k（Π_k：the proportion of markers in zero effect size）鞋喇，依據(jù)極大似然估計(jì)的原理，我們構(gòu)造參數(shù)的似然函數(shù) L_Πk 如下：（參考文章：https://www.nature.com/articles/s41467-019-12653-0#Sec18眉撵，補(bǔ)充材料的45頁）

公式11

公式12

在mcmc過程中，Π_k（k > 0）服從 dirichlet 分布，η 代表當(dāng)前輪次下拒絕-接收采樣中矾睦，接收采樣的次數(shù)（當(dāng)前輪次晦款，比方說 iter = 3，接收采樣的次數(shù)為 2枚冗，則 η = 2缓溅；iter = 10，接收采樣的次數(shù)為 6赁温，則 η = 6）

公式13

代碼部分

首先介紹下所用到的數(shù)據(jù)：

1.pheno： 一些表型信息和metadata的內(nèi)容坛怪，

pheno

包括測序樣品的id，性別束世，采樣季節(jié)等基本信息酝陈，T1可能是表型信息（連續(xù)型）
2.geno：數(shù)字基因型矩陣（n × m, n為個(gè)體數(shù)，m為snp數(shù)毁涉，本例中為 600 × 1000）沉帮，接收 0, 1, 2 或者 -1, 0, 1 這兩種形式的矩陣，但矩陣中不能含有NA

geno矩陣

3.geno.id：代表的是geno矩陣的樣品id名稱贫堰，該例子中g(shù)eno.id的長度為600

geno.id

4.map：類似于VCF文件穆壕，列名依次表示為SNP名稱，染色體編號(hào)其屏，SNP位置和基因型

map

該例子中基因型的編碼規(guī)則為A1A1編碼為2喇勋，A1A2編碼為1，A2A2編碼為0偎行，其中A1代表SNP上第一個(gè)等位基因川背，即參考等位基因（通常出現(xiàn)頻率較多）；A2代表SNP上第二個(gè)等位基因蛤袒，即替代等位基因（通常出現(xiàn)頻率較少）

第二部分熄云，輸入的介紹：

# 建立的混合線性模型
formula = T1 ~ season + bwt + (1 | loc) + (1 | dam)
#  一些表型信息和metadata的內(nèi)容
data = pheno
# 數(shù)字基因型矩陣
M = geno
# geno矩陣的樣品id名稱
M.id = geno.id
# 方法選擇
method = "BayesCpi"
# 類似于VCF文件，見 4
map = map
# 所設(shè)定的零效應(yīng)SNP和有效應(yīng)的SNP的比例（零效應(yīng)的SNP比例為0.98妙真，有效應(yīng)的SNP的比例為0.02）
Pi = c(0.98, 0.02)
# MCMC 估計(jì)參數(shù)時(shí)的最大迭代次數(shù)
niter = 20000
# nburn 與 thin 和估計(jì)參數(shù)有關(guān)
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)
# 將隨機(jī)效應(yīng)項(xiàng) loc 和 dam 提取出來
rand_term <- unlist(str_extract_all(myformula, "(?<=(\\+ |~)\\(1 \\| )[:\\w\\d]+(?=\\))"))
R <- NULL

# 將隨機(jī)效應(yīng)項(xiàng)對(duì)應(yīng)的內(nèi)容提取出來
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)
}

# 將固定效應(yīng)項(xiàng)對(duì)應(yīng)的內(nèi)容提取出來
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)

其中：

1. myformula 在該例子中代表 T1 ~ season + bwt + (1 | loc) + (1 | dam)缴允，這個(gè)式子表示 表型T1作為響應(yīng)變量，固定效應(yīng)項(xiàng)為season和bwt珍德，隨機(jī)效應(yīng)則不考慮隨機(jī)效應(yīng)項(xiàng)的斜率练般，只考慮隨機(jī)效應(yīng)項(xiàng)的截距，分組的隨機(jī)因子分別是loc和dam
   2.上述循環(huán)中锈候，將隨機(jī)效應(yīng)項(xiàng)對(duì)應(yīng)的內(nèi)容提取出來在該例子中指的是!
   pheno
   將紅框中的內(nèi)容提取出來
   3.最終的 fixed_formula 提取出來為字符串 T1 ~ season + bwt

# 對(duì)應(yīng)多張表的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.")

# 將固定效應(yīng)按條件轉(zhuǎn)換為設(shè)計(jì)矩陣
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
# 獲得隨機(jī)效應(yīng)項(xiàng)的內(nèi)容
R <- R[!yNA, , drop = FALSE]

# 提取對(duì)應(yīng)的響應(yīng)變量值
y <- data[!yNA, as.character(formula[2]), drop = TRUE]
if(!is.matrix(M)){M <- as.matrix(M); gc()}
    
# 將id對(duì)不上的 geno：數(shù)字基因型矩陣賦給 Mp
Mp <- M[yNA, , drop=FALSE]
# 將id對(duì)上的 geno：數(shù)字基因型矩陣賦給 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)

其中：

1. 固定效應(yīng)的設(shè)計(jì)矩陣 X 為
   X

   其中season為因子型變量薄料，且將Autumn作為參考物，Spring晴及，Summer和Winter都與它相比都办；而bwt為連續(xù)型變量嫡锌，直接用一列向量表示

2. y 為表型值，在該例子中為一列長度是300的向量
3. 此時(shí)經(jīng)過id對(duì)應(yīng)后的M矩陣（經(jīng)過id對(duì)應(yīng)篩選后的geno：數(shù)字基因型矩陣琳钉；300 × 1000）和表型值y（長度為300）势木，就可以做后續(xù)貝葉斯分析了

接下來重點(diǎn)解析函數(shù) Bayes()，這一部分是用cpp寫的

// 首先是定義變量
# y 值得是表型數(shù)據(jù)歌懒，數(shù)據(jù)類型為一列向量
arma::vec &y,
# X 是 geno：數(shù)字基因型矩陣啦桌，數(shù)據(jù)類型為矩陣
arma::mat &X,
# 定義模型類型，數(shù)據(jù)類型為字符串及皂，該例子為 BayesCpi
std::string model,
# 所設(shè)定的零效應(yīng)SNP和有效應(yīng)的SNP的比例（零效應(yīng)的SNP比例為0.98甫男，有效應(yīng)的SNP的比例為0.02），數(shù)據(jù)類型為一列向量
arma::vec Pi,
# 以下這些是計(jì)算時(shí)所需要的數(shù)值型變量
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代碼中對(duì)一些基本參數(shù)進(jìn)行了判斷

omp_setup(threads);

# 判斷表型數(shù)據(jù) y 中是否含有 NA值
if(y.has_nan())   throw Rcpp::exception("NAs are not allowed in y.");
# 判斷表型數(shù)據(jù) y 的長度是否與geno：數(shù)字基因型矩陣的行數(shù)相同
if(y.n_elem != X.n_rows)  throw Rcpp::exception("Number of individuals not equals.");
# 判斷模型類型碍拆，并賦予數(shù)字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;
# 判斷所設(shè)定的零效應(yīng)SNP和有效應(yīng)的SNP的比例中的元素是否小于2個(gè)
if(Pi.n_elem < 2)  throw Rcpp::exception("Pi should be a vector.");
# 判斷所設(shè)定的零效應(yīng)SNP和有效應(yīng)的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.");
# 判斷所設(shè)定的零效應(yīng)SNP和有效應(yīng)的SNP的比例若治，其數(shù)值是否介于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 的容器大小端幼，為兩個(gè)長度的向量
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 的容器大小，為兩個(gè)長度的向量
    fold_.resize(2);
}
# 判斷SNPs方差解釋的向量是否與零效應(yīng)SNP和有效應(yīng)的SNP比例向量等長
if(fold_.n_elem != Pi.n_elem){
    throw Rcpp::exception("length of Pi and fold not equals.");
}

# 獲得表型向量 y 的元素個(gè)數(shù)
int n = y.n_elem;
# 獲得 geno：數(shù)字基因型矩陣的列數(shù)
int m = X.n_cols;
# 計(jì)算表型向量 y 的方差
double vary = var(y);
double h2 = 0.5;

int n_records = (niter - nburn) / thin;

# 定義參數(shù)
double* dci;
int nc = 0;
# C_ 為 X 矩陣弧满，即固定效應(yīng)的設(shè)計(jì)矩陣
mat C_;
vec beta;
vec betasd;
mat beta_store;

# 此時(shí)的C_即為X矩陣婆跑，即固定效應(yīng)的設(shè)計(jì)矩陣，以下代碼的目的是對(duì)矩陣進(jìn)行基本判斷庭呜，
## 1.判斷C_與X矩陣的行數(shù)是否相等滑进；
## 2.C_矩陣是否含有NA；
## 3. nc值為C_矩陣（即X矩陣的列數(shù)）的列數(shù)募谎，更改向量 beta的大小為 nc 個(gè)長度的向量
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);
}

# 定義一個(gè)長度為nc的向量cpc
vec cpc;
if(nc){
    cpc.resize(nc);
    #pragma omp parallel for
    # cpc 的元素為C_矩陣（即X矩陣）第 i 列的列向量的點(diǎn)乘
    for(int i = 0; i < nc; i++){
        # C_.col(i) 相當(dāng)于第 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 是環(huán)境變量的自由度，默認(rèn)為 -1
if(dfvr.isNotNull()){
    dfr = as<double>(dfvr);
}else{
    dfr = -1;
}

# 定義s2vr值近哟，s2vr 是環(huán)境變量的標(biāo)準(zhǔn)化參數(shù)，默認(rèn)為 0
if(s2vr.isNotNull()){
    s2r = as<double>(s2vr);
}else{
    s2r = 0;
}

mat vr_store;
mat estR_store;

# 定義三個(gè)空向量鲫寄，下面有用
## sp_mat 矩陣型向量吉执；vector<string> 字符型向量
vector<sp_mat> Z;
vector<sp_mat> ZZ;
vector< vector<string> > Z_levels;

# 定義R矩陣，R矩陣是隨機(jī)效應(yīng)項(xiàng)的設(shè)計(jì)矩陣
if(R.isNotNull()){
    # 定義一個(gè)字符矩陣 R_
    CharacterMatrix R_ = as<CharacterMatrix>(R);
    # 判斷 R_ 矩陣的行數(shù)（隨機(jī)效應(yīng)項(xiàng)內(nèi)容）是否等于 X 矩陣的行數(shù)（固定效應(yīng)項(xiàng)設(shè)計(jì)矩陣）地来，即樣本量是否相等
    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_ 矩陣的列數(shù)
    nr = R_.ncol();
    # vr 向量的長度為 nr戳玫，與R_ 矩陣的列數(shù)相等
    vr.resize(nr);
    # vrtmp 向量的長度為 nr，與R_ 矩陣的列數(shù)相等
    vrtmp.resize(nr);
    # R_idx 向量的長度為 nr + 1未斑，且R_idx 向量的第一個(gè)元素為 -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() 函數(shù)的意義是對(duì)容器 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_矩陣（隨機(jī)效應(yīng)項(xiàng)的內(nèi)容）的每一列向量提取出來，賦予Ri
        CharacterVector Ri = R_(_, i);
     # makeZ() 的作用是將隨機(jī)效應(yīng)項(xiàng)的內(nèi)容轉(zhuǎn)換為稀疏的設(shè)計(jì)矩陣
     # 本例中為兩列府阀，將這兩列隨機(jī)效應(yīng)項(xiàng)的內(nèi)容轉(zhuǎn)換為稀疏的設(shè)計(jì)矩陣全部賦予 Z_info
        List Z_info = makeZ(Ri);
     #  Z_info[0] 代表稀疏設(shè)計(jì)矩陣缆镣，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;
    }
        # 建立兩個(gè)零矩陣estR，estR_store為矩陣
    estR.zeros(n_levels);
    estR_store.zeros(n_levels, n_records);
}

其中

1. 函數(shù) resize(size) 作用是更改容器大小试浙，vec a; a.resize(size)相當(dāng)于將向量a的容積改為2董瞻；fill(value) 的作用是對(duì)容器中的所有位置都填充元素 value；函數(shù) push_back(value) 作用是將一個(gè)新的元素加到vector的最后面
2. 對(duì)于矩陣 X田巴，X.col(i) 代表將 X 矩陣的第 i 個(gè)列向量取出來
   X矩陣 固定效應(yīng)設(shè)計(jì)矩陣
   
   若X.col(2)钠糊，相當(dāng)于圖中的第二列列向量（紅框），例子中的X矩陣為固定效應(yīng)的設(shè)計(jì)矩陣 X壹哺，該例子中固定效應(yīng)矩陣 X 點(diǎn)乘后的結(jié)果為
3. nc 為 C_ 矩陣或 X 矩陣的總列數(shù)抄伍，而C_或X矩陣是固定效應(yīng)項(xiàng)的設(shè)計(jì)矩陣；dfvf 是環(huán)境變量的自由度管宵，默認(rèn)為 -1截珍；s2vr 是環(huán)境變量的標(biāo)準(zhǔn)化參數(shù)，默認(rèn)為 0啄糙；R矩陣是隨機(jī)效應(yīng)項(xiàng)的內(nèi)容：
   R矩陣
   本例中笛臣，R矩陣第一列表示 loc，第二列表示 dam
4. makeZ()（見補(bǔ)充函數(shù) 1） 的作用是將隨機(jī)效應(yīng)項(xiàng)的內(nèi)容（R矩陣或R_矩陣）
   隨機(jī)效應(yīng)項(xiàng)的內(nèi)容 R 矩陣
   
   轉(zhuǎn)換為稀疏的設(shè)計(jì)矩陣
   隨機(jī)效應(yīng)項(xiàng)稀疏設(shè)計(jì)矩陣
   上圖為R矩陣或R_矩陣的第一列（也就是第一個(gè)隨機(jī)效應(yīng)變量loc隧饼，該例子有兩個(gè)loc和dam）作為例子（上圖點(diǎn)代表 0）沈堡；Z_info 這個(gè) list 有兩個(gè)元素：Z_info[0] 為隨機(jī)效應(yīng)loc的稀疏設(shè)計(jì)矩陣，Z_info[1] 為因子變量的 level
5. 由于

sp_mat z = Z_info[0];
Z.push_back(z);
ZZ.push_back(z.t() * z);

Z 代表隨機(jī)效應(yīng)稀疏的設(shè)計(jì)矩陣（loc和dam）燕雁；ZZ代表每一個(gè)隨機(jī)效應(yīng)變量的稀疏設(shè)計(jì)矩陣（loc和dam）的點(diǎn)乘：z.t() * z

接下來需要對(duì)方差-協(xié)方差矩陣進(jìn)行定義

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;

# 判斷方差-協(xié)方差矩陣
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;

# 對(duì)模型進(jìn)行判斷
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++){
    # 計(jì)算 X 矩陣每一個(gè)列向量的平方和與方差，
    # X 矩陣為固定效應(yīng)的設(shè)計(jì)矩陣
    vec Xi = X.col(i);
    xpx[i] = sum(square(Xi));
    vx[i] = var(Xi);
}

# 計(jì)算 X 矩陣列向量方差和
double sumvx = sum(vx);
int nvar0 = sum(vx == 0);

# 判斷其他參數(shù)拐格，并賦值
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);

# 輸出的日志文件僧免，分別解釋每個(gè)參數(shù)的意義
 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)的參數(shù)：

1. Number of observations的參數(shù) n 代表所觀察到有表型值的個(gè)體數(shù) int n = y.n_elem;
2. Observations with genotype的參數(shù) (n - ne) 代表所觀察到的基因型數(shù)量，默認(rèn)下的 ne = 0捏浊，因此基因型數(shù)目與所觀察到有表型值的個(gè)體數(shù)相等懂衩；int ne = 0;
3. Observations with imputed genotype的參數(shù) ne 代表的是由于缺少基因型而做的impute，默認(rèn)下的 ne = 0金踪，int ne = 0;
4. Number of covariates 的參數(shù) (nc + 1) 代表的是協(xié)方差（協(xié)變量）的數(shù)量浊洞，默認(rèn)下的 nc = 0，即int nc = 0;胡岔，相當(dāng)于協(xié)方差的數(shù)量為 1
5. Number of envir-random effects 的參數(shù) nr 代表隨機(jī)效應(yīng)（或者環(huán)境效應(yīng)）的數(shù)量法希，默認(rèn)下 nr 等于隨機(jī)效應(yīng)內(nèi)容矩陣 R 矩陣的列數(shù)（R 矩陣見下）
   R矩陣
6. Number of markers 的參數(shù) m 代表的是固定效應(yīng)矩陣的設(shè)計(jì)矩陣變量個(gè)數(shù)； int m = X.n_cols;
7. π for markers in zero effect size 代表零SNP效應(yīng)的比例 Pi[0]
8. π for markers in non-zero effect size 代表非零SNP效應(yīng)的比例 Pi[1]
9. Total number of iteration 代表的是MCMC的迭代次數(shù)靶瘸，默認(rèn) niter = 50000
10. Total number of burn-in 代表的是
11. Frequency of collecting 代表的是采樣頻率苫亦，每相隔thin輪采一次樣
12. Phenotypic var 代表的是表型值的方差毛肋，double vary = var(y);
13. Environmental var 代表的是環(huán)境變量的方差：

vec vrtmp;
vrtmp.resize(nr);
vrtmp.fill(vary * (1 - h2) / (nr + 1)); ## h2默認(rèn)等于 0.5

環(huán)境變量的方差 vrtmp 為長度為 nr （見第5條）的向量，默認(rèn)的環(huán)境變量所有的值相等（vrtmp與表型值方差vary有關(guān)）屋剑，等于vary * (1 - h2) / (nr + 1)

14. Genetic var 的參數(shù) vara_ 代表遺傳方差润匙，是一個(gè)雙精型數(shù)值：

double vara_
vara_ = ((dfvara_ - 2) / dfvara_) * vary * h2; # 默認(rèn)下h2為0.5，dfvara_ 為 4

顯然遺傳方差與表型值方差vary有關(guān)

15. Inv-Chisq gpar 的參數(shù) dfvara_ 代表遺傳方差的自由度饼丘，默認(rèn)為 4趁桃，dfvara_ = 4
16. Residual var 的參數(shù) vare_ 代表殘差方差的先驗(yàn)值，vare_ = vary * (1 - h2) / (nr + 1);肄鸽，其中 vary 代表表型值的方差卫病，h2 = 0.5，nr 為隨機(jī)效應(yīng)內(nèi)容矩陣 R 矩陣的列數(shù)（見第5條）典徘，該例子中 nr = 2
17. Inv-Chisq epar 代表殘差方差的自由度蟀苛，默認(rèn)為 -2，dfvare_ = -2;
18. Marker var 的參數(shù) varg 代表固定效應(yīng)的方差逮诲，默認(rèn)double varg = vara_ / ((1 - Pi[0]) * sumvx);帜平，其中 vara_ 代表遺傳方差（參見第14條），Pi[0] 代表零SNP效應(yīng)的比例（見第7條）梅鹦，sumvx 代表固定效應(yīng)矩陣列向量的方差之和

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

因此固定效應(yīng)的方差與遺傳方差和固定效應(yīng)的設(shè)計(jì)矩陣列向量方差有關(guān)

MCMC 估計(jì)參數(shù)部分（根據(jù)示例數(shù)據(jù)的內(nèi)容裆甩，代碼做了部分刪減）

# 進(jìn)行時(shí)間的計(jì)算
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);
# 定義一個(gè)與表型值的個(gè)體數(shù)長度相等的向量，向量元素值全為 1
vec one(n); one.fill(1);
# 定義一個(gè)減去表型均值后的表型觀察值向量齐唆，即每個(gè)表型觀測者減去其均值
vec yadj = y - mu;
double* dyadj = yadj.memptr();
double* dxi;

# 這里的niter代表總迭代輪數(shù)嗤栓，iter代表迭代到第幾輪
for(int iter = 0; iter < niter; iter++){
    // sample intercept
    # 對(duì)截距 μ 進(jìn)行采樣
    ## 每迭代一輪需要對(duì)表型值的均值重新進(jìn)行正態(tài)分布采樣
    ## 函數(shù) norm_sample() 見補(bǔ)充函數(shù) 2
    ## 這一步的作用是想通過多次采樣獲得穩(wěn)定的表型值均值
    ## mu_ 滿足均值為 sum(yadj) / n，方差為 sqrt(vare_ / n) 的正態(tài)分布
    ## 代碼 sum(yadj) / n 為 yadj（減去表型均值后的表型觀察值向量）的均值
    # 以下的 mu_ 可以理解為表型向量均值的估計(jì)值
    mu_ = - norm_sample(sum(yadj) / n, sqrt(vare_ / n));
    mu -= mu_;
    # 初始的截距值為 dyadj箍邮，初始截距 dyadj = mu_ + y - mu (dyadj, 減去表型均值后的表型觀察值向量)茉帅，即 y*
    # 本質(zhì)上此時(shí)的 dyadj 代表表型向量的估計(jì)值，詳見下面第 7 條
    daxpy_(&n, &mu_, one.memptr(), &inc, dyadj, &inc);
        
    # 對(duì)固定效應(yīng)的系數(shù)進(jìn)行MCMC采樣锭弊，參見理論的第一部分和公式 7
    for(int i = 0; i < nc; i++){
        # 這里的C_矩陣（或C矩陣）為固定效應(yīng)設(shè)計(jì)矩陣 X
        # 將 C_矩陣（或C矩陣堪澎，即固定效應(yīng)設(shè)計(jì)矩陣 X）的每一列轉(zhuǎn)換為二進(jìn)制
        dci = C_.colptr(i);
        # 向量 beta 的中所有元素的初始值為 0，這里的 oldgi 代表上一輪采樣的系數(shù) β
        oldgi = beta[i];
        # 即cpc 的元素為C_矩陣（即X矩陣）第 i 列的列向量的點(diǎn)乘
        # 即cpc[i] 代表C_矩陣（即X矩陣）第 i 列向量的點(diǎn)乘
        v = cpc[i];
        # 計(jì)算固定效應(yīng)設(shè)計(jì)矩陣 X （C 矩陣）的每一列向量與初始截距（即為表型值向量）向量的點(diǎn)積
        rhs = ddot_(&n, dci, &inc, dyadj, &inc);
        rhs += v * oldgi;
        # 這里的變量 gi 代表新采樣的 β 值
        gi = norm_sample(rhs / v, sqrt(vare_ / v));
        # 計(jì)算前后兩次采樣的 β 值的誤差味滞，gi 代表新采樣的 β 值樱蛤；oldgi 代表上一次采樣的 β 值
        # gi_ 代表前后兩次采樣的 β 值的誤差
        gi_ = oldgi - gi;
       # 計(jì)算最優(yōu)化的目標(biāo)函數(shù)，見理論的第一部分和公式 8
       # 和下面補(bǔ)充說明第 8 條
        daxpy_(&n, &gi_, dci, &inc, dyadj, &inc);
        # 新獲得的 gi 值賦予給向量 beta 的對(duì)應(yīng)元素
        beta[i] = gi;
    }
        
    # 對(duì)隨機(jī)效應(yīng)的系數(shù)進(jìn)行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;
        # 計(jì)算公式6中刹悴，隨機(jī)效應(yīng)正態(tài)分布均值的分子部分，Z 和 ZZ 詳細(xì)見下面第 9 條
        r_RHS = Z[i].t() * yadj;
        # 其中這里的 estR_tmp 相當(dāng)于前一輪采樣的隨機(jī)效應(yīng)系數(shù) rj 的估計(jì)值
        r_RHS += ZZ[i] * estR_tmp;
        # 計(jì)算公式6中攒暇，隨機(jī)效應(yīng)正態(tài)分布均值的分母部分
        # vare_ / vrtmp[i] 可以解釋為遺傳力
        r_LHS = diagvec(ZZ[i]) + vare_ / vrtmp[i];
        
        # 進(jìn)行隨機(jī)效應(yīng)系數(shù)的采樣，estR_tmp 代表待采樣的隨機(jī)效應(yīng)系數(shù) rj 的估計(jì)值
        # 參數(shù) qr 代表隨機(jī)效應(yīng)內(nèi)容 R 矩陣的列數(shù)（即一共有 qr 個(gè)隨機(jī)效應(yīng)項(xiàng)）
        # estR_tmp[qi] 代表第 qi 個(gè)隨機(jī)效應(yīng)的系數(shù) ( r )
        for(int qi = 0; qi < qr; qi++){
            # r_RHS[qi] / r_LHS[qi] 為正態(tài)分布的方差
            estR_tmp[qi] = norm_sample(r_RHS[qi] / r_LHS[qi], sqrt(vare_ / r_LHS[qi]));
        }

        # 計(jì)算最優(yōu)化的目標(biāo)函數(shù)子房，見理論的第一部分和公式 8
        vec diff = Z[i] * (estR.subvec(R_idx[i] + 1, R_idx[i + 1]) - estR_tmp);
        # 計(jì)算最優(yōu)化的目標(biāo)函數(shù)形用，詳細(xì)說明見下面第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);
        # 計(jì)算隨機(jī)效應(yīng)系數(shù) r 的方差
        vr[i] = var(estR_tmp);

        estR.subvec(R_idx[i] + 1, R_idx[i + 1]) = estR_tmp;
    }
    
   # 對(duì)Genetic marker effects的系數(shù) α 進(jìn)行采樣就轧，即基因型矩陣 M
   # 同時(shí)計(jì)算遺傳方差 Vg
   # 這里先考慮 BayesCpi 方法的參數(shù)估計(jì)
    switch(model_index){
        case 4:
            logpi = log(Pi);
            # 這里的初始 s 為：vec s = zeros(n_fold)；默認(rèn)的 n_fold = 2
            # 其中 logpi[0] 代表 logpi 的第一個(gè)元素零效應(yīng)SNP的比例
            s[0] = logpi[0];
            // sumvg = 0;
            vargi = 0;
            for(int i = 0; i < m; i++){
                if(!vx[i])   continue;
                # 這一塊的變量 X 為基因型矩陣 M田度，只不過變量定義為X
                # 切勿與固定效應(yīng)的設(shè)計(jì)矩陣 X 混淆
                # 將基因型矩陣 M 的第 i 列取出
                dxi = X.colptr(i);
                # 變量 xpx 為基因型矩陣 M 的每一個(gè)列向量的平方和所構(gòu)成的集合
                xx = xpx[i];
                # 這里的 g[i] 初始值為 vec g = zeros(m); 長度為 m (m 為 M 矩陣的列長妒御，元素全為 0)
                oldgi = g[i];
                # 這一塊詳細(xì)見下面的第 12 條
                rhs = ddot_(&n, dxi, &inc, dyadj, &inc);
                # 這里的 oldgi 相當(dāng)于上一輪采樣的 α，初始的 α = 0
                if(oldgi)   rhs += xx * oldgi; // 首次采樣時(shí)镇饺，oldgi = 0乎莉，拒絕采樣時(shí) oldgi = 0
                lhs = xx / vare_;
                logdetV = log(varg * lhs + 1);
                uhat = rhs / (xx + vare_ / varg);
                # 這一塊計(jì)算的是公式 11 的內(nèi)容，目的是構(gòu)造接收概率 acceptProb 
                s[1] = -0.5 * (logdetV - (rhs * uhat / vare_)) + logpi[1];
                # 拒絕-接收采樣奸笤，定義接收率
                acceptProb = 1 / sum(exp(s - s[0]));
                # 參考下面第 12 條
                # indistflag 作用是決定本次采樣是否被接收
                # 從均勻分布隨機(jī)取一個(gè)值與 acceptProb 比大小惋啃，acceptProb 的計(jì)算參考公式 12
                # 若從均勻分布中隨機(jī)取的一個(gè)值小于 acceptProb，indistflag 值被則賦予 0监右，否則被賦予 1
                indistflag = (uniform_sample()) < acceptProb ? 0 : 1;
                snptracker[i] = indistflag;
                # 接收采樣
                if(indistflag){
                    # 這里對(duì) genetic marker effect 的系數(shù)進(jìn)行采樣
                    v = xx + vare_ / varg;
                    # 采樣分布參考公式 9 采樣分布部分边灭，詳細(xì)介紹參考下面第 13 條
                    # 這里的 gi 相當(dāng)于本輪采樣的 α
                    gi = norm_sample(rhs / v, sqrt(vare_ / v));
                    # 定義最優(yōu)化函數(shù)參考公式 9 的最優(yōu)化函數(shù)部分
                    gi_ = oldgi - gi;
                    daxpy_(&n, &gi_, dxi, &inc, dyadj, &inc);
                    // sumvg += (gi * gi) * vx[i];
                    vargi += (gi * gi);
                # 拒絕采樣
                }else{
                # 拒絕采樣的最優(yōu)化函數(shù)
                    gi = 0;
                    if(oldgi){
                        gi_ = oldgi;
                        daxpy_(&n, &gi_, dxi, &inc, dyadj, &inc);
                    }
                }
                g[i] = gi;
            }
            # 詳細(xì)解釋查看下面第 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 的方法執(zhí)行下面語句，!fixpi = TRUE健盒，n_fold = 2
            # fold_snp_num 代表當(dāng)前輪次下拒絕-接收采樣中绒瘦，接收采樣的次數(shù)
            # 設(shè) Πk 滿足 dirichlet 分布
            if(!fixpi)  Pi = rdirichlet_sample(n_fold, (fold_snp_num + 1));
            break;
    }
        
    // collect the parameters to obtain posterior estimation
    # 獲得采樣的數(shù)據(jù)，并不是每一個(gè) iter 都采樣扣癣，而是隔 nburn 采樣一次
    if(iter >= nburn && (iter + 1 - nburn) % thin == 0){
        # 儲(chǔ)存采樣的均值（截距） μ惰帽，count 為記錄采樣的次數(shù)
        mu_store[count] = mu;
        if(!fixpi)  pi_store.col(count) = Pi;
        # 儲(chǔ)存采樣的基因型系數(shù)矩陣 α 的方差
        vara_store[count] = vara_;
        
        # 儲(chǔ)存采樣的基因型系數(shù)矩陣 α 的值
        g_store.col(count) = g;

        # 儲(chǔ)存采樣的固定效應(yīng)系數(shù)矩陣 β 的值
        if(nc){
            beta_store.col(count) = beta;
        }
        
        # 儲(chǔ)存采樣的隨機(jī)效應(yīng)系數(shù)矩陣 r 的值
        if(nr){
        # 儲(chǔ)存采樣的隨機(jī)效應(yīng)系數(shù)矩陣 r 值的方差
            vr_store.col(count) = vr;
        # 儲(chǔ)存采樣的隨機(jī)效應(yīng)系數(shù)矩陣 r 的值
            estR_store.col(count) = estR;
        }
        count++;
    }
}

上述代碼中出現(xiàn)的 oldgi 代表上一輪采樣的參數(shù)（β，α等父虑；固定效應(yīng)该酗，genetic marker effect 的系數(shù)）；oldgi 代表本輪采樣的參數(shù)（β频轿，α等垂涯；固定效應(yīng)，genetic marker effect 的系數(shù)）航邢； estR_tmp[qi] 代表第 qi 個(gè)隨機(jī)效應(yīng)的系數(shù) ( r )

其中：

1. 上述代碼中出現(xiàn)的 C_或C 代表固定效應(yīng)設(shè)計(jì)矩陣 X
2. BLAS 庫是cpp中一個(gè)進(jìn)行線性代數(shù)計(jì)算的庫耕赘，用這個(gè)庫進(jìn)行計(jì)算時(shí)首先要將矩陣或者向量轉(zhuǎn)換為二進(jìn)制：

# 將向量轉(zhuǎn)換為二進(jìn)制
double* dyadj = yadj.memptr();
dci = C_.colptr(i);
# 計(jì)算點(diǎn)積
rhs = ddot_(&n, dci, &inc, dyadj, &inc);

而函數(shù)ddot_(N, X, INCX, Y, INCY)的作用是用于計(jì)算兩個(gè)向量的點(diǎn)積，參數(shù) N 代表向量的長度膳殷；參數(shù) X 代表待計(jì)算的向量 A操骡；參數(shù) INCX 代表待計(jì)算的向量 A 元素之間的間距（默認(rèn)為1）；參數(shù) Y 代表待計(jì)算的向量 B赚窃；參數(shù) INCY 代表待計(jì)算的向量 B 元素之間的間距（默認(rèn)為1）册招；

3. 函數(shù) 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 代表系數(shù) a，const double * x 代表待輸入的向量 x辱匿，const MKL_INT incx 代表參數(shù) INCX 代表待計(jì)算的向量 x 元素之間的間距（默認(rèn)為1）键痛，double * y 代表待計(jì)算的向量 y炫彩，const MKL_INT incy 代表參數(shù) INCY 代表待計(jì)算的向量 y 元素之間的間距（默認(rèn)為1）；因此函數(shù)所作運(yùn)算為：ynew = a*x + y絮短，然后將新的 y（ynew）賦予 y
4. 在對(duì)計(jì)算隨機(jī)效應(yīng)系數(shù)的mcmc采樣過程中

r_RHS = Z[i].t() * yadj;
r_RHS += ZZ[i] * estR_tmp;

相當(dāng)于計(jì)算公式6中正態(tài)分布均值的分子部分江兢，即

5. 在對(duì)計(jì)算隨機(jī)效應(yīng)系數(shù)的mcmc采樣過程中

r_LHS = diagvec(ZZ[i]) + vare_ / vrtmp[i];

相當(dāng)于計(jì)算公式6中正態(tài)分布均值的分母部分，

6. 在對(duì)計(jì)算隨機(jī)效應(yīng)系數(shù)的mcmc采樣過程中丁频，vare_ / r_LHS[qi] 相當(dāng)于
7. 為什么進(jìn)行以下操作：

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

筆者認(rèn)為經(jīng)過上述代碼運(yùn)行后杉允，初始截距 dyadj = mu_ + y - mu（而 mu__new = mu - mu__采樣后，則mu__new表示表型向量均值的真實(shí)值與估計(jì)值之間的誤差）的目的是創(chuàng)造一個(gè)表型向量的估計(jì)值席里，即 y*

8. 對(duì)于固定效應(yīng)系數(shù)的最優(yōu)化目標(biāo)函數(shù)：daxpy_(&n, &gi_, dci, &inc, dyadj, &inc);的目的是計(jì)算最優(yōu)化的目標(biāo)函數(shù)叔磷，見理論的第一部分和公式 8，相當(dāng)于執(zhí)行了公式 8胁勺，即dyadj = dyadj + dci(oldgi - gi)世澜，其中 oldgi - gi 代表 β_j^i-1 - β_jⁱ，dci 代表 X_j署穗，dyadj 代表 y*
9. Z 和 ZZ 分別代表什么寥裂？Z 代表隨機(jī)效應(yīng)稀疏的設(shè)計(jì)矩陣的list（loc和dam），例如 Z_loc為
   是一個(gè) 300 X 50 的稀疏矩陣（300 為隨機(jī)效應(yīng)的R矩陣的行數(shù)案疲，50為隨機(jī)效應(yīng)項(xiàng) loc 的因子水平為 50）封恰，點(diǎn)代表 0 ，loc 的因子水平為 50：
   loc 的因子水平
   
   而 Z_dam 為
   是一個(gè) 300 X 150 的稀疏矩陣（300 為隨機(jī)效應(yīng)的R矩陣的行數(shù)褐啡，150為隨機(jī)效應(yīng)項(xiàng) dam 的因子水平為 150）诺舔，點(diǎn)代表 0 ，dam 的因子水平為 150：
   dam 的因子水平
   而 ZZ_loc = t(Z_loc) %*% Z_loc备畦；ZZ_dam = t(Z_dam) %*% Z_dam
10. 對(duì)于隨機(jī)效應(yīng)系數(shù)的最優(yōu)化目標(biāo)函數(shù)：

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*

11. 基因型矩陣M（為了區(qū)別于固定效應(yīng)的設(shè)計(jì)矩陣 X低飒，基因型矩陣在全文中僅稱為 M 矩陣），M 矩陣懂盐，即geno：數(shù)字基因型矩陣（n × m, n為個(gè)體數(shù)褥赊，m為snp數(shù)，本例中為 600 × 1000）莉恼，接收 0, 1, 2 或者 -1, 0, 1 這兩種形式的矩陣拌喉，但矩陣中不能含有NA，
    M矩陣
12. 有關(guān)于公式11的代碼部分解釋：

rhs = ddot_(&n, dxi, &inc, dyadj, &inc);
if(oldgi)   rhs += xx * oldgi; # xx 為 M 矩陣其中一列的平方和

上述代碼計(jì)算的 rhs

這里的 oldgi = g[i] 表示：g[i] 初始值為 vec g = zeros(m); 長度為 m（m 為 M 矩陣的列長俐银，元素全為 0）尿背；uhat = rhs / (xx + vare_ / varg);代表的是

logdetV = log(varg * lhs + 1);計(jì)算的是

而 s[1] = -0.5 * (logdetV - (rhs * uhat / vare_)) + logpi[1]; 計(jì)算的是公式11

公式11

而所有的這些計(jì)算目的是構(gòu)造接收函數(shù)（或者是接收值，因?yàn)閷?duì)參數(shù) α 的估計(jì)采用的是拒絕-接收采樣）

公式12

這個(gè)概率值就為接收函數(shù)或者接收值捶惜；一般情況下田藐，n_Π = 2，第一個(gè)元素表征具有零效應(yīng)的SNP比例，第二個(gè)元素表征具有非零效應(yīng)的SNP比例汽久。每次采樣的時(shí)候從均勻分布中隨機(jī)取一個(gè)數(shù)字 num 與上圖的 acceptProb 作比較茴晋，若 num < acceptProb 接收采樣 ，否則拒絕采樣回窘；該例子由于 Pi = c(0.98, 0.02)，因此 nk = 2;Π0 = 0.98;Π1 = 0.02市袖，Π0 = 0.98 代表具有零效應(yīng)的SNP比例啡直；Π1 = 0.02 代表具有非零效應(yīng)的SNP比例

13. 對(duì)系數(shù) α 的采樣介紹：gi = norm_sample(rhs / v, sqrt(vare_ / v));計(jì)算的采樣分布為
    接下來計(jì)算最優(yōu)化的函數(shù)

gi_ = oldgi - gi;
daxpy_(&n, &gi_, dxi, &inc, dyadj, &inc);

上述代碼則表示最優(yōu)化目標(biāo)的函數(shù)：

14. 其中 Π_k（ k > 0，Π₀ 在該例子中為 0.98）饮睬，n_fold = 2沃缘，fold_snp_num 代表當(dāng)前輪次下拒絕-接收采樣中奶段，接收采樣的次數(shù)（當(dāng)前輪次，比方說 iter = 3舷丹，接收采樣的次數(shù)為 2，則 fold_snp_num = 2蜓肆；iter = 10颜凯，接收采樣的次數(shù)為 6，則 fold_snp_num = 6）仗扬，且 Πk 滿足 dirichlet 分布
    
    η 即為fold_snp_num症概，代表當(dāng)前輪次下拒絕-接收采樣中，接收采樣的次數(shù)（當(dāng)前輪次早芭，比方說 iter = 3彼城，接收采樣的次數(shù)為 2，則 η = 2退个；iter = 10募壕，接收采樣的次數(shù)為 6，則 η = 6）

收集數(shù)據(jù)语盈，返回結(jié)果（根據(jù)示例數(shù)據(jù)做了一定刪減）

# 定義結(jié)果的 list
List results;
List MCMCsample;

# 儲(chǔ)存隨機(jī)效應(yīng)系數(shù)矩陣的方差
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;
}

# 儲(chǔ)存基因型系數(shù)矩陣的方差舱馅，取多次迭代的均值
vara_ = mean(vara_store);
double varasd = stddev(vara_store);
results["Vg"] = vara_;
MCMCsample["Vg"] = vara_store.t();

# 儲(chǔ)存均值（截距）的方差，取多次迭代的均值
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){
    # 儲(chǔ)存隨機(jī)效應(yīng)系數(shù)矩陣的方差黎烈，取多次迭代的均值
    beta = conv_to<vec>::from(mean(beta_store, 1));
    betasd = conv_to<vec>::from(stddev(beta_store, 0, 1));

        # 儲(chǔ)存隨機(jī)效應(yīng)系數(shù)矩陣 β 值
    results["beta"] = beta;
    MCMCsample["beta"] = beta_store;
}

# 儲(chǔ)存基因型系數(shù)矩陣 α 值习柠，取多次迭代的均值
g = conv_to<vec>::from(mean(g_store, 1));
results["alpha"] = g;
MCMCsample["alpha"] = g_store;

# # 儲(chǔ)存非零效應(yīng)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;

# 儲(chǔ)存隨機(jī)效應(yīng)系數(shù)矩陣 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;
}

補(bǔ)充函數(shù)

1. makeZ()：makeZ() 的作用是將隨機(jī)效應(yīng)項(xiàng)的內(nèi)容（R矩陣或R_矩陣）资溃，轉(zhuǎn)換為稀疏的設(shè)計(jì)矩陣

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);
}

2. norm_sample() 的作用是進(jìn)行正態(tài)分布的采樣：

double norm_sample(double mean, double sd){
    // return R::rnorm(mean, sd);
    return mean + sd * norm_rand();
}

補(bǔ)充的代數(shù)知識(shí)：
如果將每一個(gè) βj，rj 或 αj 看作是隨機(jī)變量（假設(shè) j 代表第 j 列列向量）：