Rを使った安定同位体比データ解析（Part 4）

2024年8月17日 14:07

はじめに

安定同位体比は、生体内の代謝や、食物網における食う–食われるの関係、生息地の移動など、わたしたちの身の回りで起きている様々な「物の流れ」を調べるツールとして用いられてきました。

近年、分析共同利用拠点や依頼分析サービスが増え、ますます多くの研究者に安定同位体比分析が使われていることと思います。
また、消費者の餌資源推定や栄養ニッチの算出など、安定同位体比データを用いた解析の幅が広がるとともに、
同位体研究において、Rを活用する機運が高まっていると感じています。

この記事は、Rの使い方を学びながら、
安定同位体比データをどうやって図示するか、どのように統計解析をおこなうか、Rのパッケージをどう動かすかを伝えることが目的です。

このページはRを使った安定同位体比データの解析シリーズの Part 4 です

Part 1 >>> データファイルの作り方 / 散布図の作り方
Part 2 >>> 箱ひげ図の作り方
Part 3 >>> 集計表の作り方

Part 4 で学べること

安定同位体比データを群間比較する方法（ノンパラメトリック法）
統計解析結果をテキストファイルとして書き出す方法

統計解析

Part 1–3では、安定同位体比データを図示しました。
図を眺めていると、グループ間で安定同位体比に有意差があるのだろうかと疑問に思うかもしれません。

この記事では、2群間でデータを比較する方法と、3群以上でデータを比較する方法を紹介します。

データセットの準備

Part 1 で用いたデータセットを使います。
ただし、全グループを比較すると計算量が増えて大変なので、一部（一次生産者=POMとPeriphyton）を抜粋して説明していきます。

dt <- read.csv("Isotope.csv", header = T)

colnames(dt)[4] <- c("delta15N")
colnames(dt)[5] <- c("delta13C")

dt$delta15N <- round(dt$delta15N, digits = 1)
dt$delta13C <- round(dt$delta13C, digits = 1)

それでは、次のシナリオを考えます。

①湖の一次生産者 [付着藻類と植物プランクトン] の間で安定同位体比に違いがあるか？
②採取年 [2017年, 2018年, 2019年] によって湖の一次生産者に違いがあるか？

グループ間における安定同位体比の違いを統計的に議論するため、①ではノンパラメトリックの2群間比較をおこない、②ではノンパラメトリックの多重比較をおこないます。

2群間比較

①湖の一次生産者 [付着藻類と植物プランクトン] の間で安定同位体比に違いがあるか？を検証します。
まずは、一次生産者の炭素安定同位体比を比較してみましょう。
今回は、ウィルコクソンの順位和検定 (Wilcoxon rank sum test) を採用します。
Rでは、"exactRankTests"パッケージでウィルコクソンの順位和検定を利用できます。

library(exactRankTests)
library(dplyr)

data.POM <- dt %>% 
  dplyr::filter(Type == "POM")

data.Peri <- dt %>% 
  dplyr::filter(Type == "Periphyton")

wilcox.exact(x=data.POM$delta13C,y=data.Peri$delta13C,paired=F)

解説

#1 library(exactRankTests) でライブラリー"exactRankTests"を呼び出します。
はじめて使う場合は、install.packages("exactRankTests")でライブラリーをダウンロードしましょう。

#2 library(dplyr) でライブラリー"dplyr"を呼び出します。
データ操作（行・列の抽出、削除など）で非常によく使われるパッケージです。

#3,4 dplyr::filter() で、一次生産者のデータセットからPOMのデータのみを取得します。

#5,6 dplyr::filter() で、一次生産者のデータセットからPeriphytonのデータのみを取得します。

#7 wilcox.exact() で、ウィルコクソンの順位和検定をおこないます。
関数のルールは、wilcox.exact(x=ベクトル1,y=ベクトル2,paired=F) です。
今回は、POMのδ13CとPeriphytonのδ13Cを比較したいので、wilcox.exact(x=data.POM$delta13C,y=data.Peri$delta13C) としました。
data.POM$delta13C で、データフレーム"data.POM"の列"delta13C"を取得します。
つまり、POMのδ15NとPeriphytonのδ15Nを比較したい場合は、wilcox.exact(x=data.POM$delta15N,y=data.Peri$delta15N) となります。

また、ウィルコクソンの順位和検定をおこなう場合は必ず paired=F を指定してください。
paired=T を指定すると、ウィルコクソンの符号順位検定（ウィルコクソンの順位和検定とは別物！）が実行されます。

さて、気になるp値は p-value < 2.2e-16 でした。
有意水準をα = 0.05と仮定した場合、δ13CはPOMとPeriphytonの間で有意差があると言えます。
※Wは統計量です。

【参考】

多重比較

次に、②採取年 [2017年, 2018年, 2019年] によって湖の一次生産者に違いがあるか？を検証します。
POMのδ13Cを用いてテストしてみましょう。
ここでは、スティール・ドゥワス検定 (Steel-Dwass test) を採用します。
Rでは、"NSM3"パッケージでスティール・ドゥワス検定を利用できます。

library(NSM3)

data.POM.year <- split(data.POM,data.POM$Year)

list.POM.d13C <- list(
  data.POM.year$'2017'$delta13C,
  data.POM.year$'2018'$delta13C,
  data.POM.year$'2019'$delta13C
)

