2020年12月20日日曜日

tidyverse版: 複数生物種の個体数集計データを縦に伸ばすor縮める(集計データ← →データフレーム変換)

 (以前の記事をtidyverseで作り直してみました。まだ旧来のRコードの方が馴染みがあるので、tidyverseに慣れている方からすると不格好かもしれません。)

多数の種が含まれる生物群集データを解析する時、エクセルなどにデータを整理していると、各種をインデックスにするか要因のひとつとして整理するか、ケースによって変わってくるかと思います。とくに各種を応答変数とするか説明変数に使用するかあたりで必要になってくるかと。完全なデータフレーム型にしておいた方が統計解析はやりやすいけれど、膨大な行数となってエクセルの限界にあっさりと到達したりもするし…(cf. エクセルのvlookup関数)。


こんな時に、一部集計データとデータフレーム型とを自在に一発変換できたら効率がよいと思い、やり方をまとめてみました。正直、入力の手間が一番省けるのは部分的には集計されたデータだったりしますよね(下記のdata1のような形)。

データの整形にはdplyrパッケージの関数などを使用します。tidyverseパッケージとして関連パッケージも一括でインストールが可能です。

require(tidyverse) # 要インストール、tidyverseのパッケージ群が一括で呼び出される。

# 例えば、こんな種毎に集計されているデータがある時
# (tidyverse標準のtibble形式で作成)
data1 <- tibble(
year = c(rep("y08", 3), rep("y09", 3)), 
site = c(1:3, 1:3), 
depth = seq(1, 5, length=6), 
sp1 = c(6:1), sp2 = c(1:6), sp3 = c(0:5))

# 出力するとこんな感じです。<chr>は文字型、<int>は整数型、<dbl>は連続変数型を示す。
> data1
# A tibble: 6 x 6
  year   site depth   sp1   sp2   sp3
  <chr> <int> <dbl> <int> <int> <int>
1 y08       1   1       6     1     0
2 y08       2   1.8     5     2     1
3 y08       3   2.6     4     3     2
4 y09       1   3.4     3     4     3
5 y09       2   4.2     2     5     4
6 y09       3   5       1     6     5



########################################################

# data1をデータフレーム型に縦に伸ばす場合
data2 <- data1 %>% pivot_longer(c(-year, -site, -depth), names_to="species", values_to="abundance")
# pivot_longerは集計された種名の表を縦長に伸ばす。種名を表す新たな種をspecies, 個体数をabundanceとする。
# 形式を保持したい変数は、c(-year, -site, -depth)のように記す(c()は必要)。ひとつだけならば-yearのように書けば良いが、複数の場合は-c(year, site, depth)ではないことに注意。一般的なRの記法と異なっており、厄介な点。
# "%>%"はパイプと呼ばれる。パイプの前のデータを後の処理に渡す意味。いちいちdata1$yearのように書かなくていい、with関数やattach関数と似ている。可読性は良くなっているかもだけれど、デバッグがしづらい気もする。

> data2 # 出力してみます
# A tibble: 18 x 5
   year   site depth species abundance
   <chr> <int> <dbl> <chr>       <int>
 1 y08       1   1   sp1             6
 2 y08       1   1   sp2             1
 3 y08       1   1   sp3             0
 4 y08       2   1.8 sp1             5
 5 y08       2   1.8 sp2             2
 6 y08       2   1.8 sp3             1
 7 y08       3   2.6 sp1             4
 8 y08       3   2.6 sp2             3
 9 y08       3   2.6 sp3             2
10 y09       1   3.4 sp1             3
11 y09       1   3.4 sp2             4
12 y09       1   3.4 sp3             3
13 y09       2   4.2 sp1             2
14 y09       2   4.2 sp2             5
15 y09       2   4.2 sp3             4
16 y09       3   5   sp1             1
17 y09       3   5   sp2             6
18 y09       3   5   sp3             5


########################################################
# もう一段階、伸ばしてみましょう。一個体あたり一行というデータセットへ。多項ロジットなど、カテゴリカルな解析に便利そうです。

data3 <- select(data2, -abundance) %>% slice(unlist(map2(1:nrow(data2), data2$abundance, rep)))
# data2のabundanceの数に応じてrepで行数を増やし、sliceで行を選択。map2はおおよそmapplyに相当、2変数を取るsapply。slice(unlist(map2(,,,,)))はもう少しスッキリ書けないものだろうか。

> data3 
# A tibble: 57 x 4
   year   site depth species
   <chr> <int> <dbl> <chr>  
 1 y08       1   1   sp1    
 2 y08       1   1   sp1    
 3 y08       1   1   sp1    
 4 y08       1   1   sp1    
 5 y08       1   1   sp1    
 6 y08       1   1   sp1    
 7 y08       1   1   sp2    
 8 y08       2   1.8 sp1    
 9 y08       2   1.8 sp1    
10 y08       2   1.8 sp1    
# … with 47 more rows




