タグ:散布図

第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.まとめ

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

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

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

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

前回、Rのggplot2を使ってグラフを作成する基本を紹介しました。

今回は散布図を作り方を紹介します。
併せてggplot2を使うのに必要なaes関数にも慣れていきたいと思います。

1.データの読み込み

今回もggplotパッケージが含まれているtidyverseパッケージを読み込みます。
#tidyverseパッケージを読み込みます。
if(!require(tidyverse)) install.packages("tidyverse", repos = "http://cran.us.r-project.org")

#既にtidyverseパッケージをインストールしている方は以下でもOK
library(tidyverse)

#データ取り込みます。今回はdatという変数にデータを入れます
url <- "https://github.com/mitti1210/myblog/blob/master/data01.csv?raw=true"
dat <- read.csv(url)


データを確認します。最初の数行を確認するのはhead関数を使います
head(dat)

歩行(0:非自立、1:自立)と氏名、性別、年齢、身長、体重、SIAS(脳卒中の評価)、BBS(バランスの評価)、MMSE(認知機能の評価)が入っています(架空のデータです)。




2.散布図の基本的な作り方

散布図を作るにはgeom_point関数を使います。

スクリーンショット 2019-05-23 0.07.49


今回はp1という変数名に横軸に身長、縦軸に体重の散布図を作成します。
またtheme〜というのはMacで日本語が□□と文字化けしないようにフォントを指定しています。
Windowsではなくて大丈夫(なはず)です。

p1 <- 
  dat %>%
  ggplot()+
    theme_gray(base_family = "HiraKakuPro-W3")+
    geom_point(aes(x = 身長, y = 体重))
p1
もしくは
p1 <- 
  ggplot()+
    theme_gray(base_family = "HiraKakuPro-W3")+
    geom_point(aes(x = 身長, y = 体重), data = dat)
p1

スクリーンショット 2019-05-23 0.30.30


基本はこれで完成です。ただもっと形を変えてみます。

点の色を赤にしてみる

色を変えるにはcolor =○○を使います。
主要な色であれば"red"や"blue"など色名を""で囲みます。
"#00BFC4"などRGBで指定することもできます。

p2 <- 
  ggplot() +
  theme_gray(base_family = "HiraKakuPro-W3")+ 
  geom_point(aes(x = 身長, y = 体重), color = "red",  data = dat)
p2

スクリーンショット 2019-05-23 0.32.14


点の色を性別で分けてみる

男女別に色分けもできます。
p3 <-
  ggplot() +
  theme_gray(base_family = "HiraKakuPro-W3")+ 
  geom_point(aes(x = 身長, y = 体重, color = 性別), data = dat)
p3

スクリーンショット 2019-05-23 0.33.39


p2とp3の違いがわかるでしょうか?

スクリーンショット 2019-05-23 0.43.27


ここでaes関数に関してもう少し詳しく説明します。


aes関数について

aes関数はデータフレームの変数名グラフに指定する場合に使います
x軸やy軸はイメージしやすいですが色・形・サイズなどに関しては①変数名の値を使いたければaes関数に入れる②変数名を使わず自分で値を指定したい場合はaes関数に入れないとなります。

スクリーンショット 2019-05-23 0.53.37


これを応用するともっとグラフを変えることができます。


点の大きさを変える。形も変える。

点の大きさはsizeで、形はshapeで指定します。

shapeは番号によって形が選べます。

スクリーンショット 2019-05-22 20.29.48


p4 <-
  ggplot() +
  theme_gray(base_family = "HiraKakuPro-W3")+ 
  geom_point(aes(x = 身長, y = 体重, color = 性別), size = 3, shape = 3, data = dat)
p4

スクリーンショット 2019-05-23 0.57.45


形を歩行の値で変更する

shapeをaes関数に入れてみます。

p5 <-
  ggplot() +
  theme_gray(base_family = "HiraKakuPro-W3")+ 
  geom_point(aes(x = 身長, y = 体重, color = 性別, shape = 歩行), size = 3, data = dat)
p5

エラー: A continuous variable can not be mapped to shape

あれっ?エラーが出ました。連続変数はshapeには使えませんよと怒られました

これは歩行が0,1で今のところ連続変数だからです。
確かに上のshapeの番号で1.5といった番号はエラーが出そうです。
なのでこのグラフの中だけ歩行をカテゴリー変数に変えます。factor関数を使います。

p5 <-
  ggplot() +
  theme_gray(base_family = "HiraKakuPro-W3")+ 
  geom_point(aes(x = 身長, y = 体重, color = 性別, shape = factor(歩行)), size = 3, data = dat)
p5

スクリーンショット 2019-05-23 0.59.20


これでエラーが消えました。ただ凡例がfactor(歩行)なんてなってしまいました。

凡例の名前だけでなく、軸やタイトルの名前を変えてみます。
名前を変えるのはlabs関数です。

p6 <-
  ggplot() +
  theme_gray(base_family = "HiraKakuPro-W3")+ 
  geom_point(aes(x = 身長, y = 体重, color = 性別, shape = factor(歩行)), size = 3, data = dat)+
  labs(x = "身長", y = "体重", color = "性別", shape = "歩行", title = "散布図", caption = "これは仮想のデータです")
p6
スクリーンショット 2019-05-23 1.01.39
これで変わりました。


横軸がカテゴリー変数の場合

横軸がカテゴリー変数でも散布図は使えます。

p7 <- 
  ggplot()+
  theme_gray(base_family = "HiraKakuPro-W3")+ 
  geom_point(aes(x = 性別, y = 身長), data = dat)
p7

スクリーンショット 2019-05-23 1.08.50


ただ重なってよくわかりません。このように縦軸や横軸にカテゴリー変数を使う場合は点の位置をずらすとわかりやすくなります。

点を重ならないようにずらすのは2つの方法があります。

①geom_point関数の中でposition = "jitter"を追加する

p8 <- 
  ggplot()+
  theme_gray(base_family = "HiraKakuPro-W3")+ 
  geom_point(aes(x = 性別, y = 身長), position = "jitter", data = dat)
p8

スクリーンショット 2019-05-23 1.09.56


②geom_jitter関数を使う
p9 <- 
  ggplot()+
  theme_gray(base_family = "HiraKakuPro-W3")+ 
  geom_jitter(aes(x = 性別, y = 身長, color = 性別), data = dat)
p9

スクリーンショット 2019-05-23 1.10.04

これで男女の分布の違いが見やすくなりました。



まとめ

今回は散布図について紹介しました。
Rを使った散布図はxとyの値を変えるだけで作れますし、colorやshapeを使うことで更に情報を付け加えることができます。

今回のデータセットは他にも変数があるのでぜひ他の散布図も作成してみてください。

↑このページのトップヘ