線形モデルの枠組み
ここでは、解析のさまざまな目的に対応できるとても汎用性の高い枠組みである、”線形モデル”について概念的な説明をおこなう。
線形モデルは、ある要因（説明変数）が別の要因（目的変数）に与える影響を検証するものである。
説明変数は複数でもよく、その場合には複数の説明変数それぞれの影響が個別に検証される。
ANOVA（分散分析）、回帰分析、ANCOVA、重回帰分析などなど、すべて線形モデルが使われている。
まず、線形モデルが統計解析全体の中でどのような位置を占めるかを説明し、次に線形モデルの考え方のエッセンスを図を使って説明する。最後に、Rによる一般線形モデルの記法を例を用いて解説する。
線形モデルの位置づけとその種類
線形モデルは「ある要因が他の要因に与える影響（効果）を検証したい」という場面で用いられる。
一見、特定の目的にしか使えないように感じられるかもしれないけれど、「群間の差を検証したい」という目的も線形モデルの枠組みで扱える（言い換えれば「”複数群”という要因の影響を検証したい」ということなので）ことを考えると、線形モデルの対応できる場面はかなり広い。
下の図は、よく行われる統計解析とその目的を雑多にまとめたものである。
線形モデルにもいろいろある…という印象を受けるかもしれないが、実のところ大きく分ければ一般線形モデル、GLM、ランダム効果入りモデルの3種類しかない（ほかにはGLSなど）。
どれを使うかは、1. 目的変数の正規性・等分散性がみたされているか？ 2. 同一個体や同一場所からの反復測定があるか？の主に2つの分かれ道で決まる（正規性とは、目的変数が正規分布かどうかということ。実践的には、データの分布が一山型で左右対称であれば良いと思う。等分散性とは、説明変数の変化に関わらず目的変数のばらつきが変わらないということ。これが満たされないと、群間（あるいは説明変数の値の違いで）でばらつきに違いが生まれる）。
なお、PCA、RDA、CCAやMANOVAは線形モデルを多変量に拡張したものである。
線形モデルの考え方
図2
各説明変数が目的変数に与える影響の、線形モデルによる検証の背景にあるエッセンスは、”変動（分散）の分解”である。
たとえば、図3のようなデータYが得られて、シンボルの形と色によって、Yの値が異なるのかどうか知りたいとする。
図3
まずはシンボルの形で2群に分け、各群の平均値を算出してみよう。
図4
左が〇の群、右が△の群であり、グレーの実線と点線が各群の平均値を示している。
確かに△の平均値は〇の平均値に比べて高くなっている。一方で、各群の点は各群の平均値の周りで大きくばらついていることもわかる。
ここで、〇群と△群の平均値のばらつき（2つしかないのでわかりづらいが）を群間変動と呼び、〇群内・△群内のばらつきを群内変動と呼ぶ。
両者は全変動との間に、以下のような関係をもつ。
全変動＝群間変動＋群内変動
つまり、全変動は群間変動と群内変動に分解することができる。
この関係に基づいた図4のような操作を”変動の分解”と呼ぶ。
図に戻ると、各群の平均値には確かに違いがみられるが、群内のばらつきが互いの平均値を大きく超えてばらついている。つまり、群間変動に対して群内変動がかなり大きい。そのため、異なるシンボルの形の間に本当は差があったとしても、それが検出しづらくなっている恐れがある。
いったん最初に戻り、今度は、1. 色の影響を調べて、2. 次に形の影響を調べる、という方法を取ることにする。
色によって差があるか調べるため、異なる色で3つに分け、各群の平均値を算出する。
図5
先ほどと同じように、全変動を（色の）群間変動と群内変動に分解することができた。赤が小さく青が大きく、緑はその中間に位置していることがわかる。
今度は、各群の群内変動そのものに注目してみよう。よく見ると、赤群でも青群でも緑群でも、〇が下に、△が上に位置する傾向にありそうだ。言い換えると、この群内変動はシンボルの形の影響を受けている可能性がある。
このような、形の影響を含む群内変動だけにもっと注目するため、色の影響は無視して考えたい。
そのために、少しトリッキーだが、各群を全平均に向かって上下方向に丸ごと平行移動させ、各群の平均の水平線が全平均の水平線に一致するように揃えてしまおう（各点ごとに、その点が属する群の平均と全平均との差分を取ったことになる）。
図6
これは、全変動からの色の群間変動の引き算に相当する。これができるのは、全変動＝群間変動＋群内変動という関係があるからである。
色の影響は取り除かれたので、すべてを一群にまとめなおそう。
図7
〇と△の差を比較するため（形の影響）、これを形でさらに2群に分け、各群の平均値を算出する。
図8
図4と比較すればわかるように、最初に全データを形で分けた時より、明らかに群内変動が小さくなっている。色の影響を取り除いたことにより、形の影響が見やすく（高い精度で検証できるように）なったのである。
図6-7と同じ操作をして、形の影響も取り除いてしまおう。
図9
山なりの密度プロットを示す変動の小さい点の集団が残された。このように、全変動からモデルに入れた要因の変動を取り除いたものを、モデルの誤差または残差と呼ぶ。線形モデルの各要因の影響の有意性は、基本的に群間変動と誤差を比べることで行われる。
変動が分解される過程をおさらいすると以下のようになる。
図10
これを改めて式で書き下すと、教科書等でよく目にするお馴染みの式になる。
このように、線形モデルは変動を分解することで各要因の影響を高い精度で検証できる枠組みである。上の説明では、色の影響を取り除くことで形の影響が見やすくなると説明したが、誤差からは形と色双方の影響が取り除かれているため、見方を変えれば色の影響もより高い精度で検証できるようになっている。このことは、適切な要因のセットをモデルに組み込むことで、互いの影響をより検証しやすくなることを意味している。（ただし、要因を増やし過ぎるとモデルが複雑になるので避けた方がよい）
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個のデータが存在する。
> summary(m1)
Call:
lm(formula = Value ~ 1, data = data)
Residuals:
Min 1Q Median 3Q Max
-21.696 -7.709 -1.001 7.931 21.817
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.4441 0.8831 16.36 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 9.674 on 119 degrees of freedom
>
> mean(data$Value)
[1] 14.44413
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関数はカッコ内で指定したデータフレームの列名によってベクトルを表すことができるようになる関数である。
> summary(m2)
Call:
lm(formula = Value ~ Symbol, data = data)
Residuals:
Min 1Q Median 3Q Max
-17.6929 -6.9881 -0.9527 6.6103 17.8139
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.441 1.141 9.152 2.02e-15 ***
Symbolt 8.005 1.613 4.962 2.37e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 8.837 on 118 degrees of freedom
Multiple R-squared: 0.1726, Adjusted R-squared: 0.1656
F-statistic: 24.62 on 1 and 118 DF, p-value: 2.371e-06
>
> with(data, tapply(Value, Symbol, mean))
c t
10.44139 18.44688
ここでも、Coefficientsのリストにある、InterceptとSymbolの行のEstimateに注目しよう。
InterceptのEstimateとマニュアルで計算したcの平均が、InterceptのEstimate+SymboltのEstimate（10.441+8.005）とbの平均が一致している。図4の平均値も参照のこと。また、Residual standard errorの値が9.674から8.837に下がっている。
> summary(m3)
Call:
lm(formula = Value ~ Color, data = data)
Residuals:
Min 1Q Median 3Q Max
-14.8119 -4.1897 0.1287 3.6635 14.4522
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.0566 0.9875 6.133 1.20e-08 ***
Colorblue 18.0040 1.3966 12.892 < 2e-16 ***
Colorgreen 7.1586 1.3966 5.126 1.18e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 6.246 on 117 degrees of freedom
Multiple R-squared: 0.5902, Adjusted R-squared: 0.5832
F-statistic: 84.26 on 2 and 117 DF, p-value: < 2.2e-16
>
> with(data, tapply(Value, Color, mean))
red blue green
6.056595 24.060584 13.215225
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に下がっている。
このように、線形モデルとは、「目的変数の平均値が説明変数の値によってどう変化するかを表現するモデル」であると言うことができる。
> summary(m4)
Call:
lm(formula = Value ~ Symbol + Color, data = data)
Residuals:
Min 1Q Median 3Q Max
-11.1367 -3.2716 -0.0953 3.6270 10.4495
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.0538 0.8712 2.357 0.0201 *
Symbolt 8.0055 0.8712 9.189 1.89e-15 ***
Colorblue 18.0040 1.0670 16.873 < 2e-16 ***
Colorgreen 7.1586 1.0670 6.709 7.54e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.772 on 116 degrees of freedom
Multiple R-squared: 0.7628, Adjusted R-squared: 0.7567
F-statistic: 124.4 on 3 and 116 DF, p-value: < 2.2e-16
>
> aggregate(Value ~ Symbol + Color, data, mean)
Symbol Color Value
1 c red 2.513413
2 t red 9.599777
3 c blue 19.974886
4 t blue 28.146282
5 c green 8.835865
6 t green 17.594584
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)
> summary(m5)
Call:
lm(formula = Value ~ Symbol * Color, data = data)
Residuals:
Min 1Q Median 3Q Max
-11.5133 -3.5515 0.0192 3.5699 10.0728
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.513 1.073 2.342 0.0209 *
Symbolt 7.086 1.518 4.668 8.35e-06 ***
Colorblue 17.461 1.518 11.503 < 2e-16 ***
Colorgreen 6.322 1.518 4.165 6.08e-05 ***
Symbolt:Colorblue 1.085 2.147 0.505 0.6142
Symbolt:Colorgreen 1.672 2.147 0.779 0.4376
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.8 on 114 degrees of freedom
Multiple R-squared: 0.7641, Adjusted R-squared: 0.7538
F-statistic: 73.87 on 5 and 114 DF, p-value: < 2.2e-16
>
> aggregate(Value ~ Symbol + Color, data, mean)
Symbol Color Value
1 c red 2.513413
2 t red 9.599777
3 c blue 19.974886
4 t blue 28.146282
5 c green 8.835865
6 t green 17.594584
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
のように、完全に一致しているはずである。