Rと統計解析


このページでは、Rでの統計解析を中心に解説します。


更新履歴:ページを作成(2011/1/29)
2011年1月29日からのアクセス数: -



GLMM

GLMについては知っているものとして解説します。


こんなときはGLMMだ!

  • 影響を検討したい変数以外に、応答変数に明らかに影響を与えている要因がある。
  • その要因に関して、要因として評価はしたくないが、無視もできない。
このような、「興味はないが応答変数に影響を与える変数」のことを、Random effectといい、これを含んだGLMをGLMM(Generalized linear mixed model)と呼びます。私たちが取るほとんどのデータは、実はRandom effectを含んでいます(林学関係でこれがないと断言できる状況なんてないと思うが)。

Random effectの例~世の中Random effectだらけ!
  • 個体差(同じ処理区なのに......)
  • ブロック間差
  • その場所に固有な何か
  • 時間(年)変動
ということで、GLMMはいまや解析の最も基本的な部分を占めつつあります。

RにおけるGLMMの実行関数

  • glmmML(): パッケージglmmMLに収録。library(glmmML)で使える。
長所: stepAIC()など既存の関数が適用しやすい。安定している。
短所: 扱える誤差構造はbinomial、poissonのみ。
     入れられるRandom effectは1個だけ
基本形: #REはRandom effect、dはデータフレーム名
> glmmML(y ~ x1 + ..., cluster=RE, family=**, d)
  • lmer(): パッケージlme4に収録。library(lme4)で使える。
長所: 扱える誤差構造がgaussian、binomial、poissonなど多い。
     Random effectは複数指定可能(しかも切片と傾きどちらも可)。
短所: 関数の仕様が特殊なため、既存の関数が適用できないことが多い。
基本形:
> lmer(y ~ x1 + ... + (1 | RE),
+      family=**, method="Laplace", d)

それぞれの解析例を示します。
> library(glmmML)
> res <- glmmML(surv ~ light, cluster=FL,
+               family=binomial, d)
> summary(res)
coef se(coef) z Pr(>|z|)
(Intercept) -2.369 0.376 -6.300 2.99e-10
light 13.111 2.126 6.166 6.99e-10
Standard deviation in mixing distribution: 1.2e-05
Std. Error: 0.2354
Residual deviance: 256.1 on 221 degrees of freedom
AIC: 262.1
GLMと大きな差はありません。Null devianceが出力されないぐらいが違いでしょうか。

> library(lme4)
> res <- lmer(surv ~ light + (1 | FL),
+             family=binomial, method="Laplace", d)
> res
Generalized linear mixed model fit using Laplace
Formula: surv ~ light + (1 | FL)
Data: d
Family: binomial(logit link)
AIC BIC logLik deviance
262.1 272.4 -128.1 256.1
Random effects:
Groups Name Variance Std.Dev.
FL (Intercept) 5e-10 2.2361e-05
number of obs: 224, groups: FL, 29
Estimated scale (compare to 1 ) 0.9923228
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.3773 0.3767 -6.311 2.77e-10 ***
light 13.1629 2.1299 6.180 6.41e-10 ***
---
...... 一部省略
こちらは大分様子が違います。まずモデルの当てはまりを示すAICなどが示されており、その下にRandom effectがどの程度のばらつきを説明しているかが表示され、その下に説明変数に関する情報が載せられています。結果のオブジェクトをfixef()に入れると説明変数が、ranef()に入れるとRandom effectの値が表示されます。

lmer()の推定結果は要注意

2010年10月13日追記。East_scrofaさんのところで色々と実験がなされています。結果を用いるときは、妥当性をよく確かめて下さい。


各種一般検定関数

こんなときは検定だ!

  • すでに検討したい要因が絞られていて、その要因の影響の差を検討したい。
  • 厳密な制御環境にあり、検討したい要因以外の影響を考慮する必要がない。
  • 指導者(あるいはEditor)に検定しろって言われて逆らえない^^;)(このケースが一番多いか?)

