[スポンサーリンク]

一般的な話題

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

[スポンサーリンク]

化学研究で役に立つデータ解析入門:エクセルでも立派な解析ができるぞ編ではエクセルでも立派な回帰分析ができることを紹介しました。しかし、回帰分析は奥が深く、エクセルの機能では、できないこと、できても手間がかかることがたくさんあります。今回は、フリーソフトウェアの統計解析向けのプログラミング言語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. 第100回有機合成シンポジウム記念特別講演会に行ってきました
  2. 【速報】2013年イグノーベル化学賞!「涙のでないタマネギ開発」…
  3. 産官学の深耕ー社会への発信+若い力への後押しー第1回CSJ化学フ…
  4. 人工DNAを複製可能な生物ができた!
  5. SNS予想で盛り上がれ!2023年ノーベル化学賞は誰の手に?
  6. IKCOC-15 ー今年の秋は京都で国際会議に参加しよう
  7. 有機合成化学協会誌2021年12月号:人工核酸・Post-com…
  8. 顕微鏡で有機分子の形が見えた!

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

注目情報

ピックアップ記事

  1. 創薬・医療系ベンチャー支援プログラム”BlockbusterTOKYO”選抜プログラムへの参加チーム募集中!
  2. 市川アリルシアナート転位 Ichikawa Allylcyanate Rearrangement
  3. 文具に凝るといふことを化学者もしてみむとてするなり⑥:実験室でも長持ち「ステンレス定規」の巻
  4. Reaxys Prize 2017ファイナリスト発表
  5. 高反応性かつ取扱い容易な一酸化炭素の代用試薬,N-ホルミルサッカリン
  6. Advanced Real‐Time Process Analytics for Multistep Synthesis in Continuous Flow
  7. EDTA:分子か,双性イオンか
  8. 有機化学を俯瞰する –古代ギリシャ哲学から分子説の誕生まで–【前編】
  9. 日本のお家芸、糖転移酵素を触媒とするための簡便糖ドナー合成法
  10. TED.comで世界最高の英語プレゼンを学ぶ

関連商品

ケムステYoutube

ケムステSlack

月別アーカイブ

2020年9月
 123456
78910111213
14151617181920
21222324252627
282930  

注目情報

最新記事

期待度⭘!サンドイッチ化合物の新顔「シクロセン」

π共役系配位子と金属が交互に配位しながら環を形成したサンドイッチ化合物の合成が達成された。嵩高い置換…

塩基が肝!シクロヘキセンのcis-1,3-カルボホウ素化反応

ニッケル触媒を用いたシクロヘキセンの位置および立体選択的なカルボホウ素化反応が開発された。用いる塩基…

中国へ行ってきました 西安・上海・北京編①

2015年(もう8年前ですね)、中国に講演旅行に行った際に記事を書きました(実は途中で断念し最後まで…

アゾ重合開始剤の特徴と選び方

ラジカル重合はビニルモノマーなどの重合に用いられる方法で、開始反応、成長反応、停止反応を素反応とする…

先端事例から深掘りする、マテリアルズ・インフォマティクスと計算科学の融合

開催日:2023/12/20 申込みはこちら■開催概要近年、少子高齢化、働き手の不足の影…

最新の電子顕微鏡法によりポリエチレン分子鎖の向きを可視化することに成功

第583回のスポットライトリサーチは、東北大学大学院 工学研究科 応用化学専攻 陣内研究室の狩野見 …

\脱炭素・サーキュラーエコノミーの実現/  マイクロ波を用いたケミカルリサイクル・金属製錬プロセスのご紹介

※本セミナーは、技術者および事業担当者向けです。脱炭素化と省エネに貢献するモノづくり技術の一つと…

【書籍】女性が科学の扉を開くとき:偏見と差別に対峙した六〇年 NSF(米国国立科学財団)長官を務めた科学者が語る

概要米国の女性科学者たちは科学界のジェンダーギャップにどのように向き合い,変えてきたのか ……

【太陽ホールディングス】新卒採用情報(2025卒)

■■求める人物像■■「大きな志と好奇心を持ちまだ見ぬ価値造像のために前進できる人…

細胞代謝学術セミナー全3回 主催:同仁化学研究所

細胞代謝研究をテーマに第一線でご活躍されている先生方をお招きし、同仁化学研究所主催の学術セミナーを全…

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

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

PAGE TOP