線形モデルの枠組み

ここでは、解析のさまざまな目的に対応できるとても汎用性の高い枠組みである、”線形モデル”について概念的な説明をおこなう。

線形モデルは、ある要因（説明変数）が別の要因（目的変数）に与える影響を検証するものである。

説明変数は複数でもよく、その場合には複数の説明変数それぞれの影響が個別に検証される。

ANOVA（分散分析）、回帰分析、ANCOVA、重回帰分析などなど、すべて線形モデルが使われている。

まず、線形モデルが統計解析全体の中でどのような位置を占めるかを説明し、次に線形モデルの考え方のエッセンスを図を使って説明する。最後に、Rによる一般線形モデルの記法を例を用いて解説する。

線形モデルの考え方

各説明変数が目的変数に与える影響の、線形モデルによる検証の背景にあるエッセンスは、”変動（分散）の分解”である。

たとえば、図3のようなデータYが得られて、シンボルの形と色によって、Yの値が異なるのかどうか知りたいとする。

まずはシンボルの形で2群に分け、各群の平均値を算出してみよう。

図に戻ると、各群の平均値には確かに違いがみられるが、群内のばらつきが互いの平均値を大きく超えてばらついている。つまり、群間変動に対して群内変動がかなり大きい。そのため、異なるシンボルの形の間に本当は差があったとしても、それが検出しづらくなっている恐れがある。

いったん最初に戻り、今度は、1. 色の影響を調べて、2. 次に形の影響を調べる、という方法を取ることにする。

色によって差があるか調べるため、異なる色で3つに分け、各群の平均値を算出する。

このような、形の影響を含む群内変動だけにもっと注目するため、色の影響は無視して考えたい。

そのために、少しトリッキーだが、各群を全平均に向かって上下方向に丸ごと平行移動させ、各群の平均の水平線が全平均の水平線に一致するように揃えてしまおう（各点ごとに、その点が属する群の平均と全平均との差分を取ったことになる）。

色の影響は取り除かれたので、すべてを一群にまとめなおそう。

〇と△の差を比較するため（形の影響）、これを形でさらに2群に分け、各群の平均値を算出する。

図6-7と同じ操作をして、形の影響も取り除いてしまおう。

変動が分解される過程をおさらいすると以下のようになる。

これを改めて式で書き下すと、教科書等でよく目にするお馴染みの式になる。

このように、線形モデルは変動を分解することで各要因の影響を高い精度で検証できる枠組みである。上の説明では、色の影響を取り除くことで形の影響が見やすくなると説明したが、誤差からは形と色双方の影響が取り除かれているため、見方を変えれば色の影響もより高い精度で検証できるようになっている。このことは、適切な要因のセットをモデルに組み込むことで、互いの影響をより検証しやすくなることを意味している。（ただし、要因を増やし過ぎるとモデルが複雑になるので避けた方がよい）

Rで作る線形モデル

最後に、Rによる線形モデルの記法・作り方を説明する。

関数と記法

線形モデルはlm関数によって作成できる。lm関数では、1番目の引数にモデル式が入り、2番目に線形モデルを作るためのデータセット（データフレーム）が入る。

データ"iris"を使って、Sepal.Lengthを目的変数にして線形モデルを作成してみる。

data(iris) #データの読み込み

##線形モデル作成

m.intercept <- lm(Sepal.Length ~ 1, iris) #Y=intercept+ε

m.Species <- lm(Sepal.Length ~ Species, iris) #Y=intercept+X+ε

m.PLength <- lm(Sepal.Length ~ Petal.Length, iris) #Y=intercept+X+ε

m.SpeciesPLength <- lm(Sepal.Length ~ Species + Petal.Length, iris) #Y=intercept+X1+X2+ε

m.SpeciesbyPLength <- lm(Sepal.Length ~ Species*Petal.Length, iris) #Y=intercept+X1+X2+X1xX2+ε

m.SpeciesbyPLength2 <- lm(Sepal.Length ~ Species + Petal.Length + Species:Petal.Length, iris) #m.SpeciesbyPLengthと全く同じモデルが作成される

##作成したモデルの概要

summary(m.intercept)

summary(m.Species)

summary(m.PLength)

summary(m.SpeciesPLength)

summary(m.SpeciesbyPLength)

summary(m.SpeciesbyPLength2)

切片（intercept）は1で表されるが、切片のみモデルでなければモデル式からは省略可能である。また、記号*と:が登場しているが、これは"交互作用"項を作るための記法である。Rでは、要因1と2の交互作用X1×X2はX1:X2と表される。Y ~ X1*X2 は Y ~ X1 + X2 + X1:X2 と全く同じである。右のように書き下すのは面倒なので、モデル式上は普通は*を使うことが多い。

交互作用の説明はこのページの最後で行う。

実践

”線形モデルの考え方”のセクションを、線形モデルを使って改めて理解し直してみよう。

#データ作成

