カテゴリ: 統計

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

第4章ではt検定から始まり色々な種類の統計やp値や信頼性の考え方などを紹介してきました。

そして内容は後半の多変量解析に入り、前回は多変量解析を行う上で大切な考え方を紹介しました。

【4-17】「先にデータを取って後で統計」がなぜ危険なのか?どうすればいいのか?を多変量解析、相関・因果関係、交絡を使って解説してみる


そしてもう1つ重要な考え方があります。それは洞察(因果関係)を見たいのか?それとも予測がしたいのか?です。

これは統計学と機械学習の違いにも通じる話になります。

1.使うデータ

ここでも前回の記事と同じデータを使います。

架空のデータです。

library(tidyverse)
set.seed(25) v1 <- rnorm(100,87,13) v2 <- rnorm(100,63,10.5) v <- c(v1,v2) a1 <- as.integer(rnorm(100,35,5)) a2 <- as.integer(rnorm(100,55,5)) age <- c(a1,a2) sex <- c(sample(c(0,1),size = 100, prob = c(0.7,0.3),replace = TRUE),sample(c(1,0),size = 100, prob = c(0.7,0.3),replace = TRUE)) t1 <- as.integer(rnorm(100,120,20)) t2 <- as.integer(rnorm(100,70,20)) time <- c(t1,t2) df <- data.frame(v,time,age,sex) df

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

v:歩行速度、1分間に歩いた距離(m)とします。

time:片脚立位の時間(s)

age:年齢

sex:性別(0:男性,1:女性)

これをみて片脚立位の時間が長くなれば歩行速度が上がると考えていいかどうかとします。




2.多変量解析はどういう仕組みか?

スクリーンショット 2020-07-17 14.34.46

スクリーンショット 2020-07-17 14.27.48



まず多変量解析(ここでは重回帰分析)は左辺(y)にアウトカム、右辺に要因(x)をおいた式となります。yとxには色々呼び方がありますが同じ意味です。

そして数学的なテクニックを使い(最小二乗法や最尤法)を使い係数を求めます。

コードと解釈の仕方は次の記事で紹介します。

fit <- lm(v ~ time + age + sex)
summary(fit)

スクリーンショット 2020-07-17 14.47.59

色々見るところはあるのですが、とりあえず黄色に注目します。

Estimateというのが係数(β)になります。(Intercept)は切片です。

1分で歩いた距離 = 0.03×片脚立位時間(秒)+ -0.84×年齢(歳) + 0.42(女性であれば) + 108.33

ロジスティック回帰もそうですが、やっていることは係数を求めているということです。


3.重回帰分析とロジスティック回帰の使い分け方

重回帰分析とロジスティック回帰の使い分けはアウトカムを見ます。

重回帰分析:アウトカムが連続した数値になっている
ロジスティック回帰:アウトカムが「自立・非自立」といった2択になっている


ちなみにアウトカムが3択以上のある場合は自分の知りたいことにあわせて2択にしてしまったほうがいいです。
3択以上の解析は複雑になります。


4.xを見たいのか?yを見たいのか?

じつはこの数式には2つの見方があります。それはxを見たいのか?yを見たいのか?です。


洞察(因果関係)に興味がある=xに興味がある

まずxを見るときの視点を説明します。

上の式を再掲します
1分で歩いた距離(m) = 0.03×片脚立位時間(秒)+ -0.84×年齢(歳) + 0.42(女性であれば) + 108.33

この式は右辺に片脚立位時間、年齢、性別という3つの変数があります。

ここで年齢と女性でないという条件で片脚立位の時間が1秒増えるとどうなるでしょうか?

答えは歩行距離が0.03m増えるということになります。

逆に片脚立位・女性でない条件で年齢が1つ増えれば歩行距離は-0.84m減ります

こうやって係数を見ることで他の条件は同じだとして知りたい要因が1つ値が変わるとアウトカムはどう変化するのかを見ることができます。そういう意味ではyの値自体には大して意味がありません。


ここで前回の記事の交絡因子が関係します。

【4-17】「先にデータを取って後で統計」がなぜ危険なのか?どうすればいいのか?を多変量解析、相関・因果関係、交絡を使って解説してみる

