H22 BioS t (i) treat1 treat2 data d1; input patno treat1 treat2; cards; ; run; 1 (i) treat = 1 treat =

【解答】固定効果分散分析１

−クロスオーバー試験の解析１−

H22

年度

BioS

継続勉強会：第３回

土居正明

【データセットの作成】

ではまず、データセットを作成します。最初に「問題１」の対応のあるt検定用（(i)のプロットは除く）のデータセット です。なお、treat1：プラセボの検査値、treat2：実薬の検査値、とします。 data d1;

 input patno treat1 treat2;  cards;  1 10 11  2 12 20  3 18 17  4 18 18  5 10 15  6 13 20  ; run; です。次に、「問題1」の(i)のプロットと、「問題２」「問題３」を全て扱うデータセットです。なお、treat = 1：プラセボ、 treat = 2：実薬です。

data d2;

 input group patno period treat y;

 label group = ’群’ patno = ’被験者番号’ period = ’時期’ ;  cards;  1 1 1 1 10  1 1 2 2 11  1 2 1 1 12  1 2 2 2 20  1 3 1 1 18  1 3 2 2 17  2 4 1 2 18  2 4 2 1 18  2 5 1 2 15  2 5 2 1 10  2 6 1 2 20  2 6 2 1 13 run; です。

問題１：対応のある

t

検定

(i)プロットのプログラムは、

proc gplot data = d2;  plot y * period = patno;  by group;

run;

です（軸やsymbolの設定等は略）。このとき、プロットは以下のようになります。

(ii)対応のあるt検定の解析プログラムは以下の通りです。

proc ttest data = d1;  paired treat2 * treat1; run;

paired-t:proc

ttest

TTEST プロシジャ 統計量 平均の 平均の 標準偏差の 標準偏差の 差 N 信頼限界の下限 平均 信頼限界の上限 信頼限界の下限 標準偏差 信頼限界の上限 標準誤差 最小値 最大値 treat2 -treat1 6 -0.686 3.3333 7.3524 2.3905 3.8297 9.3928 1.5635 -1 8 t 検定 差 自由度 t 値 Pr > | t | treat2 -treat1 5 2.13 0.0862

(iii)まずデータは 群 被験者番号 実薬投与後の値 プラセボ投与後の値 （実薬）−（プラセボ）: difi 1 1 11 10 1 1 2 20 12 8 1 3 17 18 −1 2 4 18 18 0 2 5 15 10 5 2 6 20 13 7 となります。被験者番号iの（実薬）−（プラセボ）のデータをdifiとおき、 difi ∼ N(µ, σ2) と仮定します（データは全て独立）。薬剤の効果をみる仮説は、両側で考えると H0: µ = 0 H1: µ6= 0 となります。対応のあるt検定とは、この 差のデータに対して行う1群のt検定のこと です。 では、手計算を行いながら(ii)の結果と比較していきます。平均・分散の推定値はそれぞれ bµ = 1 6(1 + 8− 1 + 0 + 5 + 7) = 3.333 bσ2₌ 1 5 © (1− 3.333)2+ (8− 3.333)2+ (−1 − 3.333)2+ (0− 3.333)2+ (5− 3.333)2+ (7− 3.333)2ª = 14.667 となります。平均の推定値は、(ii)の平均の値に一致しています。次に、標準偏差 の推定値は bσ =√bσ2₌√_{14.667 = 3.8298} となり、(ii)の標準偏差の値にほぼ一致しています。では次に、平均の95％信頼区間 です。これは第１回「t検定・分散分 析とデザイン行列」の回でも見ました通り、 " bµ + t(5, 0.025) · r bσ2 6 , bµ + t(5, 0.975) · r bσ2 6 # = " 3.333− 2.571 · r 14.667 6 , 3.333 + 2.571· r 14.667 6 # = [−0.686, 7.353] となりますので、表と大体一致しています。次に、平均の標準誤差 ですが、 p V[bµ] = s c σ2 6 = r 14.667 6 = 1.563 となり、これも表とほぼ一致しています。 では、以下検定に移ります。帰無仮説はH0: µ = 0より、検定統計量(t値)は t = bµ − 0q b σ2 6 =q3.333 14.667 6 = 2.132

となり、p値はSASのデータステップで data p value;  p =2* (1 - cdf(’t’, 2.132, 5)); run; としますとp = 0.0862となり、表と一致します。

問題２：薬剤と人を固定効果とした固定効果分散分析

(i) proc glmを用いた解析プログラムは以下の通りとなります。 proc glm data = d2;  class treat patano;

 model y = treat patno/ ss3 solution p; run; 解析結果は以下の通りです。

paierd-t : proc glm

GLMプロシジャ 従属変数: y 変動因 自由度 平方和 平均平方 F値 Pr>F Model 6 123.0000000 20.5000000 2.80 0.1396 Error 5 36.6666667 7.3333333 Corrected Total 11 159.6666667 変動因 自由度 Type III平方和 平均平方 F値 Pr>F treat 1 33.33333333 33.33333333 4.55 0.0862 patno 5 89.66666667 17.93333333 2.45 0.1744

パラメータ 推定値 標準誤差 t値 Pr>|t| Intercept 18.16666667 B 2.06827894 8.78 0.0003 treat 1 -3.33333333 B 1.56347192 -2.13 0.0862 treat 2 0.00000000 B . . . patno 1 -6.00000000 B 2.70801280 -2.22 0.0776 patno 2 -0.50000000 B 2.70801280 -0.18 0.8608 patno 3 1.00000000 B 2.70801280 0.37 0.7270 patno 4 1.50000000 B 2.70801280 0.55 0.6035 patno 5 -4.00000000 B 2.70801280 -1.48 0.1997 patno 6 0.00000000 B . . . Note：X’Xは特異行列です。正規方程式には、一般化逆行列が使用されています。 文字’B’が付けれられた推定値は、一 意的な推定値ではありません。

paierd-t : proc glm

GLMプロシジャ オブザベーション 観測値 予測値 残差 1 10.00000000 8.83333333 1.16666667 2 11.00000000 12.16666667 -1.16666667 3 12.00000000 14.33333333 -2.33333333 4 20.00000000 17.66666667 2.33333333 5 18.00000000 15.83333333 2.16666667 6 17.00000000 19.16666667 -2.16666667 7 18.00000000 19.66666667 -1.66666667 8 18.00000000 16.33333333 1.66666667 9 15.00000000 14.16666667 0.83333333 10 10.00000000 10.83333333 -0.83333333 11 20.00000000 18.16666667 1.83333333 12 13.00000000 14.83333333 -1.83333333 (ii)次に、統計モデルを考えます。薬剤の影響をdi(i = 1 :プラセボ、i = 2 :実薬)、個人の影響をpj(j = 1, 2, 3, 4, 5, 6)と すると、 yij = µ + di+ pj+ ²ij (²ij ∼ N(0, σ2)) (1)

