複数の要因や連続変数を組み込んだ解析をしたい！

　これまで私たちは、ある被説明変数に関して1つの要因or説明変数を用いて解析をする方法を用いてきた。例えば、植物の高さが施肥の有無で差があるかどうかを、検定や分散分析を用いて解析してきた。しかしながら、植物の高さというのは、施肥以外の要因の効果も大いに受ける。例えば、日当たり。野外で栽培実験をしていたら、実は施肥をした鉢としなかった鉢の半々がちょうど日陰に入ってしまうことが途中でわかったとしよう。もちろん、こんなことにならないように、実験を行う前に栽培場を確認しておくことは大切だが、やってみてから不備が発覚することも実践上はたくさんある。こうなってしまったら、もうこの実験は失敗であり、1からやり直したほうが良いのだろうか？　やり直したほうが良い、という意見もあるだろうが、うまいこと日当たりの効果を分離することも可能だ。また、場合によってはさらなる発見がここからわかることがあるかもしれない。例えば、施肥と日当たりが特定の組み合わせの時、予想よりも大きくなる、みたいなことがわかるかもしれない。このような状況下では複数の要因を使って分散分析を行ってみるのが良いだろう。上記は日当たりの効果が意図せず組み込まれてしまったが、日当たり＝日光量を人間側が操作して実験した場合でも、もちろん同様の解析を使うことができる。このように、統計（検定など）を使って現実データを解釈、予測しようという試みが統計モデリングである。特に明言しなかったが、今後紹介する解析に限らず、検定や分散分析も統計モデリングの1種である。

　また、上記の例では「施肥の有無」とか「日当たりの良し悪し」みたいに要因が数値で表されていない。このような要因（変数）をカテゴリカル変数（質的変数）と呼ぶ。植物の高さに影響を与える要因は、上記のようなカテゴリカル変数だけとは限らない。例えば、温度のように10℃、15℃、20℃……（さらに10℃と15℃の間も12.5℃のように取ることができる）というように連続的な影響を与える変数も存在する。このような要因は数値として表されており、連続変数（量的変数）と呼ぶ。今まではカテゴリカル変数を要因として組み込んだ検定を中心に扱ってきたが、統計モデリングが理解できれば温度のような連続変数の影響も解析をすることが可能となる。

　特に野外からデータを取得する生態学者は、調べたい被説明変数に対してカテゴリカル変数/連続変数、2つ以上の変数を多数組み込んで解析することが多い。ある種の植物の高さに影響を与える要因を調べるだけでも、各植物が生えていた場所の、気温（連続変数）、日射量（連続変数）、湿度（連続変数）、土壌の種類（カテゴリカル変数）、土壌窒素濃度（連続変数）、土地利用（カテゴリカル変数）、隣接している植物の種類（カテゴリカル変数）……などなど、様々な要因を組み込んで解析することができる。最終的に、これらの要因のうち植物の高さに強く影響を与えていた要因を割り出すことが可能となるのである。

　上記は説明変数側の制約を緩める過程ととらえることができる。さらに、もっと後で解説することになるが、被説明変数側の制約も緩めることができる。明示的に説明はしなかったが、これまで被説明変数はほとんどの場合、連続値（連続変数という意味ではない）を扱ってきた。例えば、植物の高さのような値は10 cmとか15 cmのみならず、そのあいだの12.5 cmとか11.345 cmとか、正の値であればあらゆる値をとりえるだろう。ほかにも、温度も連続値データだ。温度に関しては正の値はもちろん、-273℃までなら負の値も取りえる。だが、世の中はそのような値ばかりではない。例えば、種子の数は0以上かつ整数の値しかとらない。1.5個の種子というのは概念上は存在しても、実際に実験から得られる値ではないだろう。このような、とびとびの値をとるものを離散値と呼ぶ（実は一瞬だけ取り扱ったことがある）。連続値と離散値では、統計上の取り扱いも、数値としての性質も異なっており、特定の状況下では区別して解析が必要になる。統計モデリングの項は今までのt検定をはじめとした解析を、様々な状況に対応できるよう拡張していく過程を理解するものだと思って読み進んでほしい。

