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


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


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


 


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


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


今回は対応のあるt検定を紹介します



まず対応のあるt検定についてはハルさんのサイトをご参照ください。



また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-03 8.11.46


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

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

デモデータ(対応のあるt検定)


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


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


実際のコードは以下になります。
前回のコードのURL(" "の中)とdestfileのdata/以降を変更するだけでOKです。
url <- "https://haru-reha.com/wp-content/uploads/2018/04/demo-paired-t-test.xlsx"
destfile = "data/demo-paired-t-test.xlsx"

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



dataフォルダにダウンロードできました!
スクリーンショット 2019-11-03 18.54.56



4.データの読み込み

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

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

スクリーンショット 2019-11-03 18.53.56


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

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

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

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

スクリーンショット 2019-11-03 19.01.56


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

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


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


5.正規性の確認

正規性の確認は【4-1】Rでt検定を行う方法でも紹介しました

今回ハルさんのサイトではdifferenceの列の正規性をQQプロットを使って確認しています。

QQプロットの出し方はいくつかありますが、ここでは2つ紹介します。

まずはqqline関数です。
qqline(y = gait$difference)
スクリーンショット 2019-11-03 19.36.51
もしEZRと同じグラフが出したい場合はcarパッケージqqPlot関数を使います(Pが大文字なので注意!)
Packagesにcarパッケージが入っていれなければinstall.packages関数でインストール後library関数で呼び出します。

ちなみに先程のqqline関数y = でしたが、qqPlot関数x = となっています。 

install.packages("car")

library(car) qqPlot(x = gait$difference)
スクリーンショット 2019-11-03 19.43.26

こうすることでEZRを開かなくても同様のことが行なえます。
確認するのが目的であればEZRでいいと思いますが、1年後に再現しようと思ったらスクリプトに残しておくと再現性が高まります。

ヒストグラムはgeom_histgramもしくはhist関数を使います。
library(tidyverse)
ggplot(data = gait)+
  geom_histogram(aes(x = difference), bins = 5)

hist(gait$difference)
スクリーンショット 2019-11-03 19.52.12

シャピロウィルク検定は【4-1】Rでt検定を行う方法でも紹介しましたがshapiro.testでした。
shapiro.test(gait$difference)
スクリーンショット 2019-11-03 19.55.47


6.対応のあるt検定を行う

対応のあるt検定は実は【4-1】Rでt検定を行う方法でも紹介したt.test関数paired = TRUEを足すだけです。ちなみに対応のないt検定はpaired = FALSEをしていることになります。何も書かなければpaired = FALSEになります。

t検定は2種類の書き方があります。その違いはデータの形にあります。
ハルさんのサイトでEZRでは右の形しかできませんでしたがRならどちらも可能です
ちなみに左をlongデータ(縦に長い)、右をwideデータ(横に広い)といいます。
スクリーンショット 2019-11-03 20.03.44

下の図【4-1】Rでt検定を行う方法で紹介しましたが、これにpaired = TRUEを付け加えます。ただlongデータで行う場合はAの並び順とBの並び順が違うと計算結果が異なりますので注意が必要です。
スクリーンショット 2019-10-27 8.09.41

今回はwideデータなので以下のようになります。
var = は外しています。

t.test(gait$pre, gait$post, paired = TRUE)

スクリーンショット 2019-11-03 20.12.08

EZRと同じ結果になりました。

更に95%信頼区間は41.46340 - 86.86993となっており、差の平均(期待値)は64.16667となっています。

どういうことかと言うと、対応のあるt検定は2つの差(ここではpost - pre)が0であると仮定しています。そして95%信頼区間に0が含まれているとpは0.05以上と判断されます。
例えば差の平均(期待値)が20でも95%信頼区間が-30 〜50といった具合です。

スクリーンショット 2019-11-03 20.37.08
上のsampleの場合、差の平均(期待値)は20とプラスなのでpostの方がより高い値と言えそうですが、95%信頼区間に0やマイナスが含まれています。言い換えると20だったけど-10にだってなり得るし、0にだってなり得るとなります。そうなればpostの方が高い値といえません。

ちなみに今回は差(期待値)が0を仮定していたので0でしたが、例えばロジスティック回帰分析だとodd比が1であると仮定するので、そのときは1を挟むかどうかを確認することになります。

p値だけでなく信頼区間を確認するとまた発見があるかもしれません。


7.まとめ

今回は対応のあるt検定を行いました。
実際には対応のないt検定にpaired = TRUEを付け加えるだけでした。

4章を順に見ていくと重複する箇所も出てきますので検索で来られた方はサイトマップを見ていただければ別の発見があるかもしれません。




8.今回使ったコード

今回使ったコードをまとめて置いておきます。
95%信頼区間のコードも載せています。


#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/04/demo-paired-t-test.xlsx"
destfile = "data/demo-paired-t-test.xlsx"

download.file(url, destfile)

#データの読み込み
library(readxl)
gait <- read_excel("data/demo-paired-t-test.xlsx", 
                   range = "B2:D32")
View(gait)

#正規性の検定

#qqplot
qqline(y = gait$difference)

#carパッケージの
install.packages("car")
library(car)
qqPlot(x = gait$difference)

#ヒストグラム
library(tidyverse)
ggplot(data = gait)+
  geom_histogram(aes(x = difference), bins = 5)

hist(gait$difference)

#シャピロ・ウィルク検定
shapiro.test(gait$difference)

#t検定
t.test(gait$pre, gait$post, paired = TRUE)

#信頼区間
gait_ttest <- t.test(gait$pre, gait$post, paired = TRUE)

conf <- data.frame(conf.low = c(-30, round(gait_ttest$conf.int[[1]], 3)), conf.high = c(50,round(gait_ttest$conf.int[[2]],3)), mean = c(20,round(gait_ttest$estimate, 3)), x = c("sample", "post - pre"))
ggplot(data = conf) +
  geom_errorbar(aes(x = x, ymin = conf.low, ymax = conf.high), width = 0.1) +
  geom_hline(yintercept = 0, color = "red") +
  geom_text(aes(label = conf.low, x = x, y = conf.low), vjust = -1) +
  geom_text(aes(label = conf.high, x = x, y = conf.high), vjust = -1) +
  geom_text(aes(label = mean, x = x, y = mean), vjust = -1) +
  geom_point(aes(x = x, y = mean)) +
  labs(x = "", y = "")+
  coord_flip()