タグ:R

前回の記事では相関係数と相関係数の検定を行う方法を紹介しました。


ただ複数の列がある場合1つずつ相関係数を調べるのは大変な作業です。
今回は相関行列について紹介します。


1.使うデータ

今回はこちらのダミーデータを使います。

歩行が自立しているか(0:非自立、1:自立)と各種データが入っています。
スクリーンショット 2020-07-15 5.42.46



性別が男性・女性になっています。

url <- "https://raw.githubusercontent.com/mitti1210/myblog/master/data01.csv"
dat <- read.csv(url)
dat$氏名 <- NULL  #氏名の列を削除


2.GGallyパッケージのggpairs関数

相関係数を手早く確認する方法としてGGallyパッケージggpairs関数があります。
この方法のメリットは相関係数と散布図が同時に見れサクッとわかるです。

library(tidyverse)
library(GGally)
ggpairs(dat)
スクリーンショット 2020-07-15 9.30.36


左下に散布図、右上には相関係数とまとめて表示されます。
数値と数値の場合は散布図、カテゴリーと数値の場合はヒストグラムと箱ひげ図が表示されます。

もしmacで日本語が□□と文字化けしていたらもう1行追加します。

ggpairs(dat)+
  theme_gray(base_family = "HiraKakuPro-W3")
スクリーンショット 2020-07-15 9.30.47



そしてggpairsの便利な機能としてはあるカテゴリーごとに色を付けることができます。
aes_string(colour="カテゴリーの列", alpha=0.5)を追加します。
カテゴリーの列はcharactor型、もしくはfactor型である必要があります(数値型だとダメ)
alphaは色の透過の程度なので0.7あたりでも大丈夫です。
ggpairs(dat,aes_string(colour="性別", alpha=0.5))
スクリーンショット 2020-07-15 9.29.32





3.cor関数とcorrplotパッケージを使う

次にcor関数corrplotパッケージで相関行列を作成します。
この方法のメリットはp値が計算できる、pythonっぽい相関行列のグラフが作れるです。

散布図行列を求めるにはcor関数を使います。

ただ男性・女性のような文字は対応できないので0,1に置き換えます。
dat$性別 <- ifelse(dat$性別 == "男性", 0, 1)

c <- cor(dat)
c
スクリーンショット 2020-07-15 9.42.04


四捨五入する時はround関数を使います。
round(c,2)
スクリーンショット 2020-07-15 9.42.38


もしスピアマンの相関係数を求めるのであればmethod="spearman"を追加します。
c <- cor(dat, method = "spearman")

ただ数値の羅列は中々見にくいものです。そこでcorrplotパッケージのcorrplot関数を使いグラフにします。
corrplot(c, method = "color", addCoef.col = TRUE)
method="color"はグラフの見た目、addCoef.col = TRUEはグラフに数値を入れるオプションなのでそのままコピペで大丈夫です。変更するのはcのところで、ここには先程もとめた相関行列を指定します。

c <- cor(dat, method = "spearman")
cの部分です。もちろん違う名前でも大丈夫です。

もしmacで日本語文字化けの問題があれば1行追加します。
par()関数はフォントの指定ができます。

library(corrplot)
par(family="HiraKakuPro-W3") #macで日本語が文字化けするときに追加。ここでフォントの指定ができる
corrplot(c, method = "color", addCoef.col = TRUE)
スクリーンショット 2020-07-15 9.47.46
これでかなり見やすくなりました。pythonっぽいです。


検定をしてp値を求める

p値を求めるのであればcorrplotパッケージのcor.mtest()関数を使います。
cor.mtestはp値、95%信頼区間の下限、上限の3つが表示されます。スピアマンの相関係数の場合はmethod = "spearman"を追加するだけなのですが、同値(タイ)があると95%信頼区間が出ないようです。

p <- cor.mtest(dat)
p
スクリーンショット 2020-07-15 7.20.20

p$p:p値
p$lowCI:95%信頼区間の下限
p$uppCI:95%信頼区間の上限

先程の相関行列にp<0.05かそうでないかを出すこともできます。
p.mat:p値の行列(ここではp$p。上の図をよく見ると$pと書いてあります。)
sig.level = ◯◯(0.05 にしたり、0.2に設定することもあります)