となります。なお、データは全て独立とします。全員のデータを縦に並べると y11 = µ + d1+ p1+ ²11 y21 = µ + d2+ p1+ ²21 y12 = µ + d1+ p2+ ²12 y22 = µ + d2+ p2+ ²22 y13 = µ + d1+ p3+ ²13 y23 = µ + d2+ p3+ ²23 y24 = µ + d2+ p4+ ²24 y14 = µ + d1+ p4+ ²14 y25 = µ + d2+ p5+ ²25 y15 = µ + d1+ p5+ ²15 y26 = µ + d2+ p6+ ²26 y16 = µ + d1+ p6+ ²16 (2) となります。これをベクトル・行列を用いて表記します。まず、 y =                            y11 y21 y12 y22 y13 y23 y24 y14 y25 y15 y26 y16                            =                            10 11 12 20 18 17 18 18 15 10 20 13                            , X1=                            1 1 0 1 0 0 0 0 0 1 0 1 1 0 0 0 0 0 1 1 0 0 1 0 0 0 0 1 0 1 0 1 0 0 0 0 1 1 0 0 0 1 0 0 0 1 0 1 0 0 1 0 0 0 1 0 1 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0 1 0 1 0 0 0 0 1 0 1 1 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 1 1 1 0 0 0 0 0 0 1                            , β₁=                    µ d1 d2 p1 p2 p3 p4 p5 p6                    , ² =                            ²11 ²21 ²12 ²22 ²13 ²23 ²24 ²14 ²25 ²15 ²26 ²16                            とおくと、 y = X1β1+ ² (² ∼ N(0, σ 2_I 12)) (3) となります。ここで、SASは制約条件d2= 0, p6= 0を入れて考えます。このとき、(2)は y11 = µ + d1 +p1 +²11 y21 = µ +p1 +²21 y12 = µ + d1 +p2 +²12 y22 = µ +p2 +²22 y13 = µ + d1 +p3 +²13 y23 = µ +p3 +²23 y24 = µ +p4 +²24 y14 = µ + d1 +p4 +²14 y25 = µ +p5 +²25 y15 = µ + d1 +p5 +²15 y26 = µ +²26 y16 = µ + d1 +²16 (4)

となり、 X11=                            1 1 1 0 0 0 0 1 0 1 0 0 0 0 1 1 0 1 0 0 0 1 0 0 1 0 0 0 1 1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 0 0 1 0 1 1 0 0 0 1 0 1 0 0 0 0 0 1 1 1 0 0 0 0 1 1 0 0 0 0 0 0 1 1 0 0 0 0 0                            , β₁₁=               µ d1 p1 p2 p3 p4 p5               とおくと、(3)は y + X11β11+ ² (² ∼ N(0, σI12)) (5) と書き直せます。 (iii)さて、これで準備が整いましたので、SASの出力を手計算で求めることとします。ただし、計算を自然に行う順番で見 ていきます。ですので、以後の解説は、SASのアウトプットの順番とは異なりますのでご注意ください。まずは、推定値か らです。 パラメータβ₁₁の最小二乗推定値は bβ11= (X011X11)−1X011y (6) で得られます。これを少しずつ計算すると X0₁₁X11=               1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 0 0 1 0 1 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0                                          1 1 1 0 0 0 0 1 0 1 0 0 0 0 1 1 0 1 0 0 0 1 0 0 1 0 0 0 1 1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 0 0 1 0 1 1 0 0 0 1 0 1 0 0 0 0 0 1 1 1 0 0 0 0 1 1 0 0 0 0 0 0 1 1 0 0 0 0 0                            =               12 6 2 2 2 2 2 6 6 1 1 1 1 1 2 1 2 0 0 0 0 2 1 0 2 0 0 0 2 1 0 0 2 0 0 2 1 0 0 0 2 0 2 1 0 0 0 0 2              

となります。この逆行列を求めると、 (X011X11)−1 =               12 6 2 2 2 2 2 6 6 1 1 1 1 1 2 1 2 0 0 0 0 2 1 0 2 0 0 0 2 1 0 0 2 0 0 2 1 0 0 0 2 0 2 1 0 0 0 0 2               −1 =               7 12 − 1 6 − 1 2 − 1 2 − 1 2 − 1 2 − 1 2 −1 6 1 3 0 0 0 0 0 −1 2 0 1 1 2 1 2 1 2 1 2 −1 2 0 1 2 1 1 2 1 2 1 2 −1 2 0 1 2 1 2 1 1 2 1 2 −1 2 0 1 2 1 2 1 2 1 1 2 −1 2 0 1 2 1 2 1 2 1 2 1               です。また、 X011y =               1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 0 0 1 0 1 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0                                          10 11 12 20 18 17 18 18 15 10 20 13                            =               10 + 11 + 12 + 20 + 18 + 17 + 18 + 18 + 15 + 10 + 20 + 13 10 + 12 + 18 + 18 + 10 + 13 10 + 11 12 + 20 18 + 17 18 + 18 15 + 10               =               182 81 21 32 35 36 25              

となります。以上より(6)は、 bβ11= (X011X11)−1X011y =               7 12 − 1 6 − 1 2 − 1 2 − 1 2 − 1 2 − 1 2 −1 6 1 3 0 0 0 0 0 −1 2 0 1 1 2 1 2 1 2 1 2 −1 2 0 1 2 1 1 2 1 2 1 2 −1 2 0 1 2 1 2 1 1 2 1 2 −1 2 0 1 2 1 2 1 2 1 1 2 −1 2 0 1 2 1 2 1 2 1 2 1                             182 81 21 32 35 36 25               =               7 12 · 182 − 1 6 · 81 − 1 2· 21 − 1 2· 32 − 1 2· 35 − 1 2· 36 − 1 2· 25 −1 6· 182 + 1 3· 81 −1 2· 182 + 1 · 21 + 1 2· 32 + 1 2· 35 + 1 2· 36 + 1 2· 25 −1 2· 182 + 1 2· 21 + 1 · 32 + 1 2· 35 + 1 2· 36 + 1 2· 25 −1 2· 182 + 1 2· 21 + 1 2· 32 + 1 · 35 + 1 2· 36 + 1 2· 25 −1 2· 182 + 1 2· 21 + 1 2· 32 + 1 2· 35 + 1 · 36 + 1 2· 25 −1 2· 182 + 1 2· 21 + 1 2· 32 + 1 2· 35 + 1 2· 36 + 1 · 25               =               18.1667 −3.3333 −6 −0.5 1 1.5 −4               となります。これより、bβ11の各成分と見比べると、 bµ = 18.1667 b d1=−3.3333 bp1=−6, bp2=−0.5, bp3= 1, bp4= 1.5, bp5=−4 となります。制約条件より、db2= 0, bp5= 0より、全パラメータの最小二乗推定値は                                      bµ = 18.1667 ← Intercept b d1 = −3.3333 ← treat1 b d2 = 0 ← treat2 bp1 = −6 ← patno1 bp2 = −0.5 ← patno2 bp3 = 1 ← patno3 bp4 = 1.5 ← patno4 bp5 = −4 ← patno5 bp6 = 0 ← patno6 となり、それぞれSASの出力と一致しています。また、これはbβ1の各成分となります。  では次に、この推定値から予測値・残差を計算していきます。このモデルにおける 予測値 をby1とおくと、 by1= X1bβ1

