※上記の広告は60日以上更新のないWIKIに表示されています。更新することで広告が下部へ移動します。

Ⅱ型の回帰分析

2006年12月1日からこれまでの訪問者数
-
ver 2.0に合わせて修正しました(2006/12/1)


RMA回帰って?

RMA回帰とは

Reduced Major Axisを用いた回帰。SMA (Standardised Major Axis) あるいはGM (Geometric Mean) とも呼ばれる (Sokal and Rohlf, 1995)。RMAかSMAで表記されることが多いが、本稿では便宜的にRMAで表記する。

2変量の関係を検討する

→補間ではなく、回帰がよく用いられる。実際に計測されるデータはいろんな誤差を含んでいるから。
回帰:線形・非線形が存在。
線形回帰:数学的に最適解を算出できる。
非線形回帰:数学的な最適解は存在しない。

線形回帰:I型とII型が存在

I型:X軸には誤差がないことを仮定。XからYを予測することを目的とする。
II型:X軸とY軸に誤差があることを仮定。XとYの関係を検討することを目的とする。
I型とII型で目的は異なる。XとYの関係の解析をしたい時はII型回帰のほうが適切であることは以前から指摘。
→何故か?I型回帰はXからYを予測することが目的であり、その傾きは (xとyの積和) / (xの平方和) で定義されている。そのため、X軸側に誤差が含まれる場合、実際よりも傾きが緩く推定されるから。

II型の回帰

MA (Major Axis)
RMA (or SMA or GM)
Bartlett's three group Method
などがある。詳しい内容を知りたい方はSokal and Rohlf (1995)を参照されたい。

MAとRMAの違い

MA:XとYの誤差を等価に考える。具体的には、各測定点と回帰直線の距離が最短になるような回帰直線をMAとする。すなわち、MA回帰直線に対し垂直に引かれた各測定点への直線の距離の和が最小になるようにする。MAの求め方に関しては<粕谷 (1997) の回帰のところに詳細な解説がある。

RMA:XとYの誤差を等価に考える点ではMAと同じだが、MAは各軸の単位が変わると別の回帰直線になってしまう。そのため、まずデータ (x, y) の2つの変数に関し、平均0、標準偏差1に変換する (Reduced)。このデータから傾きを生成し、元のデータに当てはめる。実際にはSMAの傾きは、変換前のデータを使って、
によって簡単に算出できる。

ちなみに、MAもRMAも (xの平均, yの平均) を必ず通るから、傾きが求まれば切片も計算できる。

これらの回帰直線を処理間で比較する方法

I型回帰:ANCOVA (共分散分析)
II型回帰 (アロメトリーなどxとyに明確な独立-従属関係が考えられない場合、xとyの関係自体に興味がある場合はより適切):???

→近年までII型回帰の関数関係を比較する方法は存在しなかった。しかし、近年比較的合意が得られている方法として、RMA回帰による比較がある。

RMA回帰における回帰直線間の比較

RMA回帰間の関係を検討する上で必要な仮定

  • データの直線性。散布図を描いてみて直線的でないなら何らかの変数変換を考える必要がある。
  • 等分散性
  • 正規性←満たされなくてもロバストであるとはされているが...

SMATRの頭の中をのぞく -どのような計算をしているのか?-