因果関係(効果があるかどうか)を考える時はどれだけ重要な交絡因子が取り除けているかが説得力に繋がります。そのためにはどの因子が大事な交絡因子なのかを考える必要がありますが、それはその分野の知識によって決まります。
スクリーンショット 2020-07-17 15.10.39


よって大事なことは右辺には興味がある要因と交絡因子をできるだけ入れることが大切です。

そして知りたい要因は介入が可能であるものが基本となります。
例えば知りたい要因を年齢にしてしまうと、もし有意差が出たとしても年齢を変えることはできません。

なぜここを強調するかというと統計の手法には変数選択という自動的に変数を削除することができるオプションがあるからです。

医療統計でいちばん有名なのはステップワイズ法です。
(下図はEZRの重回帰分析の設定画面)
スクリーンショット 2020-07-17 15.20.46

交絡因子の影響を除きたいのに交絡因子自体を計算から除いてしまうとそれは本当にこの研究でいいたかったことなのか?となってしまいます。

そのためxに興味があるときは変数選択は基本行いません(強制投入といいます)。

*ただ本によっては意見が違ったりします。理学療法の研究ではステップワイズが使われていることも多く、もしかしたら査読で「ステップワイズを」と言われる可能性はあるかもしれません。その雑誌に合わせることも必要なのかもしれませんし、生物統計家のご意見も伺ってみたいところです。

以下の2サイトではステップワイズはよくないという記事になります。







予測に興味がある=yに興味がある。

次にy(左辺)を見るときの視点を紹介します。

上の式を再掲します
1分で歩いた距離(m) = 0.03×片脚立位時間(秒)+ -0.84×年齢(歳) + 0.42(女性であれば) + 108.33

これは片脚立位時間、年齢、性別がわかれば左辺の1分間で歩いた距離を予測することができるということになります。

スクリーンショット 2020-07-17 17.02.35

つまり右辺はどれだけいい予測になる因子を集められるかが大切になります。
そのため交絡因子がどうかや介入が可能かどうかは重要な問題ではありません。
ただ、勘違いしやすいのですが因子が多ければいいモデルになるわけではありません
先程と違い精度に貢献しない因子を外すこともします。

そして予測に興味がある場合に大切なのはモデルを作ることより作ったモデルが新しいデータでも通用するかを試すことで信頼性が増すということです。

これは今回のデータのみに当てはまりがよく、未知のデータでは通用しないという過学習という問題が起こっている可能性があるからです。モデルを作っただけでは過学習が起こっているかどうか判断できません。今回のモデルで使わなかったデータで予測がどれくらい的中するか再検証する必要があります。誤差だったり精度・感度・特異度はこういった再検証に用いられます。



加えて医療統計ではあまりみられませんが、取ったデータから新たな変数を作ることも可能です。

例えばアウトカムを歩行距離として、下肢の左右各関節で疼痛(VAS)を測っていたとします。ただ痛いところは人それぞれですし、変数が多くなりすぎるので変数をまとめたいとします。
  • 各関節の最大のVASを採用する
  • 各関節の平均点を採用する
  • VAS5を超えた関節の数を採用する
  • 他の点数とかけ合わせた新たな変数を作る
他にも色々あるかもしれません。そうやって精度の高い予測を作ることになります。
こういった作業は特徴量エンジニアリングといい、機械学習のテクニックになります。

ただし解釈の難しい因子を作ってしまうと解釈(説明可能性)が難しくなります
上の例だと痛い関節の数はイメージが付きますが、疼痛×筋力みたいな数値はその数値にどういった意味があるのか?と聞かれても「予数値で臨床的な意味は・・・だけど予測の精度が上がる数値」となってしまします。

精度を追求するか説明可能性を追求するかはトレードオフな関係ではありますが、医療現場では説明可能性がある因子を使った方が作った方が使ってくれやすいという印象があります。理由を説明しにくいモデルは「ホントにー?」と胡散臭く思われる可能性があります。

まずは説明しやすいモデル、その後に需要があれば難しいモデルに挑戦するような考えがいいかと思います。


5.因果関係でも予測でも入れないほうがいい変数がある

ただ因果関係・予測どちらにも言えることですが入れないほうがいい因子があります。それは多重共線性が関係します。