=                            1 1 0 1 0 0 0 0 0 1 0 1 1 0 0 0 0 0 1 1 0 0 1 0 0 0 0 1 0 1 0 1 0 0 0 0 1 1 0 0 0 1 0 0 0 1 0 1 0 0 1 0 0 0 1 0 1 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0 1 0 1 0 0 0 0 1 0 1 1 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 1 1 1 0 0 0 0 0 0 1                                               bµ b d1 b d2 bp1 bp2 bp3 bp4 bp5 bp6                    =                            1 1 0 1 0 0 0 0 0 1 0 1 1 0 0 0 0 0 1 1 0 0 1 0 0 0 0 1 0 1 0 1 0 0 0 0 1 1 0 0 0 1 0 0 0 1 0 1 0 0 1 0 0 0 1 0 1 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0 1 0 1 0 0 0 0 1 0 1 1 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 1 1 1 0 0 0 0 0 0 1                                               18.1667 −3.3333 0 −6 −0.5 1 1.5 −4 0                    =                            18.1667− 3.3333 − 6 18.1667 + 0− 6 18.1667− 3.3333 − 0.5 18.1667 + 0− 0.5 18.1667− 3.3333 + 1 18.1667 + 0 + 1 18.1667 + 0 + 1.5 18.1667− 3.3333 + 1.5 18.1667 + 0− 4 18.1667− 3.3333 − 4 18.1667 + 0 + 0 18.1667− 3.3333 + 0                           

=                            8.8334 12.1667 14.3334 17.6667 15.8334 19.1667 19.6667 16.3334 14.1667 10.8334 18.1667 14.8334                            となります。これは表の値とほぼ一致しています。次に、このモデルにおける 残差 をe1とおきますと、 e1= y− by1=                            10 11 12 20 18 17 18 18 15 10 20 13                            −                            8.8334 12.1667 14.3334 17.6667 15.8334 19.1667 19.6667 16.3334 14.1667 10.8334 18.1667 14.8334                            =                            1, 16666 −1.1667 −2.3334 2.3333 2.1666 −2.1667 −1.6667 1.6666 0.8333 −0.8334 1.8333 −1.8334                            となります。これも、表とほぼ一致しています。 (iv)では次に、分散分析表を求めます。そのためには、帰無仮説のもとでのモデルが必要となります。分散分析表における 帰無仮説・対立仮説は H0: d1= d2 かつ p1= p2= p3= p4= p5= p6（薬剤の影響も個人の影響も全くない） H1:       それ以外 です。これより、帰無仮説のもとでのモデル を考えます。このモデルは、(1)からdiとpjを除いたもの であり、 yij = µ + ²ij

となります。全データを並べると、 y11 = µ + ²11 y21 = µ + ²21 y12 = µ + ²12 y22 = µ + ²22 y13 = µ + ²13 y23 = µ + ²33 y24 = µ + ²24 y14 = µ + ²14 y25 = µ + ²25 y15 = µ + ²15 y26 = µ + ²26 y16 = µ + ²16 となります。ベクトル・行列表記しますと、 y = µ112+ ² (² ∼ N(0, σ2I12)) となります。これより、このモデルでのµの最小二乗推定値bµ0は、 bµ0= (1012112)−11012y (7) となります。ここで、 10₁₂112= ³ 1 1 1 1 1 1 1 1 1 1 1 1 ´                            1 1 1 1 1 1 1 1 1 1 1 1                            = 12 より、(10₁₂112)−1= ₁₂1 であり、 1012y = ³ 1 1 1 1 1 1 1 1 1 1 1 1 ´                            10 11 12 20 18 17 18 18 15 10 20 13                           

= 10 + 11 + 12 + 20 + 18 + 17 + 18 + 18 + 15 + 10 + 20 + 13 = 182 となりますので、(7)は bµ0= (1012112)−11012y = 1 12· 182 = 15.1667 となります。これより、このモデルにおける予測値・残差by10, e10は by10=bµ0112=                            15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667                            (8) e10= y− by10=                            10 11 12 20 18 17 18 18 15 10 20 13                            −                            15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667                            =                            −5.1667 −4.1667 −3.1667 4.8333 2.8333 1.8333 2.8333 1.8333 −0.1667 −5.1667 4.8333 −2.1333                            (9) となります。以上より、線形モデルの一般論から ky − by10k2 | {z } Corrected Total =ky − by1k2 | {z } Error +kby1− by10k2 | {z } Model となります*1_。 そして、データの数・制約条件を入れたあとのパラメータ数 は y :データ12個  → 12 by10:パラメータµのみ  → 1 by1:制約条件を入れた後のパラメータµ, d1, p1,· · · , p5 → 7

*1_{意味は、「Corrected Total：H0のモデルで説明できない部分」「Error：H1のモデルで説明できない部分」「Model：H0のモデルと H1のモデルの}

説明能力の差」です。また、残差で書き直すと、

ke10k2=ke1k2+kby1− by0k2 となります。

