Rで統計解析するための備忘録的なやつです。

それと、前回の延長戦です。
【 R 】多項式回帰

#使用するデータ
前回同様、今回もこのデータを使います。

DNase

DNaseのELISAのデータらしいです（よくわかんないですが）。

#必要なパッケージ
ggplot2を読み込んでおきます。

#ライブラリの読み込み
library(ggplot2)

#データの読み込み
データを読み込みます。

#データの読み込み
DF <- DNase

#ひとまずプロット
ggplot2でプロット。
確認用なのでシンプルに。

#とりあえずプロット
ggplot(DF, aes(x = conc, y = density)) +
  geom_point()

描画の結果がこんな感じ。

濃度依存的に密度が高くなっています。

#多項式回帰分析
Rで多項式回帰分析を行うときは前回の二通りのやり方を紹介しました。
【 R 】多項式回帰

今回も二次式と三次式を回帰式とします。

二次式

y = b_0 + b_1x_1 + b_2x_1^2
\\

三次式

y = b_0 + b_1x_1 + b_2x_1^2 + b_3x_1^3
\\

#回帰分析
LM2 <- lm(density ~ conc + I(conc^2), data = DF)
LM3 <- lm(density ~ conc + I(conc^2) + I(conc^3), data = DF)

**coefficients()**で回帰係数をチェックしてみます。

coefficients(LM2) #LM2モデルの係数を確認する
coefficients(LM3) #LM3モデルの係数を確認する

その結果

> coefficients(LM2) #LM2モデルの係数を確認する
(Intercept)        conc   I(conc^2) 
 0.12197753  0.31807970 -0.01501489 
> coefficients(LM3) #LM3モデルの係数を確認する
(Intercept)        conc   I(conc^2)   I(conc^3) 
 0.05779811  0.46110796 -0.05173350  0.00206439

#プロットに多項式回帰曲線を描画
**geom_smooth()**関数で多項式回帰曲線を加えます。

引数は二次式の場合、
method = lm, formula = y ~ x + I(x^2), se = FALSE

三次式の場合は、
method = lm, formula = y ~ x + I(x^2) + I(x^3), se = FALSE

I(x^i)をpoly(x, i)を使ったやり方にしてもできます。
信頼区間は書きたくないのでse = FALSE としました。

#プロットに回帰曲線を追加する
ggplot(DF, aes(x = conc, y = density)) +
  geom_point() +
  geom_smooth(method = lm, formula = y ~ x + I(x^2), se = FALSE) +
  geom_smooth(method = lm, formula = y ~ x + I(x^2) + I(x^3), se = FALSE)

こんな感じで描画されます。

ちなみにエクセルでこれを描くとこんな感じです。

見た目が味気ないので、色々いじった描画がこちら。

#見た目の味付け
ggplot(DF, aes(x = conc, y = density)) +
  geom_point(color = "grey", size = 2, alpha = 0.7) + #実測値をプロット
  stat_summary(fun.y = "mean", geom = "point", #平均値をプロット
               shape = "diamond", 
               size = 5, color = "blue", alpha = 0.7) +
  stat_summary(fun.ymin = function(x)mean(x) - sd(x), ##標準偏差をエラーバー
               fun.ymax = function(x)mean(x) + sd(x),
               geom = "errorbar",
               size = 0.8, alpha = 0.7, width = 0.2) +
  geom_smooth(method = lm, formula = y ~ x + I(x^2), #二次式
              linetype = "dashed", #線のタイプを破線
              se = FALSE, color = "tomato") + 
  geom_smooth(method = lm, formula = y ~ x + I(x^2) + I(x^3), 
              linetype = "longdash", #線のタイプを破線（長め）
              se = FALSE, color = "seagreen1") +
  xlab("Concentration") +
  ylab("Density") +
  theme_bw()

#最終的なスクリプト

#DNaseのELISAデータに対する多項式回帰

#Rをきれいにする
rm(list = ls())

#ライブラリの読み込み
library(ggplot2)

#データの読み込み
DF <- DNase

#とりあえずプロット
ggplot(DF, aes(x = conc, y = density)) +
  geom_point() 

#回帰分析
LM2 <- lm(density ~ conc + I(conc^2), data = DF)
LM3 <- lm(density ~ conc + I(conc^2) + I(conc^3), data = DF)

coefficients(LM2) #LM2モデルの係数を確認する
coefficients(LM3) #LM3モデルの係数を確認する

#プロットに回帰曲線を追加する
ggplot(DF, aes(x = conc, y = density)) +
  geom_point() +
  geom_smooth(method = lm, formula = y ~ x + I(x^2), se = FALSE) +
  geom_smooth(method = lm, formula = y ~ x + I(x^2) + I(x^3), se = FALSE)

#見た目の味付け
ggplot(DF, aes(x = conc, y = density)) +
  geom_point(color = "grey", size = 2, alpha = 0.7) + #実測値をプロット
  stat_summary(fun.y = "mean", geom = "point", #平均値をプロット
               shape = "diamond", 
               size = 5, color = "blue", alpha = 0.7) +
  stat_summary(fun.ymin = function(x)mean(x) - sd(x), ##標準偏差をエラーバー
               fun.ymax = function(x)mean(x) + sd(x),
               geom = "errorbar",
               size = 0.8, alpha = 0.7, width = 0.2) +
  geom_smooth(method = lm, formula = y ~ x + I(x^2), #二次式
              linetype = "dashed", #線のタイプを破線
              se = FALSE, color = "tomato") + 
  geom_smooth(method = lm, formula = y ~ x + I(x^2) + I(x^3), 
              linetype = "longdash", #線のタイプを破線（長め）
              se = FALSE, color = "seagreen1") +
  xlab("Concentration") +
  ylab("Density") +
  theme_bw()