########################################################

# data3を一つ前の状態に戻してみます(観測地点あたりで集計、といったところ)
data4 <- data3 %>% count(year, site, depth, species, name="abundance") 
# year, site, depth, species毎にデータ頻度をカウントし、カウント結果をabundance変数として追加

> data4

# A tibble: 17 x 5
   year   site depth species abundance
   <chr> <int> <dbl> <chr>       <int>
 1 y08       1   1   sp1             6
 2 y08       1   1   sp2             1
 3 y08       2   1.8 sp1             5
 4 y08       2   1.8 sp2             2
 5 y08       2   1.8 sp3             1
 6 y08       3   2.6 sp1             4
 7 y08       3   2.6 sp2             3
 8 y08       3   2.6 sp3             2
 9 y09       1   3.4 sp1             3
10 y09       1   3.4 sp2             4
11 y09       1   3.4 sp3             3
12 y09       2   4.2 sp1             2
13 y09       2   4.2 sp2             5
14 y09       2   4.2 sp3             4
15 y09       3   5   sp1             1
16 y09       3   5   sp2             6
17 y09       3   5   sp3             5

# sp3 のゼロ個体データだけは消えてしまったが、当然と言えば当然




########################################################

# data4の種数部分を集計表状に横に伸ばす(最初の状態=data1の形にする場合)

data5 <- data4 %>% pivot_wider(names_from=species, values_from=abundance, values_fill=0)
 

# pivot_widerは縦長データから横長データへ変換する、見出しとなるspeciesと対応する値abundanceを展開する。生じるNAはvalues_fill=0で置き換えた。

> data5 # 出力してみます。
# A tibble: 6 x 6
  year   site depth   sp1   sp2   sp3
  <chr> <int> <dbl> <int> <int> <int>
1 y08       1   1       6     1     0
2 y08       2   1.8     5     2     1
3 y08       3   2.6     4     3     2
4 y09       1   3.4     3     4     3
5 y09       2   4.2     2     5     4
6 y09       3   5       1     6     5



########################################################
# 出現する種のリストにひどく変動がある場合など、1つのセルの中に種名を羅列したくなるケースもあると思います。そういうデータからの変換例も追加しておきます。

data6 <- tibble(site=paste0("s", 1:10), month=rep(c(1:2), each=5), 
species=c("スズメ", "","スズメ、ヒヨドリ、シジュウカラ", "ムクドリ、スズメ", "ヒヨドリ、スズメ", "スズメ、キジバト", "スズメ、ムクドリ", "スズメ、ヒヨドリ", "", "ムクドリ"))

# 例えばこんなかんじのデータセット。空欄さえあります…
# A tibble: 10 x 3
   site  month species                         
   <chr> <int> <chr>                           
 1 s1        1 "スズメ"                        
 2 s2        1 ""                              
 3 s3        1 "スズメ、ヒヨドリ、シジュウカラ"
 4 s4        1 "ムクドリ、スズメ"              
 5 s5        1 "ヒヨドリ、スズメ"              
 6 s6        2 "スズメ、キジバト"              
 7 s7        2 "スズメ、ムクドリ"              
 8 s8        2 "スズメ、ヒヨドリ"              
 9 s9        2 ""                              
10 s10       2 "ムクドリ"                      

# 注:このような日本語データのでの入力の場合、普通に読み込もうとすると文字化けしやすいので、readxlパッケージのread_excel関数などで読み込むのがトラブルが少なくていいと思っています。csvでは機種依存文字やmac/windows間での文字コード問題が混乱のもとなので。

# こんなデータでも次のようにすれば、通常の集計データに変換できます。data1の形態へ変換してみましょう。

sptabs <- data6 %>% # data6をベースに集計作業
pull(species) %>% # 種名部分をベクトルとして抽出
map(str_split, pattern="、") %>% # "、"で文字列を分割
map_dfr(table) %>% # 要素毎に集計テーブルを作成
map_dfr(replace_na, 0) %>% # naを0で置き換え
select(-starts_with("...")) # 出現無しを削除
data7 <- select(data6, -species) %>% # 元データの見出しを抽出
mutate(sptabs) # 種集計テーブルと結合

data7 # 出力してみます(出来上がり)

# A tibble: 10 x 7
   site  month スズメ  シジュウカラ ヒヨドリ ムクドリ キジバト
   <chr> <int> <table> <table>      <table>  <table>  <table> 
 1 s1        1          1       0            0        0        0       
 2 s2        1          0       0            0        0        0       
 3 s3        1          1       1            1        0        0       
 4 s4        1          1       0            0        1        0       
 5 s5        1          1       0            1        0        0       
 6 s6        2          1       0            0        0        1       
 7 s7        2          1       0            0        1        0       
 8 s8        2          1       0            1        0        0       
 9 s9        2          0       0            0        0        0       
10 s10       2          0       0            0        1        0