より、自由度 は (CorrectedTotal)ky − by10k2 : 12− 1 = 11 (Error) ky − by1k2 : 12− 7 = 5 (Model) kby1− by10k2: 7− 1 = 6 となります。数値を求めていきますと、各 平方和 は (CorrectedTotal)ky − by10k2=ke10k2 = (−5.1667)2+ (−4.1667)2+ (−3.1667) + 4.83332+ 2.83332+ 1.83332 +2.83332+ 2.83332+ (−0.1667)2+ (−5.1667)2+ 4.83332+ (−2.1333)2 = 159.5230 (Error)ky − by1k2=ke1k2 = 1.16662+ (−1.1667)2+ (−2.3334) + 2.33332+ 2.16662+ (−2.1667)2 +(−1.6667)2+ 1.66662+ 0.83332+ (−0.8334)2+ 1.83332+ (−1.8334)2 = 36.6667 (Model)kby1− by10k2= °° °° °° °° °° °° °° °° °° °° °° °° °° °° °° °                            8.8334 12.1667 14.3334 17.6667 15.8334 19.1667 19.6667 16.3334 14.1667 10.8334 18.1667 14.8334                            −                            15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667                            °° °° °° °° °° °° °° °° °° °° °° °° °° °° °° ° 2 = °° °° °° °° °° °° °° °° °° °° °° °° °° °° °° °                            −6.3333 −3 −0.8333 2.5 0.6667 4 4.5 1.1333 −1 −4.333 3 −0.3333                            °° °° °° °° °° °° °° °° °° °° °° °° °° °° °° ° 2 = (−6.3333)2+ (−3)2+ (−0.8333)2+ 2.52+ 0.66672+ 42 +4.52+ 1.13332+ (−1)2+ (−4.333)2+ 32+ (−0.3333)2 = 122.9199 となり、表とほぼ一致します。 次に、平均平方 です。自由度で平方和を割ってやればよいですので、 (Error) 1 5ky − by1k 2₌ 1 5· 36.6667 = 7.3334 (Model) 1 6kby1− by10k 2₌1 6 · 122.9199 = 20.4867 となり、表の値とほぼ一致しています。では次にF 値 です。平均平方の比をとってやれば F = Modelの平均平方 Errorの平均平方 = 20.4867 7.3334 = 2.7936 となります。なお、分母のErrorの平均平方には 分散σ2_{の不偏推定値 という意味もあります。つまり、} bσ2 1 = 1 kek2 = 1· 36.6667 = 7.3334 (10)

となります*2_{。これは、薬剤の影響の有無をみる検定をする際にも用います}*3_。  最後にp値 です。分子の自由度が6、分母の自由度が5ですので、データステップで data p value;  p = 1 - cdf(’F’, 2.7936, 6, 5); run; としてやりますと、p = 0.1397が得られ、表の値と大体一致しています。 (v)次に薬剤の影響を見る検定です。仮説は H0d: d1= d2 H1d: d16= d2 となります。以下、帰無仮説H0dのもとでの統計モデル を考えていきましょう。モデルは(1)からdiを除いた yij = µ + pj+ ²ij (²ij ∼ N(0, σ2)) となります。データを全て並べますと、 y11 = µ + p1+ ²11 y21 = µ + p1+ ²21 y12 = µ + p2+ ²12 y22 = µ + p2+ ²22 y13 = µ + p3+ ²13 y23 = µ + p3+ ²23 y14 = µ + p4+ ²14 y24 = µ + p4+ ²24 y15 = µ + p5+ ²15 y25 = µ + p5+ ²25 y16 = µ + p6+ ²16 y26 = µ + p6+ ²26 *2_bσ2 1の下の添え字 1 は、「問題２」（1 つ目）のものという意味です。「問題３」でも別のモデルで同様に分散 σ2の推定値を考え、それはbσ22と書 きます。 *3 F = 1 5kby1− by10k2 b σ2 1 という表現で理解しておくと、薬剤の影響を検討する際に便利です。

となります。ベクトル・行列表記を考えます。まず X1d=                            1 1 0 0 0 0 0 1 1 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 0 1 1 0 0 0 0 0 1                            , β_1d=               µ p1 p2 p3 p4 p5 p6               とおきますと、 y = X1dβ1d+ ² (² ∼ N(0, σ 2_I 12)) (11) と書けます。ここで、SASは制約条件p6= 0を入れます。このとき、 X1d1=                            1 1 0 0 0 0 1 1 0 0 0 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 1 0 1 0 0 0 1 0 1 0 0 0 0 1 1 0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0                            , β_1d1=            µ p1 p2 p3 p4 p5            とおきますと、(11)は y = X1d1β1d1+ ² (² ∼ N(0, σ 2_I 12)) と書けます。これを用いて、β_1d1の最小二乗推定値を計算していきます。 bβ_1d1= (X0_1d1X1d1)−1X01d1y (12)

となりますので、少しずつ計算していきます。まず、 X0_1d1X1d1=            1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0                                       1 1 0 0 0 0 1 1 0 0 0 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 1 0 1 0 0 0 1 0 1 0 0 0 0 1 1 0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0                            =            12 2 2 2 2 2 2 2 0 0 0 0 2 0 2 0 0 0 2 0 0 2 0 0 2 0 0 0 2 0 2 0 0 0 0 2            です。これより、 (X0_1d1X1d1)−1=            12 2 2 2 2 2 2 2 0 0 0 0 2 0 2 0 0 0 2 0 0 2 0 0 2 0 0 0 2 0 2 0 0 0 0 2            −1 =            1 2 − 1 2 − 1 2 − 1 2 − 1 2 − 1 2 −1 2 1 1 2 1 2 1 2 1 2 −1 2 1 2 1 1 2 1 2 1 2 −1 2 1 2 1 2 1 1 2 1 2 −1 2 1 2 1 2 1 2 1 1 2 −1 2 1 2 1 2 1 2 1 2 1           

となります。次に、 X0_1d1y =            1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0                                       10 11 12 20 18 17 18 18 15 10 20 13                            =            10 + 11 + 12 + 20 + 18 + 17 + 18 + 18 + 15 + 10 + 20 + 13 10 + 11 12 + 20 18 + 17 18 + 18 15 + 10            =            182 21 32 35 36 25            です。これより、(12)は bβ1d1= (X01d1X1d1)−1X01d1y =            1 2 − 1 2 − 1 2 − 1 2 − 1 2 − 1 2 −1 2 1 1 2 1 2 1 2 1 2 −1 2 1 2 1 1 2 1 2 1 2 −1 2 1 2 1 2 1 1 2 1 2 −1 2 1 2 1 2 1 2 1 1 2 −1 2 1 2 1 2 1 2 1 2 1                       182 21 32 35 36 25            =            1 2· 182 − 1 2· 21 − 1 2· 32 − 1 2 · 35 − 1 2· 36 − 1 2· 25 −1 2· 182 + 1 · 21 + 1 2· 32 + 1 2· 35 + 1 2· 36 + 1 2· 25 −1 2· 182 + 1 2· 21 + 1 · 32 + 1 2· 35 + 1 2· 36 + 1 2· 25 −1 2· 182 + 1 2· 21 + 1 2· 32 + 1 · 35 + 1 2· 36 + 1 2· 25 −1 2· 182 + 1 2· 21 + 1 2· 32 + 1 2· 35 + 1 · 36 + 1 2· 25 −1 2· 182 + 1 2· 21 + 1 2· 32 + 1 2· 35 + 1 2· 36 + 1 · 25           

