はじめての ChemoSpec

1. Professor Emeritus of Chemistry & Biochemistry, DePauw University, Greencastle IN USA.

hanson@depauw.edu

## Warning in knitr::write_bib(c("knitr", "mclust", "baseline", "hyperSpec"), :
## package(s) hyperSpec not found

この vignette は、 ChemoSpec バージョン 5.3.11 にて作成されています。

ChemoSpecは，核磁気共鳴（NMR），赤外（IR），ラマン，蛍光X線（XRF）などのスペクトルデータをトップダウンで探索的に分析するための関数群である。 スペクトルのプロットや検査、ピークアライメント、階層クラスター分析（HCA）、主成分分析（PCA）、モデルベースクラスタリングなどの機能を備えています。この種の高次元データに適したロバストな手法が利用できます。 ChemoSpec は、メタボロミクスの研究のように、サンプルが治療群と対照群に分類されるような構造化された実験のために設計されています。 グラフィカルな出力は、出版品質のプロットのために一貫してフォーマットされています。 ChemoSpec は非常にユーザーフレンドリーであり、使用可能な結果を素早く得ることができるように意図されています。 典型的な操作方法を説明した vignette を用意しています。

1 はじめに

Varmuza and Filzmoser (Varmuza and Filzmoser 2009) の定義によると、計量化学 (chemometrics) とは

“… 数学的・統計的ツールによる化学データからの関連情報の抽出.”

これは、計量化学的なアプローチで扱える問題や課題の豊富さを考慮すると、適切な広義の定義である。 ここでは、変数（周波数）が多く、サンプル数が比較的少ないスペクトルデータセットに焦点を当てています。 このような多変量、高p、低nのデータセットには、アルゴリズム上の課題がありますが、これまでも知識のある人たちが取り組んできました。 特に、多変量ケモメトリックス分析の実用的かつ理論的な背景については、Varmuza/Filzmoser著(Varmuza and Filzmoser 2009) を強くお勧めします。 ここで紹介されている関数の中には、彼らや他の人たちがRコミュニティのために彼らのパッケージで提供している関数のラッパーに過ぎないものもあります。 もう一つの優れたテキストは、Ron Wehrens氏(Wehrens 2011) によるものです。

ChemoSpecは、XRF(Panchuk et al. 2018)、UV-Vis、NMR、またはMIRやNIRを含むIRデータなどの分光データのケモメトリックス分析のために開発されました。 また、ChemoSpecはクロマトグラフィーデータ（下記参照）や円偏光二色性のようなあまり一般的ではない技術でも動作します。¹ ChemoSpecの目的は、Rを初めて使うような幅広い研究者がケモメトリクスツールを容易に利用できるようにすることです。 このアプローチは完全に探索的で教師なし、つまり「トップダウン」です(Wishart 2007)。 ChemoSpecは、異なる履歴を持つサンプル、つまり、異なるクラス、カテゴリー、グループに分類されるサンプルに対応するように設計されています。 例えば、治療群と対照群、あるいは単に異なる標本（赤い花と青い花）などです。 ChemoSpecは、可能な限りユーザーフレンドリーに設計されており、豊富なエラーチェック、役立つ警告、一貫したインターフェースを備えています。 また、スタイルや注釈に一貫性があり、出版物やポスターに使用するのに適したグラフィックを作成します。 関数のドキュメント作成には細心の注意が払われていますが、この vignette は ChemoSpec によるデータ解析を学ぶための最良の出発点となります。 ChemoSpecは、分光器で通常行われる作業を複製することを意図していません。

ChemoSpec の中心となるのは Spectra オブジェクトです。 これはあなたのデータを保存し、Rで利用できるようにする場所です。 このようにデータを保存してチェックすれば、すべての分析を簡単に行うことができます。 現在、ChemoSpec にはいくつかの組み込みデータセットが同梱されていますが、ここでは SrE.IR というデータセットを使用してデモを行います。

ここでは、ChemoSpec の学習を始めるにあたり、少なくとも R の基本的な知識があり、適切なワークフローが設定されていることを前提としています。 ここで説明されている関数の詳細なヘルプについては、コンソールで ?function_name と入力してください。 また、?ChemoSpecと入力し、下部にあるインデックスリンクをクリックすると、利用可能なすべての関数が表示されます。

パッケージ、関数、関数の引数、データセットなどの R 「オブジェクト」の名前は、ファイル名 と同様に タイプライター フォントを使用しています。 また、コンソールで発行するコマンドや出力は薄いグレーの背景で表示され、用途や目的に応じて色分けされていますが、これは優れた knitr パッケージ (Xie 2021) の提供によるものです。

ところで、もしあなたがChemoSpecを試してみて、便利だと感じたり、質問や意見、提案があれば、ぜひ私に教えてください。 あなたが使っているバージョンには、すでに多くのユーザーの意見が取り入れられていますが、あなたの意見を加えてみませんか？ 考えられるバグや機能の要望は、Github issues システムを使って文書化してください。

2 ワークフローのサンプル

このサンプル探査は、典型的な ChemoSpec のワークフローを説明するために設計されています。 重要なのは、どのようにコマンドを実行するか、どのようなオプションが利用可能で一般的に使用されているか、そして分析を行う順序を説明することです。

これらのコマンドのバージョンをスクリプトファイルに入れて、作業中にソースを作成することができます。 そうすれば，簡単に変更することができますし，すべてを再現することができます。 これを行うには、空の R ドキュメントを開き、コマンドを入力します。 これを、My_First_ChemoSpec.R のような名前で保存します。 そして、その一部を切り取ってコンソールに貼り付けて実行するか、全体をソースにして実行するかのどちらかです。

source("My_First_ChemoSpec.R")

計量化学の典型的なワークフローを図に示します。 データの性質によっては、これらのステップのいくつかが無関係であったり、省略されたり、順序を変更する必要があるかもしれません。 以下のセクションではその例をご紹介します。

Figure 2.1: A typical workflow. For a given data set, some steps may be omitted and the order changed. That is part of what is meant by exploratory data analysis!

3 ChemoSpecへのデータの取り込み

生データセットを ChemoSpec にインポートする方法は2つあります。 1つは関数 files2SpectraObject で、これは生データが1つのディレクトリに別々のファイルとして存在し、各ファイルには周波数列と強度列が含まれていることを想定しています。 ヘッダ行はあってもなくてもよく、データは任意の区切り記号（通常はコンマ，タブ，セミコロン，スペース）で区切ることができます。 また、コンマやピリオドを小数点として使用することもできます。 これらのオプションにより、様々な機器で書かれた様々な規則のファイルからデータをインポートすることができます²。 また，JCAMP-DX 形式のファイルをインポートすることもできます。

2つ目の関数は matrix2SpectraObject で、データの行列を含む1つのファイルがあることを前提としています。 この行列は，1列目に周波数、残りの列に個々のサンプルの強度を持つ必要があります。 ファイルにはヘッダ行があり、そこにはサンプル名が含まれていなければなりません（ただし、周波数を示す最初のエントリは無視されます）。 この要件以外には、上述のすべての柔軟性があります。

詳細は ?files2SpectraObject のヘルプを読んでください。 また、... 引数には特に注意してください。

データファイルの名前は、クラスのメンバーシップを符号化したシステムを使ってつけるのがよいでしょう。 例えば、データセットに治療群と対照群が含まれている場合、あるいは類似したクラス/グループ情報がある場合、この情報はファイル名から得られるはずです。 引数 gr.crit は、ファイル/サンプル名の grep 処理の基礎となり、そこから各サンプルがグループに割り当てられ、同様に色が割り当てられます。 もしサンプルがグループに分類されない場合は、それも良いのですが、それでも gr.crit に何かを与えなければなりません – すべてのファイル名に共通する一つの文字列を与えれば良いのです。 もちろん、このアプローチは、楽器から出てきたファイルをどのように分析するかを考えながらファイル名をつけることを促すものであり、それは実験計画に依存します。 計画を立てることは悪いことではありません。

