Rを使った多重比較

2010年4月9日のブログ記事を貼り付け編集しました。)

GLHT関数は、General Linear Hypothesesのことでmultcompというパッケージの中に入っています。

まずRのGUIの「パッケージ」→「パッケージのインストール」→「ミラーサイトの選択(どこでもいい)」 リストから multcompというパッケージを選び、インストールします。

また「パッケージ」→「パッケージの読み込み」からリストの中にあるパッケージ[multcomp]を読み込みます。

これでglht関数が使えるようになります。

mydata<-read.csv(・・・・

のようにmydataに

要因X(control=Cかmoderately mowed=MMかheavily mowed=HM)と

目的変数Y=種数が入ったデータを読み込んだとします

~~~~~~~~~~~~~~~~~~~~~

のように入れます。

データは完全な架空データです。

まずglmを使ってポアソン回帰します。

結果を

glmmodelというオブジェクト(いれもの)に格納しろ、という命令が下の通りです。

glmmodel<-glm(Y~X,data=mydata,family=poisson)

結果を詳しく見るために

summary(glmmodel)と打ち込むと

> summary(glmmodel)

Call:

glm(formula = Y ~ X, family = poisson, data = mydata)

Deviance Residuals:

Min 1Q Median 3Q Max

    • 1.04868 -0.54871 -0.08828 0.30426 1.86508

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) 1.6487 0.1961 8.407 < 2e-16 ***

XHM -1.1787 0.4043 -2.915 0.00355 **

XMM -0.1226 0.2863 -0.428 0.66843

---

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 19.6610 on 14 degrees of freedom

Residual deviance: 8.4022 on 12 degrees of freedom

AIC: 60.157

Number of Fisher Scoring iterations: 4

こんなメッセージが出てくるはずです。

(dispersion parameterが1以下になっていてデータ不足の様相ですが、今は気にしないことにします)

対照区と重度の刈り取り区HMとの間に有意差がある、と出ています。

負の二項分布を仮定する時はMASSパッケージを

パッケージを読み込むから読み込むか、

library(MASS)

と打ち込むかでglm()の代わりにglm.nbを使いますが今は省略。

ここでglht()関数の出番です。

glht()関数はlm()線形モデル関数や、glm()GLM関数の結果オブジェクトを対象にしています。

multicomparison<-glht(glmmodel,linfct=mcp(X="Dunnet"))

とします。

今対照区VS他の2群を比べているので

Dunnet testで調整しています。

"Dunnet"のところを

"Tukey"とすると

全群を比較します。

summary(multicomparison)

とすると

> summary(multicomparison)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts

Fit: glm(formula = Y ~ X, family = poisson, data = mydata)

Linear Hypotheses:

Estimate Std. Error z value p value

HM - C == 0 -1.1787 0.4043 -2.915 0.00702 **

MM - C == 0 -0.1226 0.2863 -0.428 0.88427

---

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Adjusted p values reported)

多重性を考慮しても

HMと対照区Cの間には有意な種の減少が見られた、といえます。