corrplot(c, method = "color", addCoef.col = TRUE, p.mat = p$p, sig.level = 0.05)
スクリーンショット 2020-07-15 10.33.07

これで有意差があるものがどれか一目瞭然です。


4.corrrパッケージを使う

もう一つcorrrパッケージを紹介します。
この方法のメリットは相関が高い組み合わせ順が知りたい!というときに便利です。

相関行列を出す時はcor関数でなくcorrelate関数を使います。
使い方はcor関数と同じです。cor→correlateに変わるだけです。
スピアマンの相関係数にしたい時はmethod="spearman"を追加すればOKです。

c <- correlate(dat, method = "spearman")
c
スクリーンショット 2020-07-15 10.41.13


マイナスの値が赤字になっただけでここまでは今までと違いはありませんが、corrrの便利な所はここからになります。

まずcorrr関数は第2章で扱ったtidyverseと同じ%>%が使えます。

右上と左下は同じものなので片方を削除することができます。削除するにはshave関数を使います。
c %>% shave()
スクリーンショット 2020-07-15 10.45.16



これで左下だけになりました。

次に行列形式ではなくlong形式に変換します。streach関数を使います。
streach(na.rm=TRUE)とするとNAとなっているいらない列を消すことができます。

c %>% shave() %>% stretch(na.rm = TRUE)
スクリーンショット 2020-07-15 10.55.02



次に並べ替えます。arrange関数を使います。何も指定しなかければ相関係数が高い順から並びます。
c %>% shave() %>% stretch(na.rm = TRUE) %>% arrange()
スクリーンショット 2020-07-15 10.55.14
スクリーンショット 2020-07-15 10.55.31


ただ注意しないといけないのが相関係数は-1〜1の値を取ります。
相関がないのは0付近で、-1は1と同様に強い相関があります。
なので相関が強い順に並べるには相関係数の絶対値で並び替える事が必要になります。
相関係数の列名はrになっています。abs関数が絶対値を求める関数なのでabs(r)となります。

c %>% shave() %>% stretch(na.rm = TRUE) %>% arrange(abs(r))
スクリーンショット 2020-07-15 10.55.31


そしてこれだと相関係数が0に近い方から並んでしまったので順番を逆にします。
順番を逆にするにはdesc(並べ替えたい列)を使います。
並べ替えたい列はabs(r)でした。desc(abs(r))とします。
()が多くなってエラーが出やすくなるので注意してください。

c %>% shave() %>% stretch(na.rm = TRUE) %>% arrange(desc(abs(r)))
スクリーンショット 2020-07-15 10.55.45

相関行列を1個ずつ見てると見落としがちですが、これなら見落とすことはありません。


5.相関係数の有意差がないってどういう意味?

スクリーンショット 2020-07-15 10.33.07

ところで×がついた項目(p値が有意水準に満たさなかった)というのはどういう意味でしょうか?
これは相関係数が低いから有意差が出ないわけではありません。

相関係数の検定は「相関係数が0である」と仮説を立てます。
そしてデータから信頼区間を求め、信頼区間に0が含まれているかを調べます。
その時もし信頼区間に0が含まれていたら有意差なしとなります。

例えば上の図でいうと性別×BBSを見ると相関係数が-0.13となっています。
ただ95%信頼区間を見ると-0.27〜0.007となっています。
つまり今回は計算上-0.13となったけど実は信頼区間にマイナスの値も0もプラスの値も入っていて、符号が変わる可能性がある→正の相関があるとも負の相関があるとも言えないよね。といったニュアンスになります。もし95%信頼区間が-0.27〜-0.001であれば符号が変わらないため有意差あり(p<0.05)となります。ただ相関係数0.13というのはほとんど相関がないという意味ではありますが...

つまり相関係数の検定のp値は出した相関係数の符号が変わる可能性が高いのか低いのか?を見る検定と考え、p値が低いから相関係数が強いんだという意味ではないことに注意が必要です。相関の強さはあくまでもp値ではなく相関係数の値で評価します。



6.まとめ

今回は相関行列を紹介しました。

相関行列を出すだけならEZRでもできますが、検定やグラフの作成はEZRだけではできません。

他にもいろいろな方法はあるのですが、今回は3パターン紹介しました。

