昨年の医薬生物学の問3は行列を使って混合効果モデルを解いていく問題でした。式が意味するところとその構造を理解していれば、意外と解きやすい問題だったのではないかと思います（当日は完全スルーしましたし、当日は多分解けなかったと思いますが）。事前知識としてわかっておいた方が良い内容は以下の記事にまとめましたのでご参考ください。
平均ベクトル、分散共分散行列、ベクトルの線形変換、ガウス型線形回帰モデル【統計検定1級対策】
では見ていきましょう。
(1)
まずは混合効果モデルの分散共分散行列を求める問題です。
これはベクトルを線形変換したときに、分散共分散がどう変化するかがわかっていれば良い問題です。細かい式の変換については冒頭の記事にまとめてあるのでそちらを参照してください。さて、まず回帰式は
\(Y=X\beta+Z\gamma+\epsilon\)
ですが、Xβは確率変数ではないので定数として無視できます。
残りについてみていきますと、εはγと互いに独立であるため分散の線形性から以下のように分離できます。
\(cov(Y)=cov(Z\gamma)+cov(\epsilon)\\=ZGZ^T+R\)
あとは線形変換した場合の変形を施せば、上記のように導けます。
(2-1)
続いてはβの推定を行なっていく問題です。ただここでは変量効果をなしとしていますので、実質的には交互作用項のない単回帰分析（仮にベクトルβの行が増えれば重回帰分析）になっています。
まずは推定量の式を簡略化しろという問題です。最尤推定の式は問題文ですでに与えられているため、当てはめて計算していけば良いです。(1)の答えより
\(V=ZGZ^T+R\)
ですが、Z=0なのでVの推定量は
\(\hat V=\hat R=\hat{\sigma^2}I\)
となります。よってVの逆行列は
\(\hat V^{-1}=\frac{1}{\hat{\sigma^2}}I\)
ですので求めたい式は
\((X^TV^{-1}X)^{-1}X^TV^{-1}Y\\=(X^T\frac{1}{\hat{\sigma^2}}IX)^TX^{-1}\frac{1}{\hat{\sigma^2}}IY\\=(X^TX)^{-1}X^TY\)
となります。
ここでXで作られた行列は何だか複雑な形をしているようにみえるのですが、疑似逆行列（または一般逆行列）と呼ばれるものと同じ形をしていることが分かります。要するに右からXをかければ、単位行列になるというものです。
サンプル数と説明変数の数が合致するということはまずないので、Xは基本的に正方行列ではありませんが、\(X^TX\)とすることで正方行列になるのでこの部分の逆行列は求められるわけですね。
(2-2)
続いては先ほどの式から推定量の期待値を求めます。
ここで、計画行列Xは元々決まっている定数ですので期待値の外側に出せることに注意して式を式を変形します。
\(E[\hat\beta]\\=E[(X^TX)^{-1}X^TY]\\=(X^TX)^{-1}X^TE[Y]\\=(X^TX)^{-1}X^TE[X\beta+\epsilon]\\=(X^TX)^{-1}X^TXE[\beta]\\=\beta\)
後半の式変形ではεの期待値が0であることを用いました。計算結果から(2-1)で求めた最尤推定量が不偏性を持つことがわかります。
(2-3)
さて続いては実際の推定量の計算問題です。ここまでの式の流れに沿って行列を使って解くことができます。
問題文に記載がある通り、Xはプラセボが０、実薬を１とする計画行列なので、以下のようになります。
よって、まず\((X^TX)^{-1}\)を求めると（普通に2×2の逆行列を求めればよいので計算は省略）
となります。
さらに計算を進めて\((X^TX)^{-1}X^T\)は
となります。
あとはYのベクトルをかけて計算すると
μ=20.16666666…
β=-8.45
となります。
なんだか最後にYにかける行列が単純な構造すぎるようにみえるのですが、考えてみるとμは回帰直線を引いた際の切片になるので、ダミー変数が０であるプラセボ群の平均値を出しているだけですし、βは傾きなのでプラセボ群と実薬群の差の傾きを出せばよいだけであって、今回の問題が単回帰として非常に単純な構造をしているからだと分かります。
（別解）
やっていることは同じですが、単回帰であるため、行列を使わず各パラメータを微分することで最小二乗法で求めても同じことができます。計算過程は端折りますが、βの最尤推定値は
\(\hat\beta=\frac{\sum(x_i-\bar x)\sum(y_i-\bar y)}{\sum(x_i-\bar x)^2}\)
となります。なお\(\bar x=\frac{1}{n}\sum x_i, \bar y=\frac{1}{n}\sum y_i\)とします。
これを計算しても上記と同じ結果が得られます。
(3-1)
続いてはちゃんと変量効果を考えようというモデルの場合どうすればよいかという話です。
与えられた条件から行列Rは
なので\(ZGZ^T\)は
（以下対角にそれぞれの分散が3×3で並んだ区分行列で形成される行列）となれば良いことがわかります。Zが12×q次元ですので、この行列は12×12次元の行列です。
さてそうなるとGを以下のようにしてみると
Zはそれを行方向に３つ、列方向に３つ、成分を拡大するようにすればよいので
とおけます(3つずつ1が並び、列をずれていく)。
この辺はもうちょっと何かうまいこと説明できるような定理なりあるのかもしれませんが、私にはわかりませんでした、すみません。
(3-2)
記述問題です。ここまでで、説明してきたように[2]は変量効果を考慮しない、ただの単回帰モデルでした。明らかにこのモデルが問題なのは投与前と１週目６週目といういかにも相関がありそうなデータが全く無視されており、すべて同じものとして扱われている点です。個体間での違いが考慮されていないということですね。
これに対して変量効果モデルではその分散共分散行列をみれば分かるように、個々の個体内での分散と、投与前１週目６週目の共分散が存在しており、それぞれの相関が考慮されています。
ということでこれをまとめると公式の解答と同様になるのではないかと思われます。
非常に簡単なモデルで行列を用いた回帰モデルを勉強できる良問でした。7，8年前の問題と比べると行列を用いた問題も難しい変形や知識を迫られるというより、モデルへの理解を問う問題が多い印象です。計算を複雑にあまりできないこういった行列の問題は、うまくいけば狙いたいですね。
統計学をどんどん自分で学んで深めたい、という方へのおすすめ書籍をこちらのページの下部にまとめております。初心者向けのものから応用まで幅広く読んでいますので参考にどうぞ。
コメントを残す