pSDCFlig(list.POM.d13C, method="Asymptotic")

解説

#8 library(NSM3) でライブラリーを呼び出します。
はじめて使う場合は、install.packages("NSM3")でライブラリーをダウンロードしましょう。

#9 split() で、データセットを分割します。
関数のルールは、split(データセット, 分割基準となる変数) です。
今回は、POMのデータセットを年ごとに分割したいので、split(data.POM,data.POM$Year) としました。
そして、年ごとに分割したデータセットをオブジェクト"data.POM.year"に保管します。

#10-14 list() で、データリストを作成します。
list() は、Rでよく使われる関数です。例えば、以下のように使います。

list.A <- list(c(1,2,3),
               c(4,5,6,7,8),
               c(9,10,11,12,13,14))
list.A 

> list.A 
[[1]]
[1] 1 2 3


[[2]]
[1] 4 5 6 7 8


[[3]]

[1] 9 10 11 12 13 14

このように、複数のデータ列（ベクトル）を独立して格納しておくことができます。

今回は、最終的にスティール・ドゥワス検定をおこなうことが目標です。
スティール・ドゥワス検定の関数では、比較するデータ群がリスト化されていると便利なのです。
（後ほど解説します）

data.POM.year$'2017'$delta13C で、データセット"data.POM.year"の要素"'2017'"の列"delta13C" を取得します。

#15 pSDCFlig(list.POM.d13C, method="Asymptotic") で、スティール・ドゥワス検定をおこないます。
関数のルールは、pSDCFlig(ベクトル1,ベクトル2,ベクトル3,method="検定の種類") です。

以下に例を示します。

a <- c(1,2,3)
b <- c(4,5,6,7,8)
c <- c(9,10,11,12,13,14)
  
pSDCFlig(a,b,c, method="Asymptotic")

または、

pSDCFlig(c(1,2,3),
         c(4,5,6,7,8),
         c(9,10,11,12,13,14), 
         method="Asymptotic")

関数のなかでデータやオブジェクト（データが格納されている箱）の名前を長々と書き連ねると、コードが非常に複雑で見づらくなります。
比較するデータ群をあらかじめリスト化しておき、関数内ではそのリストを指定することで、理解しやすいコードになると思います。
また、list()はfor文でも活用しやすいので、覚えておくと便利です。

method には、Exact, Monte Carlo, Asymptotic の3種類が用意されています。今回は"Asymptotic" を使いました。
違いについては、こちらのサイトで詳しく説明されています。

全てのグループの組み合わせに対して検定結果が出力されます。
Critchlow-Fligner W Statistic が検定量、
The smallest experimentwise error rate leading to rejection is … がp値です。
有意水準をα = 0.05と仮定した場合、1–2と2–3では有意差ありですが、1-3では有意差なし、となります。

1,2,3がどのデータ群に対応しているかですが、基本的にlist()で指定した順番に並んでいます。Group sizes (各データ群のサンプルサイズ) がヒントになります。

【参考】

テキストファイルに出力する

ウィルコクソンの順位和検定やスティール・ドゥワス検定の結果をテキストファイルに保存してみましょう。

Wilcox.d13C <- wilcox.exact(x=data.POM$delta13C,y=data.Peri$delta13C,paired=F)
sink("Wilcox.d13C.txt", append = TRUE)
print (Wilcox.d13C)   
sink() 

Steel.d13C <- pSDCFlig(list.POM.d13C, method="Asymptotic")
sink("Steel.d13C.txt", append = TRUE)
print (Steel.d13C)   
sink()

解説

#16 wilcox.exact()でウィルコクソンの順位和検定をおこない、結果を"Wilcox.d13C"というオブジェクトに格納します。
（#7のコードの転用です）

#17 sink() で保存時のファイル名を指定します。
#18 print() で保存したいオブジェクトを指定します。
#19 sink() で保存します。現在のディレクトリに.txtファイルが出現しているはずです。

#20 pSDCFlig() でスティール・ドゥワス検定をおこない、結果をオブジェクト"Steel.d13C"に格納します。
（#15のコードの転用です）

#21-23 は#17-19の解説と同じです。

注意点

今回、ノンパラメトリック版の比較検定法として、ウィルコクソンの順位和検定とスティール・ドゥワス検定を紹介しました。
スティール・ドゥワス検定はウィルコクソンの順位和検定の計算方法がベースとなっています。（どちらも検定量はWなのです）
ちなみに、スティール・ドゥワス検定 pSDCFlig() で method="Exact"を指定し、2群間の比較をおこなうと、wilcox.exact()での計算と同じ結果になるはずです。
ここで注意したいのが、3群以上の比較を行う場合は、必ず多重比較、つまりスティール・ドゥワス検定を使うということです。
3群以上の差を検定する際に、2群間の検定を行う（繰り返す）と、帰無仮説が誤って棄却される確率、偽陽性が高くなるためです。

今回の内容はここまで。

一般化線形モデル、一般化線形混合モデルを用いた安定同位体比データの統計解析について、いずれ加筆するかもしれません。
（独り言：Wilcoxon rank sum test による2群間の比較結果をBonferroni法で調整するのはありなんだろうか・・・）

***************************************************
ご不明な点があれば、お気軽にご連絡ください。
***************************************************

>>> Part 5 安定同位体比分析のためのRパッケージ