[スポンサーリンク]

一般的な話題

化学研究で役に立つデータ解析入門:回帰分析の応用編

[スポンサーリンク]

化学研究で役に立つデータ解析入門:エクセルでも立派な解析ができるぞ編ではエクセルでも立派な回帰分析ができることを紹介しました。しかし、回帰分析は奥が深く、エクセルの機能では、できないこと、できても手間がかかることがたくさんあります。今回は、フリーソフトウェアの統計解析向けのプログラミング言語Rを使ってエクセルよりも簡単により深く回帰分析を使いこなす方法を紹介します。

エクセルの回帰分析ではできないこと

回帰分析は奥が深く、エクセルではできない使い方や計算されない統計量がたくさんあります。化学実験の結果を解析するにあたっては下記の機能が特に有用であり、今回はこれらに焦点を当てて解析を行っていきます。

分類変数の取り扱い:実験の条件は、温度や反応時間といった連続した数字(連続変数)だけでなく、基質の官能基や触媒の種類といった数で表すことができない条件(カテゴリー変数)もあり、これらの傾向を掴むことも重要な解析となります。しかしエクセルでカテゴリー変数を取り扱うには、各種類(例えば官能基)ごとに1か0で表現(ダミー変数)して解析する必要があり、解析前のデータの加工が必要になります。

信頼区間:信頼区間とは母集団の真の値が含まれることがかなり確信できる数値範囲のことで、この”かなり確信できる”は、ばらつきの状況によって変わり、95%(95%信頼区間)が一般的です。具体的を95%信頼区間を表すと、ある条件で100回実験を行ったら95回の結果が含まれる結果の範囲であり、予測の誤差を直感的に知るのに便利な統計量です。エクセルでも各係数における信頼区間は算出されていますが、結果を予測する際に重要なのは任意の条件における信頼区間であり、それは画像検索するとたくさんヒットするラッパ型の曲線を算出することで、あるxにおいてyがどれくらいの幅があるのかを知ることができます。エクセルの回帰分析で算出される統計量でもグラフを作ることは可能ですが、手作業でグラフを作る必要があります。

交互作用:例えばある合成において反応温度を上げれば収率は単調増加するが、同時に触媒量を増やすことで収率の増加量が加速する場合、反応温度と触媒量に正の相互作用があると言えます。逆に反応温度を上げれば収率は単調増加するが、同時に基質の量を増やすと収率の増加量が抑えられる場合。反応温度と基質量に負の相互作用があると言えます。エクセルで相互作用を調べるには単純に、条件同士を掛け合わせた新しい条件を追加して回帰分析を行えば良いのですが、相互作用ごとに行を追加すること必要で、条件が多い場合には不向きです。

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をクリック)と、解析が行われて結果が入ったデータが作られます。表を確認するときと同様にデータをダブルクリックしても解析結果を確認できますが、もう1行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 が1に程遠いので、全体としては速度超過と運転免許の発行州は相関が小さいと言えます。

回帰分析の結果、P値が低さによって右端に記号が付与される。

このようにRでは分類変数を組み込んでも問題なく解析されます。結果は値の場合とは異なり、ある基準の分類(上記ではAL)と比べた係数の大小で表されます。

信頼区間の図示

そもそもエクセルでは、回帰分析で得られた結果を実測値をプロットする機能はなく、手作業で係数を拾ってグラフを作る必要があります。しかしRを使えば5行追加するだけでグラフまで作ることができます。ここでは新たなデータセット、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])

信頼区間を表示させるには、範囲の設定が必要でそれを一行目で行います。次の行で信頼区間を計算し、それをもとに3~5行目で曲線を描きます。[,1]が回帰曲線そのもので[,2]が下限、[,3]が上限に相当します。同様にボストンの住宅価格の例を使うと下記のようになります。

回帰分析後、信頼区間を計算しグラフにプロットするコード

95%信頼区間を追加したグラフ、任意のAGEにおいてMEDVの総結果の95%が緑の間に入ることを示している。

信頼区間の違いを確認するためにZN: 25,000平方フィートを超える区画に分類される住宅地の割合=「広い家の割合」を使ってプロットしました。

ZNとMEDVの信頼区間を示したグラフ

確かに広い家の割合が多いほど町の住宅価格は上昇しますが、ばらつきも大きくなってZNが80%以上では、95%信頼区間が5000ドル以上に広がっているのがわかります。商品の有効性を示す近似曲線は、日常生活でも多々見かけますが、効果があるように見えても、ばらつきが大きく実際は有効性が疑わしい場合もあると思います。そのため統計的な信頼区間の議論はばらつきが大きい結果ほど重要であると言えます。

