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

同じ方法でspecies_biwa_data, species_biwa_data2, type_compositionにもメタデータを列として追加してみましょう。

[1-v] 整理したデータをオブジェクトファイルとしてPC上に保存する（発展的話題）

発展的話題なのでやり必要はありませんが、できあがったデータはsaveRDS()という関数を使うと、R上で保持している情報を丸々保存したまま、PC上に保存しあとで呼び出して使うことができます。

#save objects

saveRDS(species_biwa_data, "./species_biwa_data.rds")

２）比較のためのグラフ化（見える化）の基礎

[2-i] 散布図

まずはすべてのサンプルデータのばらつきを明示的にみるために散布図（scatter plot）を描いてみましょう。ここでは例として０次の有効種数（種数）esn0に注目します。ほかの有効種数については以下のサンプルコードを参考に自分でコードを書いて実行してみましょう。

#scatter plot

plot.default(species_biwa_data3$date_sampling, species_biwa_data3$esn0)

図は以下のようになるはずです。グラフの横軸はspecies_biwa_data3$date_samplingが要因（factor）として数値化されて並んでいますが、この実習はきれいなグラフを書くことを目的としていないのでとりあえずこれで良しとします。きれいなグラフが書きたかったらggplot2というキーワードでググればいろいろ情報はでてきます。

横軸の数値の見た目が違う場合は、以下のコマンドを試してください。

#scatter plot

plot.default(as.factor(species_biwa_data3$date_sampling),species_biwa_data3$esn0)

[2-ii] 箱ひげ図

もうすこしだけきれいな図は箱ひげ図（box plot）として以下のコマンドで書けます。ただし分布の情報がかなり失われてしまうため、最近はあまり使用が奨励されていません。代わりにviolin plotというのがあるのでこのサイトの解説を参考にするとよいかもしれません。箱ひげ図を描く時のポイントは、plot関数の中にY~Xの形でモデル式（model equation, formula）を書き込む必要があります。Yは描画したい値、Xは箱ひげ図で分類したい要因をデータフレームの列名で指定する必要があります。dataオプションで、データ元のデータフレームを指定します。

#box plot

plot(esn0~date_sampling,data=species_biwa_data3)

グラフは以下の左図のような感じです。横軸がちゃんと文字になっていると思います。ちなみにグラフを大きく生じしないとすべての横軸の文字は表示されないので注意です。ここで箱ひげ図が表示されなかった場合は、plot()の代わりにboxplot()を使ってみてください。

また、箱ひげ図のplotのコードは次のような書き方でも結果は一緒です。この書き方の方が、モデル式のYとXのデータが別々のベクトルやデータフレームに格納されている場合でも使えます。

plot(species_biwa_data3$esn0~species_biwa_data3$date_sampling)

箱ひげ図を描く要因を月（month_sampling）に変えると、以下のようなコマンドで、右図のようなグラフが描写できます。ちょっとだけオプションを追加して見栄えを変えました。レポートを書く際は最低限これらのx軸とy軸のラベルを指定するオプション（xlab="hogehoge", ylab="hogehoge"）を追加して作図してください。

plot(esn0~month_sampling,data=species_biwa_data3,xlab="month", ylab="species richness (esn0")

３）いろいろグラフを描いてみる

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

さて、以上、ゼロ次の有効種数(esn0)というパラメータについてのグラフ化（見える化）を解説しましたが、他の次数の有効種数(esn1, esn2)、元の個体数データ(species_biwa_data)、割合データ（species_biwa_data2）、タイプ別個体数データ(type_composition)のデータを使って、いろいろグラフ化してみましょう。たとえば、珪藻の相対個体数に月間で差はあるのかとか、特定の種（アオコを形成する種や外来種など）に注目してその相対個体数に月間で差はあるのかなど、自分でコードを書いてグラフを作ってみましょう。

[3-i] 図の保存の仕方

方法１：Plots > Export > Save as Image: これでファイルとして保存できます。

