「統計学関連なんでもあり」の過去ログ--- 048

No.22654　名義尺度の最小二乗フィッティングで追加制約がある場合の標準誤差の計算方法　　【小西】　2018/12/05(Wed) 18:18

長文失礼致します。
以下についてお分かりになる方がいましたらご教授頂けますでしょうか？

現在，交互作用のない2因子名義尺度に対する最小二乗フィッティングを行うプログラムを作成しているのですが，特殊な制約を加えた場合の標準誤差の計算方法がわからず困っております。
目標は2種類の加工装置の各個体が製品のばらつきにどの程度寄与しているかを推定するため，平均値(切片)，ばらつきがどの個体に起因しているか(回帰係数)，それら係数の信頼区間を求めることです。
例えば装置Aが3水準(A1,A2,A3)，装置Bが4水準(B1,B2,B3,B4)あるとして，AとBの組み合わせごとに製品の測定値が与えられます。ただし，いくつかの組み合わせが欠損値となることもあり得ます。

基本的にはJMPの計算結果を参考にしているのですが，測定値の偶然性や欠損値などによってデータが線形従属になった(解に自由度がある)場合，JMPでは一部の水準を除外(係数を0に固定)する計算方法を採用しています。
この方法では入力データの水準の並び順が変わっただけで除外される水準が変化してしまいますし，特定の水準を特別扱いしていることにもなるので，作成中のプログラムではこの点を改良したいと考えています。
※JMPにおけるランク落ちの対処法の詳細は
　http://www.jmp.com/japan/support/help/13/flm-sls.shtml
　の下部にある「1次従属性がある場合」のサブトピックで説明されています。

そこで，特定の水準を特別扱いせずに解を一意に決定する方法として，各個体の回帰係数(切片は除く)の2乗和を最小にする制約を設けることを考えました。
(理論的根拠をもとに考えたわけではないので，統計的に不適切だったりより良い方法があるかもしれません)
そして，等式制約付き最小二乗問題を解くことで，そのような回帰係数を求めることもできました。
しかし，信頼区間を求めるのに必要な標準誤差をこの制約下でどのように計算すればよいかがわかりません。※入力データによっては自由度が不足して標準誤差を求められないケースがあることは承知しております。
この場合の標準誤差はどのように計算すればよいでしょうか？(下記詳細もご参照ください)

■参考：JMPと同等な計算方法
・入力データ例(製品測定値と使用装置) ※線形従属になる一例

　11 A1 B1
　13 A1 B2
　12 A2 B1
　15 A2 B2
　14 A3 B3
　10 A3 B4

　1列目の測定値ベクトルをy，データ数をnとします。
　※この例では各組合せの測定値が高々1個になっていますが，実際には制限はありません。

・各因子水準をコード化(因子内の回帰係数の和は0に制約)

　　　　　A13　A23
　A1　→　　1　　0
　A2　→　　0　　1
　A3　→　 -1　 -1

　これを行列Aとします。

　　　　　B14　B24　B34
　B1　→　　1　　0　　0
　B2　→　　0　　1　　0
　B3　→　　0　　0　　1
　B4　→　 -1　 -1　 -1

　これを行列Bとします。

・コード化した入力データ

　切片　A13　A23　B14　B24　B34
　　1　　1　　0　　1　　0　　0
　　1　　1　　0　　0　　1　　0
　　1　　0　　1　　1　　0　　0
　　1　　0　　1　　0　　1　　0
　　1　 -1　 -1　　0　　0　　1
　　1　 -1　 -1　　0　 -1　 -1

・線形従属列の除外
　コード化した入力データの1列目，1〜2列目，1〜3列目，1〜4列目，1〜5列目，1〜6列目と順にランクを調べると6列目が増えたときにランク落ちしていることがわかるので，6列目を削除します。
　※より大きな入力データの場合は後続の列も同様に繰り返します。

　切片　A13　A23　B14　B24
　　1　　1　　0　　1　　0
　　1　　1　　0　　0　　1
　　1　　0　　1　　1　　0
　　1　　0　　1　　0　　1
　　1　 -1　 -1　　0　　0
　　1　 -1　 -1　　0　 -1

　これを入力データ行列X，Xのランクをrとします。

・削減した係数を復元する行列
　ブロック行列

　(1 0 0)
　(0 A 0)
　(0 0 B)

　から上記で除外した線形従属列(この例ではB34列)を同様に除外した行列をZとします。

・最小二乗問題Xa=yを解く
　解ベクトルaから切片,A1,A2,B1,B2が決まり，Zaですべての係数が求まります。
　平均二乗誤差:MSE=(y'y-a'X'Xa)/(n-r)
　標準誤差:SE=sqrt(MSE*diag(Z*inv(X'X)*Z'))
　※実際にはQR分解を利用してこれらと等価な計算を行います。

■詳細：検討中の2乗和を最小化する計算方法
※以下ではI_kはk次単位行列です。

・各因子水準のコード化(ここでは因子内の回帰係数の和を0に制約しません)
　A=I_3, B=I_4

・コード化した入力データ

　切片　A1　A2　A3　B1　B2　B3　B4
　　1　　1　 0　 0　 1　 0　 0　 0
　　1　　1　 0　 0　 0　 1　 0　 0
　　1　　0　 1　 0　 1　 0　 0　 0
　　1　　0　 1　 0　 0　 1　 0　 0
　　1　　0　 0　 1　 0　 0　 1　 0
　　1　　0　 0　 1　 0　 0　 0　 1

　これを入力データ行列X，Xのランクをrとします。

・等式制約付き最小二乗問題を解く
　Ca=0　　　　　　　※切片以外の二乗和を最小化
　s.t.　Da=0　　　　※因子内の回帰係数の和を0に制約
　　　　X'Xa=X'y　　※正規方程式で最小二乗誤差に制約

　ここでCは以下のブロック行列
　(0 I_7)

　Dは以下の行列
　(0 1 1 1 0 0 0 0)
　(0 0 0 0 1 1 1 1)

　また，以下のブロック行列をEとします。
　( D )
　(X'X)

　これにより解ベクトルaは求められますが，標準誤差SEの計算方法がわかりません。

　素人考えでJMPの計算式の形だけまねて，
　平均二乗誤差:MSE=(y'y-a'X'Xa)/(n-r)
　標準誤差:SE=sqrt(MSE*diag(pinv(X'X)))　※pinvは疑似逆行列
　としてみましたが，JMPで列の除外がない入力データの場合にJMPの結果と一致しませんでした。
　JMPの計算ではdiagの中の式(共分散行列？)にDa=0相当の制約条件が反映されているため，同様にDa=0相当およびEの零空間で二乗和を最小にする制約条件をdiagの中に反映する必要がある気がしています。