ただGGallyやcorrrはもっといろいろなことができます。

もっと詳しく知りたい方は以下のサイトがおすすめです!








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


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


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


 


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


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


今回は相関係数(ピアソン・スピアマン)を紹介します



まず相関係数についてはハルさんのサイトをご参照ください。



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

 

1.準備

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

 

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


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

次にRのコードを書くためのスクリプトファイルを作ります。
ここでは【4-15】相関係数.Rとしておきます。

スクリーンショット 2020-07-14 23.41.44


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


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

デモデータ(pearsonの相関係数)

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


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


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


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


4.データの読み込み

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

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

スクリーンショット 2019-12-08 20.31.39


データが入っているセルを確認します。
B2からC32までデータが入っています(B2:C32と表記)
スクリーンショット 2020-07-14 23.51.28




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



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

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

スクリーンショット 2020-07-14 23.53.34



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

データが正しく入っていることを確認します。
スクリーンショット 2020-07-14 23.54.38



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

5.散布図で確認

相関係数を見るときはいきなり相関係数を求めるのではなく、先に散布図を書くほうが事故が減ります。
【3-3】Rのggplot2で散布図を作るgeom_point関数で散布図を書く方法を紹介しています。
散布図を書くのはgeom_point関数、回帰直線を書くのはgeom_smooth関数です。
geom_smooth関数では(method="lm")を必ず付ける必要があります。
これは直線ですよという意味です。直線以外にも色々あるんです。
ちなみに色を付けるときはcolor=を、se=FALSEは今回はあってもなくても構いません。

library(tidyverse)
ggplot(data=borg,aes(x=`6MWD`,y=VC))+
  geom_point()+
  geom_smooth(method = "lm", color="red", se=FALSE)
注意!
列名や変数名の最初の文字が数字の時(今回でいうと6MWD)の場合、x=6MWDのようにするとエラーが出ます。こういった場合は` `(shiftキーを押しながら@マークを押す)でエラーを回避できます。でもなるべくなら最初の文字は数字じゃない方が後々便利です。


スクリーンショット 2020-07-15 0.06.04


他にも「もっとかんたんに散布図だけみたい」というときはplot関数があります。

plot(borg$`6MWD`, borg$VC)
ここでも`6MWD`と` `を忘れないようにします。
スクリーンショット 2020-07-15 0.05.51

こうみると6MWDとVCは直線上の関係にあるのが見て取れます。
こういった場合はピアソンの相関係数が使えます。
逆に外れ値が多いデータや直線でないデータでは相関係数は使えなくなります。

外れ値がない場合
スクリーンショット 2020-07-15 0.34.07

外れ値が5個あった場合
スクリーンショット 2020-07-15 0.33.58


そして二次関数のようなデータはうまく行きません。他にも「10代〜30代は良好、40〜50代は不良、60代以降は良好」みたいなデータもうまくいきません。
(相関係数は「一方が上がる」するともう一方が「上がるor下がる」ということだけで、2時間数のように途中で傾向が変わるデータには対処できない)
スクリーンショット 2020-07-15 0.36.06
上記のような例は相関係数を出すだけではわかりません。散布図で確認して初めて分かることなので、相関係数は「まずは散布図!」とイメージしておいてください。


6.相関係数の検定の実施

相関係数の検定は1行で終わりますが主に2種類あることを知る必要があります。

パラメトリック:ピアソンの相関係数
ノンパラメトリック、外れ値が多い:スピアマンの相関係数



どちらも相関係数の検定はcor.test(1つ目の列, 2つ目の列)です。

cor.test(borg$`6MWD`,borg$VC)
スクリーンショット 2020-07-15 1.02.00

もしスピアマンの相関係数を求めるときはmethod = "spearman"を追記します。
cor.test(borg$`6MWD`,borg$VC, method = "spearman")
スクリーンショット 2020-07-15 4.00.34


【4-2】RでMann-Whitney U 検定を行う方法でも紹介しましたが、ノンパラメトリックの手法では列の中に同じ数値があると警告がでます。

スクリーンショット 2020-07-15 4.07.07
このままでも大丈夫ではありますが(EZRも同じ字問題が起こっているが警告が出ないようにしている)、もし「警告をどうにかしたい!」とあれば、coinパッケージのspearman_test関数を使うという方法もあります(そこまでしている人はいるのでしょうか?)

