Microsoft Excelを用いたケモメトリクス計算(4)–主成分回帰–

J. Comput. Chem. Jpn., Vol. 10, No. 1, pp. 32–43 (2011) © 2011 Society of Computer Chemistry, Japan

1 はじめに

実験的に得られた化学データを数学的•統計的手法に より解析する「ケモメトリクス」が，近年盛んに用いら れるようになってきた．しかし，ケモメトリクスに関す る教科書類がいくつか出版 [1–4]されている一方で，大 学等の授業ではほとんど扱われていないのが実情であ る．そこで本シリーズは，より多くの学生や研究者がケ モメトリクスを用いることができるようにするために， 最も普及している計算ソフトMicrosoft Excel (以後Excel と記す)によるケモメトリクス計算法を考案することを 目的としている． 今 回 は， 主 成 分 回 帰 (principal component regression:  PCR) [1,3–5] を Excel ワークシート上で実行する．PCR は， 代 表 的 な 多 変 量 解 析 手 法 で あ る， 主 成 分 分 析 (principal component analysis: PCA) [1–4,6] と， 重 回 帰 (multiple linear regression: MLR)を結びつけた手法である． PCRではMLRの欠点であった多重共線性の問題を解決 できるという特徴がある．ただしPCAの実行にあたっ ては固有値問題を解決する必要があり，3変量以上(3行3 列以上の行列)になると固有値問題を解くことが難しく なってしまう．そこで，多変量解析の教科書などでは， あらかじめコンピュータで計算した出力を示したり，付 属のプログラムやマクロで解決したりすることが一般的 である [7–12]．本シリーズ2回目では，5行5列の行列の 固有値問題をExcelワークシート上で反復法の一つ，累 乗法によって解く方法を示した [13]．この手法を応用し て教科書に載っているPCAの例題や，教科書レベルを 超える"多"変量の実際問題に対してExcelのワークシー ト上で計算の展開•解決ができれば，学習者や研究者が PCAを理解するための強力な補助となると思われる． 本稿ではまず，前報 [14]と同様に吸光スペクトルを用 いた定量分析を想定し，スペクトル生成，PCAの実行， PCR検量モデル作成および定量を行うためのワークシー トを作成する．さらに，MLRやLambert–Beer則に基づ いた分析法 [14](LBA)との定量性の比較を行う．

Microsoft Excelを用いたケモメトリクス計算(4)

–主成分回帰–

吉村 季織

*

，福原 亘治，三ツ木 健一郎，高柳 正夫

国立大学法人東京農工大学農学府, 〒183-8509東京都府中市幸町3-5-8

*e-mail: yosimura@cc.tuat.ac.jp

(Received: October 21, 2010; Accepted for publication: December 1, 2010; Advance publication: March 23, 2011) 近年，化学データを数学的•統計的手法により解析する「ケモメトリクス」が頻繁に用いられるよう になってきた．しかし，日本の大学の化学教育の場ではほとんど取り上げられていない．ケモメトリク スや数値計算の専用ソフトウェアを使うことなく，現在最も普及しているソフトウェアのひとつである Microsoft Excel (Excel)の基本機能を用いてケモメトリクス計算を行うことができれば，多くの教育•研究 機関で役立つものと思われる．シリーズ4回目は，主成分回帰を取り扱う．検量モデル作成用試料および 未知試料のスペクトルの生成，NIPALSアルゴリズムによる主成分分析，得られた主成分スコアによる回 帰(主成分回帰)を行うためのワークシートを作成した．さらに，重回帰およびLambert–Beer則に基づい た分析法との比較を行い，主成分回帰の優位性を示した． キーワード: Microsoft Excel, ケモメトリクス, 主成分分析, 主成分回帰, NIPALSアルゴリズム

技術論文

本稿の内容は，Microsoft Windows上のExcelを想定し ている．そのため，Mac OS 上の Excelにより計算を行 う場合には，キー操作などを適宜読み替えてほしい．ま た，本稿中の図では，Excelのワークシートは列幅を初 期状態より狭い5.0 (45 ピクセル)にしてあるので，その 点に注意して計算値の桁数や最終桁の値などを見てほし い．

2 主成分回帰

2.1 概念

前報 [14]で取り扱ったMLRや，本報で取り扱うPCR は目的変量yを複数の説明変量 xmの線形結合で表される 予測モデル， 1 M m m m

y

f x

=

′ =