数ある因子同士の中で相関係数が高すぎるものが混じっていると多重共線性が発生し結果が不安定にあるという特徴があります。正確な表現ではありませんが、計算の過程で◯÷0みたいな変なことが起こってしまうイメージです。

そのため独立変数同士で相関行列を求め、相関係数が高いものがあるかを調べたりVIFというものを確認したりします。

もし多重共線性が見つかれば相関が高い因子のうちのどちらかを省くか、合成した新しい変数を採用するかになります。

相関行列の求め方は【4-16】Rで相関行列を作成するGGally, cor, clrrplot, corrr関数を紹介しますで紹介しています。



6.まとめ

重回帰分析やロジスティック回帰という回帰分析は計算方法やコードの書き方は同じでも医療統計として使う場面と機械学習として使う場面があります。知りたい要因の因果関係に興味があるのか?アウトカムの予測に興味があるのかで集めるデータや解析や検証の行い方に違いが出てきます。


知りたい要因の因果関係に興味がある→医療統計的な考え方
スクリーンショット 2020-07-17 15.10.39
→データを取り始める前に重要な交絡因子が無いかを十分に検討する
→選択バイアス、情報バイアスもできるだけ入らないよう注意する
→多重共線性がないか確認し、あれば対処する
→分析の結果から考察する


アウトカムの予測に興味がある→機械学習的な考え方
スクリーンショット 2020-07-17 17.02.35
→データを取り始める前に予測に重要な要因は何かを検討する
→多重共線性を調べる
→データを予測作成用のデータと検証用のデータに分ける
→予測作成用のデータでモデルを作る。できたモデルの精度を確認する
→検証用のデータで精度を再検証する
(モデルを作成しただけだと実は本当の精度はわからない)


医療統計ではほとんど前者になることが多い、もしくはどっちつかずになっていることが多い印象ですが、どちらに興味があるかが意識できると多変量解析の勘所が上がるのではないかなと思います。


前回・今回と長くなってしまいましたが、次は重回帰分析の実際について紹介します!


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


今回は重回帰、ロジスティック回帰分析、Cox比例ハザード回帰といった多変量解析を紹介します。


研究を行う時、多変量解析を知っているか知らないかで研究の幅や解釈の仕方が大きく変わります。


個人的にはt検定より先に多変量解析を知っておいたほうがいいと思っています。


実際にRでコードを書く前に、多変量解析でできることとできないことを紹介します。


そしてついやってしまいがちな「取れるだけデータを取って後で統計を考えよう」が危険な理由も紹介します。



1.今回使うデータ

架空のデータを使用します。

library(tidyverse)
set.seed(25) v1 <- rnorm(100,87,13) v2 <- rnorm(100,63,10.5) v <- c(v1,v2) a1 <- as.integer(rnorm(100,35,5)) a2 <- as.integer(rnorm(100,55,5)) age <- c(a1,a2) sex <- c(sample(c(0,1),size = 100, prob = c(0.7,0.3),replace = TRUE),sample(c(1,0),size = 100, prob = c(0.7,0.3),replace = TRUE)) t1 <- as.integer(rnorm(100,120,20)) t2 <- as.integer(rnorm(100,70,20)) time <- c(t1,t2) df <- data.frame(v,time,age,sex) df

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

v:歩行速度、1分間に歩いた距離(m)とします。

time:片脚立位の時間(s)

age:年齢

sex:性別(0:男性,1:女性)

これをみて片脚立位の時間が長くなれば歩行速度が上がると考えていいかどうかとします。
(実際には片脚立位訓練をした群としなかった群でRCTを行うのでしょうが、考え方の練習としてこういった形にしています)


2.多変量解析を知らなければ...

今までの知識で散布図と回帰直線、相関係数を求めてみます。

散布図と回帰直線は
【3-3】Rのggplot2で散布図を作るgeom_point関数で紹介しました。
回帰直線はgeom_smooth(method="lm")を付け加えます。
つなげる時は + です。

ggplot(df,aes(x=time,y=v))+
  geom_point()+
  geom_smooth(method = "lm")

スクリーンショット 2020-07-15 14.02.21
右肩上がりの何やら良さそうなグラフが出てきました。



相関係数と検定を行います。
相関係数に関しては
【4-15】Rで相関係数を評価する方法で紹介しています。

cor.test(df$v,df$time)

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