spearman_test(1つ目の列 ~ 2つ目の列, data=データ)という書き方になります。

library(coin)
coin::spearman_test(`6MWD`~VC, data=borg)

cor(borg$`6MWD`,borg$VC, method = "spearman")
スクリーンショット 2020-07-15 4.17.43


ただspearman_test関数は検定はしてくれますが、相関係数自体を出してくれません。
そのためcor関数で相関係数を求める必要があります。


7.まとめ

今回は相関係数の求め方を紹介しました。

今回紹介していませんが、ハルさんのサイトにもあるように相関関係と因果関係は別物という考えはとても大切です。

そして、相関係数を出す前に散布図を作るという習慣をつけることが大切です。

散布図を出した上で直線的な関係性がありそうであれば相関係数を求めてください。

データを扱う中で「外れ値を外す方法」に関して検討してみます。

apply関数の記事で以下のコメントがありました。




大変勉強になりました。
列ごとに、ファイルの値(NAがある)が2SDより高い値をNAにしたいです。
Data01 <-read_excel("raw.xlsx",1)
Data02 <-apply(Data01, 2, sd, na.rm = TRUE)
Data01[Data01 > 2*Data02]<- NA
以上のコードをして最後は列ごとにうまく出来なかったです。
ご提言があれば大変嬉しいです。

自分もそうですが、実際にやってみるとうまくいかないことってたくさんあります。
今回は質問を一緒に考えながらどう考えるか?どんな解決方法があるかを考えて、今までの理解を深めるきっかけになればと思います。

1.データの作成

コードを確認してみます。
Data01 <-read_excel("raw.xlsx",1)
これはraw.xlsxの1シート目を読み込むという意味なので、今回は仮のデータを作ります。

set.seed(2020)
a <- c(rnorm(998,10,1),NA,30)
b <- c(rnorm(998,20,1),NA,40)
c <- c(rnorm(998,30,1),NA,50)
Data01 <- data.frame(a,b,c)
head(Data01)
簡単に言うと以下のとおりです。
a : 平均10, 標準偏差1のデータ998個にNAと30を加えた1000個のデータ
b : 平均20, 標準偏差1のデータ998個にNAと40を加えた1000個のデータ
c : 平均30, 標準偏差1のデータ998個にNAと50を加えた1000個のデータ

スクリーンショット 2020-03-28 2.34.45
この赤色が2SDを超えたデータです。やりたいことはこの赤い点のデータを削除することです(これの是非については今回問いません)

2.apply関数のおさらい

次の行を確認します。
Data02 <-apply(Data01, 2, sd, na.rm = TRUE)
apply関数【2-1】Rのfor関数、apply関数を使ってまとめて標準偏差などの統計量を求める方法で紹介しました。下はイメージ図です。
スクリーンショット 2020-03-27 13.46.36
上図ではsd関数(標準偏差)を求めました。
また欠損値(NA)があるのでrm.na = TRUEをつけることでNAの値を外して計算してくれます。

スクリーンショット 2020-03-27 13.54.01
それぞれの列の標準偏差が表示されました。
ただ元データが2標準偏差を超えているかどうかを確認するにはこれだけでは足りません。


3.どう考えるか?
スクリーンショット 2020-03-28 4.00.30

それぞれのデータが2SDを超えているかどうかを判定する必要があります。
具体的には(平均-2SD) < それぞれの値 < (平均 + 2SD)の条件を満たしているかどうかとなります。

そのため以下の2つを行えばなんとかできそうです。
①a,b,cそれぞれの平均±2標準偏差の値を求める
②条件に合えばそのまま、条件にあわなければNAと変換する


 4.平均±2SDを求める

まずa,b,cの平均と標準偏差を求めます。
これには前述のapply関数かsummarize関数を使います。
まずはapply関数でやってみます。

Data01_mean <-apply(Data01, 2, mean, na.rm = TRUE)
Data01_sd <-apply(Data01, 2, sd, na.rm = TRUE)
Data01_mean
Data01_sd

low <- Data01_mean - 2*Data01_sd
high <- Data01_mean + 2*Data01_sd
low
high
スクリーンショット 2020-03-28 3.17.41