もっと簡単に、モダンなグラフを書きたい場合には、ggplot2という追加パッケージが便利です。Rではパッケージを追加することで機能を追加することかでき、具体的には、右下のウィンドウの「Packages」のタブで追加したパッケージ名にチェックを入れるだけ使えるようになります。先ほどのAGEとMEDVの値を使って信頼区間の入ったグラフを作るには、

ggplot2を使って信頼区間を図示するためのコード

の2行で下のようなグラフを作ることができます。

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で最新の解析方法を取り扱うこともできます。次回は、より深く回帰分析を使いこなす方法を紹介します。

関連書籍

関連リンク

Zeolinite

投稿者の記事一覧

ただの会社員です。某企業で化学製品の商品開発に携わっています。社内でのデータサイエンスの普及とDX促進が個人的な野望です。

関連記事

  1. 化学系面白サイトでちょっと一息つきましょう
  2. 「新反応開発:結合活性化から原子挿入まで」を聴講してみた
  3. マテリアルズインフォマティクスでリチウムイオン電池の有機電極材料…
  4. 「全国発明表彰」化学・材料系の受賞内容の紹介(令和元年度)
  5. 日本のお家芸、糖転移酵素を触媒とするための簡便糖ドナー合成法
  6. 【22卒就活イベント(東京・大阪)/修士1年 技術系職種志望者対…
  7. アルコールのカップリング、NHC塩がアルとおコール
  8. セライトのちょっとマニアックな話

コメント、感想はこちらへ

注目情報

ピックアップ記事

  1. 市販の新解熱鎮痛薬「ロキソニン」って?
  2. ケミカル・ライトの作り方
  3. 2009年10大化学ニュース
  4. ヴィドマン・ストーマー シンノリン合成 Widman-Stoermer Cinnoline Synthesis
  5. デーブナー・フォン=ミラー キノリン合成 Doebner-von Miller quinoline synthesis
  6. 生物学的等価体 Bioisostere
  7. 京大北川教授と名古屋大学松田教授のグループが”Air Liquide Essential Molecules Challenge”にて入賞
  8. 超多剤耐性結核の新しい治療法が 米国政府の承認を取得
  9. カストロ・ステファンス カップリング Castro-Stephens Coupling
  10. ノーベル化学賞:下村脩・米ボストン大名誉教授ら3博士に

関連商品

ケムステYoutube

ケムステSlack

月別アーカイブ

2020年9月
 123456
78910111213
14151617181920
21222324252627
282930  

注目情報

最新記事

十全化学株式会社ってどんな会社?

私たち十全化学は、医薬品の有効成分である原薬及び重要中間体の製造受託を担っている…

化学者と不妊治療

これは理系の夫視点で書いた、私たち夫婦の不妊治療の体験談です。ケムステ読者で不妊に悩まれている方の参…

リボフラビンを活用した光触媒製品の開発

ビタミン系光触媒ジェンタミン®は、リボフラビン(ビタミンB2)を活用した光触媒で…

紅麹を含むサプリメントで重篤な健康被害、原因物質の特定急ぐ

健康食品 (機能性表示食品) に関する重大ニュースが報じられました。血中コレステ…

ユシロ化学工業ってどんな会社?

1944年の創業から培った技術力と信頼で、こっそりセカイを変える化学屋さん。ユシロ化学の事業内容…

日本薬学会第144年会付設展示会ケムステキャンペーン

日本化学会の年会も終わりましたね。付設展示会キャンペーンもケムステイブニングミキ…

ペプチドのN末端でのピンポイント二重修飾反応を開発!

第 605回のスポットライトリサーチは、中央大学大学院 理工学研究科 応用化学専…

材料・製品開発組織における科学的考察の風土のつくりかた ー マテリアルズ・インフォマティクスを活用し最大限の成果を得るための筋の良いテーマとは ー

開催日:2024/03/27 申込みはこちら■開催概要材料開発を取り巻く競争や環境が激し…

石谷教授最終講義「人工光合成を目指して」を聴講してみた

bergです。この度は2024年3月9日(土)に東京工業大学 大岡山キャンパスにて開催された石谷教授…

リガンド効率 Ligand Efficiency

リガンド効率 (Ligand Efficacy: LE) またはリガンド効率指数…

実験器具・用品を試してみたシリーズ

スポットライトリサーチムービー

PAGE TOP