相関係数0.48と中等度の相関が出ました。p値も2.4×10の-13乗と有意差もあります。

ではこれで
片脚立位の時間が長くなれば歩行速度が上がると言えるでしょうか?

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

実は言えません。


3.別の視点で分析してみる

ではなぜ言えないのでしょうか?今度は違う視点で分析してみます。
先程の散布図に年齢別で色をつけてみます。

ggplot(df,aes(x=time,y=v))+
  geom_point(aes(col=age))

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

色が薄い方が年齢が高くなっています。
こうみると解釈が変わってきます。

今回は年齢が若くて片脚立位時間も歩行速度もいい人と、高齢で片脚立位時間も歩行速度も悪い人がいる集団なだけだった、年齢に関係なく片脚立位が長いから歩行速度が上がった人は実際にはどれくらいいるだろうか・・・?と気になってしまいます。


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

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




上記のどちらが妥当化はこのままでは判断できません。

性別に関しても同様です。

ggplot(df,aes(x=time,y=v))+
  geom_point(aes(col=as.factor(sex)))

スクリーンショット 2020-07-15 14.44.10
なんとなく男性の方が歩行速度が早いように見えます。


実は次回以降の記事で行う重回帰分析を行うと年齢と性別の影響を除くと片脚立位は有意差が出なかったという結果になります。
コードと解釈の仕方は次の記事で紹介します。

fit <- lm(v ~ time + age + sex)
summary(fit)

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


3.相関関係と因果関係の違い

(ダミーデータですが)今回の結果では片脚立位と歩行速度には相関関係はあったが、因果関係はなかったという解釈をしました。

ここで出てくるのが相関関係因果関係です。

スクリーンショット 2020-07-15 15.12.44
相関関係は原因は別としてある数値が変化するともう一方の数値が変化するという関係性を示しています。

相関があったは「効果があった」とは別物です。効果があるかどうかは何とも言えません。

先程の例でいうと
片脚立位と歩行速度には相関関係があったが、原因は片脚立位と別なところにあった(年齢など)と解釈することができます。

また
片脚立位の時間が長かった人は歩行速度が早い人が多かった歩行速度が早かった人は片脚立位の時間が長い人が多かったと表現しても違和感がありません。





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

いっぽう、因果関係は原因と結果の関係になっていて効果があったと言えるものです。
そのためAの数値を変えることでBの数値を変えることが期待できます。

方向性は一方向となります。
相関関係と違い年を取ればと歩行速度は下がりやすくなる速く歩けない人はを取りやすくなる(文字通り11ヶ月で1歳年を取るようになるという意味)では意味合いが全然違います。因果関係とはそういうことです。

そしてB→Aの因果関係はありませんが、相関があるように見えます。これを見せかけの相関といいます。

見せかけの相関が起こる例は交絡、因果の向きが逆である時、全くの偶然です。


交絡バイアス

交絡とは原因と結果両方に影響を与える因子のことです。
今回は年齢や性別が交絡因子にあたります。
交絡因子が残ったままだと因果効果(本当に効果があるのかどうか)ははっきりと言えず、そこが研究の限界の1つとなります。そのため、交絡は補正することがひつようで、重回帰分析やロジスティック回帰分析といった多変量解析は交絡を補正する方法の1つとなります。

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

全くの偶然

全くの偶然で相関関係が起きる場合があります。

海賊の数が減るにつれて、同時に地球温暖化が大きな問題となってきた。

wikipediaより



相関関係か因果関係かを見極めるポイントに関しては原因と結果の経済学に詳しい解説があります。



そして因果関係を証明するにはバイアスをどれだけ排除できるかがカギになります。




4.重回帰・ロジスティック回帰・Cox比例ハザード回帰は交絡を調整する

では重回帰分析やロジスティック回帰、Cox比例ハザード回帰は何をしてくれるのでしょうか?

これらの多変量解析は
交絡の調整を行ってくれます。

スクリーンショット 2020-07-15 21.14.16
重回帰分析で片脚立位と年齢を変数として入れると
年齢の影響を除いた上での片脚立位が歩行速度に与える影響を推定することができます。性別も一緒に入れると年齢と性別の影響を除いた上で影響を推定します。ただ他にも交絡因子はあります。例えば下肢筋力などです。

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