files2SpectraObject や matrix2SpectraObject の出力は Spectra オブジェクトです。 これは R で言うところのオブジェクトで、データだけでなく、関数の引数で提供されたデータに関するその他の情報も含んでいます。

典型的な例を挙げてみましょう。 例えば、花のエッセンシャルオイルのNMRファイルが30個入ったフォルダがあったとします。 そのうち18個は提案されたある亜種からのもので、12個は別の亜種からのものだったとします。 さらに、調査中の問題は、これらの2つの亜種の分類学に関係していると仮定しましょう。 もし、ファイル名が sspA1.csv …. sspA18.csv とか sspB1.csv …! sspB12.csv のような名前のファイルであれば、以下のコマンドでファイルを処理し、目的の Spectra オブジェクトを作成することができます。:

ssp <- files2SpectraObject(
  gr.crit = c("sspA", "sspB"),
  gr.cols = c("red", "blue"),
  freq.unit = "ppm",
  int.unit = "peak intensity",
  descrip = "Subspecies Study",
  out.file = "subsp")

これにより，files2SpectraObject はファイル名に "sspA" と "sspB" という文字列が含まれているかどうかを調べ，これらを用いてサンプルをグループに分けます。 sspA*.csv のサンプルには赤が、sspB*.csv のサンプルには青が割り当てられます（色の選択を事前に計画するためのいくつかの提案や，?colorSymbol については，Section ‾ref(sec:colsymsec) を参照してください）。 このコマンドを実行すると、ディレクトリに subsp.RData という新しいファイルが作成され、ワークスペースに ssp という Spectra オブジェクトが作成されて探索できるようになります。 後日、データを再度インポートする必要はありません。 保存したバージョンを使用して、以下のように好きな名前を付けることができます（関数 loadObject はパッケージ R.utils のものです）。

SubspeciesNMR <- loadObject("subsp.RData")

これで、すぐに使えるようになりました。

3.1 クロマトグラムを扱う

この vignette とパッケージに含まれるすべての言語は、スペクトルの解析を目的としていますが、ChemoSpecは、生データとしてのクロマトグラムでもうまく動作します。 この場合、時間が周波数の代わりになるのは当然ですが、それ以外の分析はほとんど同じです。 唯一の違いは、files2SpectraObject などでデータをインポートする際に、周波数の単位を以下のように指定することです: freq.unit = "time (minutes)"

3.2 組み込みデータセット

ChemoSpecにはいくつかの組み込みデータセットが同梱されています。 SrE.IR はこのヴィネットで使われているセットです。 これは、男性のBPHの治療によく使われるノコギリヤシ（Serenoa repens_）から抽出されたエッセンシャルオイルの14個のIRスペクトルのコレクションで構成されています。 14個のスペクトルは、異なる小売店のサンプルであり、ラベルの説明に基づいて2つのカテゴリーに分けられています：adSrE（粗悪な抽出物）とpSrE（純粋な抽出物）です。 adSrE（adulterated extract）とpSrE（pure extract）です。 adulterated extract には通常オリーブオイルが添加されていますが、これはBPHには効果がありません。 データセットには，月見草オイル（EPO）とオリーブオイル（OO）の2つのスペクトルが含まれています。 後者の2つのオイルは大部分がトリグリセリドの混合物であるのに対し，SrEのサンプルは大部分が脂肪酸である。 その結果、この2つのグループのスペクトルは異なります。 グリセリドにはエステルカルボニルの伸びがあり、O-Hの伸びはありませんが、脂肪酸には酸カルボニルの伸びとカルボン酸のOHに一致するO-Hの伸びがあります。

また，対応するNMRスペクトルのセットである「SrE.NMR」と，「CuticleIR」も含まれています。 後者は植物_Portulaca oleracea_のキューティクル（葉の表面）の一連のIRスペクトルです。 このデータは、ATRサンプリング装置に葉を静かに挟み込んで撮影したものです。 植物は2つの異なる温度で栽培され、2つの異なる遺伝子型（品種）が使用されました（古典的なG×E、遺伝子型×環境の実験）。

このヴィネットでは、サンプルのスペクトルがかなり異なっており、ほとんどのケモメトリクス手法で良好な分離が得られることから、SrE.IRデータセットを例に挙げています。 一方、「CuticleIR」のスペクトルは、より微妙な違いがあり、結果的に解析が難しくなっています。 これらのデータセットの詳細については、コンソールで ?data_set_name と入力してください。

3.3 データの要約

最初にすべきこと、これは非常に重要なことですが、データが良い状態であることを確認することです。 まず、作成したデータセットを要約し、データの範囲やその他の詳細が期待通りであることを確認します。

data(SrE.IR) # makes the data available
sumSpectra(SrE.IR)

## 
##  Serenoa repens IR quality study 
## 
##  There are 16 spectra in this set.
##  The y-axis unit is absorbance.
## 
##  The frequency scale runs from
##  399.2123 to 3999.837 wavenumber
##  There are 1868 frequency values.
##  The frequency resolution is
##  1.9286 wavenumber/point.
## 
## 
##  The spectra are divided into 4 groups: 
## 
##   group no.   color symbol alt.sym
## 1 adSrE  10 #984EA3     15       d
## 2   EPO   1 #377EB8      2       b
## 3    OO   1 #4DAF4A      3       c
## 4  pSrE   4 #E41A1C      1       a
## 
## 
## *** Note: this is an S3 object
## of class 'Spectra'

sumSpectra はいくつかの情報を提供してくれますが、そのうちのいくつかを紹介します。

4 スペクトルをプロットする

これまでのところ、すべてがうまくいっていると仮定すると、いよいよスペクトルをプロットして検査してみましょう。 すべてのスペクトルにアーチファクトやその他の潜在的な問題がないかチェックするのが良い方法です。 これを行うための関数が3つあります。

• LoopThruSpectra は、大量のスペクトルを検査する手間を省きます。 一度に1つのスペクトルを表示し、先に進む前にコンソールにリターンが入力されるのを待ちます。 詳細はヘルプページを参照してください。
• plotSpectra は一般的な使用と出版物並みのグラフィックを目的としています。
• plotSpectraJS は plotSpectra のインタラクティブバージョンです。 これはウェブページにデータを表示し、必要に応じてスペクトルのオフセットやズームを行うことができます。 非常に大きなデータセットの場合は遅いかもしれませんが、それを避ける方法についてはヘルプページを参照してください。

plotSpectra を使った基本的なプロットを図 4.1 に示します。 ここでは各カテゴリーから1つのスペクトルを選んでプロットしています。 カルボニル領域と C_sp2-H 領域が明らかに異なることがわかります。

# We'll make a fancy title here
# and re-use in other plots
myt <- expression(
  bolditalic(Serenoa)~
  bolditalic(repens)~
  bold(Extract~IR~Spectra))
plotSpectra(SrE.IR,
  main = myt,
  which = c(1, 2, 14, 16),
  yrange = c(0, 1.6),
  offset = 0.4,
  lab.pos = 2200)

Figure 4.1: Sample plot.

データセットの強度範囲や、プロットするスペクトルの数に応じて、引数 yrange, offset, amplify を手動で調整する必要がありますが、通常は数回の反復作業で済みます。 ただし、offset と amplify は関数内で乗算されるので、片方を増やすともう片方を減らさなければならないことに注意してください。 例えば、これらのスペクトルのカルボニル領域だけに焦点を当てたいとしたら、引数 xlim を追加することができる。 実際には、図 4.2 のように、スペクトルの数を少なくして、振幅を大きくして、詳細を見てみましょう．

