Rと統計の入門

このチュートリアルのゴールは「一般化線形混合モデルを用いた統計分析」の理解です。そのためには「モデル」や「統計」の理解が必要であるため、導入として以下に (i) モデルが必要な理由と (ii) 最終的に理解したい対象を考え、 (iii) 過程で必要な概念を整理していきます。

なお、このチュートリアルのコンセプトは「手を動かしながら理解する」なので、下のような作業場を使っていきます¹。アドブロックをオンにしてるとセッションが壊れる場合があるので本サイトではオフにしてください²。


    # print 関数を使って
    # "Hello LME!" と出力してみてください
    print("Hello LME!")


    print("Hello LME!")


    test_function("print")
    success_msg("Great job!")

Use the print function to say "Hello LME!"

モデル: 現象を理解するためのツール

まずは「なぜモデルが必要なのか」を考えるため、現象の理解を目的とする研究を想定します。例として「主語述語の距離と認知的負荷(読みづらさ)の関係」や「主語–述語の一致と認知的負荷の関係」などが挙がります。ここで認知的負荷を「文中の単語の読み時間」として観測できるなら、 (a) 何がどの程度影響するのかを知れば現象の理解につながります。また (b) その影響の確からしさを知れば報告の信憑性を支持できます。

上の単語の読み時間(ms)を$y$、主語–述語間の単語数を$x_1$、主語–述部の一致(0/1) を$x_2$ で表現します。ただ $y = x_1 + x_2$ にはならないことから³、 $y$と各要因$x_1, x_2$の関係を調整する $a_1$ や $a_2$ も考えます。例えば、$x_1$ が読み時間に影響しない場合は $a_1 = 0$ 、大きく効果がある場合は大きな値、減少させる場合は負の値になります。するとさっきの式は $y = a_1x_1 + a_2x_2 + b$ と書けます。

一度用語の整理をすると $y = a_1x_1 + a_2x_2 + b$ のような式は「線形モデル」と呼ばれ、観測値$y$を応答変数(response variable)、要因$x_1$や$x_2$を説明変数(explanatory variable)と呼びます。また説明変数にかける値 $a_1$や$a_2$ を傾きとよび、 $b$ を切片と呼びます。傾きや切片は $\beta$ や Estimate と呼ばれることもあり、特に傾きは説明変数のウェイト、重み、重要さを表します。

ここでもう一度 $y = a_1x_1 + a_2x_2 + b$ の式を見ると、先ほどの問に答えられそうなことがわかります。つまり (a) 説明変数(主語述語の距離など)の傾きを比較すれば影響の大きさがわかります。さらに、もし(b) 影響の確からしさがわかれば説明変数の影響を統計的な作法で主張できます。したがって、現象の理解という目的は応答変数と説明変数を結ぶモデルの作成で達成できます。これがモデルを学ぶ意義になります。

【補足】なお上は少し偏った例ですが、同様に手持ちの課題を従属変数と説明変数に落とし込めたとします⁴。画面の注視時間でもアンケートの値でもなんでもOKで、その変数間の関係を上の(a)と(b)で調べられます。

一般化とノイズ

次に「最終的に理解したい対象」を考えますが、結論を言えば変数間の直線以外の関係や応答変数が従う様々な分布、そして説明変数の構造的ノイズを扱えるモデルです。先ほどの例では反応時間を線形モデルで表現しましたが、線形モデルは説明変数と応答変数(読み時間)が直線の関係にあることや応答変数が負の値や小数の値を取れると仮定しています。

ただ、アンケート(0–5など)や文の容認性(0–1)は小数点を取れず、テストの点(0–100)も勉強時間に対し直線的に増えるとは限りません。そこで線形モデルにパーツを加えて説明変数に制限を設け、変数間の関係に直線以外も扱えるように一般化したのが「一般化線形モデル」です。これがチュートリアルの中間地点になります。

次の「構造的ノイズ」に関しては読み時間の個人差を考えましょう。読み時間では説明変数に反応しやすい人、しづらい人がいます。この個人差のような、何かに条件付けられた構造的なノイズを「ランダム効果」、全体に共通する効果を「固定効果」として分けてモデリングすることで、本当に見たい効果を見れます。これが「一般化線形混合モデル」です。このチュートリアルの最終到達目標になります。

冒頭でこのチュートリアルのゴールは一般化線形混合モデルを用いた統計分析の理解だと述べました。つまり、$y = a_1x_1 + a_2x_2 + b$ のような線形モデルを一般化させ、さらに構造的なノイズを考慮できるモデルを作り、 (a) 要因の貢献度合いと (b) 貢献度合いの確からしさ(統計) を検証するのです。そうすれば様々な現象の理解を構造的なノイズを排除して理解でき、研究の目的を達成する一助となります。

次に実際に手を動かしながら最初のステップである線形モデルの作成と統計分析を進めて行きましょう。その過程で必要な概念を整理していきます。

モデルの基本要素