因果関係(効果があるかどうか)を考える時はどれだけ重要な交絡因子が取り除けているかが説得力に繋がります。そのためにはどの因子が大事な交絡因子なのかを考える必要がありますが、それはその分野の知識によって決まります。決してデータが集まれば自動的に決まるものではありません。そのためデータを取り始める前に交絡が何かを決めておく必要があります。データを取り終えたあとで「◯◯忘れてた!」となると取り返しがつきません。

また交絡因子の数で必要な症例数が決まります(1つの交絡因子で約10症例追加)。そのため症例数が50人であれば5つの変数しか投入できないため、何例集めるのか?5つの変数をどれにするのか(重要なものから5つを選ぶ)を先に決めておきます。

相手に
確かにこの影響まで除いた上で有意差が出ているのなら確かに効果があるかもしれないな。と感じさせることができるかどうかではないでしょうか。よく学会でt検定や相関係数で効果があることが示唆されたという発表もあるのですが、「でも◯◯の影響があるんじゃない?」と思われていたら多変量解析(もしくはRCT)が必要だったのかもしれません。



5.交絡因子と中間因子

ということで交絡因子は除外する必要があります。しかし中間因子を取り除いてはいけないという注意点があります。

スクリーンショット 2020-07-16 6.52.09


中間因子は矢印の方向が交絡因子と違います。

興味のある項目(A)とアウトカム(B)の間にある因子のことです。

中間因子を多変量解析に加えてしまうとAの効果があるかどうかわからなくなるため注意が必要です。

ただ交絡因子と中間因子どっち?は統計知識ではなくその分野の知識で決まることなのでデータを取り始める前に予め決めておくことが必要です。


6.交絡以外のバイアスは取り除けない

多変量解析で取り除くことができる交絡バイアスですが、実はバイアスは交絡以外にも複数あります。

細かくは紹介しませんがここではいくつか紹介します。

選択バイアス

・日中の電話調査は電話に出ることができる一部の高齢者しか集まらないためデータに偏りが出た
・自宅退院した人を1年間かけて調査すると退院時より点数の平均が上がっていた。しかし実は死亡・再入院で悪い点数の人が減ったことが原因で平均点が上がったように見えただけであった。
・(実は効果がある研究だったが)自分の病院は重症な人が多く集まる傾向で結果が出なかった。複数の病院でいろいろな患者さんを集めてみたら結果が違ったかもしれない(逆もまた然り)。
・集めたデータに欠損値が多かった。欠損値があるデータは重症な人が多かった。



情報バイアス

・評価者が患者さんが「治療郡」「対照群」のどちらに振り分けられるか知っていて、評価時に「さぁ!評価しますよ!」だったり「あんまり無理しなくていいですからね」って無意識に言ってしまう。
・評価器具の使い方に上手・下手が出て評価者間で違ってばらつきが出てしまう。
・評価基準がスタッフによって微妙に違いばらつきが出る


こういったバイアスは重回帰分析やロジスティック回帰分析ではどうにもできません。
これらのバイアスには統計ではなく研究デザインで対応することになります。
ただ欠損値に関しては条件が整えば多重代入法という統計のテクニックで解決できる可能性があります。
また施設間での研究では各評価だけでなく、施設ごとの特徴を考慮したマルチレベル分析といったものもあります。


7.「評価を取れるだけ取って後で考えようとデータを取り始める」ことをやってはいけない理由

ここまで説明するとなぜ「評価を取れるだけ取って後で考えよう」とデータを取り始めると危ないかが説明できます。どれだけいろいろな種類・数のデータを集めても、大事な交絡が抜けていたら元も子もないからです。




8.まとめ:研究の初心者でもできること

ただこういった話を続けていると難しすぎて研究が全く始められません。
臨床で疑問に思っていることを形にすることはとても意義のあることです。
個人的には最初から完璧を目指す必要はないと思います。

ただ興味のある因子に対して以下のことに意識を向けておくだけでも最初の研究としては一歩を踏み出せるのではないかなと思います。

・交絡因子は何がありそうか?
・その中で計測できそうなものはどれか?
・なるべく欠損なくとるよう努力する
・必要なサンプルサイズ(症例数)を先に知っておく