plotSpectra(SrE.IR,
  main = myt,
  which = c(1, 2, 14, 16),
  yrange = c(0, 0.6),
  offset = 0.1,
  lab.pos = 1775,
  xlim = c(1650, 1800))

Figure 4.2: Detail of the carbonyl region.

これらのサンプルプロットでは、IRスペクトルを2つの方法で表示しています。 まず，x軸が “逆”になっています．これは，もとのスペクトルがもともと周波数の昇順で保存されていたからです（これは必ずしもそうではありません）。 これは、先ほどの例では xlim = c(1800, 1650) のように、xlim 引数を必要な順番で与えることですぐに修正できます。 第二に，これらの例の縦軸のスケールは吸光度である。 IRを構造解明に用いる場合、縦軸は通常、ピークが下を向くような%Tである。 しかし、計量化学には吸光度モードが適しています。 オリジナルのスペクトルをそのように記録し、それに慣れてください。

plotSpectra の引数 which には、プロットしたいスペクトルの整数ベクトルを指定します。 各スペクトルを行列の行に見立てて、列に強度を入れる（各列が特定の周波数値に対応する）とすると、これは行番号と考えることができます。 各行にどのサンプルが入っているのか、どうやって判断するのか気になりますよね。 これには grep コマンドを使うのが最適です。 例えば、オリーブオイルがどの行/サンプルにあるかを知りたい場合、以下のいずれかの方法で見つけることができます。

# if there are only a few spectra
# show all of the names
SrE.IR$names

##  [1] "CVS_adSrE" "ET_pSrE"   "GNC_adSrE" "LF_adSrE"  "MDB_pSrE"  "NA_pSrE"  
##  [7] "Nat_adSrE" "NP_adSrE"  "NR_pSrE"   "NSI_adSrE" "NW_adSrE"  "SN_adSrE" 
## [13] "Sol_adSrE" "SV_EPO"    "TD_adSrE"  "TJ_OO"

# if there are a lot of spectra,
# grep for the desired names
grep("OO", SrE.IR$names)

## [1] 16

5 データ前処理のオプション

データの前処理にはいくつかのオプションがあります。 主なものとしては、データを正規化するかどうか、データをビンに入れるかどうか、データをスケーリングするかどうかなどがあります。 ベースライン補正も典型的な処理で、NMRデータセットによってはアラインメントを行う必要があります。 データのスケーリングはPCAルーチンで扱われているので、Section 7.2 をご覧ください。 Engel et al. (2013) には、前処理についての良い議論があります。 Karakach, Wentzell, and Walter (2009) には ¹H NMR データのエラーソースに関する良い議論があります。

5.1 ベースラインドリフトの補正

ChemoSpec では、ベースライン (Liland and Mevik 2020) の乱れを補正するために、baseline パッケージ内の関数を使用しています。 また、baselineSpectra という関数は、必要に応じてオリジナルと補正後のベースラインを表示することができるので、手法を選択する際に便利です。 図はその典型的な例です。 メソッド rfbaseline は IR スペクトルに適しています。 retC = TRUE は補正されたスペクトルを新しい Spectra オブジェクトに入れるので、今後も使用することができます（実際に使用します）。

SrE2.IR <- baselineSpectra(SrE.IR,
  int = FALSE,
  method = "modpolyfit",
  retC = TRUE)

Figure 5.1: Correcting baseline drift.

5.2 アライメント

¹H NMR のデータでは、スペクトルを揃えることが望ましい場合があります。 これは、いくつかのタイプのプロトンに著しいシフトを引き起こす可能性のある希釈、イオン強度、またはpHの変化を部分的に補うものです。 広くてなだらかなピークを持つスペクトルにはこの問題はありません（UV-VisやIRなど）。 ChemoSpec はこの目的のために clupaSpectra 関数を提供しています。 例を ここ で見ることができます。

5.3 バケット化またはビン化

もう1つの前処理として、ビン化またはバケット化があります。 これは、周波数のグループを1つの周波数値に折りたたみ、対応する強度を合計するものです。 この処理を行う理由は2つあります。

• 大きなデータセットのコンパクト化。 これは歴史的な問題です。 Rのアルゴリズムは非常に高速で、大きなデータセットがあってもそれほど遅くはなりません。
• 上述のように、¹H NMRシフトを補正する。

この例はプロセスを示していますが、IRデータでは必要ありません。:

tmp <- binSpectra(SrE.IR, bin.ratio = 4)
sumSpectra(tmp)

## 
##  Serenoa repens IR quality study 
## 
##  There are 16 spectra in this set.
##  The y-axis unit is absorbance.
## 
##  The frequency scale runs from
##  402.1052 to 3996.945 wavenumber
##  There are 467 frequency values.
##  The frequency resolution is
##  7.71425 wavenumber/point.
## 
## 
##  The spectra are divided into 4 groups: 
## 
##   group no.   color symbol alt.sym
## 1 adSrE  10 #984EA3     15       d
## 2   EPO   1 #377EB8      2       b
## 3    OO   1 #4DAF4A      3       c
## 4  pSrE   4 #E41A1C      1       a
## 
## 
## *** Note: this is an S3 object
## of class 'Spectra'

この結果を、フルデータセットの sumSpectra と比較してみてください（Section 3.3）。 特に、ビン化処理のために周波数分解能が下がっていることに注意してください。 データセットを指定された bin.ratio で分割できるようにするために、（警告を出しながら）いくつかのポイントを落とした後、データポイントは平均周波数とグループ化された強度の合計に置き換えられます。 データの微細構造と bin.ratio によっては、重要なピークが異なる bin に分割されることがあります。 この問題を解決しようとする、より洗練されたビン化アルゴリズムが文献にありますが、現在のところ ChemoSpec には実装されていません (Anderson et al. (2008), De Meyer et al. (2008), Sousa, Magalhaes, and Castro Ferreira (2013))。 上記のようにスペクトルを揃えるのが良いでしょう。

5.4 正規化

正規化は normSpectra 関数で行います。 通常、サンプルの前処理によって濃度に差が生じたデータを正規化します。 例えば、取り扱い中に希釈された体液や、研究対象の生物の生理的状態によって変化したデータなどです。 SrE.IR データセットは、オイル抽出物をATR装置に直接置いて撮影したもので、希釈することはできないので、正規化はあまり適切ではありません。 正規化のオプションについては、ヘルプページを参照してください。 文献には、正規化の問題について多くの有益な議論が掲載されています (Craig et al. (2006) Romano, Santini, and Indovina (2000) Berg et al. (2006) Varmuza and Filzmoser (2009) Zhang et al. (2009))。

6 データセットとスペクトル領域の編集

スペクトルをプロットして検査する過程で、問題のあるスペクトルやサンプルが見つかることがあります。 おそらく、機器によるアーチファクトがあるのではないでしょうか。 また、データセットからサンプルのサブグループを削除して、結果の違いを確認することもできます。

6.1 個別サンプルの削除

特定のサンプル、あるいは特定の条件を満たすサンプルを削除するには、removeSample 関数を使います。 この関数は rem.sam という引数に基づいて grep 処理を行いますので、grep の貪欲さには注意が必要です。 例えば、TD_adSrE というサンプルにアーティファクトがあり、削除する必要があるとします。 コマンドは次のようになります。

noTD <- removeSample(SrE2.IR,
  rem.sam = c("TD_adSrE"))
sumSpectra(noTD)

