如何计算育种值的标准误及准确性（accruacy）

参考书籍：王金玉, 陈国宏. 《数量遗传与动物育种》. 东南大学出版社, 2004，P200，第四节：单性状BLUP育种值估计

育种值的标准误可以反映估算的精确性，但是它不能对不同性状或者同一性状不同的试验的BLUP的精确性进行比较，因为单位、试验不同，它仅仅能比较同一试验的准确性，这也限制了它的应用。如果我们想要对不同试验的育种值的准确性进行比较，就需要用到参数：准确性（accuracy）

文中标准差为100，遗传力为1/3，可以计算出：加性方差组分：Va=100^2*(1/3)=3333.33，误差方差组分为：V3=6666.67，因为这里遗传力已知，我们可以设置方差组分固定，进行混合线性求解。

数据

dat <- data.frame(id=c(4,5,6),sire = c(1,3,3),dam=c(2,2,4),y=c(200,170,180))


dat

asreml模型

asreml模型中，默认的是根据数据估计出方差组分，然后计算出K值，求解混合线性方程组，得到BLUE值和BLUP值。由于这里数据量较少，默认估计的方差组分误差较大。书中给出了遗传力和方差组分，因此我们需要将方差组分予以固定，然后求解。

for( i in1:3) dat[,i] <- as.factor(dat[,i]) # 转化为因子


library(asreml)ainv <- asreml.Ainverse(dat[,1:3])$ginv # 计算亲缘关系逆矩阵


Va <- 3333.33


names(Va) <- c("F") #固定加性组分


Ve <- 6666.67


names(Ve) <- c("F") #固定残差组分


mode <- asreml(y ~ 1, random=~ ped(id,init=Va), family=asreml.gaussian(dispersion = Ve), ginverse =list(id=ainv),data=dat)

方差组分

summary(mode)$varcomp

计算育种值和标准误

summary(mode,all=T)$coef.random

注意，这里的std error是标准误，它的计算方法是LHS（混合线性方程组左边的矩阵）求逆，然后得到对角线对应值da,标准误=sqrt(Ve*da)，LHS如图：

计算育种值的可靠性

根据计算公式：r=sqrt(1-dK)，这里的K由遗传力得到，K=Ve/Va=2，由于asreml中se已经知道，我们可以对公式进行转化：

r=sqrt(1-(se^2/Ve)(Ve/Va))=sqrt(1-se^2/Va)，

这个公式还可以进一步的考虑近交系数
r=sqrt(1-se^2/((1+F)Va)(Gilmour et al.,2014)

blup <-summary(mode,all=T)$coef.randomblup <- as.data.frame(blup)blup$ID <- row.names(blup)blup <- blup[,c(4,1,2)]se <- blup$`std error`blup$r <- sqrt(1-se^2/Va)blup

结果中，solution为BLUP育种值，std error是标准误，r是准确性。由此我们可以知道，asreml中BLUP的准确性计算的方法。注意，这里Va和Ve是固定的，如果是真实数据，Va要用summary中的加性效应的方差组分，这里不再赘述。

考虑近交系数的准确性的计算

F <- asreml.Ainverse(dat[,1:3])$inbreedF

blup$r_in <- sqrt(1-se^2/((1+F)*Va))blup

p <- predict(mode,"id")$predictionsstr(p)

ASReml: Thu Aug 03 10:45:47 2017     LogLik         S2      DF      wall     cpu     -9.6504   6666.6700     2  10:45:48     0.0     -9.6504   6666.6700     2  10:45:48     0.0     -9.6504   6666.6700     2  10:45:48     0.0     -9.6504   6666.6700     2  10:45:48     0.0Finished on: Thu Aug 03 10:45:48 2017LogLikelihood Converged List of 2 $ pvals:Classes 'asremlPredict' and 'data.frame':    6 obs. of  4 variables:  ..$ id             : int [1:6] 1 3 2 4 5 6  ..$ predicted.value: num [1:6] 186 180 184 188 179 ...  ..$ standard.error : num [1:6] 75.5 71.2 68.1 57.2 58.6 ...  ..$ est.status     : chr [1:6] "Estimable" "Estimable" "Estimable" "Estimable" ...  ..- attr(*, "heading")= chr "- The predictions are obtained by averaging across the hypertable\n  calculated from model terms constructed solely from factor"| __truncated__ $ avsed: Named num 65.2  ..- attr(*, "names")= chr "overall"

t <- p$pvals$standard.error

blup$r_int <- sqrt(1-t^2/(Va))blup

t

1. 75.5068648089995




2. 71.1648879012397




3. 68.0795355796378




4. 57.1852656089811




5. 58.5823629291693




6. 55.7531699473181