ただこのサイトでは研究デザインに関しては解説できていません。
多変量解析を行う上ではSTROBE声明を確認するともっと先にすすめると思います。



次回も重回帰・ロジスティック回帰を使う前に知っておきたいこととして「洞察をしたいのか?予測をしたいのか?」という話をします。リハビリの世界ではあまり使い分けられれている印象がありませんが、知っておくことで結果の解釈の幅が大きく変わります。

この話は機械学習にも通じる話となりますので、機械学習に興味がある方はぜひ読んでみてください。

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


ただ複数の列がある場合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.まとめ

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

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

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

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

【1-3】Rのソフトはどれを使えばいいか?目的別チャートを作りましたではRを使うためのソフトを紹介しました。


その中で統計を勉強してみたいという方はRstudio + EZR(Rコマンダー)の組み合わせをおすすめしました。

Rを使っている方は「EZR(Rコマンダー)しか使っていない」や「Rstudioしか使っていない」という方が多いのですが、実はRstudioからRコマンダーやEZRを簡単に呼び出すことができます。

今回はRstudioからEZRを呼び出す方法を紹介します。


1.EZRのインストールをRstudioだけで行う。

EZRのダウンロード、インストールは自治医科大学附属さいたま医療センターのページで紹介されていますが、Rstudioからでも行えます。

①RstudioのInstallボタンを押す
Rstudioの真ん中あたりにInstallボタンがあります。これはRの様々なパッケージをプログラムを書かずにダウンロード、インストールすることができます。
スクリーンショット 2020-04-24 19.34.48

②RコマンダーとEZRをインストール
次にRコマンダーとEZRをダウンロードします。

Packagesという欄にrcmdrと入力します。
そうすると候補にRcmdrとRcmdrPlugin.EZRという2つが見つかります。
スクリーンショット 2020-04-24 19.34.58

まずはRcmdrをインストールし、次にRcmdrPlugin.EZRをインストールします。

この作業は今回やってしまえば次回以降は行う必要はありません。


2.EZRを呼び出す

EZRを呼び出すには2つ手順が必要ですが慣れるとすぐにできるようになります。

①Rコマンダーを呼び出す(2種類の方法があるので好きな方でいいです)
②RコマンダーからEZRを呼び出す


①Rコマンダーを呼び出す

1つめの方法はPackagesからRcmdrをクリックして呼び出す方法です。
スクリーンショット 2020-04-24 19.35.03
2つめの方法は下のコードを入力する方法です
library(Rcmdr)

慣れるとコードを書いた方が早いですが、最初はクリックする方法の方が早いかもしれません。
(自分がそうでした)


②RコマンダーからEZRを呼び出す

EZRはRコマンダーから呼び出します。

ツール→Rcmdrのプラグインをロード…を選択
スクリーンショット 2020-04-24 19.35.13


RcmdrPlguin.EZRが選択されていることを確認してOK
スクリーンショット 2020-04-24 19.35.21

OKを押して再起動
スクリーンショット 2020-04-24 19.35.28

これで完了です。
見た目は似ていますが「標準メニュー」となっていればEZRになっています。
スクリーンショット 2020-04-24 19.35.36



3.注意点
・RstudioのPackageにRcmdrPlugin.EZRがあるのですが、クリックしても何も反応しません。上記の手順が必要です。

・EZRを使うと結果はRstudioのコンソール画面に出てきます

・1度開いたEZRを×で閉じてしまうと、同じ方法を行っても開かないことがあります。一旦Rstudioを再起動するかプロジェクトを開き直すとまたできるようになります。

・グラフはRstudioのPlotではなく別ウィンドウで開くことがあります。

この辺りはパソコンのOSでも挙動が違うかもしれません。私はmacを使っていますが、もしWindowsで違う挙動があればこの記事でも紹介しますので連絡いただけると嬉しいです。


4.まとめ
今回はRstudioとEZRを同時使いするための方法を紹介しました。

Rstudioだとどう使えばいいかわからない。でもいずれ使えるようになってみたい。という方は同時使いをがおすすめです。EZRの手軽さ、それがどのようなコードで書かれているのかを同時に体験することができると思います。

またEZRの使い方はハルさんのシロート統計学講座がおすすめです!
無料ですが、これだけでEZRの使い方はほとんど理解できると思います。






↑このページのトップヘ