複数の要因を組み込んだ解析

　一元配置分散分析の項では、この分析を3標本以上で差を明らかにするための解析として紹介した。実は、分散分析の役割はそれだけにとどまらず、もっと広範な解析に用いることができる。それが複数の要因を解析に組み込み、どの要因が被説明変数の差を生み出しているかを明らかにできる点である。こんな解析ができるようになれば、被説明変数に影響を与えうる複数の要因をあらかじめ調査しておき、それを解析に組み込むことで、重要な要因の発見に役立てることができる。特に、野外で調査した時に得られたデータの中で、始めから被説明変数にどの要因が強く影響を与えているか、わからないことが多い。野外調査から得られた知見を、分散分析（あるいは今後紹介していくことになるより一般的な解析）により解析して重要な要因を明らかにした後は、今度は実験室でその要因を操作して被説明変数に影響が出るかを明らかにすることで、より強固な知見として発表ができるだろう。

本ページでは二元配置分散分析Two-way ANOVAなどを紹介する、その前に、これらの解析を行う上での前提知識を紹介する。

分散分析は線形モデルの一種

　分散分析は線形モデルと呼ばれる解析手法の一種である。この線形モデルには、分散分析のほかに線形回帰（ここでは単回帰のこと）や多重線形回帰（いわゆる重回帰）なども含まれている。実は、説明変数がカテゴリカル変数か連続変数か、や、複数の要因を説明変数に含んでるかどうか、で呼び名が違っているだけで、その内実はほとんど差がない。分散分析の概念が理解できれば、おのずと線形モデル群の概念も理解できるようになっていく。線形モデルとはおおざっぱに言えば、「要因の効果を足し算で表したモデルを使って被説明変数を解析すること」を言う。数学の世界で、足し算で表される関係を線形である、と表現するためにこのような呼び名が定着している。例えば、植物の高さに影響を与えうる要因として、「施肥の有無」と「日当たりの良し悪し」を組みこむとしよう。この時、

「植物の高さ（被説明変数）」 ~ a（切片項） 

　　　　　　　　　　　　　　　 + b ×「施肥の有無（説明変数）」

　　　　　　　　　　　　　　　 + c ×「日当たりの良し悪し（説明変数）」

のように表せると仮定し、この時のa、b、cを推定しようとするのが、線形モデルの目標である。ここで「=」ではなく「~」を使っているのは、等式を解いているわけではなく分析をしているから、くらいの気持ちで見てほしい。a、b、cの大きさが十分大きければその要因は強く被説明変数に影響を与える=有意差を生み出す要因である。上記のような視点では紹介してしなかったが、検定統計量はこの係数にあたる要素である。このとき、上記の式の右辺を線形予測子linear predictorと呼ぶ。さて、温度のように連続変数が説明変数ならまだしも、施肥の有無のような数値で表されていないカテゴリカル変数をどうやって数値として扱うのか、疑問に思う読者もいるだろう。後述するように、解析上はカテゴリカル変数はダミー変数というものに置き換えられて解析されている。

　別に足し算で表される、ということさえ理解していれば問題ないので、このあとの由来まで理解しなくてもよいが、一応ざっくりと示しておく。足し算で表されることから、要因の影響は直線状になるような性質を示す（下図）。例えば、直線y = ax + bの説明変数xと被説明変数yの関係は線形関係である。x1のときy1、x1 + 1のときy2、x1 + 2のときy3……という対応関係があるとしよう。このとき、y2はy1とaの足し算だけで表すことができる。y3も同様だ。このあと、y4、y5となっていっても、aを足す数を変えるだけで表現できることがわかるだろう。もちろん、x1 + 0.1のときだろうが、x1 + √2だろうが、aに0.1や√2を掛けてy1に足せば、その時のyを表現できる。逆に、aを足す数を変えるだけでxとyの関係を表現できるときは、これがxとyが直線状に存在する＝線形関係にあるということである。これが、線形という呼び名の由来になっている。

