三角ダイアグラムにおける組成データの信頼領域の描画

はじめに

地質学の分野では,さまざまな組成データを三角ダイアグラムにプロットすることがあります.このとき,その組成データがグラフの中のどの程度の範囲に分布するのかを知りたいことがあります.従来は,この範囲をフリーハンドで描いたり,平均値と標準偏差から計算した六角形の領域を描くことが行われてきました.しかし,これらの描画方法は,恣意的であったり,統計的に不十分であるという問題がありました (Weltje, 2002).

三成分の組成データに対し,統計的に厳密な信頼領域を描くことはAitchison (1986)(または,Aitchison 2003)やWeltje (2002)によって構築されてきました.ここでは,Weltje (2002)の方法にしたがって,三角ダイアグラム中の信頼領域を求めてみます.信頼領域の計算とグラフの作成には,オープンソースのR (Ihaka and Gentleman, 1996; R Development Core Team, 2006)を用います.RはS 言語に基づいたデータ解析用の環境で,フリーで公開されているにもかかわらず,高機能な統計解析とグラフィックを扱えます.Unix, Macintosh, Windows のどのOS上でも動きます.Rの使い方は,ここここが非常に参考になります.

ここで紹介する方法は,データに0を含むと計算ができません.それは,三角ダイアグラム内の信頼領域を描くために,成分数が3,自由度が2であることを前提としているので,成分に0があると成分数が2や1となり次元が異なるためです.もし,データが0を含み,その0が測定誤差によるものであれば,あらかじめ0を適切な方法で置換することが必要です.

ちなみに,Rで組成データを扱うためのcompositionsというパッケージもあります.

install.packages("compositions")
library(compositions)

plot(acomp(dat))

適切なr (a scaling of the radius) を設定し,

tt <- acomp(dat); ellipses(mean(tt),var(tt),r=2,col="red")

とすると,同じようなグラフが描けます.

多項分布による信頼領域

測定がベルヌーイ試行(Bernoulli trial)を満たすとき,つまり測定に独立性や定常性が認められるとき(モード組成分析など),一つの試料に含まれる各成分の測定結果(組成)は二項分布に従います.そして,その試料の測定結果は,各成分の二項分布を合わせた多項分布となり,組成がとりうる範囲は,母組成の信頼水準$ (1-\alpha)$の信頼領域によって求めることができます($ \alpha$は有意水準).

ある均質なサンプル(地質試料)から得た複数個の試料(薄片など)のモード組成の信頼領域は,そのサンプルが取り得る組成の範囲を示すことになります.また,この方法を同じ試料に複数回適用することによって,分析誤差を見積ることができます.

三成分の多項信頼領域

計算はWeltje (2002)のAppendix A.2によります.まる写しはよくないと思いますので,要点だけにします.三成分($k=3$,自由度$ df=2$)のとき,信頼水準$(1-\alpha)$の信頼領域がとりうる範囲は,

$\displaystyle \chi^2_{(1-\alpha)} = \sum_{i=1}^{k}\frac{(n_i - Np_i)^2}{Np_i}$ (1)

となります.ここで,$ N$は総カウント数,$ n_i$は各成分のカウント数,$ p_i$は各成分の母組成の割合(未知), $ \chi^2_{(1-\alpha)}$は有意水準$\alpha$,自由度2のときの$ \chi^2$の棄却限界値です.$ \sum_{i=1}^{3} p_i = 1$$ \sum_{i=1}^{3} \hat{p}_i = 1$であることを使って,この式を解き,$ p_i$を解きます.

これをRで使うためのスクリプトがmultinomial.Rです.三角ダイアグラムのところで説明したternary.plotも必要になります.

使用例

multinomial.Rをダウンロードして,作業ディレクトリに置きます.

source("multinomial.R")

組成データの例として,Noda (2005)で測定した網走湾の現世海成砂の砂粒組成を用います.このデータをex1.datとして保存します.総カウント数Nが必要ですので,組成データにする前の生データを用意します.

