杉野 勇 (SUGINO Isamu)　お茶の水女子大学・人間発達科学専攻・応用社会学コース担当

　本文では紙幅の関係で式変形の途中などを省略したので，ここでもう少し詳しく説明するPDFを上げておこう。以下の画像の上でクリックするとPDFへのリンクが開く。

　ここでも本文では紙幅の関係で式変形の途中などを省略したので，もう少し詳しく説明するPDFを上げておく。


　　Rでは，自由度df1，df2のF分布の右端5%領域の限界値は，qf(p=.95, df1, df2)もしくはqf(p=.05, df1, df2, lower.tail=F)として求められる。グラフは，curve(df(df1, df2), from=0, to=4) などとすると描ける。以下は，グループの数mが3，総ケース数nが30を想定した数値例である。


　グラフにいくつかオプションなどで装飾を施してみよう。

　まずは使用する変数の基本的な情報（度数分布や要約統計量，ケース数）を確認しておく必要がある。欠損値NAが含まれる場合はひと手間余計にかかる事が多いので注意しよう。また，本文と異なり，データフレイムの名前はdata01としている。



　本文の箱ひげ図は次のスクリプトで描いている。結果は本文にある通りである。

　いよいよ（一要因）分散分析を行う。aov( )関数，lm( )関数，oneway.test( )関数のうち，まずは最初の二つを説明する。

■aov( )関数



　出力結果を用いて，「グループ平均はすべて等しい」というゼロ仮説の有意性検定をF分布を用いて行う事が出来る。この例ではF統計量は約3.05，有意確率は約2.9%となり，5%水準で有意となった。


■lm( )関数



　次のスクリプトを実行すると，このF分布の5%棄却域（赤）とF値（紫）の関係をカラーで図示できる。関数ごとに分けて実行すると関数の意味が分かりやすい。


　なお，各種の分散説明率（R-squared, Eta-squared, partial Eta-squared）については，二元配置分散分析の2-1 平方和の分解の末尾に例示を付け加えておく。

　一つ一つの検定の有意水準をp1，そうした検定を3つ繰り返した場合の最終的な有意水準をp3とすると，p3 = p1 + p1 + p1 - 3*p1*p1 + p1*p1*p1 である。p1を.05（5%）とすると，p3=.142625となる。そこでp1=.05/3とすると，p3は5%以内に収まる。


　世帯年収fincomeと従業上の地位jobで実際にボンフェローニ法を実行してみよう。共通の標準偏差を用いるか否かというオプションpool.sdによって結果が変わる。本文ではデフォルトのpool.sd=Tの結果を述べたが，ここでは両方実行する。

　デフォルトはプールした標準偏差を使う，つまり全ての群で等分散を仮定しているが，pool.sd=Fとすると各群での標準偏差を使用する（等分散を仮定しない）。等分散を仮定しないt検定とはつまりWelchの検定であり，上の出力の後半は，各ペアでウェルチの検定を行って計算された有意確率に検定の回数（ここでは6回）をかけた結果が出力されている。例えばjob=1とjob=2だけでWelchの検定を行うと有意確率は.005345494となる。これに全ての検定のペアの数6をかけると，.03207297となり，上記のボンフェローニ法の結果の1と2の有意確率0.032に一致する。


　本文で触れた，無調整の検定結果を一応示しておく。

　指定しうる多重比較の方法は，ボンフェローニ"bonferroni"，ホルム"holm", ホッホベルク"hochberg", ホンメル"hommel", ベンジャミニとホッホベルク"BH", ベンジャミニとイェクティエリ"BY"である。bonferroni，holm，hochberg，hommelは"family-wise error rate(FWER)"を調整するものであるが，BHとBYは"false discovery rate(FDR)"を統制するものである。以下のサイトの解説が分かり易い。


　hommel・等分散の例と，BH・不等分散の実行例を示しておく。


　Tukey法は，以下の最初の二行で実行し信頼区間を図示できるが，より見易くする為にオブジェクトの代入を活用したのがその後のスクリプトである。

　次のスクリプトの一行目がデフォルトの（等分散性を仮定しない）一要因分散分析，二行目が等分散の仮定を指定した一要因分散分析である。後者の結果は，先に紹介したaov( )やlm( )と一致している事を見比べよう。


　本文では簡単に触れているだけの等分散性の検定を実際に行ってみる。標準で存在しているバートレット検定bartlett.test( )と，追加パッケイジ"car"に入っているルヴィーン検定（レヴィン検定）の二つを行ってみよう。"car"パッケイジは，インストールした事がなければ最初にインストールし，その後も使用する際にはlibrary(car)として有効化しなければならない。また，bartlettはグループ変数がfactor型になっていなくても自動で要因型だと見做して分析するが，LeveneTestはfactor型にしなければならない。
　出力結果ではインストールの部分は省略した。

　二要因分散分析でも，（偏差）平方和SSの分解が基本であることには変わらないが，各グループのサイズが異なる場合（非釣り合い型データ，アンバランスデータ，非直交なデータ）には分解の仕方が複数存在し，標準のaov( )関数はその点で問題を生じる。まずはデータの準備をする。分割表でも分かる様に，セルのサイズ（グループのサイズ）はまるでばらばらである。


