阈性状 (threshold trait ):性状数值达到某一特定值时表现为正常,达不到则为不正常,如血压,血糖含量、生物的抗病力等,在数据方面以0、1表示,属于二元数据分布。ASReml也可以轻松应付阈性状,通过trait.mod参数选择binomial函数来分析。
注:trait.mod参数将来会替代为family参数。
简单示例:
res<-echidna(es0.file="dfm2.es0",
trait.mod=esr_binomial(link='marginal'),
fixed=lt~1+Rep,random=~nrm(TreeID),
residual=~units,
predict=c('Rep')
)
运行过程如下:
> res<-echidna(es0.file="dfm2.es0",
+ trait.mod=esr_binomial(link='marginal'),
+ fixed=lt~1+Rep,random=~nrm(TreeID),
+ residual=~units,
+ trace=T,
+ predict=c('Rep')
+ )
Running Echidna for analysis: lt
Tue Mar 16 16:52:15 2021
1 LogL= -673.01 0.9767 554 DF
2 LogL= -672.07 0.8883 554 DF
3 LogL= -672.02 0.8606 554 DF
4 LogL= -672.02 0.8591 554 DF
5 LogL= -672.02 0.8591 554 DF
Tue Mar 16 16:52:15 2021 LogL Converged
运行过程如下:
> Var(res)
Term Sigma SE Z.ratio
1 Residual 0.85911 0.10752 7.990234
2 nrm(TreeID) 0.57500 0.40652 1.414445
> pin(res,mulp=c(h2~V2/(V1+V2),Vp~V1+V2))
variance components are as following:
Term Sigma SE vcS
1 Residual 0.85911 0.10752 V1
2 nrm(TreeID) 0.57500 0.40652 V2
pin formula:
h2 ~ V2/(V1 + V2)
Vp ~ V1 + V2
Term Estimate SE
1 h2 0.401 0.196
11 Vp 1.434 0.322
> prd<-predict(res)
> names(prd)
[1] "heads" "pred" "ased"
> prd$pred
$pred1
Prediction Stnd_Error Ecode Rep
1 -0.08508 0.20736 E 1
2 -0.13336 0.19901 E 2
3 -0.12722 0.18984 E 3
4 0.04244 0.18837 E 4
5 -0.29345 0.18827 E 5
当然还可以进行批量分析,以之前例子类似,不做演示了。
https://wap.sciencenet.cn/blog-1114360-1277082.html
上一篇:
AFEchidna包的简单示例7--模型比较下一篇:
AFEchidna包免费对外开放