=            16.5 −6 −0.5 1 1.5 −4            となります。ここで、 bβ_1d1=            bµ bp1 bp2 bp3 bp4 bp5            であり、制約条件からbp6= 0となりますので、                            bµ = 16.5 bp1 = −6 bp2 = −0.5 bp3 = 1 bp4 = 1.5 bp5 = −4 bp6 = 0 となります。これが bβ1dの各成分となります。 これより、このモデルに基づく予測値by1dは by1d= X1dbβ1d =                            1 1 0 0 0 0 0 1 1 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 0 1 1 0 0 0 0 0 1                                          16.5 −6 −0.5 1 1.5 −4 0              

=                            16.5− 6 16.5− 6 16.5− 0.5 16.5− 0.5 16.5 + 1 16.5 + 1 16.5 + 1.5 16.5 + 1.5 16.5− 4 16.5− 4 16.5 + 0 16.5 + 0                            =                            10.5 10.5 16 16 17.5 17.5 18 18 12.5 12.5 16.5 16.5                            となります*4_{。TypeIIIの平方和のtreatの部分 は、要は「薬剤を入れないモデルH0d」と「入れるモデルH1d（H1と一致）」} で、説明能力がどのくらい違うか、を見ています。そこで、この2つのモデルの予測値by1とby1dの差を考えてやります。 *4_{なお、このモデルでのパラメータの最小二乗推定値や予測値を SAS で出力したい場合、プログラムは} proc glm data = d2;   class patno;

  model y = patno / solution ss3 p;

run;

具体的には、 kby1− by1dk2= °° °° °° °° °° °° °° °° °° °° °° °° °° °° °° °                            8.8334 12.1667 14.3334 17.6667 15.8334 19.1667 19.6667 16.3334 14.1667 10.8334 18.1667 14.8334                            −                            10.5 10.5 16 16 17.5 17.5 18 18 12.5 12.5 16.5 16.5                            °° °° °° °° °° °° °° °° °° °° °° °° °° °° °° ° = °° °° °° °° °° °° °° °° °° °° °° °° °° °° °° °                            −1.6666 1.6667 −1.6666 1.6667 −1.6666 1.6667 −1.6666 1.6667 −1.6666 1.6667 −1.6666 1.6667                            °° °° °° °° °° °° °° °° °° °° °° °° °° °° °° ° 2 = (−1.6666)2+ (1.6667)2+ (−1.6666)2+ (1.6667)2+ (−1.6666)2+ (1.6667)2 +(−1.6666)2+ (1.6667)2+ (−1.6666)2+ (1.6667)2+ (−1.6666)2+ (1.6667)2 = 33.3327 となり、表の数値と大体一致しています。なお、パラメータ数は by1:制約条件を入れた後のパラメータµ, d1, p1,· · · , p5  → 7 by1d:制約条件を入れた後のパラメータµ, p1,·, p5    → 6 より、自由度 は kby1− by1dk2: 7− 6 = 1 です。これより、平均平方 は 1 1 · kby1− by1dk 2 = 33.3327 となります。F 統計量 は、この平均平方を、上で求めました分散の推定値bσ₁2（全体の分散分析表のErrorの平均平方）で 割ってやって F = 1 1· kby1− by1dk 2 bσ2 1 =33.3327 7.3334 = 4.55 となります。p値 は、分母のbσ12の自由度が5であることを考慮して、 data p value;  p = 1 - cdf(’F’, 4.55, 1, 5); run; とおきますと、p = 0.8601が得られます。これはTypeIIIの表のTreatの部分のp値と大体一致しています。 (vi)では、(v)の結果と「問題１」の結果を比較していきます。まず最初にp値 ですが、直前に求めました薬剤の影響を見 る検定のp値は0.8601（表では0.8602）でした。これは、「問題１」の対応のあるt検定のp値と一致しています。 また、この際の 対応のあるt検定のt値と、TypeIIIのtreatの部分のF 値の関係 ですが、 (t値)2= 2.1322= 4.55 = (F値)

という関係が成り立っています。  次に、分散の推定値 について見ていきます。対応のあるt検定では、標準偏差（の推定値）は1.5635でしたので、 (対応のあるt検定のモデルの分散の推定値) = 1.53652= 14.6666 となります。一方、分散分析のモデルでの分散の推定値は(10)より、 bσ2 1= 1 Errorの自由度kek = 1 5 · 36.6667 = 7.3334 となりました。これより、対応のあるt検定の分散は、分散分析モデルの分散の2倍になっている ことが分かります。詳し い説明は補助資料２「対応のあるt検定について」に回しますが、要は被験者番号jの人の実薬・プラセボのデータy1j, y2j が独立で、分散がそれぞれσ2になりますので、difj = y1j− y2jの分散は V [difj] = V [y1j] + V [y2j] = σ2+ σ2= 2σ2 となるのです。  最後に、薬剤の影響 です。対応のあるt検定の表での群間差の平均3.3333は、「Treat2−Treat1」でした。一方、分散分 析の計算結果より、投与群の推定値は ・treat1:−3.3333 ・treat2: 0 となっています。これより、「treat2− treat1=3.3333」となり、一致しています。

問題３

(i)解析プログラムは proc glm data = d2;  class treat patno period;

 model y = treat patno period / solution p ss3; run; となり、出力は以下の通りです。

patno(fixed) + period

GLMプロシジャ 従属変数: y 変動因 自由度 平方和 平均平方 F値 Pr>F Model 7 124.3333333 17.7619048 2.01 0.2606 Error 4 35.3333333 8.8333333 Corrected Total 11 159.6666667

変動因 自由度 Type III平方和 平均平方 F値 Pr>F treat 1 33.33333333 33.33333333 3.77 0.1240 patno 5 89.66666667 17.93333333 2.03 0.2562 period 1 1.33333333 1.33333333 0.15 0.7174 パラメータ 推定値 標準誤差 t値 Pr>|t| Intercept 17.83333333 B 2.42670330 7.35 0.0018 treat 1 -3.33333333 B 1.71593836 -1.94 0.1240 treat 2 0.00000000 B . . . patno 1 -6.00000000 B 2.97209242 -2.02 0.1137 patno 2 -0.50000000 B 2.97209242 -0.17 0.8746 patno 3 1.00000000 B 2.97209242 0.34 0.7534 patno 4 1.50000000 B 2.97209242 0.50 0.6403 patno 5 -4.00000000 B 2.97209242 -1.35 0.2496 patno 6 0.00000000 B . . . period 1 0.66666667 B 1.71593836 0.39 0.7174 period 2 0.00000000 B . . . Note:X’Xは特異行列です。正規方程式には、一般化逆行列が使用されています。 文字’B’が付けれられた推定値は、一意 的な推定値ではありません。

patno(fixed) + period