set.seed(101)

c.red <- rnorm(20, 3, 5)

c.blue <- rnorm(20, 20, 5)

c.green <- rnorm(20, 10, 5)

t.red <- rnorm(20, 10, 5)

t.blue <- rnorm(20, 27, 5)

t.green <- rnorm(20, 17, 5)

data <- data.frame(

  Value=c(c.red, c.blue, c.green, t.red, t.blue, t.green), 

  Symbol=c(rep("c", 60), rep("t", 60)), 

  Color=c(rep("red", 20), rep("blue", 20), rep("green", 20), 

          rep("red", 20), rep("blue", 20), rep("green", 20)))

data$Symbol <- as.factor(data$Symbol)

data$Color <- as.factor(data$Color)

head(data)

以上のスクリプトをコピペすると、上記セクションのデータを作成することができる。

Symbolをシンボルの形に、Colorをシンボルの色と読み替えてほしい。

各々の群にはそれぞれ20個のデータが存在する。

まず、図3を線形モデルで表現してみよう。

図3の上部にある式の通りに、線形モデルを作成してみる。

m1 <- lm(Value ~ 1, data)

summary(m1)

mean(data$Value)

Coefficientsのリストのうち、(Intercept)の行のEstimateを、mean関数で計算した全平均と比較してみると、完全に一致している。図3の全平均とも見比べてみてほしい。また、Residual standard errorの数値は”全変動”を示している(sd(data$Value)を計算すると、これも9.674となる)。

次に、図4を表現してみよう。

図4のような、シンボルの形という要因による変動の分解は、線形モデルの説明変数にシンボルの形（Symbol）を投入することで行うことができる。

m2 <- lm(Value ~ Symbol, data)

summary(m2)

with(data, tapply(Value, Symbol, mean))

なお、tapply関数はベクトルに対してカテゴリごとに指定の関数を適用する関数であり、with関数はカッコ内で指定したデータフレームの列名によってベクトルを表すことができるようになる関数である。

ここでも、Coefficientsのリストにある、InterceptとSymbolの行のEstimateに注目しよう。

InterceptのEstimateとマニュアルで計算したcの平均が、InterceptのEstimate+SymboltのEstimate（10.441+8.005）とbの平均が一致している。図4の平均値も参照のこと。また、Residual standard errorの値が9.674から8.837に下がっている。

図5でもまったく同じことが可能である。

m3 <- lm(Value ~ Color, data)

summary(m3)

with(data, tapply(Value, Color, mean))

InterceptのEstimateはColorのredの平均に（6.566）、

InterceptのEstimate+ColorblueのEstimateはColorのblueの平均に（24.061）、

InterceptのEstimate+ColorgreenのEstimateはColorのgreenの平均に（13.215）、それぞれ一致している（図5の平均値参照）。

また、Residual standard errorの値が9.674から6.246に下がっている。

このように、線形モデルとは、「目的変数の平均値が説明変数の値によってどう変化するかを表現するモデル」であると言うことができる。

さらに、図5のモデルにShoriAを説明変数として投入し、図8-9と同じことを行う。

m4 <- lm(Value ~ Symbol + Color, data)

summary(m4)

aggregate(Value ~ Symbol + Color, data, mean)

なお、aggregate関数はtapply関数と機能は似ているが、関数を適用するカテゴリの種類を2つ以上指定できる。

Residual standard errorは6.246から4.772にさらに下がっている。

いっぽう、本来であれば、

Intercept(2.0538)がSymbol(c)Color(red)の平均(2.513413)に、

Intercept+Symbolt(10.059)がSymbol(t)Color(red)の平均(9.599777)に、

…

という感じに一致していくはずだが、モデルと実際の平均値に食い違いが見られる。

これは、このモデルがSymbolとColorとの"交互作用"を考慮していないからである。

丸・赤と三角・赤の平均値の差と、丸・青と三角・青の平均値の差と、丸・緑と三角・緑の平均値の差がすべて一致するようなことは実際にはあまりありそうにないことである。にもかかわらず、このモデルは色が赤でも青でも緑でも、形における丸と三角の差は等しく8.0055であると仮定している。

Symbolにおけるcとtの差が、Colorの水準（red,  blue, green）ごとに異なる、という効果をSymbolとColorの交互作用と呼ぶ。

最後に、交互作用も含めたモデルを作ってみる。

m5 <- lm(Value ~ Symbol*Color, data)

summary(m5)

aggregate(Value ~ Symbol + Color, data, mean)

c.redの平均       = Intercept

t.redの平均       = Intercept + Symbolt

c.blueの平均    = Intercept + Colorblue

t.blueの平均    = Intercept + Symbolt + Colorblue + Symbolt:Colorblue

c.greenの平均 = Intercept + Colorgreen

t.greenの平均 = Intercept + Symbolt + Colorgreen + Symbolt:Colorgreen

のように、完全に一致しているはずである。