R言語のサンプルデータセットを端から味わってみる、やっつけでご乱心な第３回はBODです。たまに水質汚染などのニュースで聞く生物化学的酸素要求量というもの。なんだか分かりませぬが、第３回にしてようやく時系列データを離れ、データフレームが登場しました。目出度い（何が？）
※「データのお砂場」投稿順Indexはこちら
BOD
さてRの datasets パッケージ内の大量のデータセットには、データの出所から処理のサンプル例までいろいろ情報があるものもあれば、前回のように出所が分かるだけで何やら不明なデータもあるようです。今回のBODについては処理のサンプルもある「ありがたい」方みたいです。説明は以下に。
BOD {datasets} – R Documentation
しかし、なんだか意味不明なデータならともかく、BODという名を聞いてしまいました。その名は聞いたことがある。「デカイと水質悪いんでしょ」的な理解でいましたが、実際のところはサッパリです。
BOD — Biochemical Oxygen Demand
日本語で、生物化学的酸素要求量と唱えてもサッパリです。ネットに聞くのが速いということで検索すると下記のページがヒットいたしました。IHI様、日本を代表する「重い」会社です。
上記ページより、BODの定義部分を引用させていただきましょう。
好気性バクテリアが一定時間内（20℃で5日間）に水中の有機物を酸化・分解し、浄化するのに消費される酸素の量をppm（あるいはmg/Lit.）で表したもの
そうでしたか。酸素の消費量は時間の関数なのですね。そして５日間と期間が決めてあると。
全くの蛇足ですが、IHIのH,この歳になるまで石川島播磨の「播磨」のHだと思いこんでいました。つい最近、YouTube見ていて、HはHarimaのHではなく、Heavy IndustryのH、だということを知りました。重い方のHなのね。
データを眺めてみる
まずは、データをロードして、その形式と生データを観察してみたいと思います。こんな感じ。
R言語に特有なデータフレーム形式のデータでした。しかしデータ点数すくな。６点かい。そしてTimeの単位は日みたいです。そしてdemandの単位は mg/lとな。
非線形回帰分析
データセットの説明の処理例を読むと、非線形回帰分析をやれ、ということみたいでした。３つのモデル？について何やら処理をするみたい。最初のモデル(first-order model)からして、何やら難し気な式？といって良く見ていたら
RC時定数と同じじゃん
ということに気づきました。定数にはきっとBODの流儀があるのですが、式の形は、RCフィルタとおなじじゃん、です。気を取り直して、R言語で非線形回帰分析を行うための nls関数を、処理例通りに打ち込んでみました。
coef()関数を上のように使うと、求めた定数が判明します。またnlsの結果には上記のようにいろいろ情報が含まれているようでした。
プロットしてみる
ここまで出来ると、生データに重ねて回帰曲線を描きたい、と思いました。上記の処理例には回帰直線の描き方はなかったのですが、元の式があり、係数が求まっているのだから、
- 関数を定義（しれっと実係数を書き込んでしまう）
- 元のデータの散布図を描く
- それに重ねて上で定義した関数のグラフを描く
ってな方法でいけるのでは。やってみた操作が以下に（これに一撃でたどりついた訳ではありませんよ。）
> funcFM1 <- function(x){return(19.1426*(1-exp(-exp(-0.6328215)*x))}
> plot(BOD$Time, BOD$demand)
> par(new=TRUE)
> plot(funcFM1, 1, 7, type="l", col="red", lwd=2, xlab="", ylab="")
描いたグラフがこちら。
まあ、そんなものかい。でも、抵抗とコンデンサに方形波入れる方が私には分かり易い。BODの専門家でもないので、First-Order以降の精密にするやつはパス。なんと安直な。