GLMプロシジャ オブザベーション 観測値 予測値 残差 1 10.00000000 9.16666667 0.83333333 2 11.00000000 11.83333333 -0.83333333 3 12.00000000 14.66666667 -2.66666667 4 20.00000000 17.33333333 2.66666667 5 18.00000000 16.16666667 1.83333333 6 17.00000000 18.83333333 -1.83333333 7 18.00000000 20.00000000 -2.00000000 8 18.00000000 16.00000000 2.00000000 9 15.00000000 14.50000000 0.50000000 10 10.00000000 10.50000000 -0.50000000 11 20.00000000 18.50000000 1.50000000 12 13.00000000 14.50000000 -1.50000000 (ii)では次に、統計モデルを考えます。薬剤di(i = 1 :プラセボ群、i = 2 :実薬群)、個人の効果pj(j = 1, 2, 3, 4, 5, 6)、時 期の効果をtk(k = 1, 2)とすると、 yijk = µ + di+ pj+ tk+ ²ijk (²ij ∼ N(0, σ2)) (13) となります。なお、データは全て独立とします。全員のデータを縦に並べると y111 = µ + d1+ p1+ t1+ ²111 y212 = µ + d2+ p1+ t2+ ²212 y121 = µ + d1+ p2+ t1+ ²121 y222 = µ + d2+ p2+ t2+ ²222 y131 = µ + d1+ p3+ t1+ ²131 y232 = µ + d2+ p3+ t2+ ²232 y241 = µ + d2+ p4+ t1+ ²241 y142 = µ + d1+ p4+ t2+ ²142 y251 = µ + d2+ p5+ t1+ ²251 y152 = µ + d1+ p5+ t2+ ²152 y261 = µ + d2+ p6+ t1+ ²261 y162 = µ + d1+ p6+ t2+ ²162 (14)

です。ここで、 y =                            y111 y212 y121 y222 y131 y232 y241 y142 y251 y152 y261 y162                            =                            10 11 12 20 18 17 18 18 15 10 20 13                            , X2=                            1 1 0 1 0 0 0 0 0 1 0 1 0 1 1 0 0 0 0 0 0 1 1 1 0 0 1 0 0 0 0 1 0 1 0 1 0 1 0 0 0 0 0 1 1 1 0 0 0 1 0 0 0 1 0 1 0 1 0 0 1 0 0 0 0 1 1 0 1 0 0 0 1 0 0 1 0 1 1 0 0 0 0 1 0 0 0 1 1 0 1 0 0 0 0 1 0 1 0 1 1 0 0 0 0 0 1 0 0 1 1 0 1 0 0 0 0 0 1 1 0 1 1 0 0 0 0 0 0 1 0 1                            , β₂=                         µ d1 d2 p1 p2 p3 p4 p5 p6 t1 t2                         , ² =                            ²111 ²212 ²121 ²222 ²131 ²232 ²241 ²142 ²251 ²152 ²261 ²162                            とおくと、 y = X2β2+ ² (² ∼ N(0, σ 2_I 12)) (15) となります。 ここで、SASは制約条件d2= 0, p6= 0, t2= 0を入れますので、 y111 = µ +d1 +p1 +t1 +²111 y212 = µ +p1 +²212 y121 = µ +d1 +p2 +t1 +²121 y222 = µ +p2 +²222 y131 = µ +d1 +p3 +t1 +²131 y232 = µ +p3 +²232 y241 = µ +p4 +t1 +²241 y142 = µ +d1 +p4 +²142 y251 = µ +p5 +t1 +²251 y152 = µ +d1 +p5 +²152 y261 = µ +t1 +²261 y162 = µ +d1 +²162 (16) となります。このとき、 X21=                            1 1 1 0 0 0 0 1 1 0 1 0 0 0 0 0 1 1 0 1 0 0 0 1 1 0 0 1 0 0 0 0 1 1 0 0 1 0 0 1 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0 1 1 1 0 0 0 1 0 0 1 0 0 0 0 0 1 1 1 1 0 0 0 0 1 0 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0                            , β₂₁=                 µ d1 p1 p2 p3 p4 p5 t1                

とおくと、(15)は y = X21β21+ ² (² ∼ N(0, σ 2_I 12)) となります。 (iii)ここで、β21の最小二乗推定値は bβ21= (X021X21)−1X021y (17) となります。これを少しずつ計算していきます。 X021X21=                  1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 0 0 1 0 1 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 1 0 1 0 1 0 1 0 1 0                                             1 1 1 0 0 0 0 1 1 0 1 0 0 0 0 0 1 1 0 1 0 0 0 1 1 0 0 1 0 0 0 0 1 1 0 0 1 0 0 1 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0 1 1 1 0 0 0 1 0 0 1 0 0 0 0 0 1 1 1 1 0 0 0 0 1 0 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0                            =                 12 6 2 2 2 2 2 6 6 6 1 1 1 1 1 3 2 1 2 0 0 0 0 1 2 1 0 2 0 0 0 1 2 1 0 0 2 0 0 1 2 1 0 0 0 2 0 1 2 1 0 0 0 0 2 1 6 3 1 1 1 1 1 6                 より、 (X0₂₁X21)−1 =                 12 6 2 2 2 2 2 6 6 6 1 1 1 1 1 3 2 1 2 0 0 0 0 1 2 1 0 2 0 0 0 1 2 1 0 0 2 0 0 1 2 1 0 0 0 2 0 1 2 1 0 0 0 0 2 1 6 3 1 1 1 1 1 6                 −1

=                 2 3 − 1 6 − 1 2 − 1 2 − 1 2 − 1 2 − 1 2 − 1 6 −1 6 1 3 0 0 0 0 0 0 −1 2 0 1 1 2 1 2 1 2 1 2 0 −1 2 0 1 2 1 1 2 1 2 1 2 0 −1 2 0 1 2 1 2 1 1 2 1 2 0 −1 2 0 1 2 1 2 1 2 1 1 2 0 −1 2 0 1 2 1 2 1 2 1 2 1 0 −1 6 0 0 0 0 0 0 1 3                 です。また、 X0₂₁y =                  1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 0 0 1 0 1 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 1 0 1 0 1 0 1 0 1 0                                             10 11 12 20 18 17 18 18 15 10 20 13                            =                 10 + 11 + 12 + 20 + 18 + 17 + 18 + 18 + 15 + 10 + 20 + 13 10 + 12 + 18 + 18 + 10 + 13 10 + 11 12 + 20 18 + 17 18 + 18 15 + 10 10 + 12 + 18 + 18 + 15 + 20                 =                 182 81 21 32 35 36 25 93                 となります。これより、β₂₁の最小二乗推定値(17)は bβ₂₁= (X0₂₁X21)−1X021y

