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

（下記は旧版です。新たにtidyverse版を別ページに作成しました）

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

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

（当初、使用法の面倒なreshape関数を使いやすくするための拡張関数を書いていたのですが、より使いやすい関数を発見したので差し替えました 2015.01.21）


# 例えば、こんな種毎に集計されているデータがある時、
data1 <- data.frame(
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))

# （出力するとこんな感じです）
> data1
year site depth sp1 sp2 sp3
1 y08 1  1.0     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.0     1     6     5


データの整形にはreshapeパッケージの関数、meltとcastを使用します。reshapeパッケージは事前にインストールが必要だったはず。



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

# data1をデータフレーム型に縦に伸ばす場合
require(reshape) # 要インストール
data2 <- melt(data1, id.vars=c("year", "site", "depth")) # id.varsは、いじらずに残したい変数を指定

> data2 # 出力してみます
year site depth variable value
  1 y08 1  1.0    sp1     6
  2 y08 2  1.8    sp1     5
  3 y08 3  2.6    sp1     4
  4 y09 1  3.4    sp1     3
  5 y09 2  4.2    sp1     2
  6 y09 3  5.0    sp1     1
  7 y08 1  1.0    sp2     1
  8 y08 2  1.8    sp2     2
  9 y08 3  2.6    sp2     3
10 y09 1  3.4    sp2     4
11 y09 2  4.2    sp2     5
12 y09 3  5.0    sp2     6
13 y08 1  1.0    sp3     0
14 y08 2  1.8    sp3     1
15 y08 3  2.6    sp3     2
16 y09 1  3.4    sp3     3
17 y09 2  4.2    sp3     4
18 y09 3  5.0    sp3     5



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

data3 <- data.frame(lapply(data2, function(i) rep(i, data2[,"value"])))
data3 <- data3[,-ncol(data3)]
# data2の個体数列"value"を個体数分縦に伸ばした（二行目のコードは不要な列を除去している）

> data3 # 長いので30行目より後は省略！
year site depth variable
  1 y08 1  1.0   sp1
  2 y08 1  1.0   sp1
  3 y08 1  1.0   sp1
  4 y08 1  1.0   sp1
  5 y08 1  1.0   sp1
  6 y08 1  1.0   sp1
  7 y08 2  1.8   sp1
  8 y08 2  1.8   sp1
  9 y08 2  1.8   sp1
10 y08 2  1.8   sp1
11 y08 2  1.8   sp1
12 y08 3  2.6   sp1
13 y08 3  2.6   sp1
14 y08 3  2.6   sp1
15 y08 3  2.6   sp1
16 y09 1  3.4   sp1
17 y09 1  3.4   sp1
18 y09 1  3.4   sp1
19 y09 2  4.2   sp1
20 y09 2  4.2   sp1
21 y09 3  5.0   sp1
22 y08 1  1.0   sp2
23 y08 2  1.8   sp2
24 y08 2  1.8   sp2
25 y08 3  2.6   sp2
26 y08 3  2.6   sp2
27 y08 3  2.6   sp2
28 y09 1  3.4   sp2
29 y09 1  3.4   sp2
30 y09 1  3.4   sp2




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

# data3を一つ前の状態に戻してみます（観測地点あたりで集計、といったところ）
data4 <- data.frame(xtabs(~ ., data=data3)) # data3の最後の列が集計される→Freqという列になる。xtabsの~の後のピリオドの部分は必要な変数のみでモデル式のように書いてもよい（~ year + site + depth + variable）
data4 <- data4[data4$Freq>0,] # 0が沢山出てしまうので除去する

> data4

year site depth variable Freq
    1 y08 1  1       sp1    6
    9 y08 2  1.8    sp1    5
  17 y08 3  2.6    sp1    4
  20 y09 1  3.4    sp1    3
  28 y09 2  4.2    sp1    2
  36 y09 3  5       sp1    1
  37 y08 1  1       sp2    1
  45 y08 2  1.8    sp2    2
  53 y08 3  2.6    sp2    3
  56 y09 1  3.4    sp2    4
  64 y09 2  4.2    sp2    5
  72 y09 3  5       sp2    6
  81 y08 2  1.8    sp3    1
  89 y08 3  2.6    sp3    2
  92 y09 1  3.4    sp3    3
100 y09 2  4.2    sp3    4
108 y09 3  5       sp3    5

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




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

# data4の種数部分を集計表状に横に伸ばす（最初の状態＝data1の形にする場合）
# 今度はreshapeパッケージのcast関数を使用します。

require(reshape) # 要インストール
data5 <- cast(data4, year + site + depth ~ variable, value="Freq") 
# いじらない変数を左辺に + で並べ、インデックスに用いる変数variableを ~ の右辺に設定。逆に、集計したい変数Freqはvalueの項で指定（指定しなくても、差し当たり自動で選択してくれる模様）

> data5 # 出力してみます。
year site depth sp1 sp2 sp3
  1 y08 1   1      6      1   NA
  9 y08 2   1.8   5      2   1
17 y08 3   2.6   4      3   2
20 y09 1   3.4   3      4   3
28 y09 2   4.2   2      5   4
36 y09 3   5      1      6   5

# 0個体がNAになるのが困る場合は、data5$sp3[is.na(data5$sp3)] <- 0 するなどで修正




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

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

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

# 例えばこんなかんじのデータセット。空欄さえあります…
   site month                        species
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                       ムクドリ

#（注：実際にはread.csv関数でデータシートを読み込むと思いますが、このような日本語データのでの入力の場合、普通に読み込もうとすると文字化けするので、read.csv(でーた, fileEncoding="cp932") のようにエンコーディング指定してやると良いです）。

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

break <- "、" # 区切り文字の指定
slists <- lapply(1:nrow(data6), function(s) unlist(strsplit(as.character(data6$species[s]), break))) # 出現種文字列を分解
s.level <- sort(unique(unlist(slists))) # 全出現種をリストアップしソートする
s.table <- t(sapply(1:length(slists), function(s) table(factor(slists[[s]], level=s.level)))) # サイト毎に集計
data7 <- cbind(subset(data6, select= -species), s.table) # 元データのspecies部分をs.tableに入れ替え

data7 # 出力してみます（出来上がり）

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