前回の記事(細菌コロニー形成曲線 - もうカツ丼でいいよな)は「それで実際のデータはどうやって使えばいいんだよ！！！」という点が欠けていて面白みがないのでR使ってコロニー形成曲線求めてみました。
求めるパラメータ
コロニー形成曲線は次の関係式によって表現されます。
ここでNは時間tにおける平板上のコロニー数、は平板上に最終的に(t=∞で)現れるコロニー数、λはあるコロニー未形成の細菌が単位時間内にコロニーを形成する確率、
は最初のコロニーが出現するまでの時間です。
Nとtは測定によって求められるので、、λ、
を求めるのが目的となります。これらのパラメータは平板上に最初に撒かれた細菌数や、細菌の生理状態、細菌の増殖速度などをそれぞれ反映すると言われますが、詳しいことは記事末尾の参考文献などを参照してください。
未知定数は3つなので、測定値は最低3点必要です。測定は誤差を含むのが当然と考えられます。それを考慮した上でパラメータの推定値を求める方法として線形法と非線形法の2つがあります。
線形法
最初の式を変形すると
となり、とtの間に直線関係のあることが分かります。
よって、をtに対しプロットしたときに、データ点が直線上に乗るような
を選べばいいことになりますが、そのための基準として相関係数を用います。
の値は時間とともに小さくなっていくことが予想されますから、理想的には
と時刻tの相関係数が-1になります(実際には誤差があるのでまず-1にはなりません)。
すなわち、とtの相関係数を最小化するように
を選べばいいことになります。
このような問題を解くとき、Rではoptimize関数が使用できます。実際にやってみましょう。
data <- data.frame(Time = c(1, 2, 3, 4, 6, 11), # 測定時刻(日) Count = c(17.8, 30.8, 46.2, 53.2, 60.6, 67.6)) # コロニー数 optimize( function(n.inf){ cor(data$Time, log(n.inf - data$Count)) }, c(max(data$Count), max(data$Count)*10))
データはCiNii 論文 - METHOD FOR MATHEMATICAL ANALYSIS OF BACTERIAL COUNT DATAのものを借用しました。
関数optimize()は第一引数に指定された関数を第二引数に指定された範囲内で最小化します。
は現在までのコロニー数よりいくらか大きいはずなので、現在までのカウントデータ内の最大値からその10倍までを探索範囲としました。場合によっては上限をもっと上げる必要があるかもしれません。実際には後に述べるような方法でプロットして確認するのが良いでしょう。
上のコードをRで実行すると次のような結果が返ってくるはずです。
$minimum [1] 68.74093 $objective [1] -0.9991266
の推定値は68.74093で、そのときのcor(data$Time, log(n.inf - data$Count))の値は-0.9991266ということを言っています。
なお、他のパラメータに付いてはを決定すればそこから求めることができます。
一々上のように入力するのは面倒なので関数にしておきます。ついでに他のパラメータも計算するようにしてしまいましょう。
linear.method <- function(x, y = NULL){ data <- xy.coords(x, y) n.inf <- optimize( function(n.inf){ x <- data$x # time y <- log(n.inf - data$y) # ln(n.inf - n) # plot(x, y) cor(x, y) }, c(max(data$y), max(data$y)*10))$minimum x <- log((n.inf - data$y)/n.inf) result <- c(n.inf, rev(lm(x~data$x)$coefficients)) result[2] <- -result[2] names(result) <- c("n.inf", "lambda", "Tr") result }
xy.coords()という関数を使うとベクトルの組もデータフレームも(多分他の形式も)受け取って適当な形式に変換してくれるので便利です。
実行例はこんな感じです。
> linear.method(data) n.inf lambda Tr 68.7409332 0.3799968 0.0886838
非線形法
線形法ではデータを対数変換しましたが、要するに次の式
とデータポイントとの偏差の二乗和が最小となるように3つのパラメータを設定してやればいいわけです。
Rでは関数nls()によって非線形最小二乗法を行うことができます。
data <- data.frame(Time = c(1, 2, 3, 4, 6, 11), # 測定時刻(日) Count = c(17.8, 30.8, 46.2, 53.2, 60.6, 67.6)) # コロニー数 nls(Count ~ n.inf * (1 - exp( -lambda * (Time - Tr))), data, start = c(n.inf = max(data$Count), lambda = 1, Tr = min(data$Time)))
第一引数は変数と未知定数を含む表現式オブジェクトです。ここで指定する変数は第二引数に指定するデータフレームの中に入っている必要があります。第三引数には推定値の初期値を名前付きベクトルで指定します。
実行結果は次のようになります。
Nonlinear regression model model: Count ~ n.inf * (1 - exp(-lambda * (Time - Tr))) data: data n.inf lambda Tr 68.8076 0.3830 0.2639 residual sum-of-squares: 11.01 Number of iterations to convergence: 6 Achieved convergence tolerance: 2.331e-07
線形法のときと同様に関数にしてしまいましょう。
nonlinear.method <- function(x, y = NULL){ data <- xy.coords(x, y) data <- data.frame(Time = data$x, Count = data$y) coef( nls(Count ~ n.inf * (1 - exp( -lambda * (Time - Tr))), data, start = c(n.inf = max(data$Count), lambda = 1, Tr = min(data$Time)) ) ) }
この関数の実行例はこんな感じです。
> nonlinear.method(data) n.inf lambda Tr 68.8075950 0.3829985 0.2638863
線形法と非線形法の比較
線形法ではの対数をとる関係上、
の値が小さくなる後半の測定値ほど推定値に及ぼす影響が大きく、前半の測定値は逆に影響が軽くなるという問題があります。Rなら上に示したように簡単に非線形最小二乗法が実行できるのでとりあえず非線形法でやっておけばいいと思います。
ただ折角なので両者を比較してみましょう。
データの用意
引き続きCiNii 論文 - METHOD FOR MATHEMATICAL ANALYSIS OF BACTERIAL COUNT DATAから借用します。
## データ data <- data.frame(sample = rep(1:5, c(6, 6, 5, 5, 7)), Time = c(1, 2, 3, 4, 6, 11, #sample1 1, 2, 3, 4, 6, 11, #sample2 1.708, 2.208, 2.708, 3.063, 3.729, #sample3 1.8, 2.8, 3.8, 4.8, 10.0, #sample4 1.75, 2.75, 3.08, 3.75, 4.75, 5.83, 6.83), #sample5 Count = c(17.8, 30.8, 46.2, 53.2, 60.6, 67.6, 5.4, 10.2, 15.9, 18.0, 24.8, 35.0, 104.4, 118.6, 151.4, 153.3, 174.6, 6.0, 12.2, 15.0, 16.4, 17.8, 11.6, 43.3, 48.8, 51.2, 52.0, 52.6, 56.0)) ## 個々の実験データを取り出す関数を定義 get.sample.data <- function(sample.number){ data[data$sample == sample.number,][,2:3] }
パラメータの比較
線形法と非線形法の関数は先ほど定義したものを使います。
# 線形法 t(data.frame( lapply(by(data[,2:3], data$sample, linear.method), round, 2) )) # 非線形法 t(data.frame( lapply(by(data[,2:3], data$sample, nonlinear.method), round, 2) ))
線形法の結果はこんな感じ。
| n.inf | lambda | Tr | |
|---|---|---|---|
| 68.74 | 0.38 | 0.09 | |
| 46.02 | 0.13 | 0.00 | |
| 225.96 | 0.43 | 0.15 | |
| 17.84 | 0.70 | 0.84 | |
| 56.87 | 0.63 | 0.41 |
非線形法の結果はこんな感じ。
| n.inf | lambda | Tr |
|---|---|---|
| 68.81 | 0.38 | 0.26 |
| 45.64 | 0.13 | 0.03 |
| 256.90 | 0.31 | 0.08 |
| 17.78 | 0.73 | 1.23 |
| 53.86 | 1.44 | 1.58 |
プロットして比較
推定値がどれくらい妥当なのかというのは残差とか見て判断するのだと思いますが、とりあえず今回使ったデータならプロットしてしまえば一目瞭然です。
線形法と非線形法、それぞれにより得られたパラメータから生成されるコロニー形成曲線を実際のデータの上に重ねて描いてみましょう。
サンプル1の結果はこのようになります。
ほとんど差は無いように見えますが、やはり線形法では計測初期の付近でデータポイントからの系統的なずれが目立ちます。
このこの傾向は測定後半に値がブレているようなデータでは特に顕著に現れます。例えばサンプル5のデータでは線形法の結果が明らかにマズイのがわかります。
というわけで実際にデータを処理するならやっぱり非線形法が良いようです。
ちなみにプロットはこんな感じでやりました。
## プロットで比較 gen.curve.func <- function(pars){ # c(n.inf, lambda, Tr) function(x){ pars[1] * (1 - exp( -pars[2] * (x - pars[3]))) } } linear.vs.nonlinear <- function(sample.data, main = ""){ result.linear <- linear.method(sample.data) result.nonlinear <- nonlinear.method(sample.data) plot(sample.data, main = main, xlim = c(0, max(sample.data$Time)), ylim = c(0, max(sample.data$Count)+10)) tmax <- max(sample.data$Time) + 10 time <- seq(0, tmax, by=0.1) lines(time, gen.curve.func(result.linear)(time), col = "red") lines(time, gen.curve.func(result.nonlinear)(time), col = "blue") legend("bottomright", inset = 0.05, legend = c("Linear Method", "Non-linear.method"), col = c("red", "blue"), lty = c(1, 1)) } png() for (i in 1:5){ linear.vs.nonlinear(get.sample.data(i), main = paste("Sample", i)) } dev.off()
参考文献
- CiNii 論文 - METHOD FOR MATHEMATICAL ANALYSIS OF BACTERIAL COUNT DATA : 元論文
- 微生物学の基礎 : 線形法と非線形法、FORモデルについての解説
- Rによる最適化、パラメータ推定入門 - yasuhisa's blog : 最適化についてはここ読んどけばいいと思うよ