=                 2 3 − 1 6 − 1 2 − 1 2 − 1 2 − 1 2 − 1 2 − 1 6 −1 6 1 3 0 0 0 0 0 0 −1 2 0 1 1 2 1 2 1 2 1 2 0 −1 2 0 1 2 1 1 2 1 2 1 2 0 −1 2 0 1 2 1 2 1 1 2 1 2 0 −1 2 0 1 2 1 2 1 2 1 1 2 0 −1 2 0 1 2 1 2 1 2 1 2 1 0 −1 6 0 0 0 0 0 0 1 3                                 182 81 21 32 35 36 25 93                 =                 2 3· 182 − 1 6· 81 − 1 2· 21 − 1 2· 32 − 1 2· 35 − 1 2· 36 − 1 2 · 25 − 1 6· 93 −1 6· 182 + 1 3· 81 −1 2· 182 + 1 · 21 + 1 2· 32 + 1 2· 35 + 1 2· 36 + 1 2· 25 −1 2· 182 + 1 2· 21 + 1 · 32 + 1 2· 35 + 1 2· 36 + 1 2· 25 −1 2· 182 + 1 2· 21 + 1 2· 32 + 1 · 35 + 1 2· 36 + 1 2· 25 −1 2· 182 + 1 2· 21 + 1 2· 32 + 1 2· 35 + 1 · 36 + 1 2· 25 −1 2· 182 + 1 2· 21 + 1 2· 32 + 1 2· 35 + 1 2· 36 + 1 · 25 −1 6· 182 + 1 3· 93                 =                 17.8333 −3.3333 −6 −0.5 1 1.5 −4 0.6667                 となります。これより、 bβ21=                 bµ b d1 bp1 bp2 bp3 bp4 bp5 bt1                 から、 bµ = 17.8333, bd1=−3.3333, bp1=−6, bp2=−0.5, bp3= 1, bp4= 1.5, bp5=−4, bt1= 0.6667

となります。次に、制約条件よりdb2= 0, bp6= 0, bt2= 0であることから、 bβ2=                         bµ b d1 b d2 bp1 bp2 bp3 bp4 bp5 bp6 bt1 bt2                         の各成分は                                                bµ = 17.8333 ← Intercept b d1 = −3.3333 ← treat1 b d2 = 0 ← treat2 bp1 = −6 ← patno1 bp2 = −0.5 ← patno2 bp3 = 1 ← patno3 bp4 = 1.5 ← patno4 bp5 = −4 ← patno5 bp6 = 0 ← patno6 bt1 = 0.6667 ← period1 bt2 = 0 ← period2 で与えられます。これは、表の値とほぼ一致しています。  では次に、これを用いて予測値・残差を求めます。このモデルにおける 予測値 をby2とおきますと、 by2= X2bβ2 =                            1 1 0 1 0 0 0 0 0 1 0 1 0 1 1 0 0 0 0 0 0 1 1 1 0 0 1 0 0 0 0 1 0 1 0 1 0 1 0 0 0 0 0 1 1 1 0 0 0 1 0 0 0 1 0 1 0 1 0 0 1 0 0 0 0 1 1 0 1 0 0 0 1 0 0 1 0 1 1 0 0 0 0 1 0 0 0 1 1 0 1 0 0 0 0 1 0 1 0 1 1 0 0 0 0 0 1 0 0 1 1 0 1 0 0 0 0 0 1 1 0 1 1 0 0 0 0 0 0 1 0 1                                                    17.8333 −3.3333 0 −6 −0.5 1 1.5 −4 0 0.6667 0                        

=                            17.8333− 3.3333 − 6 + 0.6667 17.8333− 6 17.8333− 3.3333 − 0.5 + 0.6667 17.8333− 0.5 17.8333− 3.3333 + 1 + 0.6667 17.8333 + 1 17.8333 + 1.5 + 0.6667 17.8333− 3.3333 + 1.5 17.8333− 4 + 0.6667 17.8333− 3.3333 − 4 17.8333 + 0.6667 17.8333− 3.3333                            =                            9.1667 11.8333 14.6667 17.3333 16.1667 18.8333 20 16.0 14.5 10.5 18.5 14.5                            となり、表の出力とほぼ一致します。次に 残差e2は、 e2= y− by2 =                            10 11 12 20 18 17 18 18 15 10 20 13                            −                            9.1667 11.8333 14.6667 17.3333 16.1667 18.8333 20 16.0 14.5 10.5 18.5 14.5                            =                            0.8333 −0.8333 −2.6667 2.6667 1.8333 −1.8333 −2 2 0.5 −0.5 1.5 −1.5                            となり、これも表とほぼ一致します。 (iv)では次に分散分析表を求めます。そのためには、帰無仮説のもとでのモデルが必要となります。今回は、固定効果に 「薬剤」「個人」「時期」の３つがありますので、分散分析表の帰無仮説・対立仮説は、 H0: d1= d2かつp1= p2= p3= p4= p5= p6かつt1= t2     （薬剤の影響も個人の影響も時期の影響も全くない） H1:それ以外 となります。これより、帰無仮説のもとでのモデル は、(13)からdi, pj, tkを除いた yijk = µ + ²ijk

となります。データを全て縦に並べると、 y111 = µ + ²111 y212 = µ + ²212 y121 = µ + ²121 y222 = µ + ²222 y131 = µ + ²131 y232 = µ + ²332 y241 = µ + ²241 y142 = µ + ²142 y251 = µ + ²251 y152 = µ + ²152 y261 = µ + ²261 y162 = µ + ²162 となります。要は、どの要因もモデルに含めない、ということですので、「問題２」の分散分析表の帰無仮説のモデルと全く 同じです。すると、予測値・残差も全く同じになります。(8)、(9)からそれぞれby10、e10とおいていまして、 by10=                            15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667                            , e10=                            −5.1667 −4.1667 −3.1667 4.8333 2.8333 1.8333 2.8333 1.8333 −0.1667 −5.1667 4.8333 −2.1333                            と書けました。以上より、線形モデルの一般論から、 ky − by10k2=ky − by2k2+kby2− by10k2 が成り立ちます。データ数・制約条件を入れたあとのパラメータ数は y :データ12個   → 12 by10:パラメータµのみ   → 1 by2:制約条件を入れた後のパラメータµ, d1, p1,· · · , p5, t1 → 8 より、自由度 は (CorrectedTotal)ky − by10k2 : 12− 1 = 11 (Error) ky − by2k2 : 12− 8 = 4 (Model) kby2− by10k2: 8− 1 = 7 となります。数値を求めていきますと、各 平方和 は (CorrectedTotal)ky − by10k2=ke10k2 = (−5.1667)2+ (−4.1667)2+ (−3.1667) + 4.83332+ 2.83332+ 1.83332 +2.83332+ 2.83332+ (−0.1667)2+ (−5.1667)2+ 4.83332+ (−2.1333)2

