タグ:分散分析

第4章では統計を扱っています。
その中で何度か分散分析のを紹介してます。

【4-10】Rで分散分析(一元配置分散分析)を行う方法
【4-12】Rで1要因反復測定分散分析(repeated-measures-ANOVA)を行う方法①
【4-12】Rで1要因反復測定分散分析(repeated-measures-ANOVA)を行う方法②

Rで分散分析を行う関数はいくつかあり、今までaov関数、lm関数→anova関数、Anova関数など紹介してきました。
ただRにはANOVA君という分散分析に特化したスクリプトがあります。

ANOVA君を使えばどんな形式の分散分析でも行え、多重比較も一気にできます。

今回はANOVA君を紹介します。



1.データの準備

今回は前回のデータを使います。記事の1~4までをご準備ください。


もし前回の記事を見るなんて面倒くさい!という方は以下のコードを実行してください。
(このコードはデータの保存場所が今までのの記事と違います)

#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-repeated-measures-anova.xlsx"
destfile <- "demo-repeated-measures-anova.xlsx" #今まで第4章の記事を読まれて実践されてる方は destfile <- "data/demo-repeated-measures-anova.xlsx"
download.file(url, destfile) #データの読み込み library(readxl) grip_rep_anova <- read_excel("data/demo-repeated-measures-anova.xlsx", range = "B2:E32") View(grip_rep_anova)

スクリーンショット 2019-11-25 22.48.23
前回は
ある運動プログラムを行い、0週目・1週目・2週目で握力を測定した、ということを想定した仮想データとなっています。更に男女の群分けも行います。

2.被験者間要因と被験者内要因と分散分析のデザインについて
分散分析を行う上で大切なのが被験者間要因と被験者内要因についてです。

スクリーンショット 2019-11-24 6.33.40
被験者間要因はいわゆる対応のないデータのことで被験者全体をA法、B法の2つに割り振るといった方法です。なので同じ人がA法とB法を行うことはありません。
被験者内要因はいわゆる対応のあるデータのことで同じ人が繰り返し測定します

そして被験者間要因は列の左側に、被験者内要因は列の右側に並べると後々分析しやすくなります。
また被験者間要因は縦に、被験者内要因は横につなげると分析を行うには都合がいいです。

仮に下のデータがあったとします。
スクリーンショット 2019-11-24 6.35.12
この並びが分散分析を行う上での基本的な形になります。

要因が3つ(sex, treatment, 繰り返し)あるのでABC
被験者間要因と被験者内要因の間にsを付ける

上記の場合はABsCデザインと言います。
また被験者間要因と被験者内要因が両方あるデザインを混合計画とも呼んだりします。

ANOVA君を使うためにはこのデザインの名付け方とデータの並べ方が鍵になります。

3.ANOVA君の導入

ANOVA君は通常のインストールでは行なえません。

ANOVA君の作者である井関龍太さんのホームページでダウンロードする必要があります。



上記記事でtxtファイルを保存します。
保存場所は作業フォルダに入れておくとすぐに使えますが、どこでも大丈夫です。
作業フォルダはgetwd関数を使えばわかります。()内は何も入れません。

getwd()


RstudioでCode→Source Fileでダウンロードしたファイルを指定すればOKです。
スクリーンショット 2019-11-28 7.54.56


4.ANOVA君を使うためのデータの準備

ANOVA君を使うには次の情報が大切です。
①何デザインか?
②デザインに併せたデータの並び順にする

①何デザインか?

今回は以下の列があります。