次に統計モデルとその基本要素である (i) data.frame (ii) formula (iii) model の概念を一つずつ抑えていきます。以下、先ほどの例で出てきた $y = a_1x_1 + a_2x_2 + b$ のような「線形モデル」を元に考えます。

data.frame

まず簡単な例として sleepstudy というデータを見てみましょう。データの概要は睡眠不足(寝不足何日目か, Days)を説明変数、反応時間(Reaction)を応答変数とし、これらの関係を被験者ごと(Subject)に見た、というものです。いかにも「Daysが上がればReactionも増えそうだな」という変数間の関係がわかります。

The average reaction time per day for subjects in a sleep deprivation study. On day 0 the subjects had their normal amount of sleep. Starting that night they were restricted to 3 hours of sleep per night. The observations represent the average reaction time on a series of tests given each day to each subject.

以降、「データを見る」といったときの「データ」は data.frame という概念を意味しているんだな、と読み替えてください。この data.frame というのは Excel を思い浮かべると分かりやすく、一行目は列名で二行目からは各列のデータを持ちます。それを R ではデータフレーム (data.frame) と呼びます。

もちろん data.frame の全てを表示することは簡単なのですが、 Rには()を右に持つ「関数」というのがあり、例として head() があります。この head() は data.frame を () の中に受け取ると最初の6件をデフォルトで返します。以下に sleepstudy という data.frame を既に読み込んであるので、 sleepstudy を head() の () に入れて動作を確認してみましょう。


    # lme4パッケージは読み込み済み
    library(lme4)


    # sleepstudy というデータが手元にある状態
    # head() の () に sleepstudy と入力して Submit
    head()


    head(sleepstudy)


    test_function("head")
    success_msg("Great job!")

Use head and summary to see the data.

sleepstudy の最初の6行を確認できたでしょうか。 1行目のReaction は249.5600(ms)、 Days は0, Subject(被験者番号)は308 となっているので、寝不足0日目の308番さんの反応時間は249msだったことがわかります。

ここまでをまとめると、

Rは data.frame の形でデータを扱う (例: sleepstudy)。
data.frame には列(変数)と値がある (例:Reactionと249.5600)。
data.frame を受け取る関数がある (例: head())。

ということがわかります⁵。

formula

次は formula という概念です。先程は $y = a_1x_1 + a_2x_2 + b$ という式を見ましたが、これを sleepstudy に当てはめるとどうなるでしょうか (Subject はいわゆる「構造的なノイズ」なので今は無視します)。

応答変数Reaction を Days(寝不足何日目か) で説明すると、 $Reaction = a \cdot Days + b$ という式になります。ここで説明変数(Days)の (a)効果と(b) 確からしさは$a$の大きさとバラツキでわかりそうです。もし $a=0$ でなさそうなら統計的に効果がある、と主張できます。これがよく論文で見る「有意差が見られた($p<.05$)」が意味するところです。

Rでは $Reaction = a \cdot Days + b$ という変数間の関係、式を Reaction ~ Days のように ~ と変数のみで表現でき、これを formula と呼びます。ここで Reaction や Days という文字列は data.frame の列名と一致していなければなりません。なので Reaction ~ a*Days + b (Rは掛け算を*で表現します) や reaction ~ days は定義できますが、後々エラーになります。

この formula はデータ内の列、つまり変数間の関係を記述しており、 Rの「model」を作るときだけでなくグラフの描画も楽になります。この model はまだ説明していないので、まずはグラフの描画を例にします。グラフの描画に使える関数は xyplot であり⁶、 formula と data.frame を受け取って変数間の関係を描画します。カンマで区切ってformulaとsleepstudy渡してましょう。


    library(lattice)  # for xyplot
    library(lme4)  # for sleapstudy data


    # 1. Reaction と Days の関係を ~ で記述し
    #    一旦 f に格納 f = `Your Answer`
    f = 

    # 2. xyplot に f と sleepstudy を渡して Submit
    xyplot(, ) # edit here


    f = Reaction ~ Days
    xyplot(f, sleepstudy)


    test_function("xyplot")
    success_msg("Great job!")

ちゃんと右肩あがりの図が描画されたでしょうか。説明が分かりづらかった場合はHintを押して答え(Solution)を確認した後、もう一度実行してみてください。

ここまでをまとめると、

Rの formula は ~ で data.frame の列(変数)同士の関係を記述できる (例: Reaction ~ Days)。
formula はグラフの作成やモデルの作成に使える(例: xyplot)。
formula の変数は data.frame 内に限定される (例: sleepstudy の記述に days は使えない)。

となります。

model

最後は model です。上の formula で $Reaction = a \cdot Days + b$ を Reaction ~ Days と表現できるものの、 Reaction ~ a*Days + b としても aやb が未知なのでエラーになると述べました。 Rの model の役割はこれら未知のパラメータ、 a(傾き)やb(切片)の推定で、これこそ我々が知りたい a の大きさや確からしさになります。

