一道數(shù)量遺傳學(xué)題

題目

設(shè)動(dòng)物個(gè)體效應(yīng)為隨機(jī)遺傳效應(yīng)（a）变隔，日糧、性別和畜舍為固定環(huán)境效應(yīng)（b）做盅，背膘厚的遺傳力為0.4，請完成以下工作：
1窘哈，建立背膘厚的線性模型
2吹榴，寫出模型的一般形式和矩陣形式
3，寫出混合線性模型方程組的各組分成分
4滚婉，獲得的估計(jì)值具有哪些特點(diǎn)
5图筹，不同日糧和性別的效應(yīng)值是多少
6，個(gè)體育種值是多少让腹，是否和表型值排序一致远剩？說明理由

處理思路

線性模型已經(jīng)很清楚：
固定因子：日糧，性別哨鸭，畜舍
隨機(jī)因子：加性效應(yīng)
觀測值：背膘厚

矩陣形式：在R語言中構(gòu)建即可
方差組分形式：因?yàn)檫z傳力為0.4民宿，可以假定加性Va=2，Ve=3像鸡，則遺傳力為：2/(2+3) = 0.4
問題4，問題5，問題6需要根據(jù)結(jié)果來解答

解決方案1：R語言

根據(jù)公式建立混合方程組只估，確定固定因子矩陣Z志群，隨機(jī)因子矩陣X，親緣關(guān)系逆矩陣solve(A)

ID <- c("A1","A2","DA","M1","C2","G","M2","CA","S","D","X")
Riliang <- c(1,1,1,1,1,1,2,2,2,2,2)
Sex <- c(2,1,2,2,1,1,2,2,1,2,1)
Sire <- c(0,0,0,0,0,"A2","A2","C2","G","A2","S")
Dam <- c(0,0,"A1","A1","A1","DA","M1",0,"M2","CA","D")
Chushe <- c(1,2,1,3,3,1,3,2,2,2,3)
mm <- c(17,20,15,30,18,12,17,16,23,19,17)
dat <- data.frame(ID,Riliang,Sex,Sire,Dam,Chushe,mm)
dat

dat

ID <- c("A1","A2","DA","M1","C2","G","M2","CA","S","D","X")
Riliang <- c(1,1,1,1,1,1,2,2,2,2,2)
Sex <- c(2,1,2,2,1,1,2,2,1,2,1)
Sire <- c(0,0,0,0,0,"A2","A2","C2","G","A2","S")
Dam <- c(0,0,"A1","A1","A1","DA","M1",0,"M2","CA","D")
Chushe <- c(1,2,1,3,3,1,3,2,2,2,3)
mm <- c(17,20,15,30,18,12,17,16,23,19,17)
dat <- data.frame(ID,Riliang,Sex,Sire,Dam,Chushe,mm)
dat

ped <- dat[,c(1,4,5)]
ainv <- asreml.Ainverse(ped)$ginv
ainv_mat <- asreml.sparse2mat(ainv)  
row.names(ainv_mat) = colnames(ainv_mat) <- attr(ainv,"rowNames")
round(ainv_mat,3)
dim(ainv_mat)

Y = matrix(dat$mm,,1);Y
u = matrix(c(1,2,1,2,1,2,3),,1);u
# x = as.matrix(dat[,c(2,3,6)])/u
X = matrix(c(1,0,0,1,1,0,0,
           1,0,1,0,0,1,0,
           1,0,0,1,1,0,0,
           1,0,0,1,0,0,1,
           1,0,1,0,0,0,1,
           1,0,1,0,1,0,0,
           0,1,0,1,0,0,1,
           0,1,0,1,0,1,0,
           0,1,1,0,0,1,0,
           0,1,0,1,0,1,0,
           0,1,1,0,0,0,1
           ),11,7,byrow = T);X
Z = diag(11);Z


tXX = t(X)%*%X;tXX
tXZ = t(X)%*%Z;tXZ

tZX = t(Z)%*%X;tZX
tZZk = t(Z)%*%Z + ainv_mat*1.5
dim(tZZk)

LHS = rbind(cbind(tXX,tXZ),cbind(tZX,tZZk))
dim(LHS)

tXY = t(X)%*%Y
tZY = t(Z)%*%Y
RHS = rbind(tXY,tZY)
dim(RHS)

library(MASS)
ab = ginv((LHS))%*%RHS
ab

R語言運(yùn)行結(jié)果

result

解決方法2： asreml處理代碼：

for(i in 1:6) dat[,i] <- as.factor(dat[,i])
moda <- asreml(mm ~ Riliang + Sex + Chushe,random = ~ ped(ID),
               ginverse = list(ID = ainv) ,data=dat,start.values = T)
t <- moda$gammas.table
t$Value <- c(2,3)
t$Constraint <- "F"
modb <- asreml(mm ~ Riliang + Sex + Chushe,random = ~ ped(ID),
               ginverse = list(ID = ainv) ,data=dat,
               G.param = t,R.param = t)
summary(modb)$varcomp
coef(modb)$fixed
dim(ab)
ab[8:18]
coef(modb)$random

asrmel運(yùn)行