ここからlow[2]、もしくはlow["b"]とすればbの平均-2SDが取り出せます。
high[1]、もしくはhigh["a"]であればaの平均+2SDとなります。

5.条件に合わないものをNAにする

次に条件に合わないものを抽出しNAにします。
データを抽出するにはデータ[条件]を使います

(aのデータ < 平均-2sd) または (平均+2SD < aのデータ)だとこうなります
Data01$a[(Data01$a < low["a"]) | (high["a"] < Data01$a)]
スクリーンショット 2020-03-28 4.08.45

条件式の記号については以下をご参照ください。今回はまたはの意味の を使いました。
スクリーンショット 2019-02-11 2.04.40
抽出した条件に <- を使うことでデータを置き換えることができます。
Data02 <- Data01
Data02$a[(Data02$a < low["a"]) | (high["a"] < Data02$a)] <- NA
Data02$b[(Data02$b < low["b"]) | (high["b"] < Data02$b)] <- NA
Data02$c[(Data02$c < low["c"]) | (high["c"] < Data02$c)] <- NA
これで完成です!
ただ他にも方法があります。どれもよく使う方法なので、もし余力があればぜひ比べてみてください。


6.tidyverseパッケージを使う(dplyrパッケージでも同じ)

第2章で紹介しているtidyverseパッケージのsummarize関数とmutate関数を使うこともできます。
tidyverseパッケージを使うときには先にlibrary関数で呼び出します。なければ先にインストールします。
#もし初めての場合はインストール
install.packages("tidyverse")
#インストール済ならlibrary関数で呼び出す
library(tidyverse)
low_high <-
  summarize(Data01,
          a_low = mean(a, na.rm = TRUE)-2*sd(a, na.rm = TRUE),
          a_high = mean(a, na.rm = TRUE)+2*sd(a, na.rm = TRUE),
          b_low = mean(b, na.rm = TRUE)-2*sd(b, na.rm = TRUE),
          b_high = mean(b, na.rm = TRUE)+2*sd(b, na.rm = TRUE),
          c_low = mean(c, na.rm = TRUE)-2*sd(c, na.rm = TRUE),
          c_high = mean(c, na.rm = TRUE)+2*sd(c, na.rm = TRUE)) 
class(low_high) #型を確認するとdata.frameとなっている
low_high <- unlist(low_high)
class(low_high) #型を確認するとnumeric(数値型)となっている
low_high
スクリーンショット 2020-03-28 3.37.43
文字は多いですが、コピペを使えば作業自体は楽になります。
ただsummarize関数の結果は表形式(data.frame型)で、このままでは数値として計算できません。
そのためunlist関数を使って数値型に戻します。
low_high[5]、もしくはlow_high["c_low"]とすればcの平均-2SDとなります。


既に第2章を読まれた方は%>%を使う事もできます。
#%>%使う場合
low_high <- summarize(Data01, a_low = mean(a, na.rm = TRUE)-2*sd(a, na.rm = TRUE), a_high = mean(a, na.rm = TRUE)+2*sd(a, na.rm = TRUE), b_low = mean(b, na.rm = TRUE)-2*sd(b, na.rm = TRUE), b_high = mean(b, na.rm = TRUE)+2*sd(b, na.rm = TRUE), c_low = mean(c, na.rm = TRUE)-2*sd(c, na.rm = TRUE), c_high = mean(c, na.rm = TRUE)+2*sd(c, na.rm = TRUE)) %>%
unlist()
low_high
スクリーンショット 2020-03-28 3.32.12

次に条件に応じて分類するifelse関数を使うこともできます。

ifelse(条件TRUEの場合FALSEの場合)

列の追加や修正を行うmutate関数を使うことでこうなります。
Data02 <- 
  Data01 %>% 
  mutate(a = ifelse((low_high["a_low"] < a) & (a < low_high["a_high"]), a, NA),
         b = ifelse((low_high["b_low"] < b) & (b < low_high["b_high"]), b, NA),
         c = ifelse((low_high["c_low"] < c) & (c < low_high["c_high"]), c, NA))
スクリーンショット 2020-03-28 4.57.19

条件式を少し変えているので確認してみてください。今回は&を使っています。

tidyverseパッケージは現在のRを使う上でとても重要になります。
ifelse,mutate関数に関してはこちらで紹介していますので難しかったという方はご参照ください。




