読者です 読者をやめる 読者になる 読者になる

StatsBeginner: 初学者の統計学習ノート

初学者が統計学、機械学習、R、Pythonの勉強の過程をメモっていくノート。

Rによる分散分析でタイプIII平方和を使う時の落とし穴

統計入門 分散分析

デフォルト設定ではダメ

 先日、RでタイプIII(タイプ3)平方和を使う方法についてエントリを書いた直後に、落とし穴があったことに気づいて、まとめたエントリを書こうと思ったんですが、勉強が進んでなくてあまりきちんと理解できておりません。しかし放置しとくのもあれなので、簡単にメモだけしておきます。


 前回のエントリでは、Rで分散分析をやるときにタイプⅢの平方和を用いたいのであれば{car}パッケージのAnova()を使うという話をしてたのですが、1つ注意しなければならないことがあります。Rのデフォルトの設定のまま使うと「間違った結果」になります!!
 結論としては、Rで

  • アンバランストなデザイン(セル間でサンプルサイズが異なる)
  • 交互作用あり
  • タイプIII平方和を使う

 という条件の分散分析を行うときは、まず最初に

options(contrasts = c("contr.sum", "contr.sum"))


 って書いとけということになります。気をつけることはこれだけです。

デザイン行列のつくりかたによって平方和が変わる

 まず、線形モデルに「デザイン行列(計画行列)」ってありますよね。変数を行列でまとめて、


{ \displaystyle
  \left( \begin{array}{ccc}
      y_1 \\
      \vdots \\
      y_n
    \end{array} \right)
=
  \left( \begin{array}{ccc}
      1        &  x_{11}  &  \ldots  &  x_{p1}  \\
      \vdots & \vdots   &  \ddots  &  \vdots  \\
      1        &  x_{1n}   &  \ldots  & x_{pn}
    \end{array} \right)
  \left( \begin{array}{ccc}
      \beta_0  \\
      \vdots \\
      \beta_p
    \end{array} \right)
+
  \left( \begin{array}{ccc}
      \epsilon_1  \\
      \vdots \\
      \epsilon_n
    \end{array} \right)
}


 というふうに書いたとして、これを、


{ \displaystyle
y = X \beta + \epsilon
}


としたとき、{ X }のところをデザイン行列といいます。解説は↓のページをみるといいと思います。要するに説明変数の係数じゃないとこを指してるわけですね。


デザイン行列について
デザイン行列(配置行列)


 分散分析の場合、線形モデルで表現すると説明変数はダミー変数になるわけですけど、そのダミー変数をどういうふうに作るかによって、平方和の計算結果が変わってきます。
 計算結果がどう変わるかという理論的な解説は、↓のPDFにとても詳しく書いてあります。

互作用に関する制約条件がTypeIIとTypeIIIの違いと考える.つまり,交互作用にセルの重みの違いを考慮した制約条件を課した結果がTypeIIであり,セルの重みを同じとした制約条件を課した結果がTypeIIIであると考える
(略)
バランスデータであればI,II,およびIIIは全て一致する.アンバランスだが要因が直交する場合は,IとIIは一致するが交互作用をモデルに含めた場合はIIIと一致しない.また,交互作用をモデルに含めない場合は,どのようなデータに対してでもIIとIIIは一致する.今回のように,アンバランスかつ要因が直交せず,モデルに交互作用まで含めた場合はI,IIおよびIIIは一般的には全て異なる.しかし,交互作用項A*Bの平方和については,交互作用まで含めたモデルの平方和と交互作用を含めないモデルの平方和の差分で定義されI,II,III全て一致する


http://www.sascom.jp/download/pdf/usergroups11_P-04.pdf


 あまり深くは考えてないのですが(だからブログに書くのを躊躇してたのですが)、結論としては要するに「交互作用項も含めて、タテに足したらゼロになるようなデザイン行列を作っとけよ!!」ってことですね。
 そうしないと、不適切な計算結果になるようです。

Rで計算する時

落とし穴に気づいたのは、自分がRで行った分散分析の結果と、他の人がSPSSで行った分散分析の結果に差があったからです。「そんな馬鹿な???」と思って調べていたら、↓のページに至りました。


myowelt: Obtaining the same ANOVA results in R as in SPSS - the difficulties with Type II and Type III sums of squares


 英語だし読むのだるいんですけど、options()関数の引数でcontrastsってのを設定しろと書いてある。
 これ、ANOVA君という分散分析の関数を公開している井関龍太氏のサイトにも分かりやすく解説が載っていました。

タイプⅢの平方和は,これらの,列を足して0になる計画行列を使った場合にのみ,適切に計算できます
他の計画行列を使った場合には,数値は計算できますが,おそらく,あまり意味のない値が出てきます。
つまり,タイプⅢ平方和では,使用する計画行列によって計算結果が変わるのです。
タイプⅠやタイプⅡでは,基本的に,このようなことはありません。


ANOVA君/平方和のタイプ/タイプⅢ平方和の計算法 - 井関龍太のページ


 要するにですね、分散分析をやりたくてタイプIIIの平方和を使うのであれば、まず始めに、

options(contrasts = c("contr.sum", "contr.sum"))


 って書いて、optionの内容を変更しておきましょうってことです。結論はこれだけ(笑)
 交互作用を考えないとか、バランストデザインになってるのであれば、結果的に気にしなくていいようですが。
 ちなみになんで2つ書くのかというと、help(options)の下の方を読むと書いてあるのですが、1個目の要素が順序なし因子に対する指定、2個目の要素が順序つき因子への指定になってるようです。

the default contrasts used in model fitting such as with aov or lm. A character vector of length two, the first giving the function to be used with unordered factors and the second the function to be used with ordered factors. By default the elements are named c("unordered", "ordered"), but the names are unused.


Rのhelp(options)より

 
 

2016.2.7追記:さらなる注意点

 Rで分析しているとき(あるいはスクリプトを書いてるとき)に、最初のほうで計画行列の指定をすると、以降もその設定が維持されてしまって、たとえばダミー変数を含む線形モデルでの回帰分析をlm()でやった場合のダミー変数の作り方にも影響しているので注意が必要です。
 2値のデータをダミー変数として投入する場合、0,1とする時が多いと思いますが、上述のような設定をしているとそうはならないので、算出される回帰係数の値が違ってきます。