## 
##  Serenoa repens IR quality study 
## 
##  There are 15 spectra in this set.
##  The y-axis unit is absorbance.
## 
##  The frequency scale runs from
##  399.2123 to 3999.837 wavenumber
##  There are 1868 frequency values.
##  The frequency resolution is
##  1.9286 wavenumber/point.
## 
## 
##  The spectra are divided into 4 groups: 
## 
##   group no.   color symbol alt.sym
## 1 adSrE   9 #984EA3     15       d
## 2   EPO   1 #377EB8      2       b
## 3    OO   1 #4DAF4A      3       c
## 4  pSrE   4 #E41A1C      1       a
## 
## 
## *** Note: this is an S3 object
## of class 'Spectra'

grep("TD_adSrE", noTD$names)

## integer(0)

sumSpectra コマンドでは、セット内のスペクトルが1つ少なくなったことを確認します。 図のように、サンプル名でgrepし直しても、サンプルが見つからないことを確認できます。 grep の最初の引数は、検索するパターンです。 そのパターンが複数の名前にマッチした場合、それらはすべて “キャッチ” されます。 例えば、“SrE” をパターンとして使用した場合、“SrE”は “adSrE” と “pSrE” に出現するので、2つのリファレンスサンプルを除くすべてのサンプルが削除されます。 これは grep 関数自体で事前に確認することができます。

SrE <- grep("SrE", SrE2.IR$names)
# show the name(s) that contain "SrE"
SrE2.IR$names[SrE]

##  [1] "CVS_adSrE" "ET_pSrE"   "GNC_adSrE" "LF_adSrE"  "MDB_pSrE"  "NA_pSrE"  
##  [7] "Nat_adSrE" "NP_adSrE"  "NR_pSrE"   "NSI_adSrE" "NW_adSrE"  "SN_adSrE" 
## [13] "Sol_adSrE" "TD_adSrE"

SrE # gives the corresponding indices

##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 15

これが “grep is greedy” の意味するところです。 このような場合、3つの選択肢があります。

• 問題のあるサンプルを手動で削除することができます（str(SrE2.IR)がその方法のアイデアを与えてくれるでしょう。
• removeSample はサンプルのインデックスも受け付けるので、上記のように grep して、実際に削除したいサンプルのインデックスをメモして、それを rem.sam で使用することができます。
• grepや正規表現に詳しい方は、より高度な検索パターンをrem.samに渡すことができます。

6.2 グループの削除

removeSample は、サンプルの名前（Spectra$names にある）を使って、個々のサンプルを識別し、Spectra オブジェクトから削除します。 また、Spectra$groupsの中の特定のグループに属するサンプルを削除する、removeGroupという関数もあります。

6.3 興味のない周波数を特定して削除する

多くのスペクトルには、分析前に除去すべき領域があります。 それは、¹H NMRの水のピークや、IRのCO₂のピークのような、情報を与えない、干渉するピークかもしれません。 また、スペクトルの中には、あまり情報を持たない領域があるかもしれません。 つまり、ノイズの多いベースラインを提供しているだけで、他には何もないのです。 例えば、IRの約1,800または1,900cm^-1から約2,500cm^-1までの領域は、大気中のCO₂の伸びと、まれに（気をつけて！）アルキンの伸びを除いて、通常はピークがない領域です。

これらの領域を見つけるのはとても簡単で、分光学の知識と検査の問題であるかもしれません。 もう1つの方法は，関数 surveySpectra を使ってスペクトル全体を調べることです。 この関数は、平均値や中央値と同様に、データセット全体の特定の周波数における強度の要約統計量（任意）を計算します。 変化の少ない領域では、平均値/中央値、上/下のサマリーラインが近くなります。 図はこのプロセスを示したものである。 もう一つの方法として、surveySpectra2 があり、これはデータを少し異なる形式で表示します。 Figure 6.2 をご覧ください。

surveySpectra(SrE2.IR,
  method = "iqr",
  main = myt,
  by.gr = FALSE)

Figure 6.1: Checking for regions of no interest.

surveySpectra2(SrE2.IR,
  method = "iqr",
  main = myt)

Figure 6.2: Checking for regions of no interest.

図6.1では、by.gr = FALSE という引数を使って、すべてのグループをまとめています。 また、全体のスペクトルを見てみました。 図6.3では、カルボニル領域だけを見ています。 黒い線は、スペクトル全体の強度の中央値です。 赤い線は四分位範囲の上限と下限を示しており、このデータセットのカルボニル領域が大きく変化していることがよくわかります。

surveySpectra(SrE2.IR,
  method = "iqr",
  main = "Detail of Carbonyl Region",
  by.gr = FALSE,
  xlim = c(1650, 1800))

Figure 6.3: Detail of carbonyl region.

最後に、surveySpectra では、データセットをグループ別に見ることができますが、これは本当はもっと便利なことです。 カルボニル領域をグループ別に見てみましょう（図6.4) 。 なお、グループのうち2つはメンバー数が少なすぎて四分位範囲を計算できないため、警告が出ており、これらは表示されていません。

surveySpectra(SrE2.IR,
  method = "iqr",
  main = "Detail of Carbonyl Region",
  by.gr = TRUE,
  xlim = c(1650, 1800))

## 
## Group EPO has 3 or fewer members
##  so your stats are not very useful...
##  This group has been dropped for display purposes!

## 
## Group OO has 3 or fewer members
##  so your stats are not very useful...
##  This group has been dropped for display purposes!

Figure 6.4: Detail of carbonyl region by group.

ここでは、1800～2500cm^-1の領域に注目してみましょう (図6.5)。

surveySpectra(SrE2.IR,
  method = "iqr",
  main = "Detail of Empty Region",
  by.gr = FALSE,
  xlim = c(1800, 2500),
  ylim = c(0.0, 0.05))

Figure 6.5: Inspection of an uninteresting spectral region.

理論的には、この領域には興味深いピークがないと予想されます。 実際、グループをプールしても、この領域の信号は非常に弱く、存在する唯一のピークは大気中のCO～2～によるものです。 この領域は主にノイズやアーティファクトなので、以下のように関数 removeFreq を使って取り除くことができます。 周波数ポイントが少なくなっていることに注意してください。³

SrE3.IR <- removeFreq(SrE2.IR,
  rem.freq = SrE2.IR$freq > 1800 &
  SrE2.IR$freq < 2500)
sumSpectra(SrE3.IR)

## 
##  Serenoa repens IR quality study 
## 
##  There are 16 spectra in this set.
##  The y-axis unit is absorbance.
## 
##  The frequency scale runs from
##  399.2123 to 3999.837 wavenumber
##  There are 1505 frequency values.
##  The frequency resolution is
##  1.9286 wavenumber/point.
## 
##  This data set is not continuous
##  along the frequency axis.
##  Here are the data chunks:
## 
##    beg.freq end.freq     size beg.indx end.indx
## 1  399.2123 1799.348 1400.136        1      727
## 2 2501.3450 3999.837 1498.492      728     1505
## 
##  The spectra are divided into 4 groups: 
## 
##   group no.   color symbol alt.sym
## 1 adSrE  10 #984EA3     15       d
## 2   EPO   1 #377EB8      2       b
## 3    OO   1 #4DAF4A      3       c
## 4  pSrE   4 #E41A1C      1       a
## 
## 
## *** Note: this is an S3 object
## of class 'Spectra'

ここで、sumSpectra がデータセットにギャップを発見したことに注目してください。 このギャップをデータで見ると、図 6.6 のようになります (sumSpectra はギャップをチェックしますが、プロットは作成しません) ; 数値結果と図の両方が提供されます。

check4Gaps(SrE3.IR$freq, SrE3.IR$data[1,])