検定の原理に関しては逐一説明しませんが、基本的には
  • 差を見ようとする集団間で差がないと仮定する。
  • 両集団からランダムサンプリングしてきて差をとった値(あるいはもう少しこねくり回した値。いわゆる検定統計量)が既存の理論分布(正規分布またはそれに類するものとすることが多い)に従うと仮定する。
  • 実際の測定値から検定統計量を計算したときに、両集団に差がないと仮定したときの検定統計量の分布においてどこに位置するかを調べ、すごくはしっこ(いわゆる95 %点よりも端。つまり珍しいゾーン)だったら、「両集団に差がないという仮定が間違っていた」と結論する。
ということです。

まずは以下の感じでオブジェクトを作っておく。
> d <- read.csv("data.csv")
> d2 <- split(d, d$syori)

t検定

t.test()を使う。
> d0 <- cbind(d2[[1]], d2[[2]])
> t.test(omosa ~ syori, data=d0, var.equal=TRUE)
var.equalをTRUEにしておけば通常のT検定、FALSEの場合(こちらがデフォルト)は、いわゆるWelchの検定となる。

U検定

wilcox.test()を使う。U検定は数学的にはWilcoxonの順位和検定と呼ばれる検定と同様であり、Rの中では、Wilcoxonの順位和検定が行われる。
> wilcox.test(omosa ~ syori, data=d0)

対応のあるt検定

t.test()の中で、paired = TRUEとすることで実行できる。

ANOVA(F検定)

aov()を使う。ただし、ANOVAは各群の母分散が等しい(等分散性)ことを仮定するので、それを確かめるために、bartlett.test()によって、等分散性の検定を事前に行っている。
> bartlett.test(omosa ~ syori, d)
> result <- aov(omosa ~ syori, d)
> summary(result) #結果の完全な取出しにはsummary()

Kruskal Wallis検定

kruskal.test()を使う。
> result <- kruskal.test(omosa ~ syori, d)

多重比較

パラメトリックな多重比較法として最も一般的に使われるTukey法の場合。
> result <- aov(omosa ~ syori, d)
> TukeyHSD(result)

ノンパラ検定の場合によく使われるStee-Dwass 検定の場合。決まった関数は用意されていないので、群馬大学の青木先生が作ったスクリプトを利用するのがよい。リンクはこちら

これを一度R上にまるごとコピペしてEnterしておくと、以下のような方法でSteel-Dwass検定が行えるようになる。多重比較したいデータのある列を最初に、カテゴリーを示す列を次に持ってくる。カテゴリー変数は数値化しておく必要があるので、as.numeric() で囲っている。
(スクリプトをR上にコピペ)
> Steel.Dwass(d$miki, as.numeric(d$syori))

ANOVA(N 元配置の分散分析)

library(car)に含まれるAnovaとlm()を使う。
library(car)
result <- Anova(lm(miki ~ omosa * syori, d), type="II")
result
式の中の*は単独項と交互作用項を含める記号。:なら交互作用項のみが投入される。

相関

cor.test()を使う。
> cor.test(d$omosa, d$miki)

回帰

lm()を使う。
> result <- lm(miki ~ omosa, d)
> summary(result)

非線形回帰

nls()を使う。用意したデータセットでは非線形回帰するようなシチュエーションがないので適当に以下のようにデータを作る。
> Height <- c(0, 0.3:12.3)
> Age <- c(0, 11.8, 26.6, 32, 35.2, 37.6, 40.2,
+          43, 46, 49.4, 53.4, 59.2, 71.4, 91)
> d <- data.frame(Height, Age)
> plot(Height ~ Age, d)
んで、
> result <- nls(Height ~ a/(1+b*exp(-c*Age)),
+               start=c(a=10, b=100, c=0.5), d)
というようする。start=c()で、初期値を決めている。うまく収束しなかった場合は、何かしらのエラーメッセージが出る。

各種クロス集計表の検定

chisq.test()やfisher.test()を使う。
> d <- matrix(c(100, 20, 70, 50), ncol=2, byrow=T)
> rownames(d) <- c(``Engineer'', ``Agriculture'')
> colnames(d) <- c(``Male'', ``Female'')
> d #データの表示
> chisq.test(d)
2£2 の表なら、fisher の正確確率検定の方がよい。
> fisher.test(d)
最終更新:2012年04月01日 16:29
添付ファイル