群集データの加工はけっこう骨の折れる作業である。
エクセル上で力技でやってしまう人もしばしば見かける。
しかし、これは絶対R上でやった方がよい。
Rでは、データの形式を加工するパッケージや関数がとても充実しており、中でもパッケージ'dplyr'や'tidyr'の使い勝手の良さは群を抜いている。
ここでは、dplyrとtidyrを用いた群集データの加工法を紹介する。
dplyrとtidyrは強力
横長データと縦長データを理解すべし
データは各行を識別（identify）できるキーの行（orキー行の組み合わせ）をつけておくべし
group_byとsummariseをマスターせよ
横長化（spread）によって未観察データを自動で埋めることができる
種Aと種Bの個体数を、サイト1~3で調査し、以下のような形式でデータを作成したとしよう。
> Community.Data
Site Species_A Species_B
1 Site1 0 10
2 Site2 1 2
3 Site3 1 3
サイト1では種Aが0個体、種Bが10個体。サイト2では種Aが1個体、種Bが2個体。サイト3では種Aが1個体、種Bが3個体観察されたことがわかる。
さて、ここで問題にしたいのはデータの内容ではなく、その形式である。
1行目にサイト1のデータ、2行目にサイト2のデータ…というようにデータが並んでいる。
一方、各種の個体数はそれぞれ別の列に並んでいる。
このようなデータ形式を、しばしばcommunity matrixと呼ぶ（2種だけなのでピンと来ないかもしれないが）。
RでPCAやNMDSなどの多変量解析を行う場合は、データをこの形式に整理する必要がある。
次に、このデータの形式を変換することを考える（何のために？というのはひとまず置いておく）。
具体的には、種の列と個体数の列を作る。
この変換のために、ここではライブラリ'tidyr'内の関数'gather'を用いる。
実際のスクリプトだが、まずは丸覚えで良い。
> library(tidyr)
> gathered.Data <- gather(Community.Data, key=Species, value=Abundance, -1)
> gathered.Data
Site Species Abundance
1 Site1 Species_A 0
2 Site2 Species_A 1
3 Site3 Species_A 1
4 Site1 Species_B 10
5 Site2 Species_B 2
6 Site3 Species_B 3
>
このようなデータ形式は、tidy dataと呼ばれる。
tidy dataの詳しい解説はここのサイトに譲るとして、この形式は、線形モデル等の統計解析を行いやすいという利点がある。
慣れない人にとっては、この形式の違いにどういう意味があるのかわかりづらいかもしれないが、まずはcommunity matrixと見比べて違いをよくよく考えてみてほしい。
最後に、このデータをcommunity matrixに戻すことを考える。
そのためには、tidyr::spreadを用いる。
> spread.Data <- spread(gathered.Data, key=Species, value=Abundance)
> spread.Data
Site Species_A Species_B
1 Site1 0 10
2 Site2 1 2
3 Site3 1 3
>
このように、spreadを使うことで、最初のデータ形式に戻すことができる。
この2つの形式はそれぞれに使いどころがあるので、両者を自由に行き来できるようになるとよい。
ライブラリdplyrとtidyrには様々な機能があるが、大きな特徴の一つに、パイプ演算子%>%がある。
パイプ演算子とはどういうものなのか、まずは使ってみる。
> Community.Data
Site Species_A Species_B
1 Site1 0 10
2 Site2 1 2
3 Site3 1 3
>
> library(tidyr)
>
> gathered.Data <-
+ Community.Data %>%
+ gather(key=Species, value=Abundance, -1)
> gathered.Data
Site Species Abundance
1 Site1 Species_A 0
2 Site2 Species_A 1
3 Site3 Species_A 1
4 Site1 Species_B 10
5 Site2 Species_B 2
6 Site3 Species_B 3
>
上では、パイプ演算子をCommunity.Dataをgathered.Dataに変換する際に使っている。
ここで、パイプ演算子はCommunity.Dataを関数gatherに渡す機能を果たしている。
パイプ演算子の便利なところは、データに対する加工を次々に行えることである。
以下の例を見てみよう。
> spread.Data <-
+ Community.Data %>%
+ gather(key=Species, value=Abundance, -1) %>%
+ spread(key=Species, value=Abundance)
> spread.Data
Site Species_A Species_B
1 Site1 0 10
2 Site2 1 2
3 Site3 1 3
>
上では、Community.Dataを関数gatherに渡し、gatherの加工を受けたものをさらに関数spreadに渡している。
パイプ演算子を使うことでスクリプトがキレイになり、後からスクリプトを見た時に、どのような加工をしたかわかりやすくなるという利点もある。
さて、以上の序論を踏まえて、群集研究におけるデータ整理の現場でしばしば直面する課題をdplyrとtidyrで解決してみよう。
種A, B, Cの個体数を、サイト1~3でそれぞれ3回調査したところ、以下のようなデータが得られた。
> Community.Data
Site Species_A Species_B Species_C
1 Site1 0 10 1
2 Site1 1 5 0
3 Site1 0 15 0
4 Site2 1 2 5
5 Site2 2 1 0
6 Site2 1 5 0
7 Site3 1 3 0
8 Site3 1 1 2
9 Site3 0 1 4
10 Site4 0 1 2
11 Site4 0 0 0
12 Site4 1 1 1
これをgatherによって縦長にしてから、spreadによって元のデータ形式に戻すことが一つ目の課題である。
ただしこの問題は少しずるくて、そもそもこのデータ形式に一つ不備がある。それを以下で明らかにしていく。
まず、序論で説明した方法で、このデータをgatherで縦長にする。
すると、以下のような縦長データが得られる。
> gathered.Data <- Community.Data %>% gather(key=Species, value=Abundance, -1)
> gathered.Data
Site Species Abundance
1 Site1 Species_A 0
2 Site1 Species_A 1
3 Site1 Species_A 0
4 Site2 Species_A 1
5 Site2 Species_A 2
6 Site2 Species_A 1
7 Site3 Species_A 1
8 Site3 Species_A 1
9 Site3 Species_A 0
10 Site4 Species_A 0
11 Site4 Species_A 0
12 Site4 Species_A 1
13 Site1 Species_B 10
14 Site1 Species_B 5
15 Site1 Species_B 15
16 Site2 Species_B 2
17 Site2 Species_B 1
18 Site2 Species_B 5
19 Site3 Species_B 3
20 Site3 Species_B 1
21 Site3 Species_B 1
22 Site4 Species_B 1
23 Site4 Species_B 0
24 Site4 Species_B 1
25 Site1 Species_C 1
26 Site1 Species_C 0
27 Site1 Species_C 0
28 Site2 Species_C 5
29 Site2 Species_C 0
30 Site2 Species_C 0
31 Site3 Species_C 0
32 Site3 Species_C 2
33 Site3 Species_C 4
34 Site4 Species_C 2
35 Site4 Species_C 0
36 Site4 Species_C 1
ここまでは上手いこと縦長になっている。
次に、これを序論のようにspreadで横長に直すことを試みる。
> gathered.Data %>% spread(key=Species, value=Abundance)
Error: Each row of output must be identified by a unique combination of keys.
Keys are shared for 36 rows:
* 1, 2, 3
* 4, 5, 6
* 7, 8, 9
* 10, 11, 12
* 13, 14, 15
* 16, 17, 18
* 19, 20, 21
* 22, 23, 24
* 25, 26, 27
* 28, 29, 30
* 31, 32, 33
* 34, 35, 36
Do you need to create unique ID with tibble::rowid_to_column()?
Call `rlang::last_error()` to see a backtrace
すると、以上のようなエラーが出てしまう。
エラーメッセージを読むと「アウトプットデータの行は、唯一のキー（の組み合わせ）によって識別できる必要あり」とある。
（rowid_to_columという関数で問題解決できるとも書いてあるが、とりあえず無視する）ここで言うアウトプットデータとは、最初のCommunity.Dataのような形式のデータのことのはず。
Community.DataにおけるキーはSite列なのだが、例えば1~3行目はどれもSite1になっており、確かに各行を唯一のキーによって識別できていない。
これを回避するためには、以下のようなデータ形式を作る必要がある。
> Community.Data
Week Site Species_A Species_B Species_C
1 1 Site1 0 10 1
2 2 Site1 1 5 0
3 3 Site1 0 15 0
4 1 Site2 1 2 5
5 2 Site2 2 1 0
6 3 Site2 1 5 0
7 1 Site3 1 3 0
8 2 Site3 1 1 2
9 3 Site3 0 1 4
10 1 Site4 0 1 2
11 2 Site4 0 0 0
12 3 Site4 1 1 1
実は、最初のデータは、種A, B, Cの個体数を、サイト1~3で”それぞれ週に1回”、計3回調査して得られたデータだった。
このようにすれば、唯一のWeek列とSite列との組み合わせにより、各行を識別することができる。
これを、もう一度gather、spreadしてみる。
> gathered.Data <- Community.Data %>% gather(key=Species, value=Abundance, -c(1:2)) #ここの最後の引数は-1ではなく、-c(1:2)にする必要がある。これは、gatherしたい行以外の行が1行目と2行目のふたつになっているから。
> gathered.Data
Week Site Species Abundance
1 1 Site1 Species_A 0
2 2 Site1 Species_A 1
3 3 Site1 Species_A 0
4 1 Site2 Species_A 1
5 2 Site2 Species_A 2
6 3 Site2 Species_A 1
7 1 Site3 Species_A 1
8 2 Site3 Species_A 1
9 3 Site3 Species_A 0
10 1 Site4 Species_A 0
11 2 Site4 Species_A 0
12 3 Site4 Species_A 1
13 1 Site1 Species_B 10
14 2 Site1 Species_B 5
15 3 Site1 Species_B 15
16 1 Site2 Species_B 2
17 2 Site2 Species_B 1
18 3 Site2 Species_B 5
19 1 Site3 Species_B 3
20 2 Site3 Species_B 1
21 3 Site3 Species_B 1
22 1 Site4 Species_B 1
23 2 Site4 Species_B 0
24 3 Site4 Species_B 1
25 1 Site1 Species_C 1
26 2 Site1 Species_C 0
27 3 Site1 Species_C 0
28 1 Site2 Species_C 5
29 2 Site2 Species_C 0
30 3 Site2 Species_C 0
31 1 Site3 Species_C 0
32 2 Site3 Species_C 2
33 3 Site3 Species_C 4
34 1 Site4 Species_C 2
35 2 Site4 Species_C 0
36 3 Site4 Species_C 1
> gathered.Data %>% spread(key=Species, value=Abundance)
Week Site Species_A Species_B Species_C
1 1 Site1 0 10 1
2 1 Site2 1 2 5
3 1 Site3 1 3 0
4 1 Site4 0 1 2
5 2 Site1 1 5 0
6 2 Site2 2 1 0
7 2 Site3 1 1 2
8 2 Site4 0 0 0
9 3 Site1 0 15 0
10 3 Site2 1 5 0
11 3 Site3 0 1 4
12 3 Site4 1 1 1
このように、キーが識別可能であれば、縦長データと横長データを行き来することができる。
ついでに、このデータを使って、層別集計をdplyrで行う方法を紹介する。
たとえば、この調査について、各サイトでの各種の累積個体数を求めたい場合があるだろう。
このような場合、dplyrではgroup_by関数とsummarise関数を用いる。
スクリプトは以下のようになる。
Community.Data.Cuml <-
Community.Data %>%
gather(key=Species, value=Abundance, -c(1:2)) %>%
group_by(Site, Species) %>%
summarise(CumulAbundance=sum(Abundance))
まず、gatherで縦長データに変換し、これをgroup_by関数に渡す。
ここでは、各サイト・各種で集計したいため、group_byではSiteとSpeciesを指定する。
ここでSpeciesも入れたいために、gatherで縦長データに変換しているわけである。
これをsummarise関数に渡し、この中では「集計後の列名=集計する関数(集計するデータ)」というような形で指定する。
すると、以下のようなデータを作ることができる。
> Community.Data.Cuml
# A tibble: 12 x 3
# Groups: Site [4]
Site Species CumulAbundance
<fct> <chr> <dbl>
1 Site1 Species_A 1
2 Site1 Species_B 30
3 Site1 Species_C 1
4 Site2 Species_A 4
5 Site2 Species_B 8
6 Site2 Species_C 5
7 Site3 Species_A 2
8 Site3 Species_B 5
9 Site3 Species_C 6
10 Site4 Species_A 1
11 Site4 Species_B 2
12 Site4 Species_C 3
なお、このような操作をするとデータシートがtibbleという形式に変換される。
これはdplyrのデータ形式で、ほとんどdata.frameと同じように使えるのだが、一部動作が異なる場合があるので注意。data.frameに戻したい場合は、as.data.frame関数で戻せる。
この例ではgather関数で縦長にしたままにしているのでアウトプットも縦長データになったが、横長で用意したい場合は以下のようにすればよい。
> Community.Data.Cuml <-
+ Community.Data %>%
+ gather(key=Species, value=Abundance, -c(1:2)) %>%
+ group_by(Site, Species) %>%
+ summarise(CumulAbundance=sum(Abundance)) %>%
+ spread(key=Species, value=CumulAbundance)
> Community.Data.Cuml
# A tibble: 4 x 4
# Groups: Site [4]
Site Species_A Species_B Species_C
<fct> <dbl> <dbl> <dbl>
1 Site1 1 30 1
2 Site2 4 8 5
3 Site3 2 5 6
4 Site4 1 2 3
累積個体数のcommunity matrixを作ることができた。
このようなデータ操作は繰り返し測定にありがちであり、覚えておくと便利だと思う。
群集調査では、調査を通してどのような種が観察されるか事前にはわからないのが普通である。
そのため、調査が進むうちに新たな種が追加されることがよくある。
その結果、以下のようなデータがしばしば出来上がる。
> gathered.Data2
Site Species Abundance
1 Site1 Species_B 10
2 Site1 Species_C 1
3 Site2 Species_A 1
4 Site2 Species_B 2
5 Site2 Species_C 5
6 Site3 Species_A 1
7 Site3 Species_B 3
8 Site4 Species_B 1
9 Site4 Species_C 2
この調査では、サイト2を調査して初めて種Aが観察され、一方サイト3では種Cは観察されず、サイト4では種Aが観察されなかった。
さて、観察されなかった種は個体数0のはずだが、それを明示するようなデータを自動的に作ることはできないだろうか？
実は、spread関数はこのような場合にも使うことができる。
spread関数にはfillという引数が用意されており、これを用いて「観察されなかった」ことを明示するデータを作ることができる。
以下を見てみよう。
> gathered.Data2 %>%
+ spread(key=Species, value=Abundance, fill=0)
Site Species_A Species_B Species_C
1 Site1 0 10 1
2 Site2 1 2 5
3 Site3 1 3 0
4 Site4 0 1 2
以上を見ればわかるように、そもそも横長データを作るためには、サイト数4×種数3=12個のデータが必要である。
そのため、gathered.Data2のような「不完全な縦長データ」を横長にするためには、どうしてもそのマスを埋める必要がある。
これを逆手に取れば、未観察データを自動的に追加することができるわけである。
fill引数は、この空白マスを埋める数字を指定する引数である。
縦長データ形式でこれを用意したい場合は、以下のようにすればよい。
> gathered.Data2.fill <-
+ gathered.Data2 %>%
+ spread(key=Species, value=Abundance, fill=0) %>%
+ gather(key=Species, value=Abundance, -1)
> gathered.Data2.fill
Site Species Abundance
1 Site1 Species_A 0
2 Site2 Species_A 1
3 Site3 Species_A 1
4 Site4 Species_A 0
5 Site1 Species_B 10
6 Site2 Species_B 2
7 Site3 Species_B 3
8 Site4 Species_B 1
9 Site1 Species_C 1
10 Site2 Species_C 5
11 Site3 Species_C 0
12 Site4 Species_C 2
見事に未観察データが0で埋まっている。
これはちょっとした小技ではあるが、群集解析では非常によく使う機会の多い技でもある。
このくらいのデータ数だったら手作業で一つずつ埋める方が早いかもしれないが、サイト数・種数ともに多くなってくると、手作業で修正するのはほぼ不可能だろう。
この方法を使えば、どんなにデータが多くても一瞬でゼロを埋めてくれる。