Figure 6.6: Identifying gaps in a data set.

##    beg.freq end.freq     size beg.indx end.indx
## 1  399.2123 1799.348 1400.136        1      727
## 2 2501.3450 3999.837 1498.492      728     1505

7 探索的データ解析

7.1 階層クラスター分析

階層クラスター分析 (Hierarchical Cluster Analysis、 以下 HCA) はクラスタリング手法の一つで、サンプル間の「距離」を計算し、デンドログラム (木のような構造; 進化学や系統学ではクラドグラムと呼ばれています)で表示します。 HCAの詳細については、他の文献を参照してください（Varmuza and Filzmoser (2009) の第6章がお勧めです）。 ChemoSpec では、サンプル間の距離を計算するためのあらゆる手法と、クラスターを識別するためのあらゆる手法を利用することができます。 典型的な例を図 7.1 に示します。

HCA <- hcaSpectra(SrE3.IR, main = myt)

Figure 7.1: Hierarchical cluster analysis.

その結果、デンドログラムができました。 縦の目盛りは、サンプル間の数値的な距離を表しています。 意外なことに，化学的に異なることがわかっている2つの基準サンプルは、他のすべてのサンプルとは別に集まっています。 驚くべきことに、様々な純粋な、あるいは混入した油の抽出物は、正確には集まっていない。 関数 hcaScores は、生のスペクトルではなく、PCAの結果を用いて同じような分析を行います。 これについては次のセクションで説明します。

7.2 主成分分析

主成分分析（以降、PCA）は、探索的データ分析の真の主力製品です。 PCAでは、グループメンバーについての仮定はありませんが、結果として得られるサンプルスコアのクラスタリングは、データを理解する上で非常に役立ちます。 PCAの理論と実践については、他でもよく取り上げられています（Varmuza and Filzmoser (2009) の第3章が良いでしょう）。 ここでは、ChemoSpec での PCA メソッドの使用方法に焦点を当てます。 簡単に言うと、PCAとは、データセットを記述するのに必要な最小数の成分を決定することであり、事実上、ノイズや冗長な情報を除去することであると考えることができます。 典型的なスペクトルを考えてみてください。 ある領域は明らかに単なるノイズです。 さらに、典型的なスペクトルのピークは、ピークが上昇し、トップアウトし、ベースラインに戻るように、かなりの数の周波数単位にわたっています。 あるピークの中のどのポイントも、ピークの強さという同じことを表しています。 さらに、あるピークのエンベロープ内の各周波数は、エンベロープ内の他のすべての周波数と相関しています（サンプルごとにピークの大きさが変わると、それらは一体となって上昇・下降します）。 PCAは、データセット内のすべてのノイズと基本的な相関関係を「過去」に見て、データセット全体を本質的な部分にまで煮詰めることができます。 残念ながら、このプロセスで発見された主成分は、通常、具体的な何かに対応するものではありません。 繰り返しになりますが、より詳細な処理を参考にしてください。

ChemoSpecで利用できるオプションの概要と、関連する機能を表 7.1 にまとめました。

Table 7.1: PCA Functions

PCA options	scaling options	function
classical PCA	no scaling, autoscaling, Pareto scaling	c_pcaSpectra
robust PCA	no scaling, median absolute deviation	r_pcaSpectra
sparse PCA	no scaling, autoscaling, Pareto scaling	s_pcaSpectra
IRLBA PCA	no scaling, autoscaling, Pareto scaling	irlba_pcaSpectra
Diagnostics		function
OD plots		pcaDiag
SD plots		pcaDiag
Choosing the correct no. of PCs		function
scree plot		plotScree
bootstrap analysis (classical PCA only)		cv_pcaSpectra
Score plots	plotting options	function
2D plots	robust or classical confidence ellipses	plotScores
3D plots
—static 3D plots		plotScores3D
—interactive 3D plots		plotScoresRGL
Loading plots		function
loadings vs frequencies		plotLoadings
loadings vs other loadings		plot2Loadings
s-plot (correlation vs covariance)		sPlotSpectra
Other		functions
HCA of PCA scores		hcaScores
ANOVA-PCA		aov_pcaSpectra

ここでは非常に多くの選択肢があります。 例を見ながら、選択肢を説明し、少なくとも言及していきましょう。 データをどのように分析するかを決めるのは、あなた自身であることを覚えておいてください。 多くの人は、さまざまな選択肢を試し、最も多くの知見を得られるものに従うでしょう。 しかし、決定権はあなたにあります。

最初のステップは、PCAを実行することです。 PCAには、古典的手法とロバストな手法の2つの主な選択肢があります。 古典的手法では、提供されたすべてのデータを使用してスコアとローディングを計算します。 ロバスト法では、データの中核部分に焦点を当てるため、一部のサンプルに重み付けがされることがあります。 この違いは重要で、データの性質によっては、2つの手法の結果が大きく異なることがあります。 この違いは、PCA法（古典的なものとロバストなものの両方）が、データセットの分散を可能な限り説明する成分を見つけようとするからです。 例えば、サンプルの取り扱いのために純粋に破損したサンプルがある場合、そのスペクトルプロファイルは他のすべてのサンプルと非常に異なる可能性があり、それは正当に外れ値と呼ばれます。 古典的なPCAでは、この1つのサンプルはデータセット全体の分散に強く寄与し、PCAスコアはそれを反映します（スコアと負荷量は外れ値に従うと言われることがあります）。 ロバストPCAでは、中央値絶対偏差のようなロバストな分散測定が用いられるので、異なる特徴を持つサンプルはそれほど大きな影響を受けません。

なお，c_pcaSpectra と r_pcaSpectra のどちらも、サンプルによる正規化は行いません。 サンプルを正規化するかどうかを決定し、正規化する場合には normSpectra を使用する必要があります。

古典的手法かロバスト手法かの選択に加え、スケーリング手法を選択する必要があります。 古典的なPCAでは、スケーリングなし、オートスケーリング、パレートスケーリングのいずれかを選択します。 古典的な分析では、データをスケーリングしない場合、大きなピークが結果に強く寄与します。 オートスケーリングを行うと、各ピークは（ノイズの「ピーク」を含めて）結果に等しく寄与します。 パレートスケーリングは、この2つの間の妥協点です。 ロバストPCAでは、スケーリングを行わないこともできますし、中央絶対偏差に基づいてスケーリングすることもできます。 中央値絶対偏差は、より極端なピークを重み付けする手段である。 測定の種類（機器）や生物学的データセットの性質に適したスケーリングオプションについては、多くの文献で推奨されています（Zhang et al. (2009) Craig et al. (2006) Romano, Santini, and Indovina (2000) Berg et al. (2006) Varmuza and Filzmoser (2009) Karakach, Wentzell, and Walter (2009)）。

ここでは、可能な限りのオプションの組み合わせを説明するための十分なスペースがありません。 図 7.2 と図 7.3 は、 スケーリングなしの古典的PCAとロバストPCAの使用方法と結果を示し、最初の2つのPCをプロットしました。 これらのプロットから、ロバストPCAと古典的手法は、プロットの全体的な外観だけでなく、各PCで説明される分散の量においても、かなり異なる結果をもたらしたことがわかります。

結果を見るためにスコアをプロットしたので、結果の2Dプロットを作成する plotScores のいくつかの機能について触れてみましょう（3Dオプションについては後ほど説明します）。 プロットの左上には、この分析の履歴を説明する注釈が表示されるので、何を見ているのかわからなくなることがないことに注意してください。 tol 引数は、サンプル名でラベル付けされるポイントの割合をコントロールします。 これは潜在的な外れ値を識別するための手段です。 ellipse 引数は、楕円を描くかどうか、どのように描くかを決定します（95%信頼区間が使用されます）。