このデータを三角ダイアグラムにプロットすると下のようになります.

source("ternary.plot")
dat <- read.table("ex1.dat", header=T, row.name=1)
postscript(file="ex1.ps", family="Helvetica")
ternary.plot(dat,label.vertices=T, cline=0, cpoint=1)
graphics.off()
ex1.png

次に,これらの測定点の多項分布による95%の信頼領域をプロットします.

ps.options(onefile=TRUE)
postscript(file="multinomial3.ps", family="Helvetica")

for (i in 1:nrow(dat)) {
  res <- multinomial(dat[i,],0.95)
  if (i == 1) {
    ternary.plot(dat[i,],label.vertices=T, cline=0, cpoint=1)
  } else {
    ternary.plot(dat[i,],label.vertices=F, cline=0, cpoint=1, addline=F)
  }
  par(new=T)
  ternary.plot(res, label.vertices=F, cline=1, cpoint=0, addline=F,cline.col="red")
  if (i != nrow(dat)) {
    par(new=T)
  }
}
graphics.off()
multinomial3.png

多変量加法ロジスティック正規分布による信頼領域

もう一つの方法は,1個の測定結果を多変量測定の1点とみなすことです (Aitchison, 1986).つまり,$n$個の標本(組成をもとめた試料の数)によって,一連の組成範囲を予測します.組成の予測範囲は,繰り返し(replicate) とみなされ,その変化パターンは母集団の成分の割合の共分散(covariance) を反映していると考えられます.もし,同一の試料から複数の試料を切り出したなら,母集団は同一と考えます.同様に,同じ地質帯に分布する同種の岩石から採取した複数の試料からも,母集団を求めることができます.

この方法では,多変量加法ロジスティック正規分布 (multivariate additive logistic normal distribution) を用います(Aitchison, 1986).この計算により,信頼水準$ (1-\alpha)$における,$ n$個の一連の測定データの2種類の信頼領域(母集団と母平均)を求めることができます.この方法で求められた信頼領域は,測定誤差を含む試料に関わる組成データのすべての変動を反映しています.また,多項分布による信頼領域では総カウント数(N) が考慮されましたが,この方法では考慮されません.

三成分の加法ロジスティック正規信頼領域

加法ロジスティック正規分布は対数比変換(Aitchison, 1986)に基づいています.この対数比変換は,組成データから定数和制約や非負性を取り除く一般的な方法です(太田・新井, 2006).対数比$y_i$は以下のように定義します.

$\displaystyle y_i = \log\left( {\frac{x_i}{x_k}} \right)= \log{x_i} - \log{x_k},\ $    where,  $\displaystyle i = 1,2,\ldots,k-1$ (2)

ここでは,$ k=3$の成分を分母にとり,$ (x_1,x_2,x_3)$$ (y_1,y_2)$に変換します.計算結果はどの成分を分母にとっても変わりません (Aitchison, 1986).