傾きの検定: 「尤度比検定」または「ならべかえ検定」のどちらか(ユーザーが選択する。デフォルトはならべかえ検定」。尤度比検定は、2つのモデルに有意差があるかどうかを、2つのモデルの対数尤度の差を-2倍した値を検定統計量とし、自由度 = (処理数) - 1のχ2分布を用いて検定する。1つのモデルは確率密度が表現できるモデルであればどのようなモデルでもよいが、2つ目は1つ目のモデルに制約を加えたものである必要がある。

しかし、RMA回帰の共通の傾きの検定では、この尤度比検定で検定すると第1種の過失を過大評価する傾向がある (らしい) ので、検定統計量をならべかえにより何度も生成し、以下の式から有意確率を算出するならべかえ検定を用意し、かつこちらをデフォルトにしている。
nはならべかえ回数、recalculatedはならべかえにより算出された検定統計量、observedはもともとのデータセットを用いて算出された検定統計量。これによって、pが0.05より小さければ処理間で傾きに差がある、と判定する。

切片の検定:通常のANOVAで検定する。ただし、各データ (x, y) を以下の要領により、(x',y') に変換し、xとyの相関を取り除いてから検定を行う。

y' = y - bx, x' = y + bx; bは共通の傾き。

SMATRの使い方

RMA回帰を処理間で比較したい貴方はコチラ!


日本でもSMATRを用いた論文が出始めている→飯島ら, 2004; Shibuya et al. 2005; Cho et al. 2005など。

特徴
  • フリーソフト
  • Windows、Mac、Linux版がある
  • 軽いのでほとんどの環境で利用できる (と思う)
  • txt形式(タブ区切り)で保存したファイルを読み込む
  • 元のファイルで一番上の列にラベルをつけ、それをソフト上で認識させることができる。
  • ソフト上で対数 (常用対数or自然対数) 変換ができる
  • グラフィック機能は一切持っていない (別に絶対必要なものではないんだが)

ソフトの入手方法

  1. http://www.bio.mq.edu.au/ecology/SMATRにアクセスし、SMATR v2.0 (smatr.exe) と SMATR user's guide (SMATR_instructions.pdf) をダウンロードする。
  2. ダウンロードしたソフトを、新しく適当に作ったフォルダの中に入れる。
  3. 以上。このソフトは単体で機能する。また、レジストリへの書き込みを行わないので、削除するときはそのままゴミ箱に入れればよい。

解析データの準備方法

Excelの例
  1. Excel上で以下のようにデータを打ち込む。各列にデータが縦に並び、1行目はデータのラベル (そのデータの名前)、2行目以降にデータを入れる(下の図参照)。
  2. 用意したデータシートにinputという名前をつけ、タブ (tab) 形式で保存する。
  3. 保存したinput.txtファイルを、先ほどのSMATRソフト本体を入れたフォルダに移動する (重要! SMATRはその本体があるフォルダ内のファイルしか認識しない)。

データ入力例

解析方法

  1. SMATRが立ち上がると、まず図1のような画面が現れる。
図1
この状態では、またデータは読み込まれていない。このソフトは、画面の指示に従いコマンド (lとかxとかキーボードの1文字が多い) を入力し、Enterすることで動作する。

最初の画面では、RMA回帰の関数関係を検討する上で必要な事項の設定を行う。
  • I - Input File: 読み込むデータファイルの名前を指定するコマンド。デフォルトはinput.txt。先ほどファイルネームをinputとしたのはこのため。必要があれば変更する。
  • O - Output File: 結果を出力するファイルネームを指定するコマンド。デフォルトはoutput.txt。(S)MATRのあるフォルダに勝手に生成される。
  • M - Method: 統計方法。デフォルトはSMA。ほかには通常の回帰(Ordinary Least Squares Regression)、主軸(Major axis)が用意されている。あえてこれらの回帰の方法の比較をしたいのでなければ変更する必要はありません。
  • T - inTercept: RMA回帰を行うときに、原点を強制的に通すかどうかの設定する。通常は必要なし。
  • E - measurement Error: 測定に伴う誤差を(RMAでも何でも)回帰に反映させるためのコマンド。これを計算するためには、同じデータを繰り返し取ってみて、その繰り返し取られたデータを用いて測定そのものにどのぐらいの誤差が含まれるか検討する。通常同じデータを何度も取れる状況にはないので、重要なことだが使うことは少ないかも。
  • C - Confidence int: 信頼区間。デフォルトは95%。
  • G - min Group size: RMA回帰時に、最低必要なデータ数。デフォルトは6。
  • R - Resample: 傾きの検定をするときに、デフォルトのならべかえ検定するか、Χ2分布を用いた尤度比検定を行うか選択。通常はならべかえ検定。
  • N - No. Iterations: 傾きの検定時に行われる並べ替えの繰り返し数。デフォルトは1000。p < 0.01以下の確率で並べ替えを行いたいときは5000とする。
  • P - Critical p-value: 有意とする確率。デフォルトは0.05。

変更する場合は、変更したいコマンドを入力しEnterする。いじるとしたらi、o、nぐらいかな?

設定が終了したら、lを入力しEnterする。すると図2のような画面に移る。

図2
この画面では、読み込んだデータの認識方法に関する設定を行う。
  • G - Grouping variable: グループ変数。いわゆる目的する処理を示している列。Gと入力してEnterすると、先ほどExcelで入力したデータファイルの1行目にあるラベルが並び、どのラベルがグループ変数なのかを聞いてくる。ラベルは数字で参照されているので、目的のラベルを示している数字を入力し、Enterする。
  • Y - y variables: Y軸にくる変数。選択法はGと同じ。
  • X - x variables: X軸にくる変数。選択法はGと同じ
  • F - Filter variable: フィルター変数。
  • T - Transformation: データ変換。デフォルトはX-Yが直線-直線 (つまり無変換)。XとYに関し、無変換のままか、対数変換するかが選択できる。わざわざ対数変換したデータを読み込ませる必要がない。これはこのソフトの優秀なところ。
  • A - H0 slope: あるひとつのRMA回帰に関し、得られる傾きがここで設定する値と有意に異なるかどうかを検定できる。2つ以上のRMA回帰間で関数関係を検討する場合はまったく無意味な設定。
  • B - H0 intercept: あるひとつのRMA回帰に関し、得られる切片がここで設定する値と有意に異なるかどうかを検定できる。2つ以上のRMA回帰間で関数関係を検討する場合はまったく無意味な設定。
  • D - Display data: 現在読み込まれているデータの値を示してくれる。それだけ。

ここでユーザは最低でもg、y、xを設定する必要がある。また、アロメトリーの場合によく用いられる対数変換をしたい場合は、tで設定すればわざわざ対数変換したデータを用意する必要がない。

例えば比較したいグループを示す、グループ変数を指定してみよう。gを押してEnterすると次の図3のような画面になる。

図3
すると、先ほど作ったinput.txtファイルで、1行目につけたラベルが表示されている。そして、各ラベルに数字がついている。smatrにどれがグループ変数か教えるには、グループ変数が書かれたラベルの数字を打ち込み、Enterする。すると図4のように、グループ変数がsmatrに認識される。

図4
この作業を、xとyに関しても行う。

ちなみに対数変換(自然、常用どちらでも)をしたい場合、tを押してEnterする。すると図5のようになる。

図5
0が常用対数、1が自然対数、2が変換しない場合である。ver 1.0であった、xかy片方だけ対数変換する項目はなくなった。実際にはほとんど使われないからだろう。

こうしてデータ認識に関する設定が終わったら、lを入力しEnterする。するとRMA回帰が各処理で生成され、処理間で傾きの検定が行われる。すると処理間で傾きの差が

有意な場合:図6

有意でない場合:図7

のような画面に移る (注:当然ですがここからの画面の数字の値はデータによって変わります)。

図6
この画面で見る必要があるのは、下から2番目のTest statistic: p = 0.001。ここで処理間の傾きに有意な違いがあれば0.05以下となり、有意な差が見られなければ0.05以上となる。

図6の場合、処理間で傾きに有意差が見られたので、どことどこの処理間に差があるのかを多重比較する。一番下の行は、処理間で傾きの多重比較を行うかどうか聞いてきているので、yを入力しEnterする。

すると処理間で傾きの多重比較が行われる。検定はここで終了。次に出てきた画面でxを入力してEnterすると終了。結果はoutput.txtに示されている。

要注意!!
(S)MATRでは、傾きの多重比較に関し、単なる対比較による有意確率しか示していない。そのため、示されている有意確率を、Bonferroni法やHolm法 (これが最も一般的。Bonferroni法は保守的すぎるのでHolmの方法が○) で調整し、family-wiseな有意確率を0.05にする必要がある。この辺の多重比較法に関しては永田・吉田 (1997) を参照。

図7
処理間で傾きに有意差が見られなかったので、共通の傾きを生成し、切片の検定に移る。共通の傾きを当てはめたときに、処理間で共通の傾きの垂直方向 (Y) と水平方向 (X) へのずれが生じているかの検定を行いますか?、と表示される。この場合はyを入力してEnterする。すると図8のような画面に移る。

図8
この画面では、切片の検定のために変換されたデータの保存、処理間で共通の傾きに対する垂直方向への差の検定、共通の傾きの水平方向での差の検定のいずれを行うか選択する。

  • S - save transformed data to disk: 処理間で共通の傾きを当てはめることが可能な場合、切片の検定を行うが、切片の検定に際し、
y' = y - bx, x' = y + bx; bは共通の傾き
によってデータを (x', y') に変換する。この変換したデータを出力するためのコマンド。
  • R - compare Residuals for shifts in elevation: 共通の傾きと垂直の方向に処理間でy'の平均値に差が見られるかを検定。
  • F - compare Fitted values for shifts along the common slope: 共通の傾きの方向に沿って処理間でx'の平均値に差が見られるかを検定。
  • E - Exit: 終了。

Rを入力しEnterした場合は、通常のANOVAにより、共通の傾きに対して垂直方向へのずれが生じているかを検定し、分散分析表が出力される(図9)。

図9
pが0.05より小さければ水平方向へのずれが有意に生じている。分散分析表の下に多重比較を行いますか?と表示されるので、pが0.05より小さければyを入力しEnterする。

すると多重比較が行われて検定は終了し、図8に戻る。

SMATR上では示されなかった結果もoutput.txtに保存されている。output.txtでは、SMATR上での設定、各処理へのRMA回帰の結果、行った統計処理の結果が順に示されている。

参考文献

  • 粕谷英一. 1998. 生物学を学ぶ人のための統計の話~きみにも出せる有意差~. 文一総合出版. 通称ピンク本 (表紙がピンク色だから).
  • Sokal, R.R., and Rohlf, E.J. 1995. Biometry: the principles and practice of statistics in biological research, 3rd edition. W.H.Freeman, San Francisco.
  • (S)MATRの作成者の一人であるWarton, I.D.さんのほーむぺーじ:http://web.maths.unsw.edu.au/~dwarton/

引用文献

  • 永田靖・吉田道弘. 1997. 統計的多重比較法の基礎. サイエンティスト社.
  • 飯島勇人・渋谷正人・斎藤秀之・高橋邦秀. 2004. 倒木上のコケの高さがエゾマツ実生の生残と成長に与える影響. 日林誌, 86 (4): 358-364.
  • Shibuya, M., Hasaba, H., Yajima, T., and Takahashi, K. 2005. Effect of thinning on allometry and needle-age distribution of trees in natural Abies stands of northern Japan. J. For. Res. 10 (1); 15-20.
  • Cho, M., Kawamura, K., and Takeda, H. 2005. Sapling architecture and growth in the co-occurring species Castanopsis cuspidata and Quercus glauca in a secondary forest in western Japan. J. For. Res. 10 (2): 143-150.