∑

   (1) を作成すること，すなわち係数 fmを求めることが目的 である (y'は yの予測値 )． MLR では説明変量の数Mが観測数 Nよりも大きくな る場合，解を得ることができない．また，M £ Nであっ ても多重共線性の問題により解が得られない場合もあ る．そこで適切な変量を選択し，これらの変量を用いて MLRを行わなくてはならない．本シリーズで取り扱っ ているスペクトルの場合などでは，Nに比べ Mが圧倒的 に大きくなる傾向がある．そのためMLRでは得られた 変量の多くが無駄になる．得られたすべての変量を活 用し，かつ多重共線性を回避できる方法が主成分回帰 (principal component regression: PCR)である． MLRとPCRの説明変量の取り扱いの概念的な差を，2 変量の場合を例に考える．変量1 (x1)と変量2 (x2)に関し て，4つの観測を得たとする．観測番号をnとすると，一 観測あたり，変量1の観測値(x1n)と変量2の観測値(x2n)が 得られる．これらの観測値はx1を横軸，x2を縦軸とすれ ば，Figure 1のように2次元上の4点として示すことがで きる (ただしFigure 1ではAppendixで述べている平均中 心化を施してある)．MLRで用いられる説明変量はx1座 標(x1n)とx2座標(x2n)であり，得られた各観測点から各軸 へ垂線を下した交点にあたる{Figure 1 (a)}． Figure 1 (b)は同じ観測点に対し，原点を通り x1軸，x2 軸とは異なる向きのp1軸とそれに直交するp2軸を引いた 場合を示している．観測点nから p1軸，p2軸へおろした 垂線の交点をそれぞれt1n, t2nとする．つまり，t1n, t2nは p1軸，p2軸を基準とした座標系の座標にあたる．そのた め，t1n, t2nを説明変量としてMLRを行うことができる． PCR では，まず PCA により pi軸方向の単位ベクトル pi (第i主成分．本稿では列ベクトルとする)を決定する． ただし piが互いに直交するという条件のもと， tinの二 乗和siが極大になるように piを求めなくてはならない． siが大きいほうから順に i = 1, 2,…とする．このようにし て得られたtinを使って目的変量 ynに対し重回帰を行う． tinをまとめた行ベクトルを [ [ S S W W W W W W W W [ [ [ [ [ [ [ [ [ [     D0/5 E3&$ Figure 1. The difference of axes and explanatory variables between MLR and PCA.

(

1

)

i

=

t

i

t

iN

t



   (2) とする． MLRでは各説明変量を軸とした直交座標系を用いる のに対し，PCRではpiという新しい直交座標系を見つ け，各点の新しい座標を用いて重回帰を行うものであ る． なおPCA, PCRの詳細については他書 [1–4]などを参考 にされたい．

2.2 計算手順

主成分は，固有値問題を解く方法で求めることが多い． 本稿ではnonlinear iterative partial least-squares (NIPALS) [4] で求める．NIPALSは固有値問題をパワー法で解くこと と実質的に同等である． 前報 [14]と同様に吸光スペクトルを例に取り上げる． L種の成分(PCAの主成分と紛らわしいので，化学種と記 すことにする)が存在しているN個の試料を，M個の波 長(l1 ~lM)で測定してスペクトルを得たとする．試料 n のスペクトルをan (列ベクトル)，濃度情報ベクトル(列 ベクトル)をcnとする．スペクトルをまとめた行列 Aを 説明変量，濃度情報をまとめた行列Cを目的変量として PCRにより検量モデルを構築する． 2.2.1 PCA pi, tiを求める． 1.  検量モデル校正用試料(以後校正用試料と記す)を測 定して得られたスペクトル行列 Acal，濃度情報行列

Ccalから平均行列 Aavg, Cavgを求める．これらを使っ

て中心平均化した行列X1, Yを算出する．

(

)

1

=

1 N

= −

avg

X

x



x

A A

, 

(

1 N

)

avg

=

= −

Y

y



y

C C

       (3) 平均中心化されたX1, Y1がPCR計算を進める上での 説明変量，目的変量となる． i = 1より始める 2.  ベクトルpiに適当な値を代入する(ただしゼロベクト ルでないこと)． 3.  tiおよび si求める T i

=

i i

t

p X

   (4) T i i i

s = t t

   (5) 4.  piを求める T T i i i i i

′ = X t

p

X t

   (6) 計算が収束していない場合は，

p

i

=

p

i

′

として3.に戻 る．収束している場合は，5.に進む． 5.  Xi+1を求める 1 i+

=

i

−

i i

X

X p t

   (7) i = i+1として2.に戻る． このようにして順次piを求めていく．M > Nならば pi は多くてi = Nまでしか求まらない．これ以上求めよう としても，Xiの各要素が0 (もしくはほぼ0)になってしま うため，もはやpiを正しく求めることができなくなる． また，化学種の数が少なく(L < N)，スペクトルにノイズ が含まれていない場合などは，i = Lまでしか求まらない． 2.2.2 PCR 検量モデルを作成する． 6.  pi, tiをまとめた行列 P, Tを作成する． 7.  Yと Tの回帰係数行列 Qを求める． Yの予測値Y′をTに回帰した式  Y QT′ =   (8) を立てる．係数Qを重回帰により算出する． T₍ T₎−1 = Q YT TT    (9) 8.  係数行列Fを求める．  _{F QP=} T₍₁₀₎ 9.  定数ベクトルf0を求める． 0

=

avg

−

avg

f

c

Fa

   (11) ゆえに，スペクトルaと濃度情報 cの関係として， 0

=

+

c Fa f

   (12) というモデル式が得られる．

Figure 2. Preparation of a worksheet (Sheet1) for generating the spectra of calibration and unknown samples. 

The ranges surrounded by dashed lines have to be filled by values following this Figure. The values in the ranges surrounded by bold  solid lines are calculated by formulas shown in Table 1.

Table 1. Procedures to construct worksheet (Sheet1) for simulating spectra.

  address formula or value etc input style, 

constraints defined name 1 C2:E4 fill this ranges with values following Figure 2. but you can  change them into optional values under constraints. > 0 (C4:E4) ≥ 0 (other cells)    2 I3:M5 Conc 3 P3:T4   4 T7 cun3max 5 T9 noise 6 I3:M4 naming only    Ccal 7 P3:T5   Cun

8 C5:E5  = 4*LN (2)/C4:E4^2 array*1  

9 B12:B32 2980 to 3020 with 2 intervals in this example. But optional.    

10 C12:E32  = C2:E2*EXP (-C5:E5*(B12:B32-C3:E3)^2)

array*1 K

11 I12:M32  = MMULT (K,Conc) Acal

12 N3:N4 input a formula " = AVERAGE(I3:M3)" into N3.   Then copy N3 and paste N4 and N12:N32  normal* 2 cavg 13 N12:N32 aavg 14 I8:M9  = Conc-cavg array*1 Y 15 I35:M55  = Acal-aavg X 16 P5:T5  = cun3max*RAND ()   17 P12:T32  = MMULT (K,Cun)+noise*(RAND ()-0.5) Aun *1 array input: when inputting an array formula in cell(s), keep pressing [Ctrl] and [Shift] key, and then hit [Enter] key.   *2 normal input: when inputting a formula in a cell, hit only [Enter] key.

2.2.3 定量 10. 濃度情報未知の試料(以後未知試料と記す)のスペク トルaunから，式(12)を用いて，濃度情報の予測値 c を計算する． un un 0

′ =

+

c

Fa

f

   (13) Nun個の未知試料のスペクトルをまとめた行列 Aunを 同時に計算する場合，f0をNun列並べた行列 F0を用意 して， un un 0

′ =

+

C

FA

F

   (14) とすることで未知試料の予測濃度行列Cを求めるこ とができる．

3 ExcelによるPCR計算

Excelの新規ブックにあらかじめ用意されているSheet1 をスペクトル生成，Sheet2をPCA計算，Sheet3はPCRモ デル作成および定量計算に用いる．

3.1 スペクトル生成

Sheet1を Figure 2に示す．各化学種のピーク波長や半 値幅は前報 [14] と同じである．14∼ 30行，37 ∼ 53行 は非表示になっているので注意してほしい．作成手順は Table 1に示した．これまでのシリーズと同様，太線で囲 まれた範囲はTable 1に従って数式を入力すること．ほと んどの場合，配列入力([Ctrl]と[Shift]を押したままの状 態で [Enter]を押す，配列数式の入力法．詳しくは既報  [15]を参照のこと)による入力になっている．数式内で 用いられる参照には，行列計算の意味の見通しを考慮し て範囲名を使っている．太破線に囲まれた範囲は制約条 件の範囲で自由に数値を設定してかまわないが，まずは Figure 2に示されている数値を入力することで計算結果 が同じ値になることを確認してほしい． Sheet1では検量モデル校正試料のスペクトルAcal，未 知試料のスペクトル Aunを生成する．校正試料の濃度 Concおよび未知試料の濃度 Cunには化学種3 (化学種とし て想定されていない汚染種)は含まれていないのものと した．平均ベクトル cavgと aavgを求める際，配列入力を しないように注意すること．平均を各行ごとに求める必 要があるので，化学種1の平均濃度をN3に求め，それを コピーしたのちN4およびN12:N32に張り付けるとよい． 校正用試料のスペクトルAcalにはノイズの設定はでき

ない(もしAcalにRAND関数を用いたノイズを含ませる

と，セルへの入力などの作業のたびにAcalの内容が変化 することになる．つまり，Sheet2で反復計算を1段階進め ようとするとAcalそのものが変化するので，反復計算が 収束せず，p1が求まらなくなる)．未知試料中の化学種 3の濃度は乱数により発生できるようにした．この最大 濃度(cun3max) はT7にて設定できる．未知試料スペクトル Aunには，T9で設定される強度のノイズを含ませること ができる．ノイズは0を中心として，設定値の幅内で等 確率に発生する．化学種1, 2, 3および校正用試料のスペ クトルを，それぞれFigure 3, 4に示した．             &6  &6  &6 

Ȝ

Figure 3. The spectra of the chemical species 1, 2 and 3. *1:  chemical  species  1,  *2:  chemical  species  2  and  *3:  chemical species 3.                 

Ȝ

Figure 4. The spectra of 5 samples for calibration.

3.2 PCA

Figure 5に PCA 計算を行うための Sheet2を示した．第

1主成分，p1を求めるための初期状態になっている．4 ∼ 20行，27 ∼ 45行は非表示になっている．作成手順は Table 2に示した． Sheet2を用いて反復計算を行い，値を収束させること でpiを求める．Excelによる反復計算法は既報 [13]にて 示してあるので，ここでは簡単に解説する．(1) Figure 5 の状態まで完成させたら，pであるI2:I22を選択してコ ピーする．(2)piの先頭であるG2を選択し，値貼り付け ( 例 :[Alt]®[E]®[S]®[V]®[Enter]) を 行 う．(3) 値 が 収 束 するまで[F4]キーを押し続ける．収束は数式バーを観察 することで判定できる． 収束したらpiと tiを，後述するようにSheet3にコピー し，結果を保存しておく．B27:F47に計算されている Xi+1をコピーし，B2:F22に値貼り付けする．piを初期化 し，再度反復計算を行う．今回の例ではp2まで求める と，X3はすべての要素がほぼ0になるため計算を終了す る．試料数は5であったのでp5まで求まる可能性があっ たが，各試料中に存在する化学種が2種であるうえ，ス ペクトルにノイズが含まれていないので，p2までしか求 まらなかった．化学種3を含めたり，スペクトルにノイ ズを設定したりするとp3以降も求まるようになる． H24のsiは計算に使われないので求めなくても構わな い．siを求めておくと，順次 piを求めていくにつれて si が小さくなっていくことがわかる．これはsiの大きい順 にi = 1, 2,…としていることの確認にもなる．今回の例で は，s1 ( = 1.463) > s2 ( = 0.886) > s3 (»0)となった． Table 2. Procedures to construct worksheet (Sheet2) of initial state for PCA.

  address formula or value input style, etc referred equation defined name

1 G2:G22 1 (G2), 0 (other cells) initial values of pi   pi 2 B2:F22  = X arrray   Xi 3 B24:F24  = MMULT (TRANSPOSE (pi),Xi)  (4) ti 4 H2:H22  = MMULT (Xi,TRANSPOSE (ti))  (6) Xt 5 I2:I22  = Xt/SQRT (SUMSQ (Xt))  (6)   6 B27:F47  = Xi-pi*ti  (7)   7 H24  = SUMSQ (ti) normal  (5)   Figure 5. Preparation of a worksheet (Sheet2) for getting principal  component vectors (pi), and score vectors (ti). This Figure shows the initial state for getting p1. The iteration calculation  is proceeded by copying I2:I22 and pasting to G2:G22 as value. Then  press [F2] key until reaching convergence.

Figure 6. Preparation of a worksheet (Sheet3) for constructing PCR model and predicting the  concentration of unknown samples.

The values of p1, p2 t1 and t2 are copied from Sheet2, after convergence on each step. The 

PCR model is constructed from 1st to 31st lows. The concentrations of unknown samples are  predicted in 32nd to 35th lows.

Table 3. Procedures to construct worksheet (Sheet3) for PCR and predicting the concentration of the unknown samples.

  address   formula or value, etc input style referred 

equation defined name 1 C2:G3  = Y   array     2 C5:G6 C5:G5 copy form B24:F24 on Sheet2 after the calculation  of the first component is converged, and paste to this  range as value.     T C6:G6 copy form B24:F24 on Sheet2 after the calculation of  the second component is converged, and paste to this  range as value. 3   B10:C30   B10:B30  copy form G2:G22 on Sheet2 after the calculation  of the first component is converged, and paste to this  range as value.     P C10:C30 copy form G2:G22 on Sheet2 after the calculation of  the second component is converged, and paste to this  range as value. 

4 D10:E11  = MMULT (TRANSPOSE (P),P) array  

5 F10:G11   = MMULT (T,TRANSPOSE (T))  array  

6 I2:J3  = MMULT (MMULT (Y,TRANSPOSE (T)),MINVERSE (MMULT 

(T,TRANSPOSE (T))))  array  (9) Q

7 J8:AD9  = MMULT (Q,TRANSPOSE (P)) array   (10) F

8 J10:J11  = cavg-MMULT (F,aavg)  array (11)  f0

3.3 PCR

PCR は Figure 6に示したように Sheet3の31行目までで 行った．作成手順はTable 3に示してある． 主成分 piとスコア tiは， pi+1とスコア ti+1を求めると 上書きされてしまうので，Sheet2にてpi, ti, が求まるご とに Sheet3へ貼り付けを行う．p1, t1が求まったら，そ れ ぞ れ Sheet3の B10:B30, C5:G5に値貼り付けする．さ らに，Sheet2にて p2, t2が求まったら，それぞれ Sheet3 の C10:C30, C6:G6に値貼り付けする．p1, p2をまとめた B10:C30を Pとする．同様に， t1, t2をまとめた C5:G6を Tとする．D10:E11には PTP，F10:G11には TTTを求めて ある．これらを通常求める必要はないが，対角行列にな ることを明確にするために求めておいた．このことは， PCRにて説明変量となるスコアtiが互いに直交している ことを確認できる．すなわち，Tを説明変量，Yを目的 変量とした重回帰は多重共線性の問題が起こらないこと を示している． 式 (9)のようにYを Tに回帰して回帰係数が Qを求め る．Pと Qが求まれば，定量の係数 Fを式(10)で求める ことができる(J8:AD9)．次いで，式(11)を用いて定数項 f0が求められる(J10:J11)．

3.4 定量計算

Sheet3 (Figure 6)の34, 35行では未知試料のスペクトル から，各化学種の濃度の推定計算を行っている．Table 3 にあるように，C34:G35に " = MMULT(F,Aun)+f0" と配列入力すればよい． 今回用いた校正用試料のスペクトル Acalと未知試料 Aunは汚染成分やノイズを含ませていなかったので，算 出されたC′_unは Sheet1のP3:T4で設した値と等しくなっ た．

4 RCR とMLRおよびLBAとの定量

性能比較

汚染やノイズのない2化学種系のデータを，2つの説明 変量(スコア)を使って解析したので，正しい結果を得る ことができた．これは，前報 [14]で示したようにMLR やLBAでも同様であった．そこで非想定成分(化学種3) による汚染や，スペクトルにノイズがある場合の定量計 算を行い，PCR, MLR, LBAの定量性能を比較した． 校正用試料中の化学種3の濃度 conc3がそれぞれ (0.03,  0.06, 0.09, 0.06, 0.03)，未知試料中の化学種3の最大濃度 c3max = 0.01, 未知試料スペクトルのノイズ強度noise = 0.1 という3つの条件を組み合わせた計算を行った．それぞ れの条件にて1000回試行を行った．ただし，乱数に関し てはあらかじめ1000回試行分の乱数を確保しておき，条 件を変えても乱数の発生内容を繰り返し再現できるよう にした．つまり，条件を変えても，j番目の試行時には 同じ乱数となる．これは，計算の再現性を得るためと， 計算法や条件による結果の差を明確にするためである． 算出された予測濃度C′_unの真の濃度 Cunに対する相対誤 差の平均値（average of relative errors : ARE）を次式， 1000 2 5 cal ,  un ,  cal ,  1 1 1

ARE

10000

ln i ln i ln i i l n

c

c

c

= = =

−

=

∑∑∑

      (15) にて求めた．ここで iは i番目の乱数のセットを用いた 試行を示す．同様の操作を MLR, LBA に関しても行っ た．PCR は第 2 主成分まで用いた場合 (PCR2)，第 3 主 成分まで用いた場合 (PCR3) に関して計算した．MLR に 関しては 2 波長 (2994, 3004) における吸光度を用いた場 合 (MLR2)，3 波長 (2994, 3000, 3004) における吸光度を 用いた場合 (MLR3) について計算した．ただし，校正試 料に汚染がない場合，2 成分系となるので，3 成分系で ある PCR3, MLR3 では検量モデルが構築できない．ま た，LBA は想定している化学種が 2 種なので，2 成分系 のみのモデル化になる．これらの結果を Table 4に示す． Table 4の太線よりも上段の数値は計算に用いた濃度や未 知試料スペクトルのノイズの条件である．下段の数値は 式 (15) で求めた ARE である． AREは値が小さいほど定量性能と言える．まず，校正 用試料に汚染種がない場合について考察する．MLRは LBA と比較して，未知試料中の汚染に強く，未知試料 スペクトルにノイズがある場合に弱いという結果は前報  [14] で示された通りである．また，PCR2と LBA で結果 が等しくなっているという特徴的な結果がみられる．こ れは，校正用試料のスペクトルが2成分系で完全にモデ ルを構築できるためであり，この二法によって作成され るモデルが等しい(言い換えると，p1とp2の線形結合で 化学種1と2の純スペクトルを再現できる)ためである． 一方，校正用試料に汚染がある場合，定量性能に大き な差は見られず，わずかにMLR2が他法に勝っていた． 3化学種系となったために，もはや PCR2と LBA の結果 は等しくならず，PCR2のほうが高性能であった．未知 試料に汚染がなく，スペクトルにノイズがある場合より も，汚染もノイズもある場合において定量性能が高かっ た．これは，汚染があることを前提としてモデルが作成

されているため，汚染がない場合にモデルが対応できて いないためと考えられる． 校正用試料に汚染がある場合は，3成分系としてモデ ル化できるため，PCR3やMLR3での定量が可能である． Table 4からも明らかなように，PCR3とMLR3はPCR2と MLR2に比べて格段に定量性能が向上している．未知試 料スペクトルにノイズがない場合ならば，正しい定量 が可能である．ノイズがある場合でも性能向上は顕著で あった．さらにPCR3はMLR3と比較してもさらに高性 能であった．また，未知試料スペクトルにノイズがある 場合は，未知試料中の汚染の有無にかかわらず同じ結果 となった(同じ結果が得られたのは，再現性のある乱数 を用いたためである)．PCR3やMLR3では3成分系を完全 にモデル化できているため，化学種3の有無も想定範囲 内である．したがって，モデル化できていないのはノイ ズ部分である．化学種3の有無にかかわらず同じノイズ が発生するため，モデルで説明できない部分(つまり誤 差)が同じになったものと考えられる．そのため，より ノイズに弱いMLR3の方がPCR3に比べて低性能となっ たといえる． MLRでは選択する説明変量を変更•追加することで， 定量性能を改善することが可能であるが，どの変量を選 ぶかは困難な問題である．これに対し，PCRでは主成 分を増やすという機械的な操作によって定量性能の向上 が期待できる．このことからも，PCRのMLRに対する 優位性が示される．しかし，PCRであっても主成分数 を増やしすぎるとノイズまでモデル化してしまい，逆 に定量性能が低下することがある．Table 4の最右列は 未知試料スペクトルだけでなく，校正用試料スペクト ルにもノイズが含まれている場合の検量性能について

示した．具体的には，Acalにあたる Sheet1の I12:M32に

" = MMULT(K,Conc)+noise*(RAND()-0.5)" と い う 数 式 を 配列入力する．Sheet2での計算が収束できるようにする ため，I12:M32をコピーし，そのままI12:M32に値貼り付 けを行う．Sheet2, 3での作業はこれまでと同様である． 今回の例では，PCAにて第4主成分まで求まった．Figure  4には，第1から4主成分を使って未知試料の濃度予測を 行ったときの AREを示してある．これまでと同様に， PCR2よりもPCR3の方がAREが小さくなった．しかし， 第4主成分まで用いたPCR4の場合，AREは減少せず増 加した．すなわち，未知試料の定量には，第3主成分ま でを用いて作成したモデルを利用した方が適切であるこ

とが分かる．PCRでは校正用試料のAcalと Ccalの関係を

モデル化しており，主成分数を増やすほどCcalをより説 明できるようになる．ただし，過剰に主成分数を増やす と，例えば今回の例のように偶然的な要素を含むノイズ までをモデル化し，校正用試料に対して過剰に適合した モデルができてしまう．そのため，このモデルによって 未知試料のAunから Cunを予測しても良い結果が得られ なくなる．そこで，いくつまでの主成分を用いることが 適当であるかの検証が必要である．検証法としては，校 正用試料とは別に濃度既知の検証用試料を用意し，主成 Table 4. Comparison of quantitative performance among PCR, MLR and LBA  conditions 

conc3 non non non fix*1 fix*1 fix*1 fix*1  cun3max 0.1 0 0.1 0.1 0 0.1 0.1 noise 0 0.01 0.01 0 0.01 0.01 0.01*7 ARE*2  PCR2*3 0.0922 0.0044 0.0923 0.0574 0.1101 0.0576 0.0565 MLR2*4 0.0877 0.0095 0.0881 0.0553 0.1089 0.0561 -LBA 0.0922 0.0044 0.0923 0.0604 0.1132 0.0605 -PCR3*5 - - - 0.0000 0.0069 0.0069 0.0088 MLR3*6 - - - 0.0000 0.0193 0.0193 -PCR4*8 - - - 0.0095  *1 conc3 is fixed, (0.03 0.06 0.09 0.06 0.03).  *2 average of relative errors calculated by eq. (15).   *3, *5 and *8 PCR using 2, 3 and 4 principal components, respectively.   *4 MLR using 2 wavelengths: 2994 and 3004.   *5 PCR using 3 principal components.   *6 MLR using 3 wavelengths: 2994, 3000 and 3004. *7 errors are included in not only Aun, but also Acal.

分数を増やしながら濃度を予測して，最適な(例えば今 回のようなAREが最小になる)主成分数を求める手法が ある．また一般にしばしば用いられている検証法とし て，検証用試料を準備することなく，校正用試料だけで 検証を行うクロスバリデーション (cross validation: CV)  [4]がある．

5 主成分 p

i

とスコアプロット

Figure 1に示されるように，主成分piは一般に説明変 量を基準としたどの軸(xm軸)とも異なる方向を示してい る．つまりすべてのxm軸座標方向の成分を含んでいる． そのため，スペクトル anの pi軸に対する座標である tin は，特定の説明変量ではなく説明変量全体を考慮した総 合的な値になる． piはスペクトル anと同じ M次元空間の中にあるベクト ルである．よってanと同様にスペクトルとして図示する ことが可能である．Figure 7はp1, p2をスペクトル様に示 したものである．p1のピーク位置は化学種1と同じ2994 であった．一方，p2では化学種2と同じ3004となってい た．すなわちp1は化学種1，p2は化学種2と強く関連のあ る軸であると予想される．しかし，p1が波長3000以降で 負になっていることからもわかるように，特定の化学種 のスペクトルそのものではないので，主成分 piは定性 的な傾向を見る指標とすべきである．今回の例では，化 学種1と化学種2が明確に区別できるような結果となっ たが，明確ではないことが多いので，主成分がどのよう な傾向を示しているかを考察することが重要になってく る． Figure 8は横軸にp1のスコア，縦軸にp2のスコアを取っ てプロットしたスコアプロットである．平均中心化を施 しているため，全試料のスペクトルの平均が原点にあた る．すなわち，各試料が平均から見てどのような傾向を 示しているかを知ることができる．今回，横軸(p1)は化 学種1，縦軸(p2)は化学種2と関係していると考えられる ことを考慮すれば，例えば試料1には化学種1が平均より 少なく，化学種2は多いことがわかる． このように主成分やスコアを用いて定性的な解釈を行 うことも可能である．

6 おわりに

PCAの実行，PCRによる検量モデル作成をExcelで行 う方法を示した．PCAは多変量解析の代表的な手法の一 つであるため，PCAを身につけていればケモメトリクス に限らず様々な場面において役立つはずである．PCAに よって求められたスコアを用いて MLRを行うPCRは， 機械的に説明変量を見出すことができ，多重共線性を起              &6  &6  S S

Ȝ

Figure 7. Comparison between the spectra of pure chemical  species and principal component vectors.

Traces  p1  and  p2  are  similar  to  the  spectra  of  chemical 

species 1 and 2, respectively, but do not coincide with them.  Especially, p1 takes negative values at over 3000 wavelength. *1: chemical species 1, *2: chemical species 2.               S S Figure 8. Score plot of calibration samples. The value near each plot is the sample number. The horizontal  axis corresponds to p1, while the vertical axis corresponds to p2.

こすことなく高性能モデルが作成できることを示した． PCR では適切な主成分数を決める必要があるだけでな く，Table 4で示されているようにMLRの方が良いモデ ルができる場合もあるなど検証が必要となる．しかし， 説明変量の数が試料数を超える場合など，多重共線性を 起こすとき，MLRでは変数選択を行わなければモデル を作成することができないのに対し，PCRは確実にモデ ルが得られるという利点がある．ケモメトリクスの学習 者には原理と共にこの点をぜひ理解してほしいと思う． • Excel ファイルに関して 本シリーズで示したExcelのブックファイルは， http://www.tuat.ac.jp/~mt2459/chemom/ にてダウンロードできるので利用してほしい． • 平均中心化 試料n (n = 1, 2, …, N)の波長lm(m = 1, 2, …, M)における 吸光度amnをまとめた行列を Aとする．同様に試料 n中 に存在する化学種l (l = 1, 2, …, L)の濃度clnをまとめた 行列をCとする．また，試料 nに関する吸光度をまと めたベクトル(スペクトル)をan，濃度をまとめたベク トルをcnとすると，

(

1

)

11 1 1 N N M MN

a

a

a

a









=

_{= }

_









A

a

a













, 

(

1

)

11 1 1 N N L LN

c

c

c

c









=

_{= }

_









C

c

c





 





   (i)

となる．A, Cの行要素の平均をそれぞれ aavgm, cavgl

avg 1

1

N m mn n

a

a

N

=

=

∑

, _avg 1

1

N l ln n

c

c

N

=

=

∑

   (ii) とする．全行をまとめて平均ベクトル aavg, cavg avg1 avg avgM

a

a









= 











a



,  avg1 avg avgL

c

c









= 











c



   (iii)

とする．anから aavg  を，cnから cavgを引いたベクトル

をそれぞれyn, xnとする．

avg

n

=

n

−

x

a a

, 

y

_n

=

c c

_n

−

_avg   (iv)

この操作を平均中心化と呼ぶ．aavg, cavgを N列まとめた

行列を Aavgと Cavg

(

)

avg

=

avg N columns avg

A

a

･･ ･･

a

_,

(

)

avg

=

avg N columns avg

C

c

･･ ･･

c

_(v) とし，平均中心化を行ったベクトルをまとめて行列 X1 と Yを定義する．

(

)

1

=

1 N

= −

avg

X

x



x

A A

,  

(

1 N

)

avg

=

= −

Y

y



y

C C

   (vi)

参考文献

  [1]  T. Hasegawa, Supekutoru Teiryô Bunseki, Kodansha  Scientific Ltd. (2005)

  [2]  T. Mitsui, Kemometorikkusu no Kiso to Ôyô, - Bunseki

Kagakuto Tahenryô Kaisekihô –, IPC Inc. (2003)

  [3]  Y.  Ozaki, A.  Uda, T. Akai,  Kagakusha no tameno

Tahenryô Kaiseki, Kodansha Scientific Ltd. (2002)

  [4]  Y. Miyashita, S. Sasaki, Kemometorikkusu - Kagaku

Patân Ninsiki to Tahenryô Kaiseki –, Kyotitsu Shuppan 

Co., Ltd. (1995)

  [5]  W. F. Massy, J. Am. Stat. Assoc, 60, 234 (1965).[Medline] 

[CrossRef]

  [6]  S. Wold, K. Esbensen, P. Geladi, Chemometrics Intell.

Lab. Syst., 2, 37 (1987). [CrossRef]

  [7]  K. Arai, S. Ishimura, Tahenryô Kaiseki no Hanashi,  Tokyo Tosho Co., Ltd. (2000)

  [8]  O.  Uchida,  Suguwakaru EXCEL niyoru Tahenryô

Kaiseki, Tokyo Tosho Co., Ltd. (2000)

  [9]  T. Yoshimuta, Y. Aoyama, Toranjisuta Gijutu Special

Gijutsuhsa no tameno Excel Katsuyô Kenkyû,  CQ 

Publishing Co., Ltd (2002)

 [10]  M. Ogura, Excel de Kantan Tahenryô Kaiseki, Kodansha  Scientific Ltd (2006)

 [11]  T. Suga, Excel de Manabu Tahenryô Kaiseki Nyûmon,  Ohmsha Ltd. (2007)

 [12]  Y. Nagata, M. Muneuchi, Tahenryô Kaiseki Nyûmon,  Saiensu–Sha Co. Ltd. (2009)

 [13]  N. Yoshimura, A. Shigetani, M. Takayanagi, J. Comput.

 [14]  N. Yoshimura, A. Shigetani, M. Takayanagi, J. Comput.

Chem. Jpn., 9, 109 (2010). [CrossRef]

 [15]  N. Yoshimura, M. Takayanagi, J. Comput. Chem. Jpn., 7,  71 (2008). [CrossRef]

Chemometrics Calculations with Microsoft Excel (4)

– Principal Component Regression –

Norio YOSHIMURA

*

, Koji FUKUHARA, Kenichiro MITSUKI and Masao TAKAYANAGI

Graduate School of Agriculture, Tokyo University of Agriculture and Technology,  

3-5-8 Saiwaicho Fuchu, Tokyo 183-8509 Japan

* yosimura@cc.tuat.ac.jp

Although chemometrics has become widely used recently for analyzing experimental chemical data, there  exist only a few instructions for the proper usage of chemometrics other than those in some introductory books.  As the fourth step of chemometrics calculations with Microsoft Excel (Excel), the principal component regression  is performed on worksheets. Three worksheets were constructed for generating the spectra of model calibration  samples and unknown samples, solving principal component analysis by NIPALS algorithm and calculating principal  component regression. The quantitative performance of principal component regression was compared with that of  multiple linear regression or the analysis based on Lambert-Beer law. Principal component regression was found to be  superior to the other two methods. Keywords: Microsoft Excel, Chemometrics, Principal component analysis, Principal component regression and  NIPALS algorithm