当有多个年份和多个地点的表型值之后，基因型只有一套，这时如果用来做GWAS的时候，就需要对表型值进行重新计算。
多年多点的表型值用于GWAS分析前，一般有以下三种方式供预处理表型：
- 求多年多点表型的均值
- 使用blup值
- 使用blue值
数据符合正太分布时使用的是lmer函数，符合偏态分布时使用glmer函数。
BLUP值和BLUE值的区别：
估计随机效应（BLUP）：用于估计随机效应的最佳线性无偏预测值（Best Linear Unbiased Predictor）
估计固定效应（BLUE）：用于估计固定效应的最佳线性无偏估计值（Best Linear Unbiased Estimator）
BLUP值预测品种将来的表现
BLUE值预测品种现在的表现。
BLUP值会对数据进行压缩。
BLUE值会对年份、地点进行矫正，得到的结果和原数据尺度一样。
多年多点无重复的模型
lmer构建的是线性回归模型，我们需要提前知道固定因子和随机因子。
lmer(FL~Taxa+(1|location)+(1|year)+(1|location:year),data=All_fiber_noNA)
上面的公式中，
- FL是指定性状的列名，是要预测的性状
- Taxa是品种的列名，是固定因子
- location是位置的列名，是随机因子
- year是年份的列名，是随机因子
- All_fiber_noNA是数据框的名称
主的效应是品种、年份、地点、还需要考虑的是品种x年份，品种x地点，年份x地点，品种x年份x地点
多年多点有重复的模型
lmer(FL~Taxa+(1|location)+(1|year)+(1|Rep)+(1|location:year)+(1|Rep:year)+(1|location:Rep)+(1|location:Rep:year)+(1|Taxa:location)+(1|Taxa:year)+(1|Taxa:Rep)+(1|Taxa:year:Rep)+(1|Taxa:year:location)+(1|Taxa:location:Rep)+(1|Taxa:year:location:Rep),data=All_fiber_noNA)
上面的公式中，
- FL是指定性状的列名，是要预测的性状
- Taxa是品种的列名，是固定因子
- location是位置的列名，是随机因子
- year是年份的列名，是随机因子
- Rep是重复，是随机因子
- All_fiber_noNA是数据框的名称
使用lme4计算blup值和blue值,这里是多年多点无重复实验
#lme4提供三个函数用于建立模型：线性 ( lmer)、广义线性 ( glmer) 和非线性 ( nlmer.)
mod1 <- lmer(FL~(1|Taxa)+(1|location)+(1|year)+ (1|Taxa:location) + (1|Taxa:year),data=dataframe)
mod2 <- lmer(FL~(1|Taxa)+(1|location)+(1|year)+ (1|Taxa:location) + (1|Taxa:year)+(1|location:year),data=dataframe)
mod3 <- lmer(FL~Taxa+(1|location)+(1|year)+ (1|year:location),data=dataframe)
mod4 <- lmer(FL~Taxa+(1|location)+(1|year)+ (1|year:location)+(1|Taxa:location),data=dataframe)
mod5 <- lmer(FL~Taxa+(1|location)+(1|year)+ (1|year:location)+(1|Taxa:location)+(1|Taxa:year),data=dataframe)
mod6 <- lmer(FL~Taxa+(1|location)+(1|year)+ (1|year:location)+(1|Taxa:location)+(1|Taxa:year)+(1:Taxa:year:location),data=dataframe)
summary(mod1)
summary(mod2)
#计算模型的方差
anova(mod1,mod2,mod3,mod4,mod5,mod6)
# refitting model(s) with ML (instead of REML)
# Data: dataframe
# Models:
# mod1: FL ~ (1 | Taxa) + (1 | location) + (1 | year) + (1 | Taxa:location) + (1 | Taxa:year)
# mod2: FL ~ (1 | Taxa) + (1 | location) + (1 | year) + (1 | Taxa:location) + (1 | Taxa:year) + (1 | location:year)
# mod3: FL ~ Taxa + (1 | location) + (1 | year) + (1 | year:location)
# mod4: FL ~ Taxa + (1 | location) + (1 | year) + (1 | year:location) + (1 | Taxa:location)
# mod5: FL ~ Taxa + (1 | location) + (1 | year) + (1 | year:location) + (1 | Taxa:location) + (1 | Taxa:year)
# mod6: FL ~ Taxa + (1 | location) + (1 | year) + (1 | year:location) + (1 | Taxa:location) + (1 | Taxa:year) + (1:Taxa:year:location)
# npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
# mod1 7 17576 17623 -8781.1 17562
# mod2 8 16248 16302 -8116.1 16232 1329.9 1 <2e-16 ***
# mod3 423 15090 17914 -7122.0 14244 1988.2 415 <2e-16 ***
# mod4 424 15092 17923 -7122.0 14244 0.0 1 1
# mod5 425 15094 17932 -7122.0 14244 0.0 1 1
# mod6 425 15094 17932 -7122.0 14244 0.0 0
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#显著的是mod2和mod3,所以选择这2个里AIC和BIC值最小的一个，
mod1和2是把Taxa当作随机因子来计算，
mod3,4,5,6是把Taxa当作固定因子来计算。
看到最终是mod2和mod3是达到了极显著水平。所以使用这2个模型进行计算blup和blue值。
评价模型优劣的指标
BIC（Bayesian Information Criterion）和AIC（Akaike Information Criterion）是模型选择中常用的两个信息准则。它们用来衡量统计模型的拟合优度和复杂度，并在比较不同模型时提供一种相对优势判断的标准。
在具体应用中，BIC和AIC的值越小越好。
计算blup值和blue值
使用mod2可以计算blup值
FL_blup <- ranef(mod2)$Taxa
colnames(FL_blup) <- "BLUP"
FL_blup <- FL_blup%>%rownames_to_column("Taxa") #这个是blup值
使用mod3计算blue值
##blue值的计算
FL_blue <- as.data.frame(fixef(mod3))
head(FL_blue)
#install.packages("lsmeans")
#植物中，需要使用lsmeans添加截距
library('lsmeans')
re <- lsmeans(mod3,"Taxa")
Blue <- data.frame(re)[,1:2] #这个是Blue值
colnames(Blue) <- c("Taxa","BLUE")
Blup_Blue <- inner_join(FL_blup,Blue,by="Taxa")
使用mod2无法计算blue值，mod3无法计算blup值。
mod2把所有的变量都当作随机因子，就没有固定因子，所以没法算固定效应值（BLUE）。
因为mod3里我把Taxa（品种）当作固定因子，所以无法算品种的随机效应值（BLUP）。
可视化一下最终的数据的分布
hist(Blup_Blue$BLUP,col = "pink")
hist(Blup_Blue$BLUE,col = "orange")
可以看出虽然mod2和mod3的模型不一样，算出的BLUP和BLUE值的分布还是一样的。
参考地址1 https://www.jianshu.com/p/f4f1b5b75830
参考地址2 https://zhuanlan.zhihu.com/p/93495710
关于lmer函数加速的问题
参考官方给的方法 https://cran.r-project.org/web/packages/lme4/vignettes/lmerperf.html
当数据比较多的时候，lmer函数会非常慢，这时就需要想办法加速了。
m0 <- lmer(y ~ service * dept + (1|s) + (1|d), InstEval,
control = lmerControl(calc.derivs = FALSE))
通过设置control = lmerControl(calc.derivs = FALSE)来加速拟合。
如果速度还是很慢，则可以试试glmmTMB 这个包。
install.packages("glmmTMB")
library(glmmTMB)
model3 <- glmmTMB(FL~Taxa+(1|location)+(1|year)+ (1|year:location),data=dataframe)
glmmTMB传参基本和lme4一致。
