/var/log/messages

Aug 25, 2016 - 5 minute read - Comments - misc

R 備忘

でびあんな配布系であれば r-base というパケジを導入すれば良い模様。あとは全画面な rxvt あたりを起動して emacs -nw で良いのかどうか。

とりあえず

を入手して目を通しはじめていたりするのですが、これ良いです。R で云々した後で Python なナニでも、と思っている次第です。備忘なメモを以降で書きちらかす方向にて。

t 検定?

  • 期待値
  • 分散
  • サンプルサイズ

というものが揃ったら t 値というものを求めてやれば有意差があるかどうかが分かる、との事。t 値は以下の式で求めることができる模様。

t 値 = 期待値の差 / 分散的なもの ÷ サンプルサイズ的なもの

これが大きければ良い、のですが大小判定のために以下の手順を踏むとのこと。

  • t 値を p 値に変換
  • t 値が大きければ p 値は小さくなる
  • p 値が 0.05 を下回れば t 値は十分大きい

で、この検定ですが以下なデータセットがあるとして

d <- c(-1, -1, 0, 0, 1, 3, 5, 6, 7)

以下で良いようです。

> > > > t.test(d)

        One Sample t-test

data:  d 
t = 2.5861, df = 9, p-value = 0.0294
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.3382386 5.0617614
sample estimates:
mean of x
      2.7

p-value は十分小さいです。

分散分析モデル

「選択肢を変えることで期待値が変わるか?」を検定したい場合に ANOVA と呼ばれる分析の方法を使う、とのこと。方法としてはデータを以下のように分解し

データ = 平均値 + 効果 + 誤差

そのうえで F 比と呼ばれる値を計算。

F 比 = 効果の大きさ / 誤差の大きさ

この値が大きければ有意差あり、と見做すとのこと。検定の方法としてはこの値から p 値を計算して、という事になるようです。

用語について

  • 自由度 : データ構造を加味したサンプルサイズ。検定によって自由度の求め方は異なるので都度調べること
  • 僅差平方和 : 分散の分子の値、誤差の僅差平方和は残差平方和あるいは RSS (Residual Sum of Squared) と呼ばれる
  • 平均平方 : 分散。F 比は効果と誤差の平均平方の比

分散分析の p 値

計算された F 比と効果の自由度、誤差の自由度を使って p 値は計算する。

> 1 - pf(16, 2, 3)

R で分散分析

以下なデータフレームが使われています。

> d2 <- data.frame(
   length = c(2, 4, 10, 12, 6, 8),
   medicine = c("A", "A", "B", "B", "C", "C")
)

lm という手続きでモデルを作る模様。

> modelANOVA <- lm(length ~ medicine, data=d2)

length が対象で対象の値を変化させている要素が medicine という形で第一引数を記載する模様。興味の対象 (応答変数) が左側で対象を変化させているもの (説明変数) が左側とのこと。

以下な形でも良いとのこと。

> modelANOVA <- lm(d2$length ~ d2$medicine)

あとは anova という関数にモデルを渡せば良い模様。

> anova(modelANOVA)
Analysis of Variance Table

Response: length
          Df Sum Sq Mean Sq F value  Pr(>F)
medicine   2     64      32      16 0.02509 *
Residuals  3      6       2

出力一部略してます。

  • Response: length は応答変数は length です、という意味
  • DF は自由度、Sum Sq は僅差平方和、Mean Sq は平均平方、F 比、p 値
  • medicine の行は効果の情報
  • Residuals の行は誤差の情報

モデルの選択

データ = 平均値 + 効果 + 誤差

という構造でデータを表現しますが、効果という書き方をしている説明変数というものが妥当かどうかは

データ = 平均値 + 誤差2

な誤差2 から誤差を引いた値が効果の大きさとなるとのことで以下の方法で確認ができるとのこと。

> model_1 <- lm(length ~ medicine, data=d2)
> model_2 <- lm(length ~ 1, data=d2)

model_2 は説明変数が固定という考え方で良いのかな。

> anova(model_1, model_2)
Analysis of Variance Table