7.関数を自作しapply関数でまとめて処理する
他の方法としては関数を自作するという方法があります。

関数名 <- function(○) {中身のコード}

今回はdrop_2sdという名前の関数を作ります(名前は何でもいいけど文字の最初に数字を使うのは×)。


a = ifelse((aの平均-2SD < a) & (a < aの平均+2SD), a, NA),
b = ifelse((bの平均-2SD < b) & (b < bの平均+2SD), b, NA),
c = ifelse((cの平均-2SD < c) & (c < cの平均+2SD), c, NA)
 ↓
function(x){
 ifelse((xの平均-2SD < x) & (x < xの平均+2SD), x, NA)
}
と共通部分をxにします。ちなみにxでなくてもcolとか別の名前をつけても大丈夫です。

drop_2sd <- function(x){
  ifelse((mean(x, na.rm = TRUE) - 2 * sd(x, na.rm = TRUE) < x) & (x < mean(x, na.rm = TRUE) + 2 * sd(x, na.rm = TRUE)), x, NA)
}
そしてapply関数で各列にdrop_2sd関数を使います。
Data02 <- apply(Data01, 2, drop_2sd)
スクリーンショット 2020-03-28 4.50.32
2行のコードで終わりました。
関数を作るのは最初は慣れないですが、使えるようになると強力です。


8.まとめ

今回は標準偏差±2SDの除去を3つの方法で紹介しました。

本やサイトで独学すると1つの方法でしか紹介してないことが多いです。

本当は自分の知っている知識でもできるはずなのに、知らない方法を偶然読んでしまいドツボにはまることもあります。そういった紆余曲折は勉強する上で必要な苦労でもありますが、いくつかの方法を見比べることでショートカットできることもあるのではと思っています。

また自分が躓いたところや他の記事でもコメントがあれば検討していきますので、ブログのコメントやtwitterで連絡いただければ幸いです。

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


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


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


 


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


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


今回はFriedman(フリードマン)検定を紹介します



まずFriedman(フリードマン)検定についてはハルさんのサイトをご参照ください。



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

 

1.準備

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

 

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


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

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

スクリーンショット 2019-12-08 20.26.22

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


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

デモデータ(Friedman検定)

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


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


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


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


4.データの読み込み

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

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

スクリーンショット 2019-12-08 20.31.39


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



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



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

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

スクリーンショット 2019-12-08 20.35.24


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

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


これでデータの取り込みは完成です。
Friedman検定は繰り返しのデータ(被験者内要因)なのでデータを横につなげています。


5.Friedman検定を行う(wideデータの場合)

ハルさんのサイトでは早速Friedman検定を行っていますのでこちらも早速行います。
今回はwideデータなのでまずwideデータの場合を紹介します。

Friedman検定はfriedman.test関数を使います。

friedman.test(cbind(繰り返しデータのある列))
スクリーンショット 2019-12-09 0.52.59

friedman.test(cbind(grip_friedman$W0,grip_friedman$W1,grip_friedman$W2))
スクリーンショット 2019-12-09 0.54.04

EZRと同じ結果が出ました。

6.Friedman検定を行う(longデータの場合)
次にlongデータでFriedman検定を行う方法も紹介します。

friedman.test(データ ~ グループ | 個人を識別するid, data = ○○)

スクリーンショット 2019-12-09 0.58.27

Friedman検定は対応のあるデータになります。

今回は30人がW0,W1,W2の3回データを取っているので3*30=90データあります。
そのため作るlongデータは上記のイメージになります。

今回はid, time, gripとしてますが自分がわかれば好きな名前を付けてOKです。
id:30人の内誰のデータかわかるためのid(氏名でも可)
time:W0, W1, W2の3つ
grip:今回測定した握力のデータ
wideデータをlongデータに変えるにはgather関数、もしくはpivot_longer関数を使います。
gather関数に関しては【2-5】Rでデータを集計するのに便利なlongデータとgather関数で紹介しています。



データの扱いはtidyverseパッケージに含まれていますので。まずlibrary関数でtidyverseパッケージを読み込みます。今回はgrip_friedman_longという名前にします。

library(tidyverse)