楕円を描かない場合は "none"、古典的に計算された信頼楕円を描く場合は"cls"、ロバストに計算された楕円を描く場合は"rob"、両者を直接比較したい場合は "both" を選択することができます。 ここで古典的、ロバストという言葉を使っているのは、PCAアルゴリズムとは何の関係もないことに注意してください — 同じ考え方ですが、PCAによって生成されたスコアの2次元配列に適用されています。 楕円の外側にある点は、外れ値の候補となる可能性が高い。

c_res <- c_pcaSpectra(SrE3.IR,
  choice = "noscale")
plotScores(SrE3.IR, c_res,
  main = myt,
  pcs = c(1,2),
  ellipse = "rob",
  tol = 0.01)

## Group EPO
##  has only 1 member (no ellipse possible)

## Group OO
##  has only 1 member (no ellipse possible)

Figure 7.2: Classical PCA scores.

r_res <- r_pcaSpectra(SrE3.IR,
  choice = "noscale")
plotScores(SrE3.IR, r_res,
  main = myt,
  pcs = c(1,2),
  ellipse = "rob",
  tol = 0.01)

## Group EPO
##  has only 1 member (no ellipse possible)

## Group OO
##  has only 1 member (no ellipse possible)

Figure 7.3: Robust PCA scores.

図 7.2 や 7.3 のようなプロットは、外れ値の可能性を示唆するものですが、ChemoSpecにはより洗練されたアプローチが含まれています。 関数 pcaDiag は、参考になる2種類のプロットを生成することができます（図 7.4 と 7.5）。 これらのプロットの意味と解釈については、Varmuza and Filzmoser, Chapter 3 (Varmuza and Filzmoser 2009) で詳しく説明されています。

diagnostics <- pcaDiag(SrE3.IR, c_res,
  pcs = 2,
  plot = "OD")

Figure 7.4: Diagnostics: orthogonal distances.

diagnostics <- pcaDiag(SrE3.IR, c_res,
  pcs = 2,
  plot = "SD")

Figure 7.5: Diagnostics: score distances.

データの内容や結果の解釈によっては、いくつかのサンプルを廃棄すべきだと判断するかもしれません。 その場合は、前述のように removeSample を使用し，PCA 分析を繰り返します。 ほとんどの人にとっての次のステップは、データを記述するのに必要なPCの数を決定することです。 これは通常、図のようなスクリープロットを用いて行います。 しかし、ChemoSpec のデフォルトは別のスタイルのスクリープロットで、私はこれの方がはるかに情報量が多いと考えています（図7.6 と 7.7）。

古典的なPCAを使用している場合は、図 7.8 ブートストラップ法で必要なPCの数を知ることもできます。 この方法は反復的で、少し時間がかかることに注意してください。 これらの結果をスクリープロットと比較すると、ブートストラップ法では、95％レベルに到達するためには、4つまたは5つのPCでは必ずしも十分ではないことを示唆しているのに対し、ScreePlot では、2つのPCで十分であることを示唆していることがわかります。

plotScree(c_res, main = myt)

Figure 7.6: Scree plot.

plotScree(c_res, style = "trad", main = myt)

Figure 7.7: Traditional style scree plot.

out <- cv_pcaSpectra(SrE3.IR,
  pcs = 5)

Figure 7.8: Bootstrap analysis for no. of principal components.

次に、スコアを3Dで表示する方法について説明します。 現在、ChemoSpec には2つのオプションがあります。 一つは lattice グラフィックスを使ったプロットで、手動で調整しなければならない静的なプロットを生成します。 もう一つは rgl パッケージに基づいたインタラクティブなプロットです。 おそらく、plotScoresRGL から始めるのが一番良いでしょう。 これはデータを調べるのに適していて、高品質で印刷することができます。 しかし、openGL グラフィックデバイスの性質上、タイトルや凡例がデータと一緒に移動するので、出版物に適したハードコピーにはならないかもしれません。 このインタラクティブなプロットは、このドキュメントでは呼び出すことができませんが、必要なコマンドは以下の通りです。

plotScoresRGL(SrE3.IR, c_res,
  main = "S. repens IR Spectra",
  leg.pos = "A",
  t.pos = "B") # not run - it's interactive!

完全な詳細については、もちろんマニュアルページの ?plotScoresRGL を参照してください。 もし同様の、そしておそらくより出版物にふさわしいプロットが必要ならば、図 7.9 のように plotScores3D を使うことができます。 この例では ellipse = FALSE を設定していますが、これは “adSrE” のデータポイントが非常に細長い楕円を形成しており、 他のポイントが非常に見にくくなるような形でプロットの限界を強いるからです。

plotScores3D(SrE3.IR, c_res,
  main = myt,
  ellipse = FALSE)

Figure 7.9: Plotting scores in 3D using plotScores3D.

PCAでは、スコアに加えて、各変数（スペクトルアプリケーションでは周波数）がスコアにどのように影響するかを示す負荷量も生成されます。 これらの負荷量を調べることは、結果を解釈する上で非常に重要です。 図 7.10 にその例を示します。 一方、PC2では、炭化水素領域に興味深いピークがあり、いくつかのピークが反対方向に作用しています。 データを実際に分析することがここでの目的ではありませんが、PC1はエステルと酸のカルボニル基に、PC2は飽和と不飽和の脂肪酸鎖（後者は C_sp2-H}$ のピークを持つ）を検出していると思われます。

plotLoadings(SrE3.IR, c_res,
  main = myt,
  loads = c(1, 2),
  ref = 1)

Figure 7.10: Loading plot.

また、関数 plot2Loadings を用いて，ある荷重を別の荷重に対してプロットすることもできます（図 7.11）。 これは、多くの変数が相関しているため（同じピークの一部であるため，図中の蛇行した線のようになっている）、 通常、分光データにはあまり役に立ちません。 しかし、プロット上の最も極端なポイントは、どのピーク（周波数）がPCのペアを区別するのに役立つのか、つまりデータのクラスタリングを推進するのに役立つのかを知ることができます。

res <- plot2Loadings(SrE3.IR, c_res,
  main = myt,
  loads = c(1, 2),
  tol = 0.002)

Figure 7.11: Plotting one loading vs. another.

しかし，より有用なアプローチとして，s-plotを使ってどの変数が最も大きな影響力を持っているかを判断することができます。 標準的な負荷量プロット（plotLoadings）は、どの周波数範囲がどの主成分に寄与しているかを示しますが、このプロットでは縦軸を自由に設定することができます。 Y軸のスケールを見ないと、主成分1などのローディングがすべて同じように寄与しているような印象を受けてしまいます。 関数 sPlotSpectra は、各頻度変数と特定のスコアの相関を、その頻度変数と同じスコアの共分散に対してプロットします。 その結果、最も影響力のある頻度変数が右上と左下の象限にあるs字型のプロットが得られます。 その例を図 7.12 詳細を図 7.13 に示します。 後者の図では、カルボニルピークの影響がよくわかります。 この方法は、Wiklund et al. (2008) で報告されています。

spt <- sPlotSpectra(SrE3.IR, c_res,
  main = myt,
  pc = 1,
  tol = 0.001)

Figure 7.12: s-Plot to identify influential frequencies.

spt <- sPlotSpectra(SrE3.IR, c_res,
  main = "Detail of s-Plot",
  pc = 1,
  tol = 0.05,
  xlim = c(-0.04, -0.01),
  ylim = c(-1.05, -0.9))

Figure 7.13: s-Plot detail.

