・お題：データが出たので、Fittingしてみたい。
・とりあえずデータセットを作成する。
>library(tidyverse)
> set.seed(1)
> X = seq(0, 10, length = 11)
> Y = (X - 1)*(X - 5)*(X - 8)
> d1 = rnorm(11, mean = 0, sd = 10)
> Y = Y + d1
> df1 = data.frame(X, Y)
> df1
・プロットする。
> p1 = ggplot(data = df1, mapping = aes(x = X, y = Y))+
+ geom_point(size = 2, color = "gray40")+
+ theme_classic()
> p1
・図示するだけなら、ggplot2のstat_smoothで何とかなりそう。設定は以下を参照。methodに線形モデルを、formulaに1次、2次、3次関数（赤、橙、緑）を、信頼区間（se）に非表示を設定した。
www.rdocumentation.org> p1+
+ stat_smooth(method = "lm", formula = y ~ poly(x, 1), color = "red", se = FALSE) +
+ stat_smooth(method = "lm", formula = y ~ poly(x, 2), color = "orange", se = FALSE) +
+ stat_smooth(method = "lm", formula = y ~ poly(x, 3), color = "green", se = FALSE)
・ちなみに、se = TRUEで信頼区間を表示するのがデフォルト。信頼区間は95%がデフォルトだけれど設定することもできる。
> p1+
+ stat_smooth(method = "lm", formula = y ~ poly(x, 3), color = "green")
・次に、もうちょっとパラメータを求めたい。lm関数を使ってみる。ちなみに、lmは線形モデルLinearModelの略っぽい。一般化線形モデル版でglm関数などもある。使い方は以下を参照。
> lm1 = lm(formula = Y ~ poly(X, 1, raw = TRUE), data = df1)
> lm2 = lm(formula = Y ~ poly(X, 2, raw = TRUE), data = df1)
> lm3 = lm(formula = Y ~ poly(X, 3, raw = TRUE), data = df1)
・図示するとこんな感じ。
> p1+
+ geom_line(mapping = aes(y = lm1$fitted.values), color = "red", linewidth = 1)+
+ geom_line(mapping = aes(y = lm2$fitted.values), color = "orange", linewidth = 1)+
+ geom_line(mapping = aes(y = lm3$fitted.values), color = "green", linewidth = 1)
・情報を見てみる。元の式を展開するとx**3-14x**2+53x-40になるので、ばらつきも考えると妥当な感じ。
> summary(lm3)
・95%信頼区間も出してみる。
> predict(lm3, interval = "confidence", level = 0.95)
・値を代入して予測したい場合は、列名にlmを作成した際のXを列名としたdataframeをpredict関数に与えるっぽい。
> new = data.frame(X = seq(0, 10, length = 100))
> new %>% head()
> fit3 = predict(lm3, new, se.fit = TRUE)
> df_fit3 = data.frame(X = new$X, Y = fit3$fit)
> df_fit3 %>% head()
・せっかくなので、図示してみる。滑らかになった。
・想定している関数をもう少し細かく設定してみる。nls関数を使う。Nonlinear Least Squaresの略っぽい。nls関数に関しては以下を参照。
・とりあえずデータセットを作る。どこかで見たような感じ。。
> Substance = rep(2**(0:10), each = 3)
> velocity = 30 * Substance / (15 + Substance) + rnorm(33)
> df2 = data.frame(S = Substance,
+ v = velocity)
> df2 %>% head()
> p2 = ggplot(data = df2,
+ mapping = aes(x = S, y = v))+
+ geom_point()
> p2
・次に、fittingしてみる。formulaのところで想定している関数を設定できる。
> nls1 = nls(data = df2,
+ formula = v ~ Vmax*S / (Km + S),
+ start = c(Vmax = 1, Km = 1),
+ trace = TRUE)
> summary(nls1)
・Vmaxが30、Kmが15を想定してデータを作成し、ばらつきを加えていたので、大体妥当な値になった気がする。
・fittingがどんな感じだったか図示してみる。
> df_pred = data.frame(S = seq(1, 1024, length = 1000))
> df_pred$v = predict(nls1, newdata = pred)
> df_pred %>% head()
> p2 +
+ geom_line(data = df_pred,
+ mapping = aes(x = S, y = v),
+ color = "red")+
+ theme_classic()
おわり
