化学研究で役に立つデータ解析入門：エクセルでも立派な解析ができるぞ編ではエクセルでも立派な回帰分析ができることを紹介しました。しかし、回帰分析は奥が深く、エクセルの機能では、できないこと、できても手間がかかることがたくさんあります。今回は、フリーソフトウェアの統計解析向けのプログラミング言語Rを使ってエクセルよりも簡単により深く回帰分析を使いこなす方法を紹介します。
エクセルの回帰分析ではできないこと
回帰分析は奥が深く、エクセルではできない使い方や計算されない統計量がたくさんあります。化学実験の結果を解析するにあたっては下記の機能が特に有用であり、今回はこれらに焦点を当てて解析を行っていきます。
分類変数の取り扱い：実験の条件は、温度や反応時間といった連続した数字（連続変数）だけでなく、基質の官能基や触媒の種類といった数で表すことができない条件（カテゴリー変数）もあり、これらの傾向を掴むことも重要な解析となります。しかしエクセルでカテゴリー変数を取り扱うには、各種類（例えば官能基）ごとに１か0で表現（ダミー変数）して解析する必要があり、解析前のデータの加工が必要になります。
信頼区間：信頼区間とは母集団の真の値が含まれることがかなり確信できる数値範囲のことで、この”かなり確信できる”は、ばらつきの状況によって変わり、95％（95％信頼区間）が一般的です。具体的を95％信頼区間を表すと、ある条件で100回実験を行ったら95回の結果が含まれる結果の範囲であり、予測の誤差を直感的に知るのに便利な統計量です。エクセルでも各係数における信頼区間は算出されていますが、結果を予測する際に重要なのは任意の条件における信頼区間であり、それは画像検索するとたくさんヒットするラッパ型の曲線を算出することで、あるｘにおいてｙがどれくらいの幅があるのかを知ることができます。エクセルの回帰分析で算出される統計量でもグラフを作ることは可能ですが、手作業でグラフを作る必要があります。
交互作用：例えばある合成において反応温度を上げれば収率は単調増加するが、同時に触媒量を増やすことで収率の増加量が加速する場合、反応温度と触媒量に正の相互作用があると言えます。逆に反応温度を上げれば収率は単調増加するが、同時に基質の量を増やすと収率の増加量が抑えられる場合。反応温度と基質量に負の相互作用があると言えます。エクセルで相互作用を調べるには単純に、条件同士を掛け合わせた新しい条件を追加して回帰分析を行えば良いのですが、相互作用ごとに行を追加すること必要で、条件が多い場合には不向きです。
Rとは
統計解析に特化したソフトウェアはいくつかありますが、ここでは、フリーのプログラミング言語Rを使っていきます。プログラミング言語といわれると、テキストエディタを使って真っ白なところにガリガリとコードを書いて何度も実行してやっとほしい結果が得られるイメージを自分は持っていますが、このRは非常に簡略的にコードを書くことが可能で、前回行ったボルトの数を数える回帰分析解析もたった2行で同様の結果を得ることができます。
前回行ったエクセルで行った解析を行うためのコード
Rstudioと呼ばれる統合開発環境もフリーで提供されていて、コードの記述支援やデータの確認も容易にできるようになっています。
Rstudioのインターフェイス
RとRstudioインストールの仕方は、ネット上にたくさん紹介されていますが、一例として下記のサイトを挙げます。
Rによる回帰分析の基礎
有用な解析に移る前にRでの回帰分析の基礎的な方法を紹介しますが、この記事では、どのようにコードを記述すれば解析ができるかについて注力し、記号の意味などは取り扱いません。コードのルールなどを知りたい場合には他のページを参照しください。また、コードは同じ機能でもいろいろな書き方があり、他のサイトでは異なるコードを使っているかもしれません。
まず、データを読み込むために右上のImport Datasetをクリックして目的のファイルを選択します。csvやtxtの場合は、From Text(base)か(readr)で、エクセルファイルの場合にはFrom Excelから読み込むことができます。エクセルで他形式のデータを読み込むときと同様に区切りの識別によって配列を確定させることができます。読み込みも x <- read.table(“data01.txt”)などを使えばスクリプトでデータを読み込むことも可能ですが、Import Datasetを使うとデータが正しく読み込まれているか確認できるので確実性があります。読み込まれたデータセットには名前（拡張子抜きのファイル名）が付けられて、ダブルクリックすると内容を確認することができます。
各ウィンドウの機能（表示されているタブの機能であり、タブを変えると異なる機能となります。）
回帰分析は、lmというコードを使います。最低限の記述で表すと、
結果ファイル <- lm( 目的変数 ~ 説明変数, data=データの名前)
で回帰分析を行うことができます。
- 結果ファイル：結果を格納するデータの名前、何でもよい
- 目的変数：結果の項目（エクセルでいうY範囲）で読み込む表のカラムの名前を入れる
- 説明変数：変数の項目（エクセルでいうX範囲）で目的変数同様に読み込む表のカラムの名前を入れる。複数の変数を取り込む場合には+で変数を足していく
- データの名前：調べるデータの名前で、読み込んだデータの名前を入力する
基本の解析はこれだけです。この一行を実行する（Sourceをクリック）と、解析が行われて結果が入ったデータが作られます。表を確認するときと同様にデータをダブルクリックしても解析結果を確認できますが、もう１行summary(データの名前)と記述するとConsoleに統計量が表示されます。
実際のコード
上記のコードを実行した時の結果、エクセルの結果と全く同じことがわかる
分類変数がある回帰分析
では、エクセルで処理できない解析方法をやってみます。分類変数を取り扱う場合も方法は同じです。ボルトを数えるデータセットには分類変数がなかったので、新たなデータセットEGViolatorsを使います。これは、ニュージャージー州の有料高速道路での車の速度超過を調べたもので、黒人の方がスピードを超過しているという疑惑を解明するために超過速度を人種や運転免許の発行州、測定地点、走行方向などとともに調べた結果です。黒人と他の人種で超過速度に差がないことが証明されていますが、今回は目的を変えて、運転免許の発行州別に速度超過の傾向に差があるのかを調べます。
実際に回帰分析を行う方法ですが、上記と全く同じで、運転免許の発行州であるLicenseを代入するだけです。
OverlimitとLicenseの関係を分析するコード
実行すると運転免許ごとに値が算出されます。これはAL（アラバマ州）の運転免許を基準に他の州の運転免許では超過速度がどう変化したを示しています。式で表すなら超過する速度＝5.67（InterceptのEstimate）+免許のEstimateとなり、例えば、表の真ん中付近のMN（ミネソタ州）の運転免許保有者は5.67+7.83=13.5 マイル/hほど統計的に有意差がある中で速度超過していることになります。発行州によっては、P値が高く有意差がないと言える場合もありますが、ミネソタ州のように有意差がある場合もあります。ただし、Adjusted R-squared（補正R2）: 0.005996 が１に程遠いので、全体としては速度超過と運転免許の発行州は相関が小さいと言えます。
回帰分析の結果、P値が低さによって右端に記号が付与される。
このようにRでは分類変数を組み込んでも問題なく解析されます。結果は値の場合とは異なり、ある基準の分類（上記ではAL）と比べた係数の大小で表されます。
信頼区間の図示
そもそもエクセルでは、回帰分析で得られた結果を実測値をプロットする機能はなく、手作業で係数を拾ってグラフを作る必要があります。しかしRを使えば５行追加するだけでグラフまで作ることができます。ここでは新たなデータセット、1970年代後半におけるボストンの住宅価格をデータセットとして使います。これは様々な条件と町の住宅価格の間に相関があるか調べたものです。
まず散布図を作成するには、
plot(Y値~X値, data=データ名)
で、lmと同様にXとYに代入したい表のカラムを入力し、データ名に取り込んだデータセットの名前を記述するだけです。例えば、AGE：1940年より前に建てられた持ち家の割合＝「古い家の割合」をXとしてMEDV：「住宅価格」（1000ドル単位）の中央値でプロットすると図が右下に表示されます。
X軸にAGE、Y軸にMEDVをプロットした散布図を作成するためのコード
グラフを作成した結果
とてもシンプルでわかりやすいですが、タイトルを入れたりプロットの色を変える場合には、コードを追加してグラフをデコレーションすることができます。
では、回帰分析を行って曲線を追加します。追加するコードは、
abline(回帰分析結果)
で、先ほどのボストンの住宅価格の例を使うと、回帰分析後にablineを入れるだけです。plotによって作られた散布図に追記する形をとるので、このコードより前に散布図が作られていないとエラーになります。
回帰分析を行ってその曲線をグラフに上書きするコード、col=’red’で曲線の色を赤色にしている。
直線を追加したグラフ
では最後に信頼区間を表示させます。
配列名<- data.frame(X軸に使うカラム名 = seq(始点, 終点, 感覚))
信頼区間名 <- predict(回帰分析結果, newdata = 配列名, interval = ‘confidence’, level = 信頼区間の％)
lines(配列名$X軸に使うカラム名, 信頼区間名[, 1])
lines(配列名$X軸に使うカラム名, 信頼区間名[, 2])
lines(配列名$X軸に使うカラム名, 信頼区間名[, 3])
信頼区間を表示させるには、範囲の設定が必要でそれを一行目で行います。次の行で信頼区間を計算し、それをもとに３～５行目で曲線を描きます。[,1]が回帰曲線そのもので[,2]が下限、[,3]が上限に相当します。同様にボストンの住宅価格の例を使うと下記のようになります。
回帰分析後、信頼区間を計算しグラフにプロットするコード
95％信頼区間を追加したグラフ、任意のAGEにおいてMEDVの総結果の95％が緑の間に入ることを示している。
信頼区間の違いを確認するためにZN： 25,000平方フィートを超える区画に分類される住宅地の割合＝「広い家の割合」を使ってプロットしました。
ZNとMEDVの信頼区間を示したグラフ
確かに広い家の割合が多いほど町の住宅価格は上昇しますが、ばらつきも大きくなってZNが80％以上では、95％信頼区間が5000ドル以上に広がっているのがわかります。商品の有効性を示す近似曲線は、日常生活でも多々見かけますが、効果があるように見えても、ばらつきが大きく実際は有効性が疑わしい場合もあると思います。そのため統計的な信頼区間の議論はばらつきが大きい結果ほど重要であると言えます。
もっと簡単に、モダンなグラフを書きたい場合には、ggplot2という追加パッケージが便利です。Rではパッケージを追加することで機能を追加することかでき、具体的には、右下のウィンドウの「Packages」のタブで追加したパッケージ名にチェックを入れるだけ使えるようになります。先ほどのAGEとMEDVの値を使って信頼区間の入ったグラフを作るには、
ggplot2を使って信頼区間を図示するためのコード
の２行で下のようなグラフを作ることができます。
ggplot2を使った作成したグラフ
もちろん、オプションのコードを増やすことでグラフの体裁を変えることもできます。
相互作用がある回帰分析
最後に相互作用がある回帰分析の方法について紹介します。相互作用を調べるときにはX値の書き方を変えるだけです。速度超過のデータセットを使って、超過速度と測定した出口（Exit）と運転免許の発行州（Licence）について関係があるかどうか回帰分析を行う場合には、result <- lm(Overlimit~Exit+License, data=Violators)と記述すると、下記のような結果が得られます。相互作用を考えなければ、超過速度と測定した出口の相関が大きいものの、特定の運転免許の発行州でも有意差がみられています。
ExitとLicenseの回帰分析の結果、補正R2は0.01034である。
そこで、相互作用を考慮した回帰分析を行います。記述をresult <- lm(Overlimit~Exit*License, data=Violatorsとするだけで、Overlimit=a×Exit+b×License+c×Exit×License+dの式の近似を行います。この結果は下記のようになり、運転免許の発行州の有意差が見られるようになった一方、Exitの有意差は見られなくなりました。相互作用の効果を見ると多くの州において負を示しています。この解釈を考えると、ExitのEstimateは正なので、番号が大きいExitほど速度超過が大きいが、多くの特定の発行州の免許保有者はその比例を強めるか弱めることになります。例えば、CO（コロラド州）の免許保有者は、平均よりも超過速度は高いが、番号が大きいExitほど、超過速度はより低くなる傾向があるということになります。全体のモデルの当てはまりを示す補正R2の相互作用を考慮することで若干向上しました。
相互作用を考慮したExitとLicenseの回帰分析の主成分の結果
相互作用を考慮したExitとLicenseの回帰分析の相互作用の結果（一部のみ掲載）
相互作用を考慮したExitとLicenseの回帰分析の結果、補正R2は0.01114である。
三つ以上の変数（例えばAとBとC）を*で接続して記述すると、A,B,C,AB,AC,AC,ABCとすべての相互作用が計算されるので特定の相互作用を除きたい場合には、考慮したい組み合わせを+で合わせます。
相互作用を考慮して当てはまりが良くなるかどうかは、データセット次第であり、当てはまりが悪くなるか条件の有意差が低くなることもあります。また相互作用を算出にはそれなりのデータ数も必要で、速度超過のデータセットでも発行州によっては相互作用の値がNAになっているのは、データ数が少なすぎて算出できなかったためです。そのため相互作用まで考慮する意味があるかどうかは、相互作用ある場合とない場合の回帰分析を比較して判断する必要があります。このように相互作用がある系の回帰分析もRでは簡単に分析することができ、単体の効果は当たり前だが、複数の条件の組み合わせによって予想外の結果が得られることはよくあり事象で、それを統計的に探して発見できることが、この回帰分析の面白さだと思います。
今回は、Rを使った回帰分析を紹介してきました。Rは世界的に認められたソフトウェアであり、FDAに薬のデータを提出する際にも使われています。また、いろいろなパッケージが開発されていて、Rで最新の解析方法を取り扱うこともできます。次回は、より深く回帰分析を使いこなす方法を紹介します。
関連書籍
関連リンク
- エクセルを利用した回帰分析の手順：図がずれているが、信頼区間の算出する式が解説されている
- 付録. M5実習4: 重回帰分析：作図から回帰分析までを手順を追って紹介している
- 有意性検定の考え方：信頼区間に関連する統計量を解説している
