2015年11月14日土曜日

回帰係数0とパラメータ0は検定が違う:黒木氏の非線形問題


偶然,ネット上で,以下のような非線形回帰の係数の統計学的検定に関する質問を見つけた。

非線形回帰の信頼性評価について(黒木)

その質問の中で,特に,気になったスレッドが,No.01788 2006/12/06(Wed) 18:09。統計解析ソフト R を使った非線形回帰分析の方法が回答された後,質問者・黒木氏が,仮説を"回帰係数(パラメータ)の値が0であるかどうか"としか置けないところが難点です,述べたことである。

黒木氏は,回帰係数(regression coefficient)とパラメータ(parameter)を同一視してしまっている,あるいは,混同している。しかも,それゆえ,値が0であるかどうかとしか置けない,と誤解してしまっている。両者は,もちろん異なるものであり,回帰分析で任意の実数の回帰係数の検定可能である

実は,統計ソフト R の非線形回帰分析などでは,厳密に言うと,係数の検定を行なっているのではなく,パラメータの検定を行なっているのである。

例えば, ab を定数として,以下のような直線を考える。

回帰直線

このとき, a は傾きを表す係数であり,かつ,パラメータである。統計ソフトでは,しばしば
帰無仮説 H0: a = 0
に対して検定が行なわれる。

しかし,ここで,p を定数として,例えば,

パラメータの一次変換

変換された回帰直線

としてみよう。

すると,パラメータは p であるが,回帰係数は p+1 である。前述のとおり, R の非線形回帰分析などでは,パラメータの検定を行なうので,
帰無仮説 H0: p = 0
に対して検定が行なわれる。

つまり,これは,回帰係数に対して
帰無仮説 H0: a = 1
という検定が行なわれることなのである。

実際に,以下のような xy データを R で検定してみよう。

x<- c(1, 2, 3, 4, 5)
y<- c(0, 3, 2, 6, 5)

line1<- lm(y~x)
(reg1<- summary(line1))  # 回帰分析

plot(x, y)
abline(line1, col="red")

まず,ここまでで,データの散布図(Scatter Diagram)と回帰直線が描ける。

分析データの散布図と回帰直線

回帰係数の検定理論は,私のホームページ
その中の6. 回帰の検定に書かれている。

例えば,傾きを a,傾きの標準誤差を Sa,としたとき,傾きが2となるか否か,つまり,
帰無仮説 H0: a = 2
に対して検定する場合,検定統計量を t とすると,以下のように表される。

回帰係数2に対しての検定統計量 t

もちろん
帰無仮説 H0: a = 0
に対して検定する場合は,

回帰係数0に対しての検定統計量 t

すなわち

回帰係数0に対しての検定統計量 t

これが,黒木氏が,仮説を"回帰係数(パラメータ)の値が0であるかどうかとしか置けないところが難点として参照したウェブページ偏回帰係数の検定が表す式である。

つまり,回帰係数の値が0であるかどうかが検定できれば,一般性を失うことなく任意の係数について検定できるのである。これは, t 統計量が母平均の差の検定に用いられるように,差の検定に使われる統計量だということからも推察される。このような理論的な部分を,少なくとも,このスレッドを書いた時点では,黒木氏は理解していなかったようである。

話を元に戻して,この式を使い,
帰無仮説 H0: a = 2
に対して, R で検定してみる。傾きとその標準誤差は,回帰分析の結果で,要素 coefficients の2行1列,2行2列のデータとして格納されているので,それを取り出して使う。

a<- reg1$coefficients[2, 1]  # 傾き
sa<- reg1$coefficients[2, 2]  # 傾きの標準誤差

df<- length(x)-2  # 自由度
(tval<- abs(a-2)/sa)  # t 値
pt(tval, df, lower=F)*2  # p 値

結果は以下の通り。

t = 1.578457
p = 0.212573


しかし,これは,前述のとおり,

変換された回帰直線

としてパラメータを検定すれば簡単である。

そのときは, nls 関数を用いる。これまた,非線形回帰のときに使われる関数だと,安易に解説している人もいるが,それだけではないのである。以下が,それを用いた検定のスクリプト。

line2<- nls(y~ (p+2)*x + b, start=c(p=1 ,b=1))
(reg2<- summary(line2))  # 回帰分析

reg2$coefficients[1, 3] # t 値
reg2$coefficients[1, 4] # p 値

実に簡単である。結果は以下の通り。

t = 1.578457
p = 0.212573


この結果は,当然だが,前述のように定義に従って計算した結果と同じである。

ここでは直線回帰(linear regression)を取り上げたが, nls 関数を使っているので,もちろん,非線形回帰(non-linear regression)において,係数が0でない場合の検定も出来る。

余談であるが,統計学における母数とは,パラメータ(parameter)のことであり,母平均とか母分散と言う時の,である。その意味で, R は係数の検定ではなく,母数の検定を行なっているのである。

統計学における母数の意味を知らない人(大学教員もか?)がかなり多い。朝日新聞の統計記事でも間違えて使われている。私の Yahoo!知恵ノート参照。

これに関連して,以下のブログも面白かった。

私の研究室トップページに,その他の統計学の話題のリスト。