[スポンサーリンク]

一般的な話題

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

[スポンサーリンク]

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

関連書籍

[amazonjs asin=”B086KQN4ZV” locale=”JP” title=”これならわかる化学のための統計手法: 正しいデータの扱い方”] [amazonjs asin=”4758120951″ locale=”JP” title=”Rをはじめよう生命科学のためのRStudio入門”]

関連リンク

Avatar photo

Zeolinite

投稿者の記事一覧

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

関連記事

  1. 留学せずに英語をマスターできるかやってみた(1年目)
  2. 軽くて強いだけじゃないナノマテリアル —セルロースナノファイバー…
  3. MEDCHEM NEWSと提携しました
  4. スルホキシイミンを用いた一級アミン合成法
  5. 2010年ノーベル化学賞予想―海外版
  6. 【5月開催】第八回 マツモトファインケミカル技術セミナー 有機金…
  7. マンチニールの不思議な話 ~ウィリアム・ダンピアの記録から~
  8. Actinophyllic Acidの全合成

注目情報

ピックアップ記事

  1. アニリン版クメン法
  2. ジアゾカップリング diazocoupling
  3. スコット・ミラー Scott J. Miller
  4. 「新規高活性アルコール酸化触媒 nor-AZADOの有用性」 第1回 Wako 有機合成セミナー オンデマンド配信を開始! 富士フイルム和光純薬
  5. 蛍光と光増感能がコントロールできる有機ビスマス化合物
  6. 眼精疲労、糖尿病の合併症に効くブルーベリー
  7. 分子マシンー2016年ノーベル化学賞より
  8. [(オキシド)フェニル(トリフルオロメチル)-λ4-スルファニリデン]ジメチルアンモニウムテトラフルオロボラート:[(Oxido)phenyl(trifluoromethyl)-lambda4-sulfanylidene]dimethylammonium Tetrafluoroborate
  9. 相撲と化学の意外な関係(?)
  10. ピラーアレーン

関連商品

ケムステYoutube

ケムステSlack

月別アーカイブ

2020年9月
 123456
78910111213
14151617181920
21222324252627
282930  

注目情報

最新記事

2024年ノーベル化学賞ケムステ予想当選者発表!

大変長らくお待たせしました! 2024年ノーベル化学賞予想の結果発表です!2…

“試薬の安全な取り扱い”講習動画 のご紹介

日常の試験・研究活動でご使用いただいている試薬は、取り扱い方を誤ると重大な事故や被害を引き起こす原因…

ヤーン·テラー効果 Jahn–Teller effects

縮退した電子状態にある非線形の分子は通常不安定で、分子の対称性を落とすことで縮退を解いた構造が安定で…

鉄、助けてっ(Fe)!アルデヒドのエナンチオ選択的α-アミド化

鉄とキラルなエナミンの協働触媒を用いたアルデヒドのエナンチオ選択的α-アミド化が開発された。可視光照…

4種のエステルが密集したテルペノイド:ユーフォルビアロイドAの世界初の全合成

第637回のスポットライトリサーチは、東京大学大学院薬学系研究科・天然物合成化学教室(井上将行教授主…

そこのB2N3、不対電子いらない?

ヘテロ原子のみから成る環(完全ヘテロ原子環)のπ非局在型ラジカル種の合成が達成された。ジボラトリアゾ…

経済産業省ってどんなところ? ~製造産業局・素材産業課・革新素材室における研究開発専門職について~

我が国の化学産業を維持・発展させていくためには、様々なルール作りや投資配分を行政レベルから考え、実施…

第51回ケムステVシンポ「光化学最前線2025」を開催します!

こんにちは、Spectol21です! 年末ですが、来年2025年二発目のケムステVシンポ、その名…

ケムステV年末ライブ2024を開催します!

2024年も残り一週間を切りました! 年末といえば、そう、ケムステV年末ライブ2024!! …

世界初の金属反応剤の単離!高いE選択性を示すWeinrebアミド型Horner–Wadsworth–Emmons反応の開発

第636回のスポットライトリサーチは、東京理科大学 理学部第一部(椎名研究室)の村田貴嗣 助教と博士…

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

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