（飞哥专属二维码，报名成功后私信我领取大礼包）
更多知识学习，欢迎关注我们2024年11月23-24日的ASReml全基因组选择培训
培训内容及安排
培训费用
培训福利涵盖：
2天GS实战培训+1年的视频课程回放+ASReml软件授权+售后技术支持
培训优惠及报名
优惠详情：
1、10月31日前可享受早鸟价，优惠200元；
2、转发朋友圈集赞满36个，优惠200元；
3、3人团报免一人费用；
4、企业用户及多年购买用户优惠力度更大，详情可咨询工作人员；
扫码报名：
微信扫一扫或者长按二维码提交报名信息。
联系我们
Nancy邓老师（15611308826）
Meggi孙老师（18201280067）
电话：010-62680244
邮箱：China@vsni.co.uk
网址（CHN）：www.vsnc.com.cn
下面介绍一下asreml如何进行固定因子的显著性分析：
混合线性模型LSmeans如何进行多重比较
asremlPlus包应用介绍
LSmeans是最小二乘均值，在混合线性模型的框架下，它是固定因子的预测值，我们下面看一下，如何使用ASReml软件，对固定因子计算LSmeans以及多重比较。
01
载入软件包和数据
今天我们介绍将要用到的是asreml和asremlPlus，asreml包可以通过下方链接，申请试用或者购买（软件试用链接：http://www.vsnc.com.cn/request-trial-version-cn/）或者私聊客服小编获取，asremlPlus可以在CRAN中安装：install.packages("asremlPlus")
library(asreml)library(asremlPlus)data(oats)str(oats)head(oats)
数据介绍，oats是燕麦裂区试验设计的产量数据。
Yield of oats with four fertilizer treatments. The yield of oats from a split-plot field trial with 3 varieties and 4 levels of nitrogen fertilizer. The experiment was a split-plot with 6 blocks of 3 main plots (varieties), each split into 4 sub-plots (nitrogen).
02
构建模型
这里，我们简化一下模型，将Nitrogen和Blocks作为固定因子：
mod = asreml(yield ~ Nitrogen, random = ~ Blocks, residual = ~ idv(units), data = oats)summary(mod)$varcompwald(mod)
运行结果：
> mod = asreml(yield ~ Nitrogen, random = ~ Blocks, residual = ~ idv(units), data = oats)ASReml Version 4.2 25/10/2024 11:00:46LogLik Sigma2 DF wall1 -371.0342 1.0 68 11:00:462 -328.8292 1.0 68 11:00:463 -283.1032 1.0 68 11:00:464 -252.7651 1.0 68 11:00:465 -238.0402 1.0 68 11:00:466 -234.8200 1.0 68 11:00:467 -234.5410 1.0 68 11:00:468 -234.5376 1.0 68 11:00:46> summary(mod)$varcompcomponent std.error z.ratio bound %chBlocks 243.3047 167.34898 1.453876 P 0units!units 255.3995 45.49174 5.614194 P 0units!R 1.0000 NA NA F 0> wald(mod)Wald tests for fixed effects.Response: yieldTerms added sequentially; adjusted for those above.Df Sum of Sq Wald statistic Pr(Chisq)(Intercept) 1 245.174 245.174 < 2.2e-16 ***Nitrogen 3 78.405 78.405 < 2.2e-16 ***residual (MS) 1.000---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
03
asreml计算LSmeans
asreml计算lsmeans，需要用predict函数，方法如下：
lsmeans = predict(mod,classify = "Nitrogen")$pvalslsmeans
结果：
> lsmeansNotes:- The predictions are obtained by averaging across the hypertablecalculated from model terms constructed solely from factors inthe averaging and classify sets.- Use 'average' to move ignored factors into the averaging set.- The ignored set: BlocksNitrogen predicted.value std.error status1 0_cwt 79.38889 7.398624 Estimable2 0.2_cwt 98.88889 7.398624 Estimable3 0.4_cwt 114.22222 7.398624 Estimable4 0.6_cwt 123.38889 7.398624 Estimable
结果可以看出，Nitrogen的四个水平的最小二乘估计分别是：79.38, 98.88， 114.22， 123.38，
怎么判断他们之间有没有达到显著性水平呢？我们这里使用LSD的方法。
首先先计算SED，代码可以设置为：sed=T：
predict(mod,classify = "Nitrogen",sed=T)结果：
$pvalsNotes:- The predictions are obtained by averaging across the hypertablecalculated from model terms constructed solely from factors inthe averaging and classify sets.- Use 'average' to move ignored factors into the averaging set.- The ignored set: BlocksNitrogen predicted.value std.error status1 0_cwt 79.38889 7.398624 Estimable2 0.2_cwt 98.88889 7.398624 Estimable3 0.4_cwt 114.22222 7.398624 Estimable4 0.6_cwt 123.38889 7.398624 Estimable$avsedmin mean max5.327074 5.327074 5.327074$sed4 x 4 Matrix of class "dspMatrix"[,1] [,2] [,3] [,4][1,] NA 5.327074 5.327074 5.327074[2,] 5.327074 NA 5.327074 5.327074[3,] 5.327074 5.327074 NA 5.327074[4,] 5.327074 5.327074 5.327074 NA
结果可以看到，SED最小值为5.32， 平均SED为5.32， 最大SED为5.32，这是因为数据是平衡的，所以两两之间的SED都是5.32。
然后，我们根据刚才wald检验计算出自由度的大小，注意，这里denDF参数：
> wald(mod,denDF="default")ASReml Version 4.2 25/10/2024 11:10:51LogLik Sigma2 DF wall1 -234.5376 1.0 68 11:10:512 -234.5376 1.0 68 11:10:513 -234.5376 1.0 68 11:10:51$WaldWald tests for fixed effects.Response: yieldDf denDF F.inc Pr(Intercept) 1 5 245.10 1.93182e-05Nitrogen 3 63 26.13 0.00000e+00$stratumVariancesdf Variance Blocks units!unitsBlocks 5 3175.0556 12 1units!units 63 255.3995 0 1
结果可以看出，Nitrogen的denDF为63，
所以，我们可以计算tvalue值：
tvalue = qt((1- 0.05/2), 63,lower.tail = T)tvalue
结果为：1.998
> tvalue = qt((1- 0.05/2), 63,lower.tail = T)> tvalue[1] 1.998341
于是，我们就可以计算两两之间的LSD值为：
> sed = predict(mod,classify = "Nitrogen",sed=T)ASReml Version 4.2 25/10/2024 11:13:54LogLik Sigma2 DF wall1 -234.5376 1.0 68 11:13:542 -234.5376 1.0 68 11:13:543 -234.5376 1.0 68 11:13:54> sed$sed * tvalue4 x 4 Matrix of class "dspMatrix"[,1] [,2] [,3] [,4][1,] NA 10.64531 10.64531 10.64531[2,] 10.64531 NA 10.64531 10.64531[3,] 10.64531 10.64531 NA 10.64531[4,] 10.64531 10.64531 10.64531 NA
都是10.645，所以我们就可以进行多重比较了：如果两者之差大于10.645，就是0.05水平的显著。如果我们想计算0.01水平的显著，tvalue计算时选择0.01即可。
Nitrogen predicted.value Lable1 0_cwt 79.38889 a2 0.2_cwt 98.88889 b3 0.4_cwt 114.22222 c4 0.6_cwt 123.38889 c
04
使用asremlPlus计算LSMEANS以及多重比较
代码如下：
# 用asremlPlus计算LSD值和两两之间的P值current.asrt <- as.asrtests(mod)Var.diffs <- predictPlus(mod, classify="Nitrogen",wald.tab = current.asrt$wald.tab,tables = "none")
Var.diffs
结果如下：
> # 用asremlPlus计算LSD值和两两之间的P值> current.asrt <- as.asrtests(mod)> Var.diffs <- predictPlus(mod, classify="Nitrogen",+ wald.tab = current.asrt$wald.tab,+ tables = "none")> Var.diffs#### Predictions for yield from NitrogenNotes:- The predictions are obtained by averaging across the hypertablecalculated from model terms constructed solely from factors inthe averaging and classify sets.- Use 'average' to move ignored factors into the averaging set.- The ignored set: BlocksNitrogen predicted.value standard.error upper.Confidence.limit lower.Confidence.limit est.status1 0_cwt 79.38889 7.398624 94.17386 64.60392 Estimable2 0.2_cwt 98.88889 7.398624 113.67386 84.10392 Estimable3 0.4_cwt 114.22222 7.398624 129.00719 99.43725 Estimable4 0.6_cwt 123.38889 7.398624 138.17386 108.60392 EstimableLSD valuesminimum LSD = 10.64531mean LSD = 10.64531maximum LSD = 10.64531(sed range / mean sed = 0 )Variance matrix of the predicted valuesNULLAll pairwise differences between predicted values0_cwt 0.2_cwt 0.4_cwt 0.6_cwt0_cwt 0.00 -19.50 -34.833 -44.0000.2_cwt 19.50 0.00 -15.333 -24.5000.4_cwt 34.83 15.33 0.000 -9.1670.6_cwt 44.00 24.50 9.167 0.000p values for all pairwise differences between predicted values0_cwt 0.2_cwt 0.4_cwt 0.6_cwt0_cwt 0.001 0.000 0.000.2_cwt 0.001 0.005 0.000.4_cwt 0.000 0.005 0.090.6_cwt 0.000 0.000 0.090Standard errors of differences between predicted values0_cwt 0.2_cwt 0.4_cwt 0.6_cwt0_cwt 5.327 5.327 5.3270.2_cwt 5.327 5.327 5.3270.4_cwt 5.327 5.327 5.3270.6_cwt 5.327 5.327 5.327
可以看出来，两两之间的P值结果直接可以计算出来：
0和0.2，以及0.4之间都达到显著水平。0.4和0.6之间没有达到显著性水平。
05
结论
asreml软件，可以很容易计算LSmeans，但是计算多重比较时，需要根据wald检验计算出denDF自由度，根据sed=T参数计算两两之间的SED，最后再计算两两之间的LSD，然后再进行显著性标识。虽然可以计算，但是稍微有点繁琐。
asremlPlus软件，通过predictPlus函数，可以一键输出LSmeans以及两两之间的显著性，非常方便，而且两者结果完全一致。
所以，如果水平数比较多时，推荐使用asremlPlus进行LSmeans的计算，以及LSmeans之间的多重比较。