多変量正規分布の信頼領域は,ホテリングの$ T^2$ (Hotelling's $ T^2$) 分布(母平均が分かっていて,母分散が未知の1変量正規分布)を用います.

$\displaystyle T^2 = n [ \overline{\bm{Y}} - \bm{Y} ]^T \bm{S}^{-1} [ \overline{\bm{Y}} - \bm{Y} ]$ (3)

ここで,$ n$はサンプル数,$ \overline{\bm{Y}}$は標本平均の列ベクトル,$ {\bm{Y}}$ は母平均の列ベクトル(未知),$\bm{S}^{-1}$は標本共分散行列の逆行列となります.

信頼水準$ (1-\alpha)$のときの母平均の信頼領域は,下の式で描かれます.

$\displaystyle [ \overline{\bm{Y}} - \bm{Y} ]^T \bm{S}^{-1} [ \overline{\bm{Y}} - \bm{Y} ] = \frac{m(n-1)}{n(n-m)} F_{(1-\alpha;df_1;df_2)}$ (4)

ここで,$ m$は変数の数(対数比変換した三成分の場合は2),自由度は$ df_1=m$$ df_2=n-m$となります.また,母集団の信頼領域は,下の式のようになります.

$\displaystyle [ \overline{\bm{Y}} - \bm{Y} ]^T \bm{S}^{-1} [ \overline{\bm{Y}} - \bm{Y} ] = \frac{m(n-1)(n+1)}{n(n-m)} F_{(1-\alpha;df_1;df_2)}$ (5)

使用例

logistic.Rをダウンロードして,作業ディレクトリに置きます.
source("logistic.R")

先程のex1.datを用います.この計算で用いるデータには,組成データ(合計を1または100に再計算したデータ)を用いてもかまいません.オプションcrに信頼水準を,methodに母平均 (mean) または母集団 (entire) を 指定します.

logistic(dat,cr=0.95,method=c("entire","mean"))

ex1.datのデータの信頼水準95%の母平均の領域は,以下のように求められます.

res <- logistic(dat,0.95,method="mean")

結果の返り値を適当なオブジェクトに付値します.ここではresに付値していますが,resにはsとYとxyzの3つの返り値を持っています.sは対数比変換したデータ,Yはsの信頼領域,xyzはYを三成分に戻したときの信頼領域です.

まず,Yをプロットしてみます.ex1.datの95%の母平均(赤)と母集団(青)の信頼領域をプロットします.

res.mean <- logistic(dat,0.95,method="mean")
res.entire <- logistic(dat,0.95,method="entire")

ps.options(onefile=TRUE)
postscript(file="logistic.Y.ps", family="Helvetica")
plot(res.mean$s,pch=20,xlim=c(-6,2),ylim=c(-6,6))
par(new=T)
plot(res.mean$Y,type="l",col="red",xlim=c(-6,2),ylim=c(-6,6),xlab="",ylab="",axes=F)
par(new=T)
plot(res.entire$Y,type="l",col="blue",xlim=c(-6,2),ylim=c(-6,6),xlab="",ylab="",axes=F)
graphics.off()
logistic.Y.png

次に,上で描いた楕円を三成分の三角ダイアグラムに投影します.

ps.options(onefile=TRUE)
postscript(file="logistic.xyz.ps", family="Helvetica")
ternary.plot(dat,label.vertices=T, cline=0, cpoint=1)
par(new=T)
ternary.plot(res.mean$xyz,label.vertices=F,cline=1,cpoint=0,cline.col="red",addline=F)
par(new=T)
ternary.plot(res.entire$xyz,label.vertices=F,cline=1,cpoint=0,cline.col="blue",addline=F)
graphics.off()
logistic.xyz.png

参考文献

Aitchison, J. (1986) The Statistical Analysis of Compositional Data
London, Chapman & Hall.

Aitchison, J. (2003) The Statistical Analysis of Compositional Data
Caldwell, New Jersey, The Blackburn Press.

Ihaka, R. and Gentleman, R. (1996) R: a language for data analysis and graphics.
Journal of Computational and Graphical Statistics vol. 5, no. 3. p. 299–314.

Noda, A. (2005) Texture and petrology of modern river, beach and shelf sands in a volcanic back-arc setting, northeastern Japan.
The Island Arc vol. 14, p. 687–707, (abstract).

太田 亨・新井 宏嘉 (2006) 組成データ解析の問題点と解決方法.
地質学雑誌 vol. 112, p. 173–187, (abstract).

R Development Core Team (2006) R: A Language and Environment for Statistical Computing.
R Foundation for Statistical Computing Vienna, Austria.ISBN 3-900051-07-0, (html).

Weltje, G. J. (2002) Quantitative analysis of detrital modes: statistically rigorous confidence regions in ternary diagrams and their use in sedimentary petrology.
Earth-Science Reviews vol. 57, p. 211–253, (abstract).

Back