被験者間要因:sex(A
被験者内要因:W0,W1,W2の繰り返し(B


【4-12】Rで1要因反復測定分散分析(repeated-measures-ANOVA)を行う方法①ではBのみで分散分析を行いました。この場合は被験者間要因がありませんのでsAデザインとなります。

【4-12】Rで1要因反復測定分散分析(repeated-measures-ANOVA)を行う方法②では被験者間要因がA被験者内要因がBなのでAsBデザインとなります。

分散分析には一元配置とか二元配置とか反復測定とか混合モデルとかいろいろあるのですが、なれない方はまず上記のデザインが付けられることを優先してもいいと思います。その後書籍で二元配置とか出てきた時にそのデータをみながら「あぁ、これは○○デザインのことだな」と理解できればいいと思います。

②デザインに併せたデータの並び順にする
そしてANOVA君はデザインの順にデータを並べることが必須です。

これはtidyverseパッケージのselect関数を使えばできます。
スクリーンショット 2019-11-25 22.48.23

今はW0,W1,W2,sexの順に並んでいます。

もしsAデザインならsexを消し、W0,W1,W2のみの列にします。
select関数は【2-4】Rで指定した列や行だけを取り出すselect関数、slice関数、filter関数を紹介しますで紹介しています。

今回はgrip_anovakun_1_wideという変数名を付けました。なんでも大丈夫です。
まずtidyverseパッケージを読み込んでから行います。

library(tidyverse)

grip_anovakun_1_wide <- grip_rep_anova %>% select(W0, W1, W2)

スクリーンショット 2019-12-02 23.30.22

もしAsBデザインならsex,W0,W1,W2の並び順にする必要があります。

grip_anovakun_2_wide <-
  grip_rep_anova %>% 
  select(sex, W0, W1, W2)

スクリーンショット 2019-12-02 23.32.04

これでデータができました。

5.ANOVA君を使う

ANOVA君で分散分析を行うにはanovakun関数を使います。

anovakun(データ, "デザイン", 各要素のグループの数)

ポイントはAやBの要素の数を入れることです。
スクリーンショット 2019-12-02 23.49.02
anovakun(grip_anovakun_1_wide, "sA", 3)
anovakun(grip_anovakun_2_wide, "AsB", 2, 3) anovakun(grip_anovakun_2_wide, "AsB", sex = c("F", "M"), time = c("W0", "W1", "W2"))

6.結果の見方
今回の結果をそれぞれ表示します。

①sAデザインの場合
anovakun(grip_anovakun_1_wide, "sA", 3)

スクリーンショット 2019-12-02 23.56.16

sAデザインと書いています。
AがW0,W1,W2だったので、a1,a2,a3と表記されます。
まず、それぞれの平均と標準偏差が書いてあります。

スクリーンショット 2019-12-03 0.03.43
つぎに反復測定分散分析に必要な球面性の検定を行っています。
球面性の検定に関しては【4-12】Rで1要因反復測定分散分析(repeated-measures-ANOVA)を行う方法①をご参照ください。

スクリーンショット 2019-12-03 0.03.48
そして分散分析表です。
sAデザインは要素が1つしかないので交互作用はありません。

スクリーンショット 2019-12-03 0.03.54
有意差があれば多重比較を行います。
多重比較をは何も指定しないとShafferの方法になります。
もしHolmの方法を行いたい!という場合は holm = TRUEを付け加えます。
anovakun(grip_anovakun_1_wide, "sA", 3, holm = TRUE)

②AsBデザインの場合
anovakun(grip_anovakun_2_wide, "AsB", sex = c("F", "M"), time = c("W0", "W1", "W2"))
スクリーンショット 2019-12-03 0.31.54
最初に各群の平均・標準偏差が書かれています。
スクリーンショット 2019-12-03 0.32.01
反復測定分散分析にひつような球面仮定性についてみています。


スクリーンショット 2019-12-03 0.33.01
次に分散分析表がでています。
今回みたいのはsex x timeの部分です。
timeで有意差がでているので時間経過の影響はあるようですが、その傾向は男女で違うようです。

スクリーンショット 2019-12-03 0.33.08
timeで有意差が出ていたので男女合わせたtimeで多重比較の結果が出ています。

スクリーンショット 2019-12-03 0.33.19
次はsex x timeの交互作用についてです。
今回みたいのは下のtime at Ftime at Mです。
女性・男性それぞれ時系列に差があるようです。

sex at W0 ~ W2はW0においてFとMに差があるか?という意味ですが今回知りたいことではありません。

スクリーンショット 2019-12-03 0.33.29
男女ともに時系列で有意差が見られたのでそれぞれの多重比較の結果が出ています。

どちらもW0 < W1 < W2で大きいようです。


7.まとめ

今回はANOVA君について紹介しました。

Rで分散分析を行うならANOVA君はあまりにも有名なので是非試してみてください!




第4章は統計を扱います。


今回も「シロート統計学」のハルさんとコラボレーションすることとなりました。


シロート統計学はEZRを使った統計分析をわかりやすく解説されています。


 


第4章はシロート統計学で使われていたEZRをRで行うとどうなるのか?といった視点で進めていきます。


今回使うデータもハルさんのサイトと同じものを使わせでいただく事になりました。それぞれ見比べることで参考にしてみてください!


まずはハルさんのサイトをご覧ください。
ここで行っている解析をRで実施する方法を紹介します。


1.データの準備

今回は前回のデータを使います。記事の1~4までをご準備ください。


もし前回の記事を見るなんて面倒くさい!という方は以下のコードを実行してください。
(このコードはデータの保存場所が今までのの記事と違います)

#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-repeated-measures-anova.xlsx"
destfile <- "demo-repeated-measures-anova.xlsx" #今まで第4章の記事を読まれて実践されてる方は destfile <- "data/demo-repeated-measures-anova.xlsx"
download.file(url, destfile) #データの読み込み library(readxl) grip_rep_anova <- read_excel("data/demo-repeated-measures-anova.xlsx", range = "B2:E32") View(grip_rep_anova)

スクリーンショット 2019-11-25 22.48.23
前回は
ある運動プログラムを行い、0週目・1週目・2週目で握力を測定した、ということを想定した仮想データとなっています。今回は握力の経時的検討に加えて男女の群分けも行います。

2.被験者間要因と被験者内要因と分散分析のデザインについて
分散分析を行う上で大切なのがvについてです。
スクリーンショット 2019-11-24 6.33.40
被験者間要因はいわゆる対応のないデータのことで被験者全体をA法、B法の2つに割り振るといった方法です。なので同じ人がA法とB法を行うことはありません。
被験者内要因はいわゆる対応のあるデータのことで同じ人が繰り返し測定します

そして被験者間要因は列の左側に、被験者内要因は列の右側に並べると後々分析しやすくなります。
また被験者間要因は縦に、被験者内要因は横につなげるとEZRで分析を行うには都合がいいです。

仮に下のデータがあったとします。
スクリーンショット 2019-11-24 6.35.12
この並びが分散分析を行う上での基本的な形になります。

要因が3つ(sex, treatment, 繰り返し)あるのでABC
被験者間要因と被験者内要因の間にsを付ける

上記の場合はABsCデザインと言います。
また被験者間要因と被験者内要因が両方あるデザインを混合計画とも呼んだりします。


今回は以下のようになります。

被験者間要因:sex(A
被験者内要因:W0,W1,W2の繰り返し(B


今回のExcelではsexの列が1番右にありますが被験者間要因を左、被験者内要因を右に書くのが習わしのようなのでAsBデザインということになります。

被験者内要因(改善の傾向)に性別が影響しているかどうかを見たいということになります。

3.反復測定分散分析(repeated measures ANOVA)を行う

前回の記事で正規分布の確認を行いましたので、早速反復測定分散分析を行います。

反復測定分散分析を行うには主に2つの方法があります。

1.carパッケージのAnova関数を使う(EZRもcarパッケージを使ってる)
2.ANOVA君を使う

ANOVA君は井関龍太さんが作成した分散分析に特化したスクリプトです。
使い方に特徴がありますが使えるようになると慣れると簡単でどの方法よりも詳細に分析できます。

ここではまず1の方法を紹介します。

反復測定分散分析を行うには次の方法になります。




①carパッケージを読み込む
ここは前回の記事と全く同じです。
まずはcarパッケージを読み込みます。
carパッケージはEZRをインストールする時に一緒にインストールされています。
もしインストールしていなければ先にインストールします。

install.packages("car")


その後library関数で呼び出します。

library(car)


②被験者要因を指定するdata.frameを準備する
ここは前回の内容と全く同じです。
今回のデータではsexの列は「MとFの2つのデータがある」とパソコンは読み取れますがW0,W1,W2はwideデータとなっているため「この3列が繰り返しのデータである」ことがパソコンは読み取れません。それがわかるためのものを2つ作ります。

まず繰り返しの順番がわかるtimeというベクトルを作ります。
このデータはfactor型である必要があるためfactor関数を使います。
そしてtimeのデータをdata.frame型に変更したidataを作ります。

 スクリーンショット 2019-11-24 10.11.33
そして大事なのがlevels = です。levels = を指定しないとfacorは五十音順になります。
今回は問題ないのですがもしpre,mid,postのように列名を付けるとlevels = を使わないとmid→post→preの順番になってしまいます。
factor関数に関しては【1-10】Rでよく使われる型について説明します。
 

③lm関数でモデルを作る
スクリーンショット 2019-11-25 23.38.40
次に分析するためのモデルを作ります。ここが前回の記事と唯一違うところです。

モデルはlm関数を使います。model(名前は何でもいい)という変数名に結果を入れます。

被験者内要因は3列になっているのでcbind関数でつなげます。
被験者間要因は~の右に入れます。もし被験者間要因が複数あれば * でつなぎます(交互作用を見たい場合)。もし交互作用を見ない場合は + でつなぎます。

そして必ず大切なのがcontrasts = の箇所です。
数学的には反復測定分散分析はタイプⅢ平方和というものを計算します。
詳しい話は井関龍太さんの記事にありますが、contrasts = を設定しないと正しいタイプⅢ平方和を計算してくれないため、この設定が必要となります。
使い方はcontrasts = list(〇 = contr.sum, △ = contr.sum, □ = contr.sum)のようにlist内で指定します。


④Anova関数で分散分析を実行する
スクリーンショット 2019-11-25 23.40.16
ここは前回の記事と全く同じです。
反復測定分散分析を行う時はcarパッケージのAnova関数を使います。Aは必ず大文字です。
res(resultの略、何でも良い)という変数名に結果を入れます。

Anova(lm関数で作ったモデル, idata = ◯◯, idesign = ~ ◯◯, type = "Ⅲ")

lm関数で作ったモデル:先程作ったmodel
idata:先程つくったidata
idesign:先程作ったtime(timeの前に を付けます)
type:"Ⅲ"または3とする(反復測定分散分析の場合は必ずⅢ)

スクリーンショット 2019-11-24 12.08.02

⑤summary関数で結果を表示する

summary(Anova関数の結果, multivariate = FALSE)

multivariate = FALSEと入れるとEZRと同じ結果が表示されます。

①〜⑤をまとめると以下になります

library(car)
time <- factor(c("W0", "W1", "W2"), levels = c("W0","W1", "W2"))
idata <- as.data.frame(time)
model <- lm(cbind(W0,W1,W2) ~ sex, 
            data = grip_rep_anova, 
            contrasts = list(sex = contr.sum))
res <- Anova(model, idata = idata, idesign = ~time, type = 3)
summary(res, multivariate = FALSE)


行が長くて大変ですが実は色の部分以外はコピペ可能です。
(被験者間要因があればまた変わります。次回記事で紹介します)
スクリーンショット 2019-11-25 23.51.29

結果はまずMauchlyの球面仮定を確認し、有意差がなければ上、あれば下をチェックします。

sex→握力は性別により変化する
time→握力は時間経過により変化する
sex:time→交互作用(時間経過に伴う握力の改善の程度は男女で異なるか)

スクリーンショット 2019-11-26 0.16.54
結果はEZRと同じになりました。

4.まとめ
今回は2要因以上ある反復測定分散分析を行いました。
前回の記事との違いは被験者間要因があるかどうかでした。
被験者間要因があるとcontrasts = 設定が必要になりますので注意が必要です。

そして前回・今回はEZRと同じ方法のcarパッケージのAnova関数を使いました。
ただどうしてもコードが長くなってしまいます。
ただRにはのANOVA君という分散分析に特化した非常に有名なスクリプトがあります。


次回はANOVA君を使った分散分析を紹介します。

5.今回使ったコード


#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-repeated-measures-anova.xlsx"
destfile <- "demo-repeated-measures-anova.xlsx" #今まで第4章の記事を読まれて実践されてる方は destfile <- "data/demo-repeated-measures-anova.xlsx"
download.file(url, destfile) #データの読み込み library(readxl) grip_rep_anova <- read_excel("data/demo-repeated-measures-anova.xlsx", range = "B2:E32") View(grip_rep_anova)

#sexの列もいれた反復測定分散分析 library(car) time <- factor(c("W0", "W1", "W2"), levels = c("W0","W1", "W2")) idata <- as.data.frame(time) model <- lm(cbind(W0,W1,W2) ~ sex, data = grip_rep_anova, contrasts = list(sex = contr.sum)) res <- Anova(model, idata = idata, idesign = ~time, type = 3) summary(res, multivariate = FALSE)





(2019年11月24日修正:8の②のfactor関数を修正しています)

第4章は統計を扱います。


今回も「シロート統計学」のハルさんとコラボレーションすることとなりました。


シロート統計学はEZRを使った統計分析をわかりやすく解説されています。


 


第4章はシロート統計学で使われていたEZRをRで行うとどうなるのか?といった視点で進めていきます。


今回使うデータもハルさんのサイトと同じものを使わせでいただく事になりました。それぞれ見比べることで参考にしてみてください!


今回は繰り返しのある一元配置分散分析を紹介します



まず繰り返しのある一元配置分散分析(ANOVA)についてはハルさんのサイトをご参照ください。



また1.準備〜4.データの読み込みまでは【4-1】Rでt検定を行う方法と全く同じ流れになります。
もし1〜4まででわからない部分があれば確認してみてください。

 

1.準備

第4章は毎回ExcelデータをダウンロードしてRを使うのでプロジェクトで管理して行うことを勧めています。

 

ここではR練習というプロジェクトを作り、Excelファイルを入れるためのdataフォルダを作っています。
これを前提に次から進めていきます。
スクリーンショット 2019-10-20 7.54.14


2.スクリプトファイルの作成

次にRのコードを書くためのスクリプトファイルを作ります。

スクリーンショット 2019-11-23 22.05.34



3.データのダウンロード


【1-13】Rで読み込みやすいExcelファイルの作り方今回もハルさんのサイトのデータを使わせていただきます。
実は分散分析を行う上で覚えておくと便利なExcelデータの作り方(列の並べ方)があります。
詳細は【1-13】Rで読み込みやすいExcelファイルの作り方で紹介しています。

デモデータ(反復測定分散分析)

この章ではRを使ってダウンロードしています。


download.file(url = “ファイルのURL”,
        destfile = “保存したい場所/ファイル名”)
urlはデモデータで右クリック → リンクのアドレスをコピー
destfileは保存場所と保存のファイル名を指定します。


実際のコードは以下になります。
前回のコードのURL(" "の中)とdestfileのdata/以降を変更するだけでOKです。
#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-repeated-measures-anova.xlsx"
destfile = "data/demo-repeated-measures-anova.xlsx"
download.file(url, destfile)
スクリーンショット 2019-11-23 22.14.17


dataフォルダにダウンロードできたことを確認します。


4.データの読み込み

データを読み込みます。
今回は【4-0】第4章を進めていく上での準備で行った方法で進めます。

View Fileでデータを確認します。

スクリーンショット 2019-11-23 22.15.57

データが入っているセルを確認します。
B2からC22までデータが入っています(B2:E32と表記)
スクリーンショット 2019-11-23 22.20.35


次にImport Datasetでデータを取り込みます。
スクリーンショット 2019-11-23 22.15.57


Import画面ではName, Sheet,Rangeを指定します。

Name:ハルさんのサイトではgripでしたがgrip_rep_anovaとします(大文字・小文字は別物とされます)
Sheet:このExcelは1つしかデータがないのでDefaultのままでOK
Range:先ほど確認したB2:E32(実は今回のように表以外に余計なものがなければ空欄でもOK)

スクリーンショット 2019-11-23 22.24.40

Importボタンを押す前に右にあるコードをコピーしスクリプトファイルに貼り付けることも忘れずに行います。
library(readxl)
grip_rep_anova <- read_excel("data/demo-repeated-measures-anova.xlsx", 
                             range = "B2:E32")
View(grip_rep_anova)

データが正しく入っていることを確認します。
スクリーンショット 2019-11-23 22.25.10


これでデータの取り込みは完成です。

5.被験者間要因と被験者内要因と分散分析のデザインについて
分散分析を行う上で大切なのが被験者間要因と被験者内要因についてです。
スクリーンショット 2019-11-24 6.33.40
被験者間要因はいわゆる対応のないデータのことで被験者全体をA法、B法の2つに割り振るといった方法です。なので同じ人がA法とB法を行うことはありません。
被験者内要因はいわゆる対応のあるデータのことで同じ人が繰り返し測定します

そして被験者間要因は列の左側に、被験者内要因は列の右側に並べると後々分析しやすくなります。
また被験者間要因は縦に、被験者内要因は横につなげるとEZRで分析を行うには都合がいいです。

仮に下のデータがあったとします。
スクリーンショット 2019-11-24 6.35.12
この並びが分散分析を行う上での基本的な形になります。
要因が3つ(sex, treatment, 繰り返し)あるのでABC
被験者間要因と被験者内要因の間にsを付ける
上記の場合はABsCデザインと言います。

今回はsexの列を使いませんので繰り返しの1要因のみで分散分析を実施します。
(sexの列を使った分析は次の記事で紹介します)

被験者間要因がありませんので、今回はsAデザインということになります。
繰り返しの群間に差があるかどうかを見たいということです。


6.正規分布かどうかを確認
分散分析を行うにはデータの分布が正規分布であることが必要です。
ハルさんのサイトに沿って正規分布の確認を行います。

正規分布の確認は【4-10】Rで分散分析(一元配置分散分析)を行う方法と基本的には同じですが、今回は繰り返しの部分がwideデータになっています。longデータであれば今までの記事で紹介した方法でまとめて出すことができますが、今回は1つずつ確認します。

今回ヒストグラムはgeom_histgram関数を使いますのでtidyverseパッケージを読み込みます。
EZRで正規性の確認をする時に併せて作成されるヒストグラムは棒が5本なのでbins = 5とすることにします。

library(tidyverse)
ggplot(data = grip_rep_anova) +
  geom_histogram(aes(x = W0), bins = 5) 

ggplot(data = grip_rep_anova) +
  geom_histogram(aes(x = W1), bins = 5) 

ggplot(data = grip_rep_anova) +
  geom_histogram(aes(x = W2), bins = 5) 
スクリーンショット 2019-11-23 22.42.57


7.シャピロ・ウィルク検定
次にハルさんのサイトでは正規性の検定を行っています。
正規性の検討はshapriro.test関数3つのカテゴリーそれぞれに行う必要があります。
ここでもwideデータなので1つずつ行います。

shapiro.test((grip_rep_anova$W0))
shapiro.test((grip_rep_anova$W1))
shapiro.test((grip_rep_anova$W2))
スクリーンショット 2019-11-23 22.44.56


すべて正規分布であることが棄却されませんでした。


8.反復測定分散分析(repeated measures ANOVA)を行う

そして反復測定分散分析を行います。

反復測定分散分析を行うには主に2つの方法があります。

1.carパッケージのAnova関数を使う(EZRもcarパッケージを使ってる)
2.ANOVA君を使う

ANOVA君は井関龍太さんが作成した分散分析に特化したスクリプトです。
使い方に特徴がありますが使えるようになると慣れると簡単でどの方法よりも詳細に分析できます。

ここではまず1の方法を紹介します。

反復測定分散分析を行うには次の方法になります。




①carパッケージを読み込む

まずはcarパッケージを読み込みます。
carパッケージはEZRをインストールする時に一緒にインストールされています。
もしインストールしていなければ先にインストールします。

install.packages("car")

その後library関数で呼び出します。
library(car)

②被験者要因を指定するdata.frameを準備する

今回のデータではsexの列は「MとFの2つのデータがある」とパソコンは読み取れますがW0,W1,W2はwideデータとなっているため「この3列が繰り返しのデータである」ことがパソコンは読み取れません。それがわかるためのものを2つ作ります。

まず繰り返しの順番がわかるtimeというベクトルを作ります。
このデータはfactor型である必要があるためfactor関数を使います。
そしてtimeのデータをdata.frame型に変更したidataを作ります。

 スクリーンショット 2019-11-24 10.11.33
そして大事なのがlevels = です。levels = を指定しないとfacorは五十音順になります。
今回は問題ないのですがもしpre,mid,postのように列名を付けるとlevels = を使わないとmid→post→preの順番になってしまいます。
factor関数に関しては【1-10】Rでよく使われる型について説明します。
 

③lm関数でモデルを作る

次に分析するためのモデルを作ります。
モデルはlm関数を使います。model(名前は何でもいい)という変数名に結果を入れます。
被験者内要因は3列になっているのでcbind関数でつなげます。
被験者間要因は今回無いので1とします。もしあると別の設定も必要ですが次の記事で紹介します。

スクリーンショット 2019-11-24 11.55.22
④Anova関数で分散分析を実行する



反復測定分散分析を行う時はcarパッケージのAnova関数を使います。Aは必ず大文字です。
res(resultの略、何でも良い)という変数名に結果を入れます。

Anova(lm関数で作ったモデル, idata = ◯◯, idesign = ~ ◯◯, type = "Ⅲ")

lm関数で作ったモデル:先程作ったmodel
idata:先程つくったidata
idesign:先程作ったtime(timeの前に を付けます)
type:"Ⅲ"または3とする(反復測定分散分析の場合は必ずⅢ)
スクリーンショット 2019-11-24 12.08.02

⑤summary関数で結果を表示する

summary(Anova関数の結果, multivariate = FALSE)

multivariate = FALSEと入れるとEZRと同じ結果が表示されます。

①〜⑤をまとめると以下になります
library(car)
time<- factor(c("W0", "W1", "W2"), levels = c("W0","W1", "W2")) idata <- as.data.frame(time) model <- lm(cbind(W0,W1,W2) ~ 1, data = grip_rep_anova) res <- Anova(model, idata = idata, idesign = ~time, type = 3 ) summary(res, multivariate = FALSE)
スクリーンショット 2019-11-24 21.02.57
行が長くて大変ですが実は色の部分以外はコピペ可能です。
(被験者間要因があればまた変わります。次回記事で紹介します)

スクリーンショット 2019-11-24 12.51.33
このあたりはハルさんのサイトと同じです。
EZRと同じ結果になります。

9.多重比較を行う

EZRの多重比較はEZR独自の関数を使っています。
そのため今回はpairwise.t.test関数を使います。
(被験者間要因があると使えないので注意!)
ただpairwise.t.test関数はlongデータで行う必要があります。

wideデータをlongデータに変えるにはgather関数、またはpivot_longer関数を使います。
wideデータとlongデータの切り替えは【2-5】Rでデータを集計するのに便利なlongデータとgather関数で紹介しています。

pairwise.t.test(数値, グループ, p.adjust.method = "◯◯")
〇〇はここでは3つの方法を紹介します。

bonferroni
holm
hochberg


grip_rep_anova_long <- 
  grip_rep_anova %>% 
  gather(key = "key", value = "value", W0:W2)

pairwise.t.test(grip_rep_anova_long$value, grip_rep_anova_long$key, paired = TRUE, p.adjust.method = "bonferroni")
pairwise.t.test(grip_rep_anova_long$value, grip_rep_anova_long$key, paired = TRUE, p.adjust.method = "holm")
pairwise.t.test(grip_rep_anova_long$value, grip_rep_anova_long$key, paired = TRUE, p.adjust.method = "hochberg")

スクリーンショット 2019-11-24 13.17.18


10.まとめ

今回は反復測定分散分析を紹介しました。
今までの手法に比べると手続きが必要ですが、慣れるとコピペでもできるようになります。
次回は被験者間要因も含めた混合計画反復測定分散分析を紹介します。


11.今回使ったコード

今回使ったコードは以下になります。
#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-repeated-measures-anova.xlsx"
destfile <- "data/demo-repeated-measures-anova.xlsx"

download.file(url, destfile)

#データの読み込み
library(readxl)
grip_rep_anova <- read_excel("data/demo-repeated-measures-anova.xlsx", 
                             range = "B2:E32")
View(grip_rep_anova)

#正規分布の確認
library(tidyverse)
ggplot(data = grip_rep_anova) +
  geom_histogram(aes(x = W0), bins = 5) 

ggplot(data = grip_rep_anova) +
  geom_histogram(aes(x = W1), bins = 5) 

ggplot(data = grip_rep_anova) +
  geom_histogram(aes(x = W2), bins = 5) 

shapiro.test((grip_rep_anova$W0))
shapiro.test((grip_rep_anova$W1))
shapiro.test((grip_rep_anova$W2))


#carパッケージのインストールしたことなければインストール
install.packages("car")

#反復測定分散分析を行う
library(car)
time<- factor(c("W0", "W1", "W2"), levels = c("W0","W1", "W2"))
idata <- as.data.frame(time)
model <- lm(cbind(W0,W1,W2) ~ 1, data = grip_rep_anova)
res <- Anova(model, idata = idata, idesign = ~time, type = 3)
summary(res, multivariate = FALSE)

#多重比較
grip_rep_anova_long <- 
  grip_rep_anova %>% 
  gather(key = "key", value = "value", W0:W2)

pairwise.t.test(grip_rep_anova_long$value, grip_rep_anova_long$key, paired = TRUE, p.adjust.method = "bonferroni")
pairwise.t.test(grip_rep_anova_long$value, grip_rep_anova_long$key, paired = TRUE, p.adjust.method = "holm")
pairwise.t.test(grip_rep_anova_long$value, grip_rep_anova_long$key, paired = TRUE, p.adjust.method = "hochberg")



第4章は統計を扱います。


今回も「シロート統計学」のハルさんとコラボレーションすることとなりました。


シロート統計学はEZRを使った統計分析をわかりやすく解説されています。


 


第4章はシロート統計学で使われていたEZRをRで行うとどうなるのか?といった視点で進めていきます。


今回使うデータもハルさんのサイトと同じものを使わせでいただく事になりました。それぞれ見比べることで参考にしてみてください!


今回は一元配置分散分析を紹介します



まず一元配置分散分析(ANOVA)についてはハルさんのサイトをご参照ください。



また1.準備〜4.データの読み込みまでは【4-1】Rでt検定を行う方法と全く同じ流れになります。
もし1〜4まででわからない部分があれば確認してみてください。

 

1.準備

第4章は毎回ExcelデータをダウンロードしてRを使うのでプロジェクトで管理して行うことを勧めています。

 

ここではR練習というプロジェクトを作り、Excelファイルを入れるためのdataフォルダを作っています。
これを前提に次から進めていきます。
スクリーンショット 2019-10-20 7.54.14


2.スクリプトファイルの作成

次にRのコードを書くためのスクリプトファイルを作ります。

スクリーンショット 2019-11-08 14.32.40


3.データのダウンロード

今回もハルさんのサイトのデータを使わせていただきます。

デモデータ(ANOVA)

この章ではRを使ってダウンロードしています。


download.file(url = “ファイルのURL”,
        destfile = “保存したい場所/ファイル名”)
urlはデモデータで右クリック → リンクのアドレスをコピー
destfileは保存場所と保存のファイル名を指定します。


実際のコードは以下になります。
前回のコードのURL(" "の中)とdestfileのdata/以降を変更するだけでOKです。
#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-anova.xlsx"
destfile = "data/demo-anova.xlsx"

download.file(url, destfile)
スクリーンショット 2019-11-08 14.30.34

dataフォルダにダウンロードできたことを確認します。


4.データの読み込み

データを読み込みます。
今回は【4-0】第4章を進めていく上での準備で行った方法で進めます。

View Fileでデータを確認します。

スクリーンショット 2019-11-08 14.24.20



データが入っているセルを確認します。
B2からC22までデータが入っています(B2:C62と表記)
スクリーンショット 2019-11-08 14.26.26


次にImport Datasetでデータを取り込みます。
スクリーンショット 2019-11-07 23.37.11

Import画面ではName, Sheet,Rangeを指定します。

Name:ハルさんのサイトではgripでしたがgrip_anovaとします(大文字・小文字は別物とされます)
Sheet:このExcelは1つしかデータがないのでDefaultのままでOK
Range:先ほど確認したB2:C62

スクリーンショット 2019-11-07 23.40.57



Importボタンを押す前に右にあるコードをコピーしスクリプトファイルに貼り付けることも忘れずに行います。
library(readxl)
grip_anova <- read_excel("data/demo-anova.xlsx", 
                         range = "B2:C62")
View(grip_anova)

データが正しく入っていることを確認します。
スクリーンショット 2019-11-08 14.36.37

これでデータの取り込みは完成です。

5.正規分布の確認

ハルさんのサイトに沿って正規分布の確認を行います。

正規分布の確認は【4-1】Rでt検定を行う方法で紹介しています。
以下のコードがイメージできない場合はまずこちらの記事をご参照ください。



まずヒストグラムで分布を確認します。
1つずつグラフを作るのは手間がかかるので、まとめてグラフを作ります。

ヒストグラムはgeom_histgram関数を使います。
グラフが重なっていると見にくいのでfacet.gridでカテゴリーごとにグラフを分けます。
facet.gridに関しては【3-4】Rのggplot2でヒストグラムを作るgeom_histogram関数で紹介しています。

library(tidyverse)
ggplot(data = grip_anova) +
  geom_histogram(aes(x = grip, fill = category), bins = 5) +
  facet_grid(category ~ .)
スクリーンショット 2019-11-08 15.27.34

6.シャピロ・ウィルク検定
次にハルさんのサイトでは正規性の検定を行っています。
正規性の検討はshapriro.test関数3つのカテゴリーそれぞれに行う必要があります。

今回はsplit関数map関数を使いまとめて行ってみます。
下の図は【4-1】Rでt検定を行う方法で使った方法ですが今回のコードと見比べてみてください。
コードがほとんど同じです。コードを使うと他でも役に立つことがわかります。

スクリーンショット 2019-10-26 22.42.55

grip_anova %>% 
  split(.$category) %>% 
  map(~shapiro.test(.$grip))
スクリーンショット 2019-11-08 16.19.14

まとめて結果が出ました。すべて正規分布であることが棄却されませんでした。


7.等分散性の確認

次に等分散性の確認を行います。

3群以上で等分散性を検定するにはBartlett検定を使います。

bartlett.test(数値の列 ~ グループ)

data = を使えば$は必要ありません。どちらでも大丈夫です。
bartlett.test(grip_anova$grip ~ grip_anova$category)

bartlett.test(grip ~ category, data = grip_anova)



8.分散分析(ANOVA)を行う

それでは本題の一元配置分散分析と多重比較を行います。

実は一元配置分散分析と重回帰分析(t検定も)は同じ分析方法です(詳しく知りたい方は一般線形モデルで検索してみてください)。出てきた結果を分散分析の形式で表示するか重回帰分析の結果で表示するかの違いだけです。
そのため一元配置分散分析を行うには分散分析のaov関数重回帰分析のlm関数で行います。
ただその後に行う多重比較ではlm関数が使えないものもあるので基本はaov関数で大丈夫です。

どちらの関数もまずモデルを作りサマリーを表示するというのが一般的な流れです。

全体的な流れは以下の通りです。
sumamry関数だけ結果の表示が変わるので注意が必要です。
スクリーンショット 2019-11-11 20.05.02


aov関数の場合

モデル名:model_grip_anova(好きな名前でOK)
従属変数:grip
独立変数:category
データ名:grip_anova
model_grip_anova <- aov(grip ~ category, data = grip_anova)
summary(model_grip_anova)
スクリーンショット 2019-11-11 20.09.17


lm関数の場合

モデル名:model_grip_lm(好きな名前でOK)
従属変数:grip
目的変数:category
データ名:grip_anova
model_grip_lm <- lm(grip ~ category, data = grip_anova)
summary.aov(model_grip_lm)
スクリーンショット 2019-11-11 20.10.45



結果はどちらもEZRと同じになりました。


9.多重比較を行う

ハルさんのサイトでは多重比較でTukeyを紹介していますが他も紹介します。
Tukeyだけ関数が違います。
スクリーンショット 2019-11-11 22.18.55


Tukeyの場合
Tukeyの場合は先程のaov関数の結果を使います。
lm関数の結果は使えませんので注意が必要です。
TukeyHSD(model_grip_anova)
スクリーンショット 2019-11-11 22.26.49


Bonferroniの場合
pairwise.t.test関数はモデルの結果を使いません。
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "bonferroni")
スクリーンショット 2019-11-11 22.24.31



Holmの場合
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "holm")
スクリーンショット 2019-11-11 22.24.40


Hochbergの場合
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "hochberg")
スクリーンショット 2019-11-11 22.24.46



10.まとめ
今回は一元配置分散分析を紹介しました。
今回は割愛しましたが、一元配置分散分析を行うときはlongデータである必要があります。

Excelでデータ収集するときに注意するか、wideデータの場合はlongデータに直してから分析を行ってください。

longデータとwideデータについてはで【2-5】Rでデータを集計するのに便利なlongデータとgather関数
紹介しています。

次回はKruskal-Wallis検定を行います。

(追記)グラフの作成
EZRでは検定と同時にグラフを作成しますが、Rでは直接作成します。
【4-1】t検定でグラフを使ったときのコードがそのまま使えます。


グラフを作るために必要な要約のコードもそのまま使います。
#データの要約
grip_anova_summary <- 
  grip_anova %>% 
  group_by(category) %>% 
  summarize(平均 = mean(grip),
              標準偏差 = sd(grip),
              '0%' = quantile(grip, 0),
              '25%' = quantile(grip, 0.25),
              '50%' = quantile(grip, 0.5),
              '75%' = quantile(grip, 0.75),
              '100%' = quantile(grip, 1),
              n = n())
grip_anova_summary
#グラフ化 ggplot(data = grip_anova_summary)+ geom_bar(aes(x = category, y = 平均), color = "black", fill = "gray", stat = "identity") + geom_errorbar(aes(x = category, ymin = 平均 - 標準偏差, ymax = 平均 + 標準偏差), width = 0.1) + labs(y = "grip") + theme_classic(base_size = 15) + scale_y_continuous(expand = c(0,0), limits = c(0,50))
スクリーンショット 2019-11-15 20.13.23
コードの違いはわかったでしょうか?
t検定で使ったデータであるgripを今回のデータであるgrip_anovaに変更しただけです。
それ以外のコードの中身は全く一緒です!

11.今回のコード
今回使ったコードのまとめになります。
#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-anova.xlsx"
destfile = "data/demo-anova.xlsx"

download.file(url, destfile)

library(readxl)
grip_anova <- read_excel("data/demo-anova.xlsx", 
                         range = "B2:C62")
View(grip_anova)


#正規分布の確認
library(tidyverse)
ggplot(data = grip_anova) +
  geom_histogram(aes(x = grip, fill = category), bins = 5) +
  facet_grid(category ~ .)


grip_anova %>% 
  split(.$category) %>% 
  map(~shapiro.test(.$grip))

#Bartlett検定

bartlett.test(grip_anova$grip ~ grip_anova$category)

bartlett.test(grip ~ category, data = grip_anova)

#一元配置分散分析
#aov関数の場合
model_grip_anova <- aov(grip ~ category, data = grip_anova)
summary(model_grip_anova)

#lm関数の場合
model_grip_lm <- lm(grip ~ category, data = grip_anova)
summary.aov(model_grip_lm)

#多重比較
#Tukey
TukeyHSD(model_grip_lm)

#Bonferroni
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "bonferroni")

#Holm
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "holm")

#Hochberg
pairwise.t.test(grip_anova$grip, grip_anova$category, p.adjust.method = "hochberg")

#データの要約 grip_anova_summary <- grip_anova %>% group_by(category) %>% summarize(平均 = mean(grip), 標準偏差 = sd(grip), '0%' = quantile(grip, 0), '25%' = quantile(grip, 0.25), '50%' = quantile(grip, 0.5), '75%' = quantile(grip, 0.75), '100%' = quantile(grip, 1), n = n()) grip_anova_summary
#グラフ化 ggplot(data = grip_anova_summary)+ geom_bar(aes(x = category, y = 平均), color = "black", fill = "gray", stat = "identity") + geom_errorbar(aes(x = category, ymin = 平均 - 標準偏差, ymax = 平均 + 標準偏差), width = 0.1) + labs(y = "grip") + theme_classic(base_size = 15) + scale_y_continuous(expand = c(0,0), limits = c(0,50))


全記事の目次はこちらをご参照ください






↑このページのトップヘ