　一方で、2次関数y = ax^2 + bx + cは線形関係が満たされていない。実際、x1 + 1やx1 + 2を代入してみればわかるが、y2やy3はy1とaやbの足し算だけでなく、x1という変数とaの掛け算された項が足されている。つまり、x1の値によってxとyの増減の関係性が変わるということである（だからこそ曲がり具合が変わる線=曲線となる）。この場合は、線形とは呼ばない（工夫をすれば線形の関係に帰着できるがここではそこまでは踏み込まない）。

カテゴリカル変数を数値として扱う術：ダミー変数

　カテゴリカル変数は、数と名前がついているが、データを入力するときの取り扱い上は数値ではない。例えば、施肥の有無とか、日当たりの良し悪しは、下記左のような形で記録がとられているだろう。ここでは、施肥が有/無をYes/No、日当たりの良し/悪しをStrong/Weakで表している。植物の高さが数値として記録がとられているのと対照的だ。このようなカテゴリカル変数を数値に変換して解析をできるようにしたものがダミー変数である。ややけったいな存在だが、R上で線形モデルを使って解析しようとするとき、Rがうまいことダミー変数に変換をしてくれるので、実はスクリプトの表記上は特に違いがない。ただし、どのようにダミー変数に変換されているかを理解することで、分散分析も線形回帰などと類似した考え方に基づくことがわかる。

　ダミー変数の実体は、Yes/NoやStrong/Weakを0/1に変換した数値である。それゆえ、上記左の例は上記右のように変換されて解析されることになる。ここで、R上で特に指定をしない場合、0/1への変換は辞書並びに応じて変換される。つまり、辞書的にはNo→Yes、Strong→Weakの順に出てくるのでNoとStrongは0、YesとWeakは1に変換される。もうこうなってしまえば、説明変数が0か1しかない線形回帰とイメージはほとんど変わらない（ただし、ダミー変数と連続変数の自由度が異なるので完全に同じわけではないが）。このように、0/1に変換されることで、分散分析が足し算的に解析をしていることがより顕著にわかるようになる。

　「植物の高さ」~ a + b ×「施肥の有無」+ c × 「日当たりの良し悪し」と表されるとき、施肥および日当たりの要因に0を代入してみる。つまり、施肥はなく、日当たりが良いときの植物の高さの推定値がaとなる。これは丁度、一次関数などにおいてx = 0のとき、切片を示している。それゆえ、統計上もすべての変数に0を代入したときの推定値を切片項interceptと呼んでいる。また、a～cはまとめて（回帰）係数coefficientと呼んでる。下図を見て理解を深めてほしい。ほかに、施肥があり(1)、日当たりが良い(0)ときの植物の高さの推定値はa + bとなる。先ほど推定したaを引けばbが求まる。また、施肥がなく(0)、日当たりが悪い(1)とき、a + cが推定値で、同様にcも求まるだろう。そして、施肥があり(1)、日当たりが悪い(1)とき、a + b + cが推定値となっている。まさにこれが、各要因の効果が足し算で表されている状態といえよう。

　ところで、「3水準以上になったらどうなるんだろう」と勘の鋭い読者は思ったかもしれない。実際、一元配置分散分析の項では3標本以上の比較をしたではないか。0/1の二者択一だったら3水準以上を扱えないはずだ。こういうときが、ダミー変数と連続変数の取り扱いの違いが如実に表れる瞬間である。例えば、下記のように薬剤X、薬剤Y、薬剤Zの3水準存在し、植物の高さに与える影響を考えよう。この時、上の例にならうなら、

「植物の高さ」 ~ a + b × 「薬剤」

という風にして、aとbを推定したいのだ。しかし、ダミー変数の性質上、これでは薬剤に3水準あってうまく扱えない。

とりあえず、デフォルトのRの設定ではXが辞書において最も早いので、Xを基準に考える。つまり、薬剤Ｘを切片として扱うということである。そして、3水準のダミー変数の扱いでは、薬剤Yを処理した時は切片に薬剤Yの係数が足される、薬剤Zを処理した時は切片に薬剤Zの係数が足される、という風に考える。下図を参考にしてほしい。

