タグ:factor

(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")



今までRでデータの変数名を変更したり条件でグループ化したりしました。







そして前回の記事では集計して平均や標準偏差などの要約を出したり、グラフを作るためにgather関数を使ってlongデータを作りました。

【2-5】Rでデータを集計するのに便利なtidyデータとgather関数



今回は棒グラフや折れ線グラフ作成に必要な平均や標準偏差などの統計量を求めます。



データは前回と同じデータを使います。

FIM.xlsx

スクリーンショット 2019-02-19 22.25.41


ダウンロードした後、プロジェクトの指定フォルダにファイルを移動させておけば、以下のコマンドで前回最後の場面まで進みます。

library(tidyverse)
library(readxl)
fim <- read_excel("FIM.xlsx", sheet = "入院時")
fim_long <- fim %>% 
  gather(食事:FIM合計, key = 項目, value = 点数, factor_key = TRUE) 
head(fim_long)


1,groop_by関数でグループ化したい項目を指定する

ますグループ化するためにはgroop_by関数を使います。

group_by(データ,列名)

%>%を使うとデータの部分は省略できます。

fim_long %>% 
  group_by(項目) 

スクリーンショット 2019-02-21 21.42.39



何も変わってないように見えますが、薄い文字のところに# Groups:項目 [21]とあります。


ちなみにグループの種類が複数でも可能です。

fim_long %>% 
  group_by(項目,性別) 

スクリーンショット 2019-02-21 21.45.43



2.列名の順番について

下の2つのコードはどういった違いがあるのでしょうか?

スクリーンショット 2019-02-22 0.15.27

表の並び順が違うだけでその後グラフを作るときには影響はないのですがイメージとしては上記のようになります。





3.統計量を出すsummarize関数


平均などの統計量などを出すにはsummarize関数を使います。

summarize(名前1 = 関数1, 名前2 = 関数2)

summarize関数の前にgroup_by関数を使っていると、グループごとの集計が出てきます。

fim_summarize <-  fim_long %>% 
  group_by(項目,性別) %>% 
  summarize(平均 = mean(点数), 
            標準偏差 = sd(点数), 
            最小値 = min(点数), 
            最大値 = max(点数))
fim_summarize

スクリーンショット 2019-02-22 0.49.09




もし保存したい場合はwrite.csv関数を使います。

wite.csv(保存する変数名, "ファイル名.csv")
ファイル名には" "と.csvを入れます

write.csv(fim_summarize,"FIM集計")
スクリーンショット 2019-02-22 0.57.54
右のfilesビューにFIM集計.csvができました。

csvファイルをExcelで読み込む時は【1-11】Rで医療統計で必要なtable1を作るtableoneパッケージについて紹介しますをご参照ください。

ポイントとしてはファイルの出力先を$A$1ではなく$B$1にしてください。
$A$1だとなぜかエラーが出ます。
スクリーンショット 2019-02-22 1.04.57

スクリーンショット 2019-02-22 1.07.22

Excelにするとわかるのですが、小数点10桁まで表示されます。
もし小数点第一位までの表示にしたい時はround関数を使います。

round(数値, x)
たとえばxが1だと小数点第二位を四捨五入して小数点第一位まで表示します。

fim_summarize <- fim_long %>%
group_by(項目,性別) %>%
summarize(平均 = round(mean(点数), 1),
   標準偏差 = round(sd(点数), 1),
   最小値 = min(点数),
   最大値 = max(点数))
fim_summarize
スクリーンショット 2019-02-22 1.18.28



4.summarize関数を使うときの注意点。

group_by関数の順番に影響する

2でも紹介しましたが、groop_by関数で指定したグループが複数の場合、summarize関数で表示される順番はgroop_by関数の影響を受けます。

スクリーンショット 2019-02-22 0.15.27

あとでExcelで並べ替えるのはただ面倒です。もしcsvで保存をする時は何を示したいかをあらかじめイメージしておくことが重要になります。



2.欠損地があるとNAとなる

平均を求めるmean関数などはどこか1つでも欠損値(空欄)があると結果はNAとなります。

氏名 <- c("A", "B", "C", "D")
年齢 <- c(55, 63, 67, 71)
test_1回目 <- c(1,2,3,NA)
test_2回目 <- c(5:8) test_3回目 <- c(9:12)
data <- data_frame(氏名, 年齢, test_1回目, test_2回目, test_3回目)
data data %>%
gather(3:5, key = 回数, value = 点数) %>%
group_by(回数) %>%
summarize(平均 = mean(点数))
スクリーンショット 2019-02-22 1.53.38


もし結果にNAが出た時はまずはそもそものデータの入れ損ねがないか確認をし対応します。

それでも欠損値がある時はNAを省いて計算する欠損値を統計の技術を使って代入するといった方法があります。

もし欠損値を省いて平均を出す場合はna.rm = TRUEを付け加えます。

data %>%
gather(3:5, key = 回数, value = 点数) %>%
group_by(回数) %>%
summarize(平均 = mean(点数, na.rm = TRUE))
スクリーンショット 2019-02-22 2.06.51


ただ実は欠損値を省いた方がいいのかどうかという問題があります。ただ業務で傾向を確認したいなどであれば欠損値を省いてもいいと思いますが、きちんと出さないと行けない場面ではこれはこれできちんと勉強する必要があります。欠測データに関しては医療統計の本では紹介されていない事が多く専門書が必要かもしれません。


欠測データ処理: Rによる単一代入法と多重代入法 (統計学One Point)



summarize関数に入れられる関数は単一の値が出るものに限る。

summarize関数で使える統計量はmean関数,sd関数,median関数など単一の値になります。

ただ、関数の中には最小値と最大値を一度に出してくれるrange関数など複数の値を出すものがあります。

test_2回目 <- c(5:8)
range(test_2回目)
スクリーンショット 2019-02-22 2.27.00


range関数をsummarize関数に入れようとするとエラーが出ます。

data %>% gather(3:5, key = 回数, value = 点数) %>% group_by(回数) %>% summarize(平均 = mean(点数, na.rm = TRUE), 範囲 = range(点数))

スクリーンショット 2019-02-22 2.36.38

「1つの値しか入れられないのに2つ入ってるよ!」と怒られています。

スクリーンショット 2019-02-22 2.41.49

特にあるのが四分位範囲です。
四分位範囲はquantile関数を使って以下のように一度に値を出すことができます。
quantile(fim$年齢, c(0.1, 0.25, 0.5, 0.75, 0.9))
スクリーンショット 2019-02-22 2.45.59


しかしsummarize関数ではまとめて入れられないので1つずつ入れる必要があります。
fim_summarize <- fim_long %>% group_by(項目,性別) %>% summarize(平均 = mean(点数), 標準偏差 = sd(点数), 最小値 = min(点数), percent_25 = quantile(点数, 0.25), 中央値 = median(点数), percent_75 = quantile(点数, 0.75), 最大値 = max(点数)) fim_summarize
スクリーンショット 2019-02-22 3.01.32



そもそもの列のfactorの順番があってるのか?

gather関数のkey列に関してはfactor_key=TRUEで五十音順ではなく、元の順番に戻すことができます。

しかし他の列でfactorの順番が合っていない可能性もあります。

もしほかの列でfactorの順番を合わせるにはmutate関数とfactor関数を組み合わせて使うことができます。

今回は女性と男性を入れ替えてみます。
fim_summarize <- fim_long %>%
mutate(性別 = factor(性別, levels = c("男性","女性"))) %>%
group_by(項目,性別) %>%
summarize(平均 = mean(点数),
     標準偏差 = sd(点数),
      最小値 = min(点数),
   percent_25 = quantile(点数, 0.25),
   中央値 = median(点数),
   percent_75 = quantile(点数, 0.75),
      最大値 = max(点数))
fim_summarize
スクリーンショット 2019-02-22 3.01.32


factor関数の使い方は【1-10】Rでよく使われる型について説明しますをご参照ください


まとめ

今回はgroup_by関数とsummarize関数と、実際に集計を出し保存するところまで紹介しました。

これでひとまず第2章で行う予定だった「データを集計しやすいように形を整え集計する」が終わりです。

スクリーンショット 2019-02-08 10.33.40


ただ今回は「できるだけExcelの時点でデータは編集しやすい形式で保存している」ことを前提に話を進めています。

もっと形が整っていないデータの扱いやここでは説明できなかった項目も多くあります。

もし「もっと知りたい」「ここの情報では足りない!」ということであれば下記のサイトなどもご参照ください。


データハンドリング入門
https://kazutan.github.io/kazutanR/hands_on_170730/index.html




↑このページのトップヘ