Model 1: length ~ medicine
Model 2: length ~ 1
  Res.Df RSS Df Sum of Sq  F  Pr(>F)
1      3   6
2      5  70 -2       -64 16 0.02509 *

上の書き方がモデル選択における分散分析の記法とのこと。

  • Res.DF は誤差の分散の自由度
  • RSS は誤差の偏差平方和 (残差平方和)

うーん、上の出力、見方がイマイチよくわからんな。Model 1 を使うことで誤差は少なくなる (効果は大きくなる?) から 1 の方が良いでしょ、という事なのかな。

回帰分析

「ある値が変化することによって、対象がどれだけ変化するか」をモデル化する手法とのこと。

変化の対象? としては

  • 分散分析は選択肢
  • 回帰分析は数値

になる模様。あと線形回帰、ということで

y = ax + b + 誤差

みたいな式も出てきていたり。傾きと切片というパラメータでモデルの推定を、とのこと。そしてパラメータの推定には最小二乗法という手法が出てきました。残差平方和 (予測値と実測値の差を二乗したもの) を最小とするパラメータをモデルに使う手法とのこと。

実際

R で optim という関数を使う模様。残差平方和を計算する関数を定義して optim に渡せば良いらしい、ということでテキストを追い掛けてみます。

以下なデータフレームを作っておきます。

d3 <- data.frame(
   beer = c(4, 3, 5, 5, 8),
   temperature = c(1, 2, 3, 4, 5)
)

R は plot という関数がある模様。

plot (beer ~ temperature, xlim=c(0,6), ylim=c(2,8), data=d3)

で、ax + b な関数作って云々、は抜きにしていきなり lm で何とかなる模様。

modelKaiki <- lm(beer ~ temperature, data=d3)
modelNull <- lm(beer ~ 1, data=d3)

anova(modelKaiki, modelNull)

anova の出力は以下でした。

> > > anova(modelKaiki, modelNull)                                                                                                                                                                                

Analysis of Variance Table

Model 1: beer ~ temperature
Model 2: beer ~ 1
  Res.Df RSS Df Sum of Sq   F  Pr(>F)
1      3   4
2      4  14 -1       -10 7.5 0.07142 .

p 値は 0.07 なので「有意な影響を与えているとは言えない」という結論とのこと。

正規線形モデル

分散分析モデルと回帰モデルを一つにまとめる、とのこと。あと、用語のフォローあり。

  • 分散分析で用いた選択肢型のデータをカテゴリ変数と呼ぶ
  • 回帰分析で用いた数値型のデータを連続変数と呼ぶ

正規線形モデルにおける応答変数は必ず連続変数とのこと。説明変数はどちらでも良いし、複数ある場合は混在しても良いとのこと。また、この章では

  • 応答変数として魚の体長
  • 説明変数として薬の有無、餌の量、という二種類

を使う模様。

解析の流れ

以下とのこと。

  • データの読み込み
  • データの図示 (大切)
  • 統計モデルの作成
  • 検定とモデル選択
  • 選択したモデルの結果を用いた予測

解析

あら? データを読み込んで

d4 <- read.csv("data/data/data0_linearModel.csv")

図示して

pairs(d4)

モデルを作って

lmModel <- lm(length ~ food + medicine, data=d4)

検定して

anova(lmModel)

あるいはモデルに対して summary 適用

summary(lmModel)


Call:
lm(formula = length ~ food + medicine, data = d4)

Residuals:
     Min       1Q   Median       3Q      Max
-25.6371  -7.8463   0.7147   8.6351  20.4050

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 45.66005    2.79274   16.35   <2e-16 ***
food         1.27075    0.04109   30.92   <2e-16 ***
medicineyes 50.93025    2.12848   23.93   <2e-16 ***
---
Signif. codes:  (ry

Residual standard error: 10.58 on 97 degrees of freedom
Multiple R-squared: 0.9346,     Adjusted R-squared: 0.9333
F-statistic: 693.1 on 2 and 97 DF,  p-value: < 2.2e-16 

これ、Intercept が切片で food および medicineyes が係数なのか。で、係数が導きだせれば予測ができると。

一旦終わり

続きも面白そうですが一旦ここで止めます。