（cf. Rいろはのトップページ http://nhkuma.blogspot.jp/p/r_5.html）
# 注：作成者は統計の専門家ではありません。内容には十分な注意を払っているつもりですが、限界があることを理解した上で参照してください。感想・苦情・間違いのご指摘、歓迎します。
（注：使用環境によっては、スペースやクオーテーションマーク"などが化けてコードがエラーになるかもしれません）
###### "古典"統計 or 統計モデリング ######
# 古典統計（ネイマン・ピアソン流？）のメリットは強力な客観性ですが、とにかく「有意差：P < 0.05」さえ出ればよいといった"悪用"が常態化しています（この論文などを参照）。
# この有意確率Pは真の確率を示すものではなく、また帰無仮説と対立仮説が１対１でない限りは、帰無仮説が棄却されたからと言って、対立仮説が採択されることを意味しません。
# Rは古典統計に"冷たく"できているように感じます。例えば、安易な使用すると痛い目をみるanova()関数など（Anova()関数を推奨、後述）。
# 予測を目的とする場合には、GLMなど統計モデリングと情報量規準（AICなど）、あるいはマシンラーニングを用いた方がよいでしょう。
# 応答変数の予測に興味があれば統計モデリング、outputに影響を与えている説明要因の方に興味があれば古典統計を選択するのがよいかと考えています。
# 統計モデリングに進むには→（cf. RいろはのGLM編）
# しかし、古典統計の使用が根強い分野では、たとえ間違っていても"〜ANOVA"という名前さえ付いていれば受け入れられるが"〜モデル"とかいう耳慣れない方法は門前払いというのはやむを得ないことだ思います…"長いものには巻かれろ"というのも無難な選択肢だと思います。
# 本項はそういった"現状"のために最小限必要な手法をまとめておきました。
#### 今回はRに内蔵のデータを使用します ###
# データCO2を呼び出す
data(CO2)
############ ANOVAや回帰を使用する上での前提条件を満たししているかチェックする ######
# 前提条件のチェック→メインの解析、の手順がすでに多重比較にあたるという議論があります。前提条件の検定さえやっておけば、という思考停止にならないためにも、最近は図示をして確かめるようにしたほうが良い気がしています。回帰分析の項に例示したようにplot(lm(モデル))のQ-Qプロットで判断したと記述するのがよい方法だと考えています。
############ 正規性の検定（例：Shapiro-Wilk test）######
# 応答変数（= 従属変数、目的変数）が正規性を満たすかをチェック
# もっとも市民権のありそうな kolmogorov-smirnov 検定
# CO2のuptake、をドルマークで繋ぐことで表現できる
ks.test(CO2$uptake, "pnorm", mean=mean(CO2$uptake), sd=sd(CO2$uptake))
# あまりにもめんどくさい…
# この関数をコピペして利用すると、データ名の記述は一回で済ませられますね。
KS.test <- function(data) ks.test(data, "pnorm", mean=mean(data), sd=sd(data))
KS.test(CO2$uptake) # 使用例
# 代替法：Rで推奨している？、論文でも近年使用例が増えつつあるShapiro-Wilk検定
shapiro.test(CO2$uptake)
# Shapiro-Wilk normality test
#
# data: CO2$uptake
# W = 0.941, p-value = 0.0007908
# 正規分布に対し有意差有り→正規性を満たしていないので、変数変換をする
CO2$uptake.L <- log(CO2$uptake + 1) # 対数変換した列をCO2に追加
CO2$uptake.R <- sqrt(CO2$uptake + 0.5) # 平方根変換した列をCO2に追加
# 再度トライする
shapiro.test(CO2$uptake.L)
shapiro.test(CO2$uptake.R)
# 結局、むしろ正規性が悪化したという珍しいケース…
# 「改善を試みたがダメだった」と、次善の策として元のまま行くしかないだろう
# ちなみに、なるべく図示して、どういうデータなのか視覚的にチェックしてみよう：
par(mfrow=c(1,3)) # 複数の図を行数、列数を指定してまとめるコード
plot(CO2$uptake)
plot(CO2$uptake.L)
plot(CO2$uptake.R)
# 片山状に下がっていっている感じ、たしかに正規分布ではない…
############ 等分散性の検定（例：Levene's test）######
# 残差に対しANOVAや回帰をする→残差に差や傾向がないことをチェックという検定。ただし等分散性のチェックは特にやらないことも多いです。そんなに影響はしません。
library(car) # car というパッケージに入ってるので呼び出す
# 注：初めから入っているパッケージではないので、インストールする必要があります
# （cf. 機能パッケージの追加方法）
leveneTest(uptake ~ Plant, data=CO2)
# CO2というデータのuptakeをPlantの群で分割して、という意味
# Levene's Test for Homogeneity of Variance (center = median)
# Df F value Pr(>F)
# group 11 0.5877 0.8328
# 72
# 今度は有意差無し→分散が異なるとは言えなかった→等分散性は否定されなかった（舌噛みそう…）
############ 回帰分析 ######
# 量的変化に対する量的反応は、ANOVA+多重比較ではなく回帰分析
# （xの変数は不等間隔でも等間隔でもOK）
reg1 <- lm(uptake ~ conc, data=CO2)
summary(reg1)
# 切片や回帰係数（Estimateの部分）、
# それらがゼロと異なっているか否かの検定結果、
# R^2値、自由度調整R^2値などを出力
# Call:
# lm(formula = uptake ~ conc, data = CO2)
#
# Residuals:
# Min 1Q Median 3Q Max
# -22.831 -7.729 1.483 7.748 16.394
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 19.500290 1.853080 10.523 < 2e-16 ***
# conc 0.017731 0.003529 5.024 2.91e-06 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 9.514 on 82 degrees of freedom
# Multiple R-squared: 0.2354, Adjusted R-squared: 0.2261
# F-statistic: 25.25 on 1 and 82 DF, p-value: 2.906e-06
library(car) # carは事前にインストールが必要なパッケージ
Anova(reg1) # 検定にかけてみる
# 回帰係数の分散分析結果が得られる
# Anova Table (Type II tests)
#
# Response: uptake
# Sum Sq Df F value Pr(>F)
# conc 2285 1 25.245 2.906e-06 ***
# Residuals 7422 82
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
par(mfrow=c(2,2)) # 4枚のグラフが出てくるが、これを並べて表示するためのコード
plot(reg1) # lmモデルをplotすると回帰診断の結果が得られる。
・左上：残差（Residuals）と当てはめ値（Fitted values）、当てはめ値に対して特に傾向がなければOK。この例では少し山なりになっているのが気になる。
・左下：標準化残差との関係、標準化することで残差の傾向の有無がハッキリする。標準化すると左上での傾向は見られなくなった。たぶんOKだろう。
・右上：Q-Qプロット（quantile versus quantile plot）、残差が正規分布に従っているならば直線的に増加する。中央付近はたいてい直線に乗ってくるが、ズレは上下の端付近に生じる。この例は判断に悩むところだが、多少正規性を満たしていなくても大きな問題ではない。
・右下：影響の強い点の存在を示す図。この例では見えていないが、Cock's distanceの点線の外側にはみ出る点があると外れ値の影響が大きいことを示す。
############ 1 way-ANOVA ######
library(car) # ←上で一度呼び出していれば不要なコードです（以下でも同様）
Anova(lm(uptake ~ Plant, data=CO2)) # 上記の回帰分析と中身が同じ構造！
# Anova Table (Type II tests)
#
# Response: uptake
# Sum Sq Df F value Pr(>F)
# Plant 4862.2 11 6.569 1.842e-07 ***
# Residuals 4844.8 72
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
############ 多重比較…の前に、ANOVA＋多重比較の乱用に対する批判 ######
# 私はポジティブに多重比較を推奨することはできません。そもそも、仮説があれば事前対比をすべきです（でも、逆に事前対比の方が市民権がないような…）。
# ANOVAは質的に異なる群間の違いを解析する検定です（対照vs処理A群、B群…。あるいは特定の生物種間の性質や反応の違いに着目する必要がある場合のような）。
# よく、X軸に段階的に何かの濃度や環境要因を取り、それに伴うYの変化（量的変化に対する量的反応）をANOVA＋多重比較で解析している例を見掛けますが、データの情報量の無駄遣いです。間違いとまでは言わないですが。
# ANOVA＋多重比較ではよほどの差がない限り隣合うペアの間に有意差は出ません。データ分布の山の裾野が触れるか触れないかぐらいかまたはそれ以上の差が必要です。単に全体としての増加減少傾向を見たい場合は、そんな極端な差があるかどうかを検定したいわけではないはずです。
# 増加減少傾向を示したいのであれば回帰分析でやりましょう。こちらは回帰直線の傾きがゼロと有意に異なっているか否かを検定した上で、傾きが増加なのか減少なのかを回帰係数の正負で示すことができます。多重比較で差が全然出ない場合でも、回帰分析では問題なく傾きが有意となるケースが多々あります。
############ post-hoc test（多重比較）######
# 多重比較：同じデータについて、何回も検定を繰り返すとそのうち有意になる確率は、繰り返す度に増大する。例えば、コイントスを１回やって表が出なくとも、もう１回やったら表が出そう、というようなこと。同様に、検定をもう１回行う場合には、p値の判定の厳しさを倍にしなければならない。つまり、0.05/2 で0.025。３グループあるときは、多重比較を３回する必要があるので、0.05/2/2/2 = 0.00625を判定基準としなければならない。というのが基本的な多重比較の考え方（Bonferroni法）。
# 多重比較については、こちらのページが分かりやすいです
http://www.med.osaka-u.ac.jp/pub/kid/clinicaljournalclub1.html
# Bonferroni法はあまりに厳しすぎるというので（検出力が落ちすぎる）、より緩めたのが下記に述べるような方法。
# Holm法（sequential-Bonferroni法）
# ☆手計算でも簡単にできるし、ベースとなる統計手法を選ばないので適用範囲が広い方法
p.adjust(probabilities, "holm") # 検定手法を選ばない、出力されている確率値を用いて計算する関数
# 使い方：probabilitiesの部分に比較したい確率を入れる、c(0.01, 0.023, 0.031)のような確率値のベクトルを適用
# 返り値は、修正後のより厳しくなった確率値。
pairwise.t.test(CO2$uptake, CO2$Plant, method="holm") # 分散分析の多重比較の場合
# 群のペア間での差異の有意確率を対戦行列で出力する（P値は修正後の値で出力される）
# TukeyHSD法
TukeyHSD(aov(uptake ~ Plant, data=CO2))
# こちらはペアで有意確率（p adj）を算出、差異の推定と95%信頼区間も出してくれる
# ANOVAとともに使用（パラメトリックな方法）
# cf. multcompでTukey法
library(multcomp) # multcompは要インストール
summary(confint(glht(lm(uptake ~ Plant, data=CO2), linfct=mcp(Plant="Tukey"))))
# linfct=mcp(変数名="Tukey")、ここでは変数名がPlantである
# linfct=mcp(変数名="Tukey")、ここでは変数名がPlantである
# Dunnett法
# 最初の群（普通は対照群）に対する差異を検定する多重比較
library(multcomp) # multcompは要インストールsummary(confint(glht(lm(uptake ~ Plant, data=CO2), linfct=mcp(Plant="Dunnett"))))
# glhtはかなり高機能な関数の模様、help(glht)で説明を見てみるとよいです
############ 多元ANOVA ######
library(car) # carは要インストール
Anova(lm(uptake ~ Type*Treatment, data=CO2))
# "*"：単効果の項 Type, Treatment、さらに交互作用の項 Type:Treatment、の全ての組み合わせを表せる
# Anova Table (Type II tests)
#
# Response: uptake
# Sum Sq Df F value Pr(>F)
# Type 3365.5 1 52.5086 2.378e-10 ***
# Treatment 988.1 1 15.4164 0.0001817 ***
# Type:Treatment 225.7 1 3.5218 0.0642128 .
# Residuals 5127.6 80
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# このAnova関数は多元の時、標準ではType II 平方和で検定する
# →単効果を優先する平方和であり、交互作用は特殊ケースであるという考え方
# 交互作用も等評価する場合は下記のようにtype="III"と明記するとType III 平方和で検定できる
options(contrasts = c("contr.sum", "contr.sum")) # Type IIIの場合はこのコードも通しておく
Anova(lm(uptake ~ Type*Treatment, data=CO2), type="III")
# Anova Table (Type III tests)
#
# Response: uptake
# Sum Sq Df F value Pr(>F)
# (Intercept) 26217.3 1 409.0389 < 2.2e-16 ***
# Type 924.0 1 14.4165 0.0002839 ***
# Treatment 134.6 1 2.1007 0.1511413
# Type:Treatment 225.7 1 3.5218 0.0642128 .
# Residuals 5127.6 80
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# Type II の時とは異なり、Treatmentが有意でなくなった
# cf. 多元ANOVAは多くとも3元で留めましょう。多すぎると摩訶不思議な交互作用が出ます。CCAやRDAなどの多変量解析が向いていると思います。
# 注：他に、全部小文字の"anova"関数があるが、こちらはType I 平方和に基づく計算をする。多元ANOVAでN 数が不揃いなときには、一般的なType II, Type III 平方和を利用した解析と異なる結果となるので、使用に注意が必要です。
# Type I はRコードで先に入れた要因から順に優先的に有意にする（Y ~ A + B + A:B なら、A, B, A:Bの順に有意になりやすくなる）。なので恣意的な解析結果となってしまう。
# 断りなく"anova"の方を解説している初心者向けテキストもよく見掛けるが注意されたし。そもそも N数が揃っていない解析がANOVAのポリシーに反しているのかもしれませんが…
############ 多元ANOVAの後で多重比較をする場合 ######
# まず、そうすることが本当に必要なのか、研究の目的と仮説を振り返って見直してください。
# 交互作用が有意でなければ、各要因それぞれにつき1wayANOVAと同様に多重比較する
# （今回は交互作用が有意でなかったが、有意だった場合のやり方も示しておく↓）
# まず、両変数をドッキングさせた新変数を作る
CO2$TyTr <- factor(paste(CO2$Type, CO2$Treatment, sep=""))
#（ただの文字列になってしまったので、要因として認識させる："factor"という関数）
TukeyHSD(aov(uptake ~ TyTr, data=CO2)) # 新変数でTukeyHSD法を実行する
################ repeated-measures ANOVA #############
# 同一サンプルに対し時間的に条件を変化させつつ複数回測定し、グループ間の差を検定したい場合に使用されてきた方法。言い換えると、全体が実験操作など同一条件変化を経験した時のグループによる応答の違いを見たい場合。単に実験前と後との比較のような2段階比較ならば対応のあるt-検定で十分、この方法は2段階より多い場合に使用。
Rで用意された関数を知りませんが検索すると自作関数が見つかります。単にサンプルIDを要因のひとつに加えた2way-ANOVAをrepeated-measures ANOVAとして紹介してる論文も見ましたが、本当にそれでいいのかどうか確信が持てません。
# 必要な仮定などややこしいので、私ならGLMMで対処します。
# repeated measures ANOVA周辺については、こちらのテキストが詳しいです（簡単ではないが）
＊追記：repeated-measures ANOVAと下記の対応のあるANOVAはどうやら同物異名らしい。一方で交互作用のないtwo-way ANOVAと同じコードを適用している例もあり、混乱があるようにも思える。いずれにしても反復が多ければ検出力が落ちていくはずで、利用可能な場合にはやはりGLMMを用いるのがよいと思う。
################ 対応のある ANOVA #############
# 同一サンプルに対し複数回測定し、全測定のデータを用いてグループ間の差を検定したい場合。単に同一サンプルに対する複数回観測。
# 一回あたりの反応の誤差や計測の誤差が大きいと思われる場合には有効かもしれません。
# 通常は、ある一時点のデータのみを用いるようにし、「統計解析には○○時間後の値を使用」すると言い切り、普通のANOVAなどでよいと思います。
data(iris) # データirisを呼び出すのに使用（前に実行していれば不要）
iris$id <- rep(1:15,each=10) # サンプル番号を追加（各10回測定）
# 同じSpeciesを繰り返し使用すると見なす。
summary(aov(Sepal.Length ~ Species + Error(id/Species), data=iris))
# 一番下の、Error: Within のところが"見たい"検定結果。
# Error: id
# Df Sum Sq Mean Sq
# Species 1 52.83 52.83
#
# Error: id:Species
# Df Sum Sq Mean Sq
# Species 2 6.679 3.34
#
# Error: Within
# Df Sum Sq Mean Sq F value Pr(>F)
# Species 2 4.60 2.3007 8.705 0.00027 ***
# Residuals 144 38.06 0.2643
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# さらにこのあとに多重比較というのは、TukeyHSD()やglht()はコードが通らなかった。seuquential-Bonferroni法を関数を使わずに手動で実行するしかないかもしれません。
# aovはサンプル数が均等であることを前提としている（Type I 平方和）。不揃いの場合は不可と考えた方がよい。
# やはり、第5部で紹介するGLMMを用いた方がスッキリする。サンプルにid付けをし、glmer(.... + (1|id), ...) のように使用すると、idによるバラツキをノイズとして切り分けることができる。
################ nested ANOVA #############
# メインの要因の効果が、バックグラウンドの条件によっても影響されるかどうかを検定したい場合に使用。
# 多元ANOVAとの違いは、バックグラウンドの要因の単効果による影響は見ない点です。
# つまり、メインの要因の単効果と、メインの要因とバックグラウンドの要因との交互作用、の影響を見るのがnested ANOVAです。
#（以下の例では repeated measures…のところのirisとidを使用します）
model.nA <- lm(Sepal.Length ~ Species/id, data=iris)
# スラッシュを使うのがnestedの指定方法
# これはSpecies + Species:id と記述するのと同じことです…メインの要因であるSpecies単効果と、Speciesとidとの交互作用を見ていることになります。
library(car) # car は事前にインストール要
Anova(model.nA) # 今回は単効果はSpeciesだけなので、anova()によるtype I平方和でも構わない。
# Anova Table (Type II tests)
#
# Response: Sepal.Length
# Sum Sq Df F value Pr(>F)
# Species 63.212 2 119.5812 <2e-16 ***
# Species:id 0.896 3 1.1301 0.339
# Residuals 38.060 144
# これでもしSpecies:id に有意差があった場合「Speciesによる違いはidによっても異なる」ことになる。
# repeated-measures ANOVAは、階層的な誤差を用いてANOVAをする。nested ANOVAは交互作用として階層の要因の影響を評価する。という違いと言えばよいのだろうか…。
# repeated-measures ANOVAは、階層的な誤差を用いてANOVAをする。nested ANOVAは交互作用として階層の要因の影響を評価する。という違いと言えばよいのだろうか…。
# 階層が増えるとか、要因が増えるとか、ちょっと構造が複雑になると、モデル式の書き方がややこしくなるし、計算自体がコケることもあります。また、それ以前に、そんな複雑な計算がまともに出来ているとは思えない。
# そんな時は、やはりGLMM。とくにバックグラウンドの要因が複数あるような場合（時間と場所、または場所が沢山あるなど）にはGLMMでないとまともな計算はできないと思います。また、複数の階層にしたい場合は、WinBUGSなどによるMCMC推定をオススメします。
# nested ANOVAは、一段階まではglhtによる多重比較が"計算可能"のようです。それ以上は手動でsequential-Bonferroni法を使用するしかないかもしれません。
# ちなみに私は、複雑な構造のnested ANOVAを頑張って実行しようとして試行錯誤した結果、どうやってもこんな複雑なデータは検定ベースでは解析できない！と断念して、統計モデリングの世界に入った経緯があります…。
# ちなみに私は、複雑な構造のnested ANOVAを頑張って実行しようとして試行錯誤した結果、どうやってもこんな複雑なデータは検定ベースでは解析できない！と断念して、統計モデリングの世界に入った経緯があります…。
################ cf. 共分散分析 ANCOVA #############
################ cf. 重回帰分析 #############
# 重回帰は、前提条件のチェックなどけっこう厄介なので、何でもかんでも変数を放り込んで"レンジでチン"的な安易な適用は危険だと思います。
################ cf. 多変量分散分析 MANOVA #############
# 1つのサンプルに対し、複数項目の応答（呼吸と光合成など）があり、説明変数に対する反応をトータルとして見る場合（変数間の相関（共分散構造）をも組み込んだ解析をしてくれる）。個別の応答タイプに対し個々の検定をするのではなく、検定を一括で１回で済ますことで、第一種の過誤を抑える狙いがある。また逆に、複数変数を同時に解析することにより、単一変数では検出できない効果を検出できる場合もある。
# manova()関数を使用。
# 前提条件は基本的にはANOVAと同様：
# ・応答変数同士の等分散性、個々の変数の正規性が必要（ただし正規性に対しては頑健）
# ・サンプル数が少ないと検定力が落ちる←サンプル数が少ないからMANOVAで複数の応答変数をトータルで見て"ゆ〜い"にしようというスケベ心には答えてくれない
# 関連（相関）のありそうな複数の指標を用いてグループ間の差を見たいケースなどに有効。全体の傾向を見るのに適している。
# MANOVA以外のやり方としては、主成分分析で応答変数を減らした後、成分ごとに解析を行うのも手だろう。共分散構造を用いて変数をまとめるので、主成分と元の変数との相関関係もちゃんと分かります。
# PerMANOVAもオススメ、library(vegan)のadonis関数で提供されています。生物群集の構造（各種の個体数）を応答変数とした解析が主目的のようですが、method="euclid"を用いれば連続変数にも適用できると思います。
################ cf. 多変量分散分析 MANOVA #############
# 1つのサンプルに対し、複数項目の応答（呼吸と光合成など）があり、説明変数に対する反応をトータルとして見る場合（変数間の相関（共分散構造）をも組み込んだ解析をしてくれる）。個別の応答タイプに対し個々の検定をするのではなく、検定を一括で１回で済ますことで、第一種の過誤を抑える狙いがある。また逆に、複数変数を同時に解析することにより、単一変数では検出できない効果を検出できる場合もある。
# manova()関数を使用。
# 前提条件は基本的にはANOVAと同様：
# ・応答変数同士の等分散性、個々の変数の正規性が必要（ただし正規性に対しては頑健）
# ・サンプル数が少ないと検定力が落ちる←サンプル数が少ないからMANOVAで複数の応答変数をトータルで見て"ゆ〜い"にしようというスケベ心には答えてくれない
# 関連（相関）のありそうな複数の指標を用いてグループ間の差を見たいケースなどに有効。全体の傾向を見るのに適している。
# MANOVA以外のやり方としては、主成分分析で応答変数を減らした後、成分ごとに解析を行うのも手だろう。共分散構造を用いて変数をまとめるので、主成分と元の変数との相関関係もちゃんと分かります。
# PerMANOVAもオススメ、library(vegan)のadonis関数で提供されています。生物群集の構造（各種の個体数）を応答変数とした解析が主目的のようですが、method="euclid"を用いれば連続変数にも適用できると思います。
0 件のコメント:
コメントを投稿