範 例 八： 單性狀混合模式
程式檔名稱 : LINEAR8.htm
PROC IML;
sigma_a2 =1; * 累加遺傳變方;
sigma_e2 =2; * 剩餘變方;
y={ 2.3, 2.2, 1.7, 1.8, 1.9, 1.9,
1.8, 1.7, 1.9, 1.7};
* y: 觀測值,五隻豬的 左、右半邊屠體背脂;
* X: 固定效應指標變數矩陣;
* 觀測值1-6 固定效應 1,餵飼方法 1;
* 觀測值7-10 固定效應 2;餵飼方法 2;
X={1 0,
1 0,
1 0,
1 0,
1 0,
1 0,
0 1,
0 1,
0 1,
0 1};
* Z: 逢機效應指標變數矩陣;
Z={ 1 0 0 0 0,
1 0 0 0 0,
0 1 0 0 0,
0 1 0 0 0,
0 0 1 0 0,
0 0 1 0 0,
0 0 0 1 0,
0 0 0 1 0,
0 0 0 0 1,
0 0 0 0 1};
* A: 5 隻豬累加性親屬係數矩陣;
A={ 1 0 0.5 0 0.25 ,
0 1 0.5 0.5 0.25 ,
0.5 0.5 1 0.25 0.5 ,
0 0.5 0.25 1 0.125,
0.25 0.25 0.5 0.125 1 };
G = A * sigma_a2;
R = I(10) * sigma_e2;
PRINT y,X,Z;
PRINT G, R;
* 以條件期望值法(conditional expectation) 求得最佳無偏預測值值, BLUP;
V = Z*G*Z` + R; * V：變方矩陣;
IV = INV(V);
b = INV (X`*IV*X) * (X`*IV*y); * 以最小平方法求得固定效應 b 矩陣的估計值;
a = G*Z`*IV*(y-X*b); * 逢機效應a 矩陣的估計值;
e = R *IV*(y-X*b); * 殘差矩陣;
PRINT b, a, e;
* 以混合模式方程式組( Mixed model equation) 求得最佳無偏預測值值, BLUP;
IR = INV(R);
XRX = X`*IR*X;
XRZ = X`*IR*Z;
ZRX = XRZ`;
ZRZG = Z`*IR*Z + INV(G);
* . 混合模式方程式.� . . ;
C1 = XRX || XRZ; * C1 . XRX XRZ . b . = . XRy ;
C2 = ZRX || ZRZG; * C2 . ZRX ZRZ+G . a . . ZRy ;
C = C1 // C2; * ......................................;
IC = INV(C); * C 矩陣為 C1,C2的結合. . . RHS ;
* 於 XRY 內 R= INV(R) . . . 右矩陣;
XRy = X`*IR*y; * ..................................... ;
ZRy = Z`*IR*y; * IC為 C的逆矩陣 . . 方程式求解;
* IC=INV(C)==> . SOL . = .IC*RHS ;
RHS = XRY//ZRy; *.......................................;
* . 單性狀簡化方程式組 . . ;
SOL = IC * RHS; * . XX XZ . b . . Xy ;
* . ZX ZZ+afa . a . . Zy ;
b = SOL(|1:2 ,|);
a = SOL(|3:7 ,|);
e = y - X*b - Z*a;
PRINT XRX, XRZ, ZRX, ZRZG, XRy, ZRy;
PRINT RHS, SOL, a, b;
輸出 8-1 混合模式原死始資料中的 Y 矩陣、X 矩陣及 Z 矩陣
Y X Z
2.3 1 0 1 0 0 0 0
2.2 1 0 1 0 0 0 0
1.7 1 0 0 1 0 0 0
1.8 1 0 0 1 0 0 0
1.9 1 0 0 0 1 0 0
1.9 1 0 0 0 1 0 0
1.8 0 1 0 0 0 1 0
1.7 0 1 0 0 0 1 0
1.9 0 1 0 0 0 0 1
1.7 0 1 0 0 0 0 1
輸出 8-2 混合模式中的 G 矩陣 及 R 矩陣
G R
1 0 0.5 0 0.25 2 0 0 0 0 0 0 0 0 0
0 1 0.5 0.5 0.25 0 2 0 0 0 0 0 0 0 0
0.5 0.5 1 0.25 0.5 0 0 2 0 0 0 0 0 0 0
0 0.5 0.25 1 0.125 0 0 0 2 0 0 0 0 0 0
0.25 0.25 0.5 0.125 1 0 0 0 0 2 0 0 0 0 0
0 0 0 0 0 2 0 0 0 0
0 0 0 0 0 0 2 0 0 0
0 0 0 0 0 0 0 2 0 0
0 0 0 0 0 0 0 0 2 0
0 0 0 0 0 0 0 0 0 2
輸出 8-3 最小平方法求得固定效應 b 向量及由條件期望值法求得
最佳無偏預測值 a 與剩餘的殘差 e
B
1.975
1.8125
A
0.125
-0.125
-0.025
-0.0625
-0.0125
E
0.2
0.1
-0.15
-0.05
-0.05
-0.05
0.05
-0.05
0.1
-0.1
輸出 8-4 混合模式方程式組中的左邊C 矩陣的次矩陣(XRX、XRZ、ZRX 及 ZRZ+G)
XRX XRZ
3 0 1 1 1 0 0
0 2 0 0 0 1 1
ZRX ZRZG
1 0 2.5 0.5 -1 0 0
1 0 0.5 2.833 -1 -0.667 0
1 0 -1 -1 3.333 0 -0.667
0 1 0 -0.667 0 2.333 0
0 1 0 0 -0.667 0 2.333
輸出 8-5 混合模式方程式組中的右邊的次矩陣(RHS 與 RHS 內的 XRY 及 ZRY)
XRY RHS
5.9 5.9
3.55 3.55
2.25
ZRY 1.75
2.25 1.9
1.75 1.75
1.9 1.8
1.75
1.8
輸出 8-6 混合模式方程式組求解後與前述方法獲得相同的解答
SOL B
1.975 1.975
1.8125 A 1.8125
0.125 0.125
-0.125 -0.125
-0.025 -0.025
-0.0625 -0.0625
-0.0125 -0.0125