最後に、PCAとHCAのアイデアをブレンドすることができます。 PCAはデータセットのノイズを除去するので（重要なPCを選択した後）、PCAのスコアに対してHCAを実行することができま。 なぜならスコアはクリーンアップされたデータを表すからです。 SrE.IR データセットを使った結果は、生のスペクトルに対してHCAを行うのと変わらないので、説明は省略しますが、コマンドは次のようになります。:

hcaScores(SrE3.IR,  c_res,
  scores = c(1:5),
  main = myt)

7.3 ANOVA-PCA

Harrington et al. (Harrington et al. 2005) (および他の数名 – Pinto et al. (2008)) は、伝統的なANOVAとPCAを組み合わせた手法を実証しました。 標準的なPCAはクラス・メンバーシップには盲目ですが、一般的には既知のクラス・メンバーシップを使用してスコア・プロットのポイントを着色します。 ANOVA-PCAは、クラス・メンバーシップを使用して，元の中心データ行列をサブマトリックスに分割します。 各サブマトリックスは、特定の因子に対応し、サブマトリックスの行は、因子の各レベルの平均スペクトルで置き換えられます。 元のデータセットは、これらの部分行列の合計に残差誤差を加えたものと考えられます。 残留誤差を各部分行列に戻してから、PCAを実行します。 これを概念的に示したのが、図7.14 と 7.15。

Figure 7.14: aovPCA breaks the data into a series of submatrices.

Figure 7.15: Submatrices are composed of rows which are averages of each factor level.

ChemoSpec において ANOVA-PCA は、 aov_pcaSpectra、aovPCAscores、 aovPCAloadings といった関数によって実装されています。 ここでの考え方は、もしある因子が有意であれば、PC1対PC2のプロットにおいて、PC1に沿った分離があるということです。 SrE.IR "のデータセットには、ANOVA-PCAを実行するのに十分なグループとレベルがありません。 しかし、aov_pcaSpectra のヘルプページには、CuticleIR データセットを用いた例があり、どのように分析を行うかを説明しています。 また、もう一つの便利な関数である splitSpectraGroups も紹介されています。 この関数を使うと、既存のグループ指定を新しい指定に分割することができます。 詳しくは ?aov_pcaSpectra をご覧ください。

7.4 mclust を使ったモデルベースのクラスタリング

PCAとHCAは教師なしの手法で、基礎となるモデルを想定していません。 HCAはスペクトルのペア間の距離を計算し、デンドログラムが完成するまで反復的にグループ化します。 PCAは、分散を最大化する成分を探し出します。 PCAでは、しばしば（そしてChemoSpec でも）グループメンバーシップでコード化されたサンプルを表示しますが、この情報は実際にはPCAでは使用されません。 サンプルグループの分類と発見されたクラスターの間の見かけ上の対応は、計算の観点からは偶然のものですが、 もちろん、これは見つけたいと思っているものです。