方法２：Plots > Export > Copy Plot to Clipboard > Copy Plot: これでクリップボードに画像が一時保存（コピー）されるのでワードを開いて貼り付けることができます。

方法３：Plots > Zoom > 右クリック：これでファイルで保存したり、クリップボードにコピーしたりできます。

４）グラフ上の見た目の差は偶然の産物ではないのか？ ーー仮説検定の方法

[4-i] 統計方法の復習

解説ファイルでt検定、分散分析、線形回帰モデルを復習しておくこと。

さらに勉強したいひとは、キーワードは一般線形モデル(general linear model)と分散分析(Analysis of Variance, ANOVA)です。

以下の解説サイトがお勧めです。

橋本洸哉さんのRで統計解析

[4-ii] サンプル月日間、サンプル月間で差があるものは何か？

仮説検定は卒論では必須のプロセスです。これまでの週の教材や課題でもすでに取り組んでいると思いますが、ここでも分散分析によってサンプル間に差があるか検定してみましょう。たとえば、上でグラフ化したところ、ゼロ次の有効種数(esn0)には月間に差があるように見えます。これを分散分析で確認してみましょう。

#test

oneway.test(species_biwa_data3$esn0~species_biwa_data3$month_sampling, var.equal=F)

分析の結果は以下のようにコンソールに表示されるはずです。

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

３）のところでいろいろグラフを作ってみたと思いますが、それに対応した分散分析をしてみましょう。

[4-iii] 月間や月日間でパターンに差があるとすればそれはなぜか？

[4-ii] のところで差が見られたパターンがあるとすれば、それはなぜでしょうか？

琵琶湖の栄養塩の量とか風速とか天気（光の量とか）とか水温とかの違いかもしれません。だとすると、それらの環境パラメータをどこかから見つけてきて、散布図を描いたり線形回帰モデルを適用して検定したりできるはずです。

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

たとえば、水温なら「琵琶湖水温 on Twitter まとめ」というサイトのデータをとりあえず信じて、プランクトンデータを取った年月日の水温を調べましょう。

他の環境要因で説明したい！ というひとはいろいろグルると公開データ出てくるはずです。

いずれにせよ、[1-iii]でメタデータを作ったときと同じように、サンプル数（標本数）17と同じ長さのベクトルを用意する必要があります。たとえば、以下のような感じです（仮想例）。20180924の水温が25.5℃、20181004の水温が20.3℃、（以下省略）なら同じサンプリング年月日の繰り返しの数だけ同じ数値を用意する必要があります。

#add environmental information

temperature_biwa <- c(25.5, 25.5, 20.3, 20.3, 20.3,(あと12個の数値))

ここで重要なのはc()で値を指定してベクトルを作るとき、それぞれの数値を""で囲ってはいけません。""で囲った文字はたとえ数字であっても数値ではなく文字として保存されてしまいます。実は、値がcharacterなのかfactorなのかnumericなのかは、実は微妙で頻出のエラーな上に、Rのバージョンアップによってもデフォルトの読み込み設定が変わるので要注意です。

この部分については、最初にエクセルファイルやcsvファイルに保存してから、read.csv()関数を使ってR上に読み込むこともできます。そのときもheaderの有無や、値が数値なのか文字なのか、などに気を付ける必要があります。

たとえば、ゼロ次の有効種数と水温の関係は以下のような散布図で表すことができます。

plot(species_biwa_data3$esn0~temperature_biwa)

上の図から明らかに、水温とゼロ次の有効種数の間には線形な関係はなさそうですが、以下のようなコードで線形回帰モデルによる検定と元データと線形回帰直線の重ね合わせのグラフが描けます。

model_esn0<-lm(species_biwa_data3$esn0~temperature_biwa)

summary(model_esn0)

plot(species_biwa_data3$esn0~temperature_biwa, xlab="water temperature", ylab="esn0")

abline(model_esn0)