Rによる珪藻群集の分析

大塚泰介(滋賀県立琵琶湖博物館)

 環境技術特集「珪藻類と環境」内の記事「これから珪藻群集の環境指標性を研究する人のために」で紹介した統計解析方法を、フリーソフト R を用いて実行する手順と、サンプルデータを掲載する。

下準備

 以下のリンクにあるサンプルデータ(カンマ区切りテキストファイル)およびRのプログラム(テキストファイル)をダウンロードする。

サンプルデータとプログラム

ダウンロードした sample_program.zip を解凍し、全てのファイルを1つのフォルダ内に(サブフォルダを置かずに)格納する。このフォルダが、R で分析を行う際の作業ディレクトリになる。以下、D ドライブの Documents の下に"珪藻群集の分析"というフォルダを置いたものとする。

 R がコンピュータにインストールされていない場合には、R をインストールする。方法については、例えば以下のサイトを参照。
http://www2.accsnet.ne.jp/~miwa/ja/R/setupReasy/

 R を起動し、画面上部のメニューから、パッケージ → パッケージのインストール を選択する。
 ダウンロードサイトを選び、必要なパッケージを一覧から選んでインストールする。
 今回の分析で必要なパッケージは、"vegan" と "rioja" である。

 作業ディレクトリを、RGui の、ファイル→ディレクトリの変更 により設定。R-console で以下のように入力してもよい。

> setwd("D:/Documents/珪藻群集の分析") # 作業ディレクトリの指定

 サンプルデータを読み込む。

> env <- read.csv("environment.csv", header=TRUE, row.names=1) # 環境データの読み込み > spe <- read.csv("species.csv", header=TRUE, row.names=1) # 種組成データ(殻数)の読み込み > spe2 <- read.csv("species_relative.csv", header=TRUE, row.names=1) # 種組成データ(相対度数)の読み込み > spe_new <- read.csv("species_testdata.csv", header=TRUE, row.names=1) # GLR-ML用のテストデータ(相対度数)の読み込み

 読み込んだデータを確認するには fix 関数を用いる。例えば R-console で
fix(env)
と入力すると、"env" として読み込まれた "environment.csv" のデータがデータエディタ内に表示される。

制約付き序列化

 パッケージ "vegan" を読み込む。

> library(vegan) # パッケージ "vegan" の読み込み

 DCA を実行する。ここでは DCA そのものの結果ではなく、環境勾配の長さを得ることが目的である。なお、全サンプルの5%以下にしか出現しない種の重みを下げている。

> dca <- decorana(downweight(spe, fraction=5)) # 殻数データに対するDCA の実行 summary(dca) # DCA の結果の概要を表示

 計算結果の最初の方を見ると、

Call:
decorana(veg = downweight(spe, fraction = 5)) 

Detrended correspondence analysis with 26 segments.
Rescaling of axes with 4 iterations.
Downweighting of rare species from fraction 1/5.

                  DCA1   DCA2    DCA3    DCA4
Eigenvalues     0.5171 0.2470 0.10646 0.13482
Decorana values 0.5268 0.1779 0.05234 0.02198
Axis lengths    3.6010 2.3215 1.98153 1.62675

このうち、DCA1 の Axis length に注目する。この値が3以上であれば CCAが、1.5以下であれば RDA がお勧めである。1.5~3の範囲であればどちらでもよい。
 今回は DCA1 の Axis length が 3.601 なので、まずは CCA を試みる。

 CCA を実行する。ここでも、全サンプルの5%以下にしか出現しない種の重みを下げている。また、説明変数に対する変数選択を行っている。

cca.0 <- cca(downweight(spe, fraction=5) ~ 1, env) # NULLモデル(説明変数なしでモデル作成;全サンプルの5%以下にしか出現しない種の重みを下げている) cca.all <- cca(downweight(spe, fraction=5) ~ ., env) # FULLモデル(すべての説明変数を用いてモデル作成) cca.select<- ordistep (cca.0, scope = formula (cca.all), direction = 'both', perm.max = 999) # 変数増減法でモデル選択 

 最終的に、"env" に含まれていた説明変数のうち、v3 と v1 のみが有意な説明変数として選択された。

Step: downweight(spe, fraction = 5) ~ v3 + v1 Df AIC F Pr(>F) - v1 1 198.03 3.4034 0.005 ** - v3 1 203.06 8.9462 0.005 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Df AIC F Pr(>F) + v4 1 196.65 1.6223 0.130 + v7 1 197.29 1.0383 0.365 + v2 1 197.41 0.9363 0.500 + v8 1 197.55 0.8104 0.610 + v6 1 197.71 0.6644 0.800 + v5 1 197.81 0.5739 0.865 

結果をトリプロットで図示する。

> plot(cca.select, type = "text") # CCA の結果をトリプロットで図示

すると、以下の図が出力される。

 続いて RDA を実行する。ここでは、平方根変換した種の比率のデータを用いる。また、CCA の場合と同様に,説明変数に対する変数選択を行っている。

spe_sqrt <- sqrt(spe2) rda.0 <- rda(spe_sqrt ~ 1, env) # NULLモデル(説明変数なしでモデル作成) rda.all <- rda(spe_sqrt ~ ., env) # FULLモデル(すべての説明変数を用いてモデル作成) rda.select <- ordistep (rda.0, scope = formula (rda.all), direction = 'both', perm.max = 999) # 変数増減法でモデル選択 plot(rda.select, type = "text") # RDA の結果をトリプロットで図示 

 RDA でも CCA と同様、説明変数として v3 と v1 のみが選択されている。また、種のスコアとサンプルのスコアのスケールが揃っていないため、種のスコアについては見づらくなっている。

GLR-ML

 まず、パッケージ"rioja"を読み込む。

library(rioja) # パッケージ "rioja" の読み込み

 RDA の第一軸との相関が最も大きかった v3 に対して、GLR を実行する。

fit <- MLRC(spe2, env$v3) # 各種に対して GLR を一括して実行 summary(fit) # 結果の概要を表示 cf <- coef(fit) # GLR の係数を "cf" に格納 write.csv(cf, "coefficient.csv", quote=FALSE, row.names=TRUE) # CSVファイルへの出力 

 CSV ファイルを Microsoft Excel であけると、次のように出力されている。

 最適点 optimum = -b1/b2/2, トレランス tolerance = √-1/b2/2,A = exp (b0-b1*b1/b2/4) と計算することにより、ガウス関数のパラメータが算出される。

 二乗項の係数 b2 が正になっている種では、トレランスは算出されず、最適点の代わりに珪藻の相対度数が極小になる点が算出されるため、注意が必要である。

 10サンプルのテストデータに対して、v3 上での推定点を予測する。

> predict(fit, spe_new) # 推定点の予測

 R-console 上に、以下の通り計算結果が出力される。

$fit
        MLRC
X1  4.612957
X2  5.908878
X3  4.928422
X4  5.701026
X5  2.570330
X6  5.232751
X7  5.031613
X8  5.801649
X9  6.871378
X10 4.237271