エクセルシートやcsvファイルにまとめた数値データはエクセル上で関数を使うなどして前処理してしまうことが多いですが、エクセル上の操作は記録の残らないマウス操作に依存するため、再現性がまったくありません。再現性を確保するため、可能な限りR上で前処理を行うことを強くお勧めします。
以下の内容は基本操作のごく一部です。ほとんどすべての情報は、以下のサイトとその内容をまとめた書籍に掲載されていますので自習することができます。
まずは準備として、この実習用に、Qドライブ直下に以下の二つのファイルをダウンロードしましょう。さらに、以下の解析に出てくる解析を実際に行るRスクリプトもQドライブ直下に作りましょう。(1)フォルダを作成し(たとえばphytoという名前にしてみる)、(2)このdataというフォルダの中に以下の二つのファイルをダウンロードしましょう(右クリックではダウンロードできないので、まずは左クリックでファイルをウェブブラウザにそのまま表示し、左上角のダウンロードボタンを押してダウンロードしてください)。さらに、以下の解説に出てくる解析を実際に行うRスクリプトを同じフォルダ(phyto)の中に作りましょう。たとえば、biwa_phyto.Rみたいな名前にしたらよいと思います。
上の二つのデータファイルを、それぞれ異なるデータフレームに保存してみます。
#load the data sets
data_biwa20180924<-read.csv("20180924_1_biwa.csv", header=TRUE)
data_biwa20181025<-read.csv("20181025_2_biwa.csv", header=TRUE)
ここらへんで分からなくなる人が多いと思いますが、data_biwa20180924, data_biwa20181025というのは、Rの中でデータを加工したり計算に使ったりするためにデータを保存しておくハコ(オブジェクト)で、データフレームという形式のものです。data_biwa20181025というデータフレームの名前を使う(スクリプトに書き込む)ことによって、そのデータフレーム内のデータにアクセスできるようになります。それに対し20180924_1_biwa.csv、20181025_2_biwa.csvはPC上に保存されているファイルで、そのファイルの名前です。PC上からならどこからでもアクセスできるファイルと、Rのスクリプト上だけでアクセス可能なデータフレームを混同しないように注意してください。ファイルにアクセスするにはファイル名を指定する一方、R上のオブジェクトの一つのタイプであるデータフレームにアクセスにするにはデータフレームの名前を指定します。
データフレームの中身をチェックしてみると、植物プランクトンのタイプ(green, diatomなど)と種名(species)と個体数(群体数)(abundance = 細胞数もしくは群体数)の一覧であることが分かります。
View(data_biwa20180924)
#Check the part of the data set
data_biwa20180924$type
RstudioのRエディターでは、data_biwa20180924$
まで書き込むとそれ以降の候補を自動表示してくれるので便利です。
あるいは、以下のようにすると各行(1行目)、各列(一列目)、特定の行-列(3行め、2列め)の値にそれぞれアクセスできます。
#access to subdata
data_biwa20180924[1, ]
data_biwa20180924[, 1]
data_biwa20180924[3,2]
次に、二つのデータフレームを結合したいと思います。というのは、***0924, ***1025のファイルはそれぞれ2018年9月24日と10月25日の植物プランクトン組成であり、ゆくゆく二つの組成を比較したりするときに一つのデータフレームに統合することが必要になってくるからです。
一連の作業は以下のようなスクリプトで行ないます。最初の二行は、各月のデータフレームの各項目に名前を付けます。重要なのは3つ目の項目に、サンプルごとに別々の名前(この例では、cells0924とcells1025) をつける一方、最初の二つは共通の名前を付けることです(typeとspecies)。その次にmerge()という関数を使ってデータを結合し、data_biwa_combinedという名前の新たなデータフレームに結果を格納します。
#combine
colnames(data_biwa20180924)<-c("type","species","cells0924")
colnames(data_biwa20181025)<-c("type","species","cells1025")
data_biwa_combined <- merge(data_biwa20180924, data_biwa20181025, all=T)
新しいデータフレーム中身の中身を看てみると、9月と10月で共通していないspeciesの部分には、欠損データ(NA)が挿入されていることが分かります。
View(data_biwa_combined)
NAは、数字ではないので後々の統計解析の時のエラーの下になるので、観察されなかった=個体数 0ということで、値をゼロにします。細かいことを言うと、is.na()という関数はインプットのオブジェクト(データフレームやベクトル等)の各要素がNAかどうかを判定して、NAなら真(TRUE)、NAでないなら偽(FALSE)を返す関数です。これをデータフレーム名[]の中に埋め込むと、TRUEの要素だけにアクセスし、ゼロを代入します(<- 0)。これによって値がNAだった要素がすべてゼロに置き換わります。
#convert all NA to zero
data_biwa_combined[is.na(data_biwa_combined)]<-0
さて、3つ以上のデータを結合したい時は、merge()関数を何回も使う必要があります。なぜならこの関数は二つのデータしか統合できないからです。仮想的に、data1, data2, data3という3つのデータを上と同じ方法で結合させたいなら以下のスクリプトみたいな感じで書きます。具体的なスクリプトは状況に合わせて自分で考えて書きましょう。
combined <- merge(data1, data2, all=T)
combined <- merge(combined, data3, all=T)
以下のデータファイルを全て読み込んで、結合すれば、データの前処理の第一段階は終了です。
ただし種名のスペルミス等があるかもしれません。一文字違いの種名とかがあった場合は、以下の藻類データベースで検索して確かめましょう。
ここでは、自分たちの計数した植物プランクトンのリストを1)のところで例として用意したものと同じようにPC上のファイルとして整理します。最終的には、1)で結合したデータフレームにさらに追加することを念頭に置いて、表のフォーマットに気を付けることが重要です。ファイルを作成しているワークフローは以下の通り5ステップとします。
[1] 植物プランクトンの集計用紙からエクセルファイルへ:各個人でMS Excelを開いて自分自身の計数結果をまとめてエクセル形式で保存する。
[2] エクセルファイルからcsvファイルへ変換:MS Excel上で、エクセルファイルをRに読み込ませるためにcsv形式のファイルに変換する。
[3] グループ内で各speciesごとに合計値を計算:Rを使ってグループ内では観測された種それぞれの合計数を計算する。
[4] すべてでのグループの集計結果と1)のデータを結合:Rを使ってグループ間のデータも結合する。
[5] [4]で作ったデータフレームを加工:異なるデータを比較可能なようにデータフレームに情報を追加する
以下はこの5ステップに関する補足とRスクリプトの解説となります。
このようなファイルを、ひとつの研究プロジェクトでは、複数のサンプルごと、ひとつのサンプルの培養時間ごとに数多く作ることになります。ファイルの準備で重要なことは以下の2つです。
[2.1] フォルダ名・ファイル名の一貫性
たとえば、"サンプル日時-処理名orサンプル地点-測定時間.csv"などのような感じです。
[2.2] 日本語フォントおよびスペース禁止
ファイル名は数字とアルファベットだけ使いましょう。またスペースをファイル名の中に含めるのは絶対にやめましょう。単語を区切りたいときはハイフン(-)やアンダーバー(_)を使いましょう。
[2.3] ファイルの中身のフォーマットの一貫性
1)のデータの前処理1のところでmerge関数で簡単に複数のファイルの情報を統合できましたが、これが可能なのは、元の各ファイルの中身の形式(フォーマット)に一貫性があるためです。この一貫性を最初から意識してエクセルファイルやcsvファイルを用意することが必要です。
このステップでは、グループ内のデータを一つの観測値として統合するために、各人がカウントしたデータの合計を計算します。その仮想例として、1)のところで作ったdata_biwa_combinedというデータフレームを加工してみます。本当は、2018年9月24日の結果と10月25日の結果は、独立した別個の観測結果のため、各種の個体数を合計することには生態学的な意義はありませんが、計算練習のために二つの観測結果の個体数を合計してみましょう。そのためには合計したい列(cells0924とcells1025)を呼び出して、単純に足し算を行い(+)、その結果を新しいオブジェクト(sum)に格納します。
#When calculating the total
sum <- data_biwa_combined$cells0924 + data_biwa_combined$cells1025
cbind.data.frameという関数を使って、足し算結果(sum)を、新しい列(column)としてdata_biwa_combinedに追加し、新しいデータフレームdata_biwa_combined2
として保存します。
data_biwa_combined2 <- cbind.data.frame(data_biwa_combined, sum)
追加した結果はView(data_biwa_combined2)
とすれば以下のようになっているはずです。
さらに、個別のデータ(cells0924とcells1025)を残しておく必要が無いなら、data_biwa_combined2から3列目と4列目を削除します。データフレームの行と列の指定は、[行, 列]の順であることに注意してください。マイナス記号(ー)をつけるとその行や列を削除することができます。
#remove uncessary columns
data_biwa_combined2 <- data_biwa_combined2[,c(-3,-4)]
さて、上の解説はあくまで計算の仮想例ですので、この実習ではまずグループ内で各メンバーが計数結果をcsvファイルにまとめてメンバー間で共有して、Rに読み込み、それらをグループ内で集計する(各種ごとに合計数を計算する)必要があります。全体のワークフローを図にすると以下のようになり、グループ内で集計する部分は青い点線で囲った部分の作業になります。この部分のスクリプトはうえで解説した関数を組み合わせれば書くことができるのでやってみましょう。
[3] の部分が完了していると、各グループでそれぞれひとつ、植物プランクトンの集計結果をまとめたデータフレームが完成しているはずです。これをグループ間で共用し、過去のデータと合わせて結合させるのがこの部分の目標になります。
個人間でRのオブジェクトを共有するには、二つの方法があります。一つは、write.csv()という関数でPC上(つまりRの外)にcsvとしてデータを書き出し、個人間でそれらのcsvファイルを共有する方法です。これは多くの場面で有効なのですが、Rのオブジェクト→csvファイルへの書き出し→csvファイルの読み込み、という流れで回り道なので、もう一つの方法をここでは使います。saveRDSという関数をsaveRDS(書き出したいRオブジェクト名、”書き出しファイル名.rds”)の形式で使うとPC上(Rの外)にRのオブジェクトをそのまま保存できます。たとえば、以下の感じです。
#save the result
saveRDS(data_biwa_combined, "biwa_data2018.rds")
逆にPC上に保存されたオブジェクトファイルをRに読み込むには、readRDS関数を用いて、以下のようにします。
data_biwa_2018 <- readRDS("biwa_data2018.rds")
以上の解説の内容をよく吟味して、上の「データ整理のワークフロー」の残りの部分、今回のデータの別のグループの集計結果および過去のデータとの結合をR上で実行してみましょう。
ここでは、上のほうで作ったデータフレームdata_biwa_combinedを加工する方法を解説します。実際に解析する際にはこのデータフレームよりもたくさんのデータを結合したファイルを使うことになると思うので適宜読み替えてください。スクリプトがすでに変更されていてこのデータフレームがどこに行ったか分からない人は改めでダウンロードし、以下のスクリプトで読み込んでください。
ダウンロード:data_biwa_combined.rds
data_biwa_combined <- readRDS("data_biwa_combined.rds")
[5.1] Typeごとに集計しTypeごとのデータだけにまとめる(発展的内容)
aggregateという関数を使うと指定した列の情報を用いてデータを集計することができます。今回のプランクトンのデータには、各種にType (緑藻=green、珪藻=diatom、藍藻=cyano、鞭毛等を持ち移動性のもの=motile)が割り当てられています。このTypeごとに合計個体数を計算することが可能です。以下のスクリプトで引数の一つ目xはどの列のデータを集計に使うかを指定し、次の引数byは何を基準に集計するかを指定し、最後の引数FUNは集計の仕方の関数(FUNction)を指定しています。
#data summarization is possible
aggregate(x=data_biwa_combined[c("cells0924", "cells1025")], by=list(data_biwa_combined$type),FUN=sum)
その結果は次の表にコンソール(console)に表示されるはずです。
Group.1 cells0924 cells1025
1 cyano 62 7
2 diatom 56 55
3 green 31 638
4 motile 0 20
この集計結果はTypeレベルに注目した統計解析用に、新しいデータフレームに保存して列名を新たに付けておきましょう。
type_aggregate_data <- aggregate(x=data_biwa_combined[c("cells0924", "cells1025")], by=list(data_biwa_combined$type),FUN=sum)
colnames(type_aggregate_data) <- c("type", "cells0924", "cells1025")
[5.2] 種ごとのデータだけまとめる
種レベルに注目して統計解析をするために、Typeの情報は取り除いてしまいましょう。
#remove type data
species_biwa_data <- data_biwa_combined[,c(-1)]
[5.3] データの整形と情報の追加
いままでのところで、種のリストのデータフレーム(species_biwa_data)とTypeのリストのデータフレーム(type_aggregate_data)が完成したはずです。これに少し修正を加えます。まず、データフレームの列の一部となっている種名のリストを行(row)の名前にコピーします。
#copy species name to row names
rownames(species_biwa_data) <- species_biwa_data$species
これにより種名をデータフレームの要素として維持する必要はなくなったので、種名が保存されている一列目は削除します。
species_biwa_data <- species_biwa_data[,-1]
さらに、このデータフレームでは表の縦方向(行方向)に種が並んでおり、表の横方向(列方向)に異なるサンプルの計数結果が並んでいる形になっています。ここでは、今後の統計解析の準備として、この縦と横の並びを逆転されます。表の横と縦を並び替える作業は、線形代数の用語で転置(transpose)といい、R上では簡単にt()という関数でできます。ついでに単に転置するとデータフレームでなくなってしまうのでas.data.frameという関数も同時に使います。
#transpose
species_biwa_data <- as.data.frame(t(species_biwa_data))
エラー処理メモ
ここでひと段落なのですが、上のようなデータをコピーしたり縦横を変えたりしていると、データの種類(数字numericか文字characterか要因factorか、など)が勝手に変わってしまうことがあります。一番注意しないといけないのは、見た目は数字のままなのにcharacterかfactorとしてRに認識されてしまう場合です。こうすると、数字に対して行う計算や統計などがうまくいきません。したがって、このようなエラーを防ぐためにclass関数を使ってデータの種類をチェックします。class(species_biwa_data[1,1])
の結果がnumericでないと困ります。
#check class
class(species_biwa_data)
class(species_biwa_data[1,1])