= 159.5230 (Error)ky − by2k2=ke2k2 = 0.83332+ (−0.8333)2+ (−2.6667)2+ 2.66672+ 1.83332+ (−1.8333)2 +(−2)2+ 22+ 0.52+ (−0.5)2+ 1.52+ (−1.5)2 = 35.3333 (Model)kby2− by10k2= °° °° °° °° °° °° °° °° °° °° °° °° °° °° °° °                            9.1667 11.8333 14.6667 17.3333 16.1667 18.8333 20 16.0 14.5 10.5 18.5 14.5                            −                            15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667 15.1667                            °° °° °° °° °° °° °° °° °° °° °° °° °° °° °° ° 2 = °° °° °° °° °° °° °° °° °° °° °° °° °° °° °° °                            −6 −3.3334 −0.5 2.1666 1 3.6666 4.8333 0.8333 −0.6667 −4.6667 3.3333 −0.6667                            °° °° °° °° °° °° °° °° °° °° °° °° °° °° °° ° 2 = (−6)2+ (−3.3334)2+ (−0.5)2+ 2.16662+ 12+ 3.66662 4.83332+ 0.83332+ (−0.6667)2+ (−4.6667)2+ 3.33332+ (−0.6667)2 = 124.331 となり、表の出力と大体一致します。  次に、平均平方 です。自由度で平方和を割ってやればよいですので、 (Error) 1 4ky − by2k 2₌1 4 · 35.3333 = 8.8333 (Model) 1 7kby2− by10k 2₌ 1 7· 124.331 = 17.7616 となり、表とほぼ一致します。では、次にF値 です。平均平方の比をとってやれば F = Modelの平均平方 Errorの平均平方 = 17.7616 8.8333 = 2.0108 となります。なお、問題2と同様、分母のErrorの平均平方には分散σ2の不偏推定値という意味もあります。つまり、 bσ2 2= 1 Errorの自由度ke2k 2₌ 1 4 · 35.3333 = 8.8333 (18) となります。ここでも、F 統計量をF = 17kby2−by10k 2 b σ2 2 と捉えておくと、薬剤の検討の際に役立ちます。  最後にp値 です。分子の自由度が7、分母の自由度が4ですので、データステップで、 data p value;  p = 1 - cdf(’F’, 2.0108, 7, 4); run; としますと、p = 0.2606が得られ、表の値とほぼ一致します。

(v)次に、薬剤の影響をみる検定です。仮説は、「問題２」と同じく H0d: d1= d2 H1d: d16= d2 となります。以下、帰無仮説H0dのもとでの統計モデル を考えていきます。なお、最初のモデル自体が「問題２」と「問題 ３」で異なりますので、H0dのもとでのモデルは、「問題２」とは異なっています。 まず、モデルは(13)からdiを除いた yijk= µ + pj+ tk+ ²ijk となります。データを全て並べますと、 y111 = µ + p1+ t1+ ²11 y212 = µ + p1+ t2+ ²21 y121 = µ + p2+ t1+ ²12 y222 = µ + p2+ t2+ ²22 y131 = µ + p3+ t1+ ²13 y232 = µ + p3+ t2+ ²23 y141 = µ + p4+ t1+ ²14 y242 = µ + p4+ t2+ ²24 y151 = µ + p5+ t1+ ²15 y252 = µ + p5+ t2+ ²25 y161 = µ + p6+ t1+ ²16 y262 = µ + p6+ t2+ ²26 です。ベクトル・行列表記を考えますと、まず X2d=                            1 1 0 0 0 0 0 1 0 1 1 0 0 0 0 0 0 1 1 0 1 0 0 0 0 1 0 1 0 1 0 0 0 0 0 1 1 0 0 1 0 0 0 1 0 1 0 0 1 0 0 0 0 1 1 0 0 0 1 0 0 1 0 1 0 0 0 1 0 0 0 1 1 0 0 0 0 1 0 1 0 1 0 0 0 0 1 0 0 1 1 0 0 0 0 0 1 1 0 1 0 0 0 0 0 1 0 1                            , β_2d=                    µ p1 p2 p3 p4 p5 p6 t1 t2                    とおきますと、 y = X2dβ2d+ ² (² ∼ N(0, σ 2_I 12)) (19)

と書けます。ここで、SASは制約条件p6= 0, t2= 0を入れます。そして、 X2d1=                            1 1 0 0 0 0 1 1 1 0 0 0 0 0 1 0 1 0 0 0 1 1 0 1 0 0 0 0 1 0 0 1 0 0 1 1 0 0 1 0 0 0 1 0 0 0 1 0 1 1 0 0 0 1 0 0 1 0 0 0 0 1 1 1 0 0 0 0 1 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0                            , β_2d1=               µ p1 p2 p3 p4 p5 t1               とおきますと、(19)は y = X2d1β2d1+ ² (² ∼ N(0, σ 2_I 12)) と書けます。これを用いて、β_2d1の最小二乗推定値を計算していきます。 β_2d1= (X0_2d1X2d1)−1X02d1y (20) となりますので、少しずつ計算していきます。まず、 X0_2d1X2d1=               1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 1 0 1 0 1 0 1 0 1 0                                          1 1 0 0 0 0 1 1 1 0 0 0 0 0 1 0 1 0 0 0 1 1 0 1 0 0 0 0 1 0 0 1 0 0 1 1 0 0 1 0 0 0 1 0 0 0 1 0 1 1 0 0 0 1 0 0 1 0 0 0 0 1 1 1 0 0 0 0 1 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0                            =               12 2 2 2 2 2 6 2 2 0 0 0 0 1 2 0 2 0 0 0 1 2 0 0 2 0 0 1 2 0 0 0 2 0 1 2 0 0 0 0 2 1 6 1 1 1 1 1 6              

です。これより、 (X0_2d1X2d1)−1=               12 2 2 2 2 2 6 2 2 0 0 0 0 1 2 0 2 0 0 0 1 2 0 0 2 0 0 1 2 0 0 0 2 0 1 2 0 0 0 0 2 1 6 1 1 1 1 1 6               −1 =               7 12 − 1 2 − 1 2 − 1 2 − 1 2 − 1 2 − 1 6 −1 2 1 1 2 1 2 1 2 1 2 0 −1 2 1 2 1 1 2 1 2 1 2 0 −1 2 1 2 1 2 1 1 2 1 2 0 −1 2 1 2 1 2 1 2 1 1 2 0 −1 2 1 2 1 2 1 2 1 2 1 0 −1 6 0 0 0 0 0 1 3               となります。次に、 X0_2d1y =               1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 1 0 1 0 1 0 1 0 1 0                                          10 11 12 20 18 17 18 18 15 10 20 13                            =               10 + 11 + 12 + 20 + 18 + 17 + 18 + 18 + 15 + 10 + 20 + 13 10 + 11 12 + 20 18 + 17 18 + 18 15 + 10 10 + 12 + 18 + 18 + 15 + 20               =               182 21 32 35 36 25 93              