#gather関数の場合 grip_friedman_long <- grip_friedman %>% rowid_to_column("id") %>% gather(key = time, value = grip, W0:W2, factor_key = TRUE)

#pivot_longer関数の場合 grip_friedman_long <- grip_friedman %>% rowid_to_column("id") %>% pivot_longer(cols = c(W0,W1,W2), names_to = "time", values_to = "grip", names_ptypes = list(time = factor()))

ポイントはもともとのExcelデータにidにあたる列がないことです。
そこでrowid_to_column関数を使います。""の中身が列名になります。
スクリーンショット 2019-12-09 1.14.03

これで準備が整いました。

longデータでFriedan検定を行う場合は改めてこうなります。

friedman.test(データ ~ グループ | 個人を識別するid, data = ○○)
friedman.test(grip ~ time | id, data = grip_friedman_long)
スクリーンショット 2019-12-09 0.54.04
もちろん結果は同じです。

7.多重比較を行う
Friedman検定で多重比較を行いたいのですがwideデータだとEZR専用のパッケージでしかできません。
そのためRで行うにはlongデータで行う必要があります。

longデータで多重比較を行うにはpairwise.wilcox.test関数を使います。簡単に言うと全ての組み合わせでMann-Whitney U検定を行いその結果を元にp値の調整を行います
どの方法を使うかはp.adjust.method = で指定します。
またタイのエラーを回避するためにexact = FALSEを入れています。
data = 

pairwise.wilcox.test(目的変数factor(グループ),  p.adjust.method = "◯◯", paired = TRUE, exact = FALSE)

#Bonferroniの方法
pairwise.wilcox.test(grip_friedman_long$grip, grip_friedman_long$time, p.adjust.method = "bonferroni", paired = TRUE, exact = FALSE)
#Holmの方法
pairwise.wilcox.test(grip_friedman_long$grip, grip_friedman_long$time, p.adjust.method = "holm", paired = TRUE, exact = FALSE)
スクリーンショット 2019-12-09 1.39.44
EZRとほぼ同じ結果になりました。

8.どうしてもwideデータのままなら手計算が必要

wideデータのまま多重比較を行う方法がEZRのパッケージしかありません。

そのためどうしてもwideデータのまま行うにはwilcox.test関数、もしくはcoinパッケージwilcox_test関数を使い、1つずつ組み合わせを行い、BonferroniHolmHochbergそれぞれの方法で自分でp値を直接計算することになります。

この方法は統計ERさんの記事に詳しく紹介されています。
更にBonferroniHolmHochbergって何をしてるの?という疑問も解決するかもしれません。
仕組みがわかればp値に2とか3とかかけるだけなので意外と単純です。
 


8.まとめ

今回はFriedman検定を紹介しました。

もしRで行う場合はwideデータで分析するかlongで分析するかをあらかじめイメージしていた方がその後の解析もスムーズかもしれません。

次は相関係数について紹介します!

9.今回使ったコード
#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-friedman-rank-sum-test.xlsx"
destfile = "data/demo-friedman-rank-sum-test.xlsx"

download.file(url, destfile)

library(readxl)
grip_friedman <- read_excel("data/demo-friedman-rank-sum-test.xlsx")
View(grip_friedman)

friedman.test(cbind(grip_friedman$W0,grip_friedman$W1,grip_friedman$W2))

#pivot_longer関数の場合
grip_friedman_long <-
  grip_friedman %>% 
  rowid_to_column("id") %>% 
  pivot_longer(cols = c(W0,W1,W2), names_to = "time", values_to = "grip", names_ptypes = list(time = factor())) 

#gather関数の場合
grip_friedman_long <-
  grip_friedman %>% 
  rowid_to_column("id") %>% 
  gather(key = time, value = grip, W0:W2, factor_key = TRUE)

friedman.test(grip ~ time | id, data = grip_friedman_long)

pairwise.wilcox.test(grip_friedman_long$grip, grip_friedman_long$time, p.adjust.method = "bonferroni", paired = TRUE, exact = FALSE)
pairwise.wilcox.test(grip_friedman_long$grip, grip_friedman_long$time, p.adjust.method = "holm", paired = TRUE, exact = FALSE)
 

記事一覧はこちらになります。


第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君はあまりにも有名なので是非試してみてください!




↑このページのトップヘ