　次に，それぞれの一要因分散分析と，性別，従業上の地位の順に変数を並べた二要因分散分析の結果を表示する。一要因の場合とは，各独立変数（主効果）の平方和も有意確率も異なる。


　ここで，独立変数の並べ順を変えてみよう。変数の並べ順を変えただけで，それぞれの独立変数の偏差平方和も有意確率もすっかり変わってしまった


　ここで，"car"パッケイジのAnova( )関数を利用する。パッケイジのインストールには上の1-5を参照。model32とmodel33は既に上で作成した分散分析のモデルである。この，投入順を変えた二つのモデルで，Anovaのタイプ?平方和の計算による分散分析を行った結果は，平方和も有意確率も全て一致していて，投入順に左右されない事が確認出来る。因みに，Anova( )関数でタイプ?の平方和を用いる事も出来る（type=3というオプションを付ける）。

model32，model33自体は以下の通りであり，Anova( )関数を使わなければtype?の計算は行われない。


　各グループのケース数や年収の平均といった基本的な情報を求めるには以下の様にすれば良い。



　ここで，効果サイズ(effect size)の指標である分散説明率について紹介しておこう。
分散分析では伝統的にη²（イータの二乗）と云う単語・統計量が使用され，回帰分析ではR²（決定係数）が用いられてきた。R²は通常，モデル全体の説明率を示すのに用いられ，η²は個々の要因の説明率を示すのに用いられる。（重回帰分析では偏回帰係数のt検定は分散説明率の増加分についてのF検定と同値であるので，実際には回帰分析でも個々の要因の説明率を示すのにも使用可能である。）
　世帯収入を，性別と従業上の地位で説明するモデルを，aov( )とlm( )で実行し，anova( )やAnova( )，heplotsパッケイジ中のetasq( )関数に与えて出力結果を確認する。


　以下に，偏イータの二乗とイータの二乗の計算の仕方を実行例とともに示す。解読してみよう。定義上必ず偏イータの二乗の方が大きくなる。

　複数の要因（独立変数）が複合的に従属変数に影響を及ぼす事を，計量分析では「交互作用intraction」と呼ぶが，質的比較分析(QCA, Qualitative Comparative Analysis)でいう「結合因果conjunctural causation」と同種の考えであると言ってよいだろう。
以下，本文のスクリプトとその結果を，より装飾を施した形で紹介してゆく。



平方和の分解方法によって結果が異なる事を確認しておこう。
　最後に，交互作用プロットの描き方を紹介する。本文のと類似のグラフが以下のスクリプトで描ける。

　性別による幸福度の平均の差の分析結果を再掲しておこう。


　これを，一要因分散分析で行ってみる。この例では等分散を仮定出来る為，関数は3種類あるので全て実行して見せよう。独立変数であるx（元はdata01$sex）はfactor型にすべきだが，二値変数の場合には結果に本質的な違いはない。

　t検定の結果とF検定の結果をそれぞれオブジェクトに格納し，並べて表示させよう。更に，オブジェクトにどんな情報が格納されているのかも一覧を確認しよう。

　t分布の自由度とF分布の残差（級内）の自由度，有意確率同士が対応している事が分かる。
t統計量とF統計量の関係は，次の様に確認してみよう。

　t分布はマイナス領域とプラス領域がある。二乗すれば同じ値になる点が正負の二つある。F分布は全て非負である。どの様に対応しているだろうか。下の様に確認すると，tが負の時と正の時の確率を合わせたものが，（tの二乗に等しい）Fの上側確率に等しい事が分かる。tの両側検定はFの上側検定に対応しているのである。


　等分散を仮定出来ずにウェルチの検定を行った，男女での世帯年収差はどうであろうか。

これもやはり，自由度，有意確率が一致し，t統計量の2乗とF統計量が一致している。anova(aov( ))関数では等分散性が仮定出来ない分散分析を行う事は出来ない。

1) の答え： 問題文で要求されていないコマンドも含んだスクリプトを，二つに分けて結果と共に掲載する。



2) の答え： まず，anova(aov( ))では独立変数の投入順が結果に影響する事を確認しておく（タイプ?平方和）。
　"car"パッケイジのAnova( )関数（タイプ?平方和）ではそうした違いが生じない事を次に確認し，交互作用を含んだモデルも実行する。そして最後に交互作用プロットを描く。




3) の答え： 交互作用含むモデルを取り上げ，aov( )，anova( )，Anova( )の結果も比較してみる。


　model1の"coefficients"を表示させると以下の通りである。

　カテゴリごとの平均は以下の通りである。

これをよく見比べると，model1$coefficients の(Intercept)は，高校―男性カテゴリの平均個人年収に一致している。また，「x2女性」の-300.9495465は，高校―男性と高校―女性の差にぴったり一致している。「x短大」の61.3189720は高校―男性と短大―男性の差にぴったり一致する。最終学歴においては「高校」が，性別においては「男性」が「基準カテゴリ」となっており，そこからの差が"coefficients"として計算されているのである。これは「ダミー変数を用いて，交互作用項も投入した重回帰分析」と同じである。
　ちょっと難しいかも知れないが，この coefficients を次の様に合計すると，各カテゴリの平均が求められる。いずれも最後の列の Sum が，上で調べた各カテゴリの平均年収に一致している事を確認しよう。


　以下の様にも計算出来る。この方がダミー変数を投入した重回帰分析の考え方をよく表現している。