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

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

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

２）有効種数計算の準備

[2-i] ループ構造

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

「繰り返し・ループ（for loop）」

[2-ii] ライブラリのインストールと読み込み

有効種数の計算では、各種の相対個体数(p_i)を使うため、群集生態学用のライブラリ（vegan）で提供されている便利な関数（decostand）を使って、個体数データを相対個体数データへ一括変換します。そのための準備としてライブラリveganのインストールします。以下のスクリプトはveganが既にインストールされていないかをチェックしてインストールされていないときだけインストールをする命令になっています。これでうまくいかない場合はGUI操作も試してください（ここ）。それでもうまく書ない場合は、ここを飛ばして[2-iii]に進んでください。

######Calculate Effective Species Number###########################

#install library only when it has not been installed

pkg.name <- "vegan"

if(!require(pkg.name, character.only=TRUE)){

  install.packages(pkg.name)

}

ちなみに以下の図のようなメッセージが出た場合は、すでにインストールされかつ読み込まれているライブラリーを再度インストールしようとしていることを意味しますので、キャンセルを押してください。

上の部分がうまくいかない場合は、ライブラリーのインストールと読み込みを諦めて、[2-iii]のvegan()を使わない場合に進んでください。

ライブラリーの読み込みはlibrary(ライブラリ名)という命令文で書きます。

#load library

library(vegan)

[2-iii] 相対個体数計算

veganライブラリーとdecostand()関数を使わない場合

ライブラリーのインストールや読み込みがうまくいかない場合は以下のコードで同じ結果になるはずですので、そのまま進んでください。

species_biwa_data2 <- species_biwa_data/rowSums(species_biwa_data)

veganライブラリーのインストールがうまくいった場合

個体数データを相対個体数データへ一括変換します。decostand(元データ、手法、行／列どちらを基準に変換するか)という形でパラメータを指定します。method="total"は相対個体数への変換で、MARGIN=1はデータフレームの行（横）方向を基準にした変換です（各サンプルデータが横方向に並んでいることに注意）。

#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,])

[3-i] 0次の有効種数の計算のトライアル

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

#trial to calculate 0th-order Effective Species Number

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

[3-ii] 1次の有効種数の計算のトライアル

次に１次の有効種数の計算には相対個体数(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となって計算がうまくいきません。したがって、View(temp)で中身を確かめると"NaN"という値が入っている要素があることに気づくでしょう。

#trial to calculate 1st order Effective Species Number

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

View(temp)

そこでNaNになる要素はゼロに変換するという処理を施します。そのためには以下のコードを実行します。

#convert NAN to 0

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

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

exp(-sum(temp))

[3-iii ] 2次の有効種数の計算のトライアル

[自分で考えてコードを書く]

同様に, [3-ii]の内容を参考に、２次の有効種数も自分で計算してみましょう。まずは上の有効種数の公式を見直してください。

それでも、Rスクリプトをどう書いてよいかわからない人は、エクセルファイル（effective_species_number.xlsx）のD～F列のセルに書き込まれている関数式を頼りにどのような計算をすればよいか（→相対個体数の二乗を各種ごとに計算し、その総和を取って最後にその逆数を計算する）を確認してください。以下の（コメントを抜いて）三行のスクリプトの？部分を適切なコードで置き換えれば計算できます。

ただし、スクリプト中の日本語コメントは文字化けするのでコピペしないこと！！

#trial to calculate 2nd order Effective Species Number

temp <- ?*?   #相対個体数を用いた計算

temp2 <- sum(?) #何かの総和を取る

1/temp2  #逆数を取る（temp2を分母に持ってくる）

４）有効種数の自動計算

[4-i ] 0次、1次の有効種数の自動計算

この計算を各行（各サンプル）について繰り返せば、すべてのサンプルに関する０、１、２次の有効種数が計算できますが、これをいちいちコードを書いていると面倒くさいです。そこで以下のようにfor loopを用いて自動化します。最初の三行は有効種数の計算結果を格納するためのベクトルの用意です。

for loopの中のlength(species_biwa_data2[,1])はデータフレームspecies_biwa_data2の行数を確認するためのもので、行数分（サンプル数分）だけループが繰り返されます。

このfor loopのポイントは、{}の中でkの値が1からlength(species_biwa_data2[,1])=17まで一つずつ大きくなっていき、ループが回るごとに異なるサンプルに対して有効種数を計算し、その計算結果を用意したベクトル（たとえばesn0）のk番目の要素に格納していく点である。

ただし、スクリプト中の日本語コメントは文字化けするのでコピペしないこと！！

#Auto calculation for effective species number (esn)

esn0 <- vector()  #0次の有効種数格納のためのベクトル

esn1 <- vector()  #1次の有効種数格納のためのベクトル

esn2 <- vector()  #2次の有効種数格納のためのベクトル

length(species_biwa_data2[,1]) #サンプル数のチェック

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

  #0th-order Effective Species Number

  esn0[k] <- sum(species_biwa_data2[k,] > 0) #ｋ番目のサンプルの0次の有効種数

  #1st-order Effective Species Number

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

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

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

  #2nd order Effective Species Number

}

View(esn0)#計算結果の確認

View(esn1)#計算結果の確認

[4-ii] 2次の有効種数の自動計算

[自分で考えてコードを書く]

うえの[4-i]の太文字の部分をよく理解したら、[3-iii]の部分で1行目のサンプルから2次の有効種数をするために自分で作ったスクリプトを微調整して（k行目のサンプルに対して計算できるようにして）上のfor loopの中に入れれば、2次の有効種数ベクトルesn2の計算が完了するはずです。

最後にView(ens2)で2次の有効種数ベクトルesn2の中身を確かめたら終了です。