こうすることで、ダミー変数の制約の下、薬剤Xの時の植物の高さはa、Yの時はa + b1、Zの時はa + b2という風に薬剤Xの効果を基準として薬剤YとZの効果を足し算で表すことができる。これは4水準以上になろうと変わらず、基準となる水準の推定値（切片）をa、残りの水準の係数をそれぞれb1、b2、……とするとき、残りの水準の効果はa + b1、a + b2、……として推定するのである。ということは、一元配置分散分析の項では説明しなかったが、内部的には各水準ごとの係数が計算されている=多重比較に似た数値は算出されている、のである。これを利用して、Tukey-Kramer法で多重比較している。

交互作用という考え方：ちょっとした非線形っぽい効果を取り扱う

　ところで、「植物の高さ」と「施肥の有無」「日当たりの良し悪し」の関係に戻ると、

「植物の高さ（被説明変数）」 ~ a（切片項） 

　　　　　　　　　　　　　　　 + b ×「施肥の有無（説明変数）」

　　　　　　　　　　　　　　　 + c ×「日当たりの良し悪し（説明変数）」

この足し算関係では取り扱えない現象が出てきてしまう。例えば、下記の図の右のように、施肥があり(1)、日当たりが悪い(1)とき、a + b + cが推定値であるはずなのに、それよりも低い（高くてもよい）場合である。こんな場合は、上記のモデル式ではうまく推定ができない。

そこでモデルの式を少し改変して扱えるようにしてみる。

「植物の高さ（被説明変数）」 ~ a （切片項）

　　　　　　　　　　+ b ×「施肥の有無（説明変数）」 

　　　　　　　　　　+ c ×「日当たりの良し悪し（説明変数）」

　　　　　　　　　　+ d ×「施肥の有無（説明変数）」 ×「日当たりの良し悪し（説明変数）」

モデルの最後に「施肥の有無」と「日当たりの良し悪し」を掛け算した項を追加した。この項を追加し、その係数dを推定することで実は、上図の右のような状況にうまく対応できる。確認してみると、「施肥の有無」と「日当たりの良し悪し」には0/1を代入するのだから、dが残るのは「施肥の有無」と「日当たりの良し悪し」の両方に1が代入されたときだけである。これにより、施肥があり(1)、日当たりが悪い(1)ときの推定値はa + b + cではなく、a + b + c + dとなり植物の高さが低くなることを数式でうまく表現できている。このような、2つ以上の要因の組み合わせで初めて生じる効果を交互作用interactionと呼び、交互作用を表している項を交互作用項と呼んでいる。上の例なら、「施肥の有無」×「日当たりの良し悪し」が交互作用項ということだ。一方で、交互作用項以外の項は主効果（単純主効果）と呼んでいる。また、dが0に近い値を示すときは、上図左のような交互作用がない状態にうまく帰着できている。このように複数の要因を解析するときには交互作用の存在に注意を払う必要がある。

　交互作用の存在は、グラフを描けば直感的に認識できる。例えば、下図のように平均値を結んだ時の傾きが平行からかけ離れていれば交互作用が影響している可能性がある。

交互作用が存在するとき、特定の条件では主効果の効果が変わるということを意味する。例えば、上記の施肥や日当たりは、施肥があることと日当たりが悪くなることは、単独では植物の高さに正の効果を及ぼした。しかし、組み合わさるとむしろこの正の効果は打ち消されて、むしろ負の効果になった。もし、交互作用がなければ2つの正の効果が組み合わさってさらに高くなるはずだったのだ。このようなことを指して、教科書的には「主効果に意味はない」と表現されることがある。ただし、この言葉をうのみにするのではなく、前述通り、特定の条件で主効果の効果が変わる、と認識しよう。

　交互作用がないとき、論文などではそれぞれの主効果がどれくらい影響が強いかなどを議論すればよい。一方で、交互作用があるときは、それぞれの主効果の組み合わせに応じた議論が必要になることを覚えておこう。例えば、上記の例なら施肥、日当たりの単独の効果を説明した後、施肥あり＆日当たり悪いときにより高くなるのではなく、低くなってしまった生物学的原因を考察することになるだろう。