１）古典的種多様性指標と、「現代的多様性指標＝有効種数」の定義・意味・仮想例

以下のファイルにまとめてあるので復習すること（実習のプリントの一部）。

effective_species_number.pdf

より詳しい内容は以下の英文総説を参照のこと。

Anne Chao et al. (2014) Unifying Species Diversity, Phylogenetic Diversity, Functional Diversity, and Related Similarity and Differentiation Measures Through Hill Numbers. Annual Review of Ecology, Evolution, and Systematics Volume 45, pp 297-324

２）有効種数計算過程のRスクリプトへの実装

有効種数の計算においては、各種ごとの計算値の総和（Σ）を求める必要があるため、R上ではループ構造によって繰り返し計算をするか、もしくはsum()関数を使う必要があります。いずれにしてもループ構造の計算が含まれることになるので、まずはプログラミング構造について、特にfor loopについて復習する必要があります（ここ）。

有効種数計算のための前処理

次に、有効種数計算過程を学ぶためのサンプルオブジェクトを現在使っているRスクリプトと同じフォルダ（つまり、Qドライブ直下）にダウンロードし保存します。

species_biwa_data.rds

以下のようにRスクリプトに読み込みます。

#load the result

species_biwa_data <- readRDS("species_biwa_data.rds")

View関数を使って中身を見れば（あるいはRstudioのGlobal Environmentのウィンドウからオブジェクトをクリックすれば）、２０１８年の９月と１０月のプランクトン個体数リストであることがわかります。

相対個体数の計算

有効種数の計算では、各種の相対個体数(p_i)を使うため、群集生態学用のライブラリ（vegan）で提供されている便利な関数（decostand）を使って、個体数データを相対個体数データへ一括変換します。decostand(元データ、手法、行／列どちらを基準に変換するか)という形でパラメータを指定します。method="total"は相対個体数への変換で、MARGIN=1はデータフレームの行（横）方向を基準にした変換です（各サンプルデータが横方向に並んでいることに注意）。

#load library

library(vegan)

#convert to the proportion

species_biwa_data2 <- decostand(species_biwa_data, method="total", MARGIN=1)

相対個体数への変換の方向が正しく指定されているか不安だったら、以下のコードで確認できます。species_biwa_data2の1行目の要素の総和（sum）を計算するという意味です。上の計算ですでに相対個体数に変換されているはずなので総和は１になるはずです。

sum(species_biwa_data2[1,])

それでは次に有効種数を計算してみましょう。有効種数の公式は以下の通りです。

有効種数の計算のトライアル（試行）

まず0次の有効種数、すなわちただの種数Sを計算してみましょう。一行目のサンプルcells0924no1についてやってみると以下のコードでできます。これは総和を計算する関数sumと条件式species_biwa_data2[1,]>0の組み合わせです。この条件式が満たされる（TRUE）要素数の合計を計算しています。

sum(species_biwa_data2[1,]>0) #species richness

次に１次の有効種数の計算には相対個体数(p_i)に関する自然対数ln()を計算する必要があります。Rでは自然対数はlog()という関数で計算できます。したがって、species_biwa_data2[1,]*log(species_biwa_data2[1,])のように計算すればOKです。R上ではベクトル同士や行列同士の計算は内積になるのではなく要素ごとの掛け算になるので便利です。ここで注意したいのは、数学的にはx*ln(x)はx→0の極限でゼロに収束するのに対し、数値計算ではうまくそのような計算はできないため、p_i = 0 (種iが存在しない)時, p_i*log(p_i)は０*-Inf = NaNとなって計算がうまくいきません。そこでNaNになる要素はゼロに変換するという処理を施します。そのためには以下のコードを実行します。

temp <- log(species_biwa_data2[1,])*species_biwa_data2[1,]

temp[is.nan(as.matrix(temp))]<-0

最後にこのtempにマイナスを掛けて指数関数の中に入れれば、１次の有効種数exp(-sum(p_i*ln(p_i)))が計算できます。

exp(-sum(temp))

同様にして、２次の有効種数も自分で計算してみましょう。

有効種数の自動計算

この計算を各行（各サンプル）について繰り返せば、すべてのサンプルに関する０、１、２次の有効種数が計算できますが、これをいちいちコードを書いていると面倒くさいです。そこで以下のようにfor loopを用いて自動化します。最初の三行は有効種数の計算結果を格納するためのベクトルの用意です。for loopの中のlength(species_biwa_data2[,1])はデータフレームspecies_biwa_data2の行数を確認するためのもので、行数分（サンプル数分）だけループが繰り返されます。

esn0 <- vector()

esn1 <- vector()

esn2 <- vector()

for(k in 1:length(species_biwa_data2[,1])) {

  temp <- log(species_biwa_data2[k,])*species_biwa_data2[k,]

  temp[is.nan(as.matrix(temp))] <- 0

  esn1[k] <- exp(-sum(temp))

  esn0[k] <- sum(species_biwa_data2[k,] > 0) 

}

ここまでで有効種数の計算は完了です。自分でやるときは、サンプルデータには2019年5月の各班のものをすべて結合し、2次の有効種数まで計算して以上の処理を完成させましょう。