Rでは model(formula, data.frame) のように入力することでmodelがパラメータを推定します。つまりmodel はheadと同様に関数であり、 formula と data.frame を受け取ります。そして formula の Reaction ~ Days には a や b が隠れていて、それらを調整して上手いこと Reaction ~ a*Days + b が成り立つようなパラメータを model は推定します。

なお model には色々な種類があり線形モデルの lm, 一般化線形モデルの glm、そして一般化線形混合モデルの glmer と抽象化されていきます。ここでは最も単純な線形モデルlmを使って挙動を確認してみましょう。 f = の右に Reaction と Days の関係を ~ で記述すると、 f に記述の内容が格納されます。


    library(lme4)  # for sleapstudy data


    # 1. Reaction と Days の関係を ~ で記述
    f = 
    # 2. lm に f と sleepstudy を渡して Submit
    lm(, ) # edit here


    f = Reaction ~ Days
    lm(f, sleepstudy)


    test_function("lm")
    success_msg("Great job!")

上のコードを走らせた結果の Coefficients(係数)を見ると、 Intercept(切片)が 251.41、 Days(傾き)が 10.47になることがわかります。先ほど以下を述べたのですが、

そして formula の Reaction ~ Days には a や b が隠れていて、それらを調整して上手いこと Reaction ~ a*Days + b が成り立つようなパラメータを model は推定します。

傾きaが10.47、切片bが251.41の時に、 Reaction ~ a*Days + b が上手いこと成り立つという事になります。つまり $Reaction = 10.47 \cdot Days + 251.41$ が最も尤もらしい式となり、意味は「基本の反応時間は251.41ms(つまり0.2秒ほど)で寝不足が1日増えるごとに10.47msずつ反応が遅くなる」と解釈できます。

もちろん分野にもよるのですが、ここで「その差って統計的に信頼できるの？」という疑問が出てきます。そこで最後に summary() という関数に上の model を渡します。関数の埋め込み(summary(lm(model, data)))を避けるために、また一旦 m に model を格納しましょう。


    library(lme4)  # for sleapstudy data


    # 1. Reaction と Days の関係を ~ で記述し
    f = 
    # 2. lm に f と sleepstudy を渡して Submit
    m = lm(, )
    # 3. m を summary に渡す
    summary()


    f = Reaction ~ Days
    m = lm(f, sleepstudy)
    summary(m)


    test_function("summary")
    success_msg("Great job!")

情報が多く出てきたと思いますが、大切なのは説明変数($Days$) と同じ行にある Pr(>|t|) の値です。分野にもよりますが、これが 0.05 より小さければ「その傾きは統計的に意味のあるものなんだよ」と安全に主張できます。なお 9.89e-15 は 9.89の.を左に15ずらしたものです(0.000…)。

ここまでをまとめると、

R の model は formula と data.frame を受け取って、formulaの未知のパラメータを推定する (例: 説明変数に対する傾きaや切片b)。
modelが求めた「要因に対する傾きや切片」を見れば要因と現象の関係がわかる (例: $Reaction = 10.47 \cdot Days + 251.41$)。
傾きの確からしさは summary() 関数で分かる (例: Days の Pr(>|t|) は 9.89e-15)。

となります。傾きの値や確からしさの求め方が気になる方はそれぞれデータ解析のための統計モデリング入門の p.24 と p.51 に関連する記述があります。

まとめ

まとめると、私達が分析時に行うことは

data.frame の用意 (例:調査や実験で応答変数(Reaction)と説明変数(Days)のデータを集めて整形)
data.frame の変数間の関係を formula で定義(例: Reaction ~ Days)
model に formula と data 与えてパラメータを推定、解釈

となります。データを用意して関係を定義しモデルを作成することで、モデルのパラメータを推定します。そして推定したパラメータの大きさと確からしさを根拠に現象を理解し主張を統計的にサポートします。

今回扱った lm は線形モデルでした。ただ「Reaction じゃなくて “正答か否か” みたいな 0/1 の値は直線でモデルできないじゃないか」や「Reactionに対するDaysの効果や、そもそも Reaction のベースの時間って人によるじゃないか」といったツッコミもありえます。それぞれの問に対応できるのが次のユニットの一般化線形モデル、最後のユニットの一般化線形混合モデルになります。

脚注

以前はゼミの勉強会などで RとRStudioをインストールして、パッケージをインストールして、あぁバージョンが合わない、そもそも文字化けする、という悲しみがあったのですが、最近はウェブで実行できます。 ↩
検証してくれたOさんに感謝します。 ↩
「5(単語数)+1(一致)」としても読み時間にはならない。 ↩
この直接観測できない対象を観測できる値で代用する作業を「Linking Hypothesis を置く」なんて言ったりします。 ↩
どうしても実体が気になるひとは R Console の方に sleepstudy とだけ打ち込んでみてください。データを全件表示しますが、これはデータが少ないことを知っている時以外は行わないでください。 ↩
後々 ggplot というリッチなツールを使う事にもなると思いますが、今は分かりやすい xyplot を使います。 ↩