mclust はモデルベースのクラスタリングパッケージで，異なるアプローチをとっています (Fraley, Raftery, and Scrucca 2020) 。 mclust はデータセットの中にグループが存在し、それらのグループが多変量正規分布していることを仮定しています。 `mclust は、反復アプローチを用いてデータセット内の様々な可能性のあるグループをサンプリングし、 ベイズ情報量規準(BIC)を用いてどのグループがデータ分布に最も適合するかを決定します。 mclust は、特定の制約に従ったグループを探します。 例えば、ある制約は、見つかったすべてのグループがデータポイントの分布が球形でなければならないというものです。 別の制約は楕円形の分布を許容するものです. 詳細はScruccaとRafteryの論文(Scrucca et al. 2017)を参照してください。 しかし、基本的な考え方は、mclustがデータセットの中のグループを探しに行き、見つけたグループをあなたが知っているグループと比較することができるということです。

ChemoSpec には mclust の機能を拡張するいくつかの関数が含まれています。 mclust では，まずBICを用いて，どのモデルがデータに最もフィットするかを判断し，その結果を図 7.16 に示します． 次に図 7.17 は mclust がデータから見つけたグループを示しています。 図 7.2 のスコアプロットと図 7.17 のmclustの結果を視覚的に比較すると面白いです。 mclustでは、2つの外れ値を残りのデータの一部と一緒にグループ化しているように見えます。 次に、mclustは見つけたグループに真のグループをマッピングします。 誤差のあるポイントはX-アウトされます。 これらの結果は図 7.18 で見ることができます。 このプロットから、mclust はこのデータセットではあまりうまくいっていないことがわかります。 一般的に、mclustのエラーの概念を使用することには非常に注意が必要です。 発見されたグループをアルゴリズム的に「真実」にマッピングすることは非常に困難です。 私はますます mclust の “truth” オプションを使わない方に傾いています。

model <- mclustSpectra(SrE3.IR, c_res,
  plot = "BIC",
  main = myt)

Figure 7.16: mclust chooses an optimal model.

model <- mclustSpectra(SrE3.IR, c_res,
  plot = "proj",
  main = myt)

Figure 7.17: mclust’s thoughts on the matter.

model <- mclustSpectra(SrE3.IR, c_res,
  plot = "errors",
  main = myt,
  truth = SrE3.IR$groups)

Figure 7.18: Comparing mclust results to the TRUTH.

また，mclust3dSpectra を用いて，同様の分析を3次元で行うこともできます。 この関数は、グループを見つけるために mclust を使用しますが、信頼楕円を描くためにmclust ではない関数を使用します。 この関数は rgl グラフィックを使用するので，ここではデモンストレーションを行うことはできませんが、コマンドは次のようになります．

# not run - it's interactive!
mclust3dSpectra(SrE3.IR, c_res)

8 付録

8.1 ここでは説明されていない関数

詳細はヘルプファイルを参照してください。

• splitSpectraGroups: この関数の使用例は ?aov_pcaSpectra にあります。
• hypTestScores: PCA スコアに対して anova を実行します。PCAのスコアを分析します。
• hmapSpectra: 系列化されたヒートマップをプロットします。系列化されたヒートマップをプロットします。
• evalClusters: 様々なクラスタリングオプションを比較します。様々なクラスタリングオプションを比較します。
• sgfSpectra: Savitzky-Golay フィルタを適用する。
• plotSpectraDist: スペクトルをプロットします。各スペクトルとリファレンススペクトルとの距離をプロットする。

8.2 色と記号のオプション

ChemoSpecでは、R に知られている任意の色名/フォーマットを使用することができます。 データをインポートする際、必要に応じて ChemoSpec は自動的に色を選択してくれます。 Spectra オブジェクトの現在の配色は、sumSpectra を用いて決定したり、conColScheme を用いて変更することができます。 色の問題についての詳しい説明は ?colorSymbol にあります。

色に加えて、Spectra オブジェクトはシンボルのリストと代替シンボルも含みます。 これらは白黒でプロットする場合や、色覚異常者がプロットを見る場合に便利です。 代替記号は単純に小文字ですが、これは plotScoresRGL やその他の rgl-graphics 駆動の関数で必要となるもので、従来の記号ではプロットできません。 図 8.1 に内蔵されているオプションの一部を示しますが、前述のように、好きなものを選択することができます．

Figure 8.1: Color and symbol suggestions.

8.3 関連パッケージ

ChemoSpec と同じタスクのいくつかを実行する他のパッケージがいくつか存在します。 ChemoSpec に最も近い機能を持つパッケージは、私の友人であり共同研究者でもある Claudia Belietes によって書かれた hyperSpec です (これらのパッケージは同時期に独立して開発されましたが、最近では私も hyperSpec の開発に貢献しています) (???) 。 また、SpectraオブジェクトとhyperSpecオブジェクトを相互変換するように設計されたパッケージもあり、 これによりパッケージ間でのデータの移動がより簡単になりました (McManus 2018)。 Rのエコシステムでは多くの開発が行われており、新しいパッケージが次々と登場しています（時には消えてしまうことも！）。 これらをチェックするのに良い場所は、多くのオープンソース分光法プロジェクトに関する情報を集めた FOSS4Spectroscopy のウェブサイトです。

8.4 謝辞

ChemoSpec の開発は、私が2007年から2008年のサバティカル期間中に始まり、フィッシャー・フェローシップの授与に大きく助けられました。 これらのプログラムは、DePauw のファカルティ・ディベロップメント委員会によってコーディネートされています。 この委員会と、これらのプログラムを最初に作成した人たちに大変感謝しています。 私の学生研究者の一人であるケリー・サマーズは、2009年の夏、予備調査の一環として、「CuticleIR」に掲載されているデータを取得しました。 また、Peter Filzmoser教授には、彼の chemometrics パッケージのアルゴリズムに関する私の質問に答えていただきました。 最後に、多くの方々がバグを指摘し、新機能を提案し、修正プログラムを提供してくださいました。 そのすべてに感謝しています。 それぞれ、コード、NEWSファイル、Github issues trackerのいずれかに記載されています。

引用文献

Anderson, P. E., N. V. Reo, N. J. DelRaso, T. E. Doom, and M. L. Raymer. 2008. “Gaussian Binning: A New Kernel-Based Method for Processing NMR Spectroscopic Data for Metabolomics.” Metabolomics 4 (3): 261–72.

Berg, R. A. van den, H. C. J. Hoefsloot, J. A. Westerhuis, A. K. Smilde, and M. J. van der Werf. 2006. “Centering, Scaling, and Transformations: Improving the Biological Information Content of Metabolomics Data.” BMC Genomics 7: 15.

Craig, A., O. Cloareo, E. Holmes, J. K. Nicholson, and J. C. Lindon. 2006. “Scaling and Normalization Effects in NMR Spectroscopic Metabonomic Data Sets.” Analytical Chemistry 78 (7): 2262–7.

De Meyer, T., D. Sinnaeve, B. Van Gasse, E. Tsiporkova, E. R. Rietzschel, M. L. De Buyzere, T. C. Gillebert, S. Bekaert, J. C. Martins, and W. Van Criekinge. 2008. “NMR-Based Characterization of Metabolic Alterations in Hypertension Using an Adaptive, Intelligent Binning Algorithm.” Analytical Chemistry 80 (10): 3783–90.

Engel, Jasper, Jan Gerretzen, Ewa Szymanska, Jeroen J. Jansen, Gerard Downey, Lionel Blanchet, and Lutgarde M. C. Buydens. 2013. “Breaking with Trends in Pre-Processing?” TRAC-Trends in Analytical Chemistry 50: 96–106. https://doi.org/10.1016/j.trac.2013.04.015.

Fraley, Chris, Adrian E. Raftery, and Luca Scrucca. 2020. Mclust: Gaussian Mixture Modelling for Model-Based Clustering, Classification, and Density Estimation. https://mclust-org.github.io/mclust/.

Harrington, PD, NE Vieira, J Espinoza, JK Nien, R Romero, and AL Yergey. 2005. “Analysis of variance-principal component analysis: A soft tool for proteomic discovery.” Analytica Chimica Acta 544 (1-2): 118–27. https://doi.org/10.1016/j.aca.2005.02.042.

Karakach, Tobias K., Peter D. Wentzell, and John A. Walter. 2009. “Characterization of the measurement error structure in 1D H-1 NMR data for metabolomics studies.” Analytica Chimica Acta 636 (2): 163–74. https://doi.org/10.1016/j.aca.2009.01.048.

Liland, Kristian Hovde, and Bjørn-Helge Mevik. 2020. Baseline: Baseline Correction of Spectra. https://CRAN.R-project.org/package=baseline.

McManus, Conor. 2018. HyperChemoBridge: Bridge Between hyperSpec & Chemospec. https://github.com/Chathurga/HyperChemoBridge.

Panchuk, Vitaly, Irina Yaroshenko, Andrey Legin, Valentin Semenov, and Dmitry Kirsanov. 2018. “Application of Chemometric Methods to Xrf-Data – a Tutorial Review.” Analytica Chimica Acta 1040: 19–32. https://doi.org/https://doi.org/10.1016/j.aca.2018.05.023.

Pinto, Rui Climaco, Veronique Bosc, H. Nocairi, Antonio S. Barros, and Douglas N. Rutledge. 2008. “Using ANOVA-PCA for discriminant analysis: Application to the study ofmid-infrared spectra of carraghenan gels as a function of concentration and temperature.” Analytica Chimica Acta 629 (1-2): 47–55. https://doi.org/10.1016/j.aca.2008.09.024.

Romano, R., M. T. Santini, and P. L. Indovina. 2000. “A Time-Domain Algorithm for NMR Spectral Normalization.” Journal of Magnetic Resonance 146 (1): 89–99.

Scrucca, Luca, Michael Fop, Thomas Brendan Murphy, and Adrian E. Raftery. 2017. “mclust 5: Clustering, Classification and Density Estimation Using Gaussian Finite Mixture Models.” The R Journal 8 (1): 205–33. https://journal.r-project.org/archive/2017/RJ-2017-008/RJ-2017-008.pdf.

Sousa, S. A. A., Alvicler Magalhaes, and Marcia Miguel Castro Ferreira. 2013. “Optimized Bucketing for NMR Spectra: Three Case Studies.” Chemometrics and Intelligent Laboratory Systems 122: 93–102. https://doi.org/10.1016/j.chemolab.2013.01.006.

Varmuza, K., and P. Filzmoser. 2009. Introduction to Multivariate Statistical Analysis in Chemometrics. CRC Press.

Wehrens, R. 2011. Chemometrics with R: Multivariate Data Analysis in the Natural Sciences and Life Sciences. Springer.

Wiklund, Susanne, Erik Johansson, Lina Sjostrom, Ewa J. Mellerowicz, Ulf Edlund, John P. Shockcor, Johan Gottfries, Thomas Moritz, and Johan Trygg. 2008. “Visualization of Gc/Tof-Ms-Based Metabolomics Data for Identification of Biochemically Interesting Compounds Using Opls Class Models.” Analytical Chemistry 80 (1): 115–22. https://doi.org/10.1021/ac0713510.

Wishart, D. S. 2007. “Current Progress in Computational Metabolomics.” Briefings in Bioinformatics 8 (5): 279–93.

Xie, Yihui. 2021. Knitr: A General-Purpose Package for Dynamic Report Generation in R. https://yihui.org/knitr/.

Zhang, S. C., C. Zheng, I. R. Lanza, K. S. Nair, D. Raftery, and O. Vitek. 2009. “Interdependence of Signal Processing and Analysis of Urine H-1 NMR Spectra for Metabolic Profiling.” Analytical Chemistry 81 (15): 6080–8.

---

1. ChemoSpec は質量分析データセット (MS) のために開発されたものではなく、テストもされていません。 概要は Chemometrics and Computational Physics Task View を参照してください。↩︎

2. 私の経験では、csvファイルのセパレータにコンマがあるとは限らず、また、当然ながら小数点以下のマークに関する慣習は世界各地で少しずつ異なります。 また、ある国に設置された機器がその国の慣習に従うとは限りません。↩︎

3. removeFreq には、除去する周波数を記述する数式も入力できます。 例えば、res <- removeFreq(SrE.IR, rem.freq = low ~ 800) とすると、最小値から800までのすべての周波数を取り除きます。 詳細は ?removeFreq を参照してください。↩︎