疫学統計セミナー

(1)

疫学統計セミナー

疫学と統計の基礎からロジスティック回帰

担当：　茅野光範

グローバルアグロメディシン研究センター獣医学研究部門

メール

: kayano@、内線：5521

第 3 回：統計ソフト R の基礎と応用

H29.1.24

1

セミナー資料：　

h3p://www.obihiro.ac.jp/~kayano/epi-stat/

(2)

このセミナーについて

内容：　疫学と統計を復習し、交絡因子とその調整方法、

　　　　　ロジスティック回帰等を紹介する　　

目標：　交絡因子調整の検定やロジスティック回帰を理解し、

　　　　　

R

等で実行できるようになる！

ポイント：疾病の規定要因（リスク因子）を正しく同定する

日時（予定）：　毎月下旬月曜

or

火曜の午後

5

時から

1.5

時間程度

スケジュール（予定）：　全４回

　　第

1

回

(11/28)

：　疫学と統計の基礎

　　第

2

回

(12/20)

：　交絡因子とその調整方法

　　第

3

回

(1/24)

：　統計ソフト

R

の基礎と応用

　　第

4

回

(2/21?)

：　ロジスティック回帰（仮）＋

α

？

2

(3)

1 回目の目標と内容

目標：コホート研究（追跡、

dynamic/ﬁxed

）と症例対照研究（

case/

control

）において、暴露が疾病に関与しているかどうかを

検証（検定）する。　（ただし、交絡は無視する）

内容：

• 

はじめに

　　　疫学とは何か、有名な疫学研究、トピック、リスク因子の同定

• 

研究方法（研究デザイン）と疾病のタイミング　　　　コホート研究（

follow-up

研究）、症例対照研究

• 

疫学で用いられる指標と統計的推測

　　　罹患率（

incidence raNo

）、有病率（

prevalence

）　　　リスク比、オッズ比、カイ二乗検定、信頼区間

3

復習

(4)

2 回目の内容：交絡因子とその調整方法

• 

復習

　　　目標・内容、　リスク因子の同定と交絡因子の影響　　　研究方法に応じたデータレイアウト

　　　暴露効果の指標と推測

• 

交絡の調整方法

1

　研究方法で調整する　（考え方の紹介）

　　　

(1)

因子範囲の制限　　　　

(2)

マッチング

• 

交絡の調整方法

2

　解析方法で調整する

　　　

(1)

層化解析（

StraNﬁcaNon

）：マンテル・ヘンツェル検定　　　

(2)

回帰分析　

←

　第４回目にやります

4

（交絡の影響は無視）

復習

(5)

今日の内容：統計ソフト R の基礎と応用

導入

•  R

の紹介

R

の基礎

•  R

のメイン画面と電卓としての利用

•  R

における「変数」と「代入」

•  R

における保存と読み込み

• 

記述統計学（平均や分散などを求める、図を描く）

R

の応用

• 

推測統計学：検定をする

　　　　　　　　（

t

検定、カイ二乗検定、フィッシャーの正確確率検定、

　　　　　　　　　マンテルヘンツェル検定）

5

(6)

R の紹介

•  R とは？

•  R の関連ソフト

•  R の資料

6

(7)

R とは？

• 

フリーの統計解析ソフト

利点

• 

複雑な統計の計算を、たった数行のコマンドで実行できる！

• 

統計で使う様々な関数が組み込まれている（パッケージが豊富！）

• 

パッケージの追加も簡単（フリー！）

　　最新の論文の方法をコマンド１つで実行出来ることもある

• 

グラフィックスが豊富

• 

反復計算、行列計算が出来る

欠点

• 

計算が速いとはいえない（が、工夫してそこそこ速くすることも出来る）

• 

データが大きすぎる場合にはフリーズすることもある　　（

n>1,000,000

など）

7

(8)

紹介：　 R の関連ソフト・アドイン

•  R Studio

：　

R

を見やすく、操作しやすくしたもの

•  R Commander

•  EzR

•  Excel

で

R

：　（

Excel

のバージョンが限られている？）

他のプログラミング言語との互換

•  Rcpp

：

R

と

C++

言語（高速！）と橋渡し

　　　　　　時間がかかる計算を

C++

で実行して、

　　　　　　結果を

R

に入れることも可能

　　　　　　例：　

R

でやると

5

時間かかる計算が、

　　　　　　　　　

C++

でやると

3

分で終わることもある、、！

8

（

100

万回の検定等）

(9)

R の資料

一般

•  R

の本：

　　『

The R Tips

』（第３版）　

舟尾暢男 (編集)

　　『統計処理ポケットリファレンス』（

Excel

＆

R

対応）：　涌井良幸・涌井貞美　　　『リトルブック：

R

による生物医学統計』

　　

h3p://www.ec.kansai-u.ac.jp/user/arakit/documents/lbrbs.pdf

•  R Tips h3p://cse.naro.aﬀrc.go.jp/takezawa/r-Nps/r.html

• 

統計ソフトウェア

R

についての

Tips

　　　

h3p://minato.sip21c.org/swNps/R.html

•  Google

で検索　（結構使える）　例：

”R

　

t

検定

”

で検索学内

• 

増田先生のページ：

R

の紹介や基本操作、検定など

h3p://www.obihiro.ac.jp/~masuday/resources/stat/index.html

• 

茅野のページ：バイオインフォマティクスのセミナー資料（

2012

年５月、６月）　　

h3p://www.obihiro.ac.jp/~kayano/bioinfo/

　

9

(10)

統計ソフト R の基礎

1.   R のメイン画面と電卓としての利用 2.  R における「変数」と「代入」

3.   データの保存・読み込み

4.   記述統計学 1: 平均や分散を求める 5.   記述統計学 2: 作図

10

2-1

　ベクトルを扱う

2-2

　行列を扱う

2-3

　配列を扱う（後述）

(11)

1. R のメイン画面と電卓としての利用

11

ここに命令文（コマンド）を入力する

R は電卓として使える

などなど

> 1+2*3

　 ⏎

[1] 7

> sin(pi/2)

^⏎

[1] 1

> 2/3

　　 ⏎

[1] 0.6666667

> 10^3

　　 ^⏎

[1] 1000

足し算とかけ算割り算

累乗

平方根

三角関数

> sqrt(2)

　　 ⏎

[1] 1.414214

コンソール画面

(12)

いろいろな演算

演算　　　演算子　例加算　　　　　

+

　　　

a+b

減算　　　　　　

-

　　　

a-b

乗算　　　　　　* 　　　

a*b

除算　　　　　

/

　　　

a/b

累乗　　　　　

^

　　　

a^b

剰余　　　　

%%

　　　

a%%b

平方根　　　　

sqrt

　　　

sqrt(a)

対数　　　　　　

log

　　　

log(a)

常用対数　　

log10

　　

log10(10) e

のベキ乗　

exp()

　　

exp(a)

サイン　　　　　

sin()

　　

sin(pi/2)

コサイン　　　

cos()

　　

cos(0)

タンジェント　　

tan()

　

tan(0)

演算　演算子　例より大きい　

>

　

a>b

より小さい　

<

　　

a

＜

b

以上　　　

>=

　

a>=b

以下　　　

<=

　　

a<=b

等しい　　　

==

　　

a==b

等しくない　

!= a!=b

12

(13)

2. R における「変数」と「代入」

変数：　数値やベクトル、行列（データ行列）、配列に名前をつけて、

R

の中に保存したもの。

R

では“オブジェクト”と呼ばれる。

13

例：変数

x

に

2

を、変数

y

に

3

を代入し、

　　

x

と

y

の中身を確認する。

> x<-2

> y<-3

> x [1] 2

> y [1] 3

> x=2

> y=3

> x [1] 2

> y [1] 3

代入の演算子：　

“<-” or “=“

> (x<-2) [1] 2

> (y<-3) [1] 3

方法

1

方法

2

方法

3

> assign(“x”,2)

　

> assign(“y”,3)

方法

4:

変数の名前：

複数文字でも数字を入れてもいいただし、、

• 

大文字と小文字は区別する　　

“x”

と

”X”

は違う

• 

先頭に数字を入れない　　

“x1”

は

OK

、

”1x”

はダメ

• 

以下の名前はつけられない。既に　　使われていて、かつ、上書き出来ない

Inf, NA, NaN, NULL, TRUE, FALSE, for, in, if, else, break, while, next, repeat, funcNon

• 

上書き可能な変数もある

pi

など

(14)

2-1 　ベクトルを扱う　

> a <- c(1,2,3,4,5)

> a

[1] 1 2 3 4 5

a

に

(1, 2, 3, 4, 5)

を代入し、確認する

> a <- 1:5

> a

[1] 1 2 3 4 5

> a <- seq(1,5)

> a

[1] 1 2 3 4 5

> a <- c(1,2,3,4,5)

> b <- c(2,3,4,5,6)

> a+b

[1] 3 5 7 9 11

a

に

(1, 2, 3, 4, 5)

を、

b

に

(2, 3, 4, 5, 6)

を

代入し、足し合わせる

a

の要素を確認する・置き換える

> a[1]

[1] 1

> b[2]

[1] 3

> a[1:3]

[1] 1 2 3

> a[1] <- 6

> a

[1] 6 2 3 4 5

> a[1:2] <- -1

> a

[1] -1 -1 3 4 5

> a[1:2] <- 6:7

> a

[1] 6 7 3 4 5

14

> rep(1,5) [1] 1 1 1 1 1

同じ要素なら、

rep

が使える

(15)

2-2 　行列を扱う

15

> x<-matrix(c(1,2,3,4,5,6),nrow=2,ncol=3)

> x

[,1] [,2] [,3]

[1,] 1 3 5 [2,] 2 4 6

> x<-matrix(c(1,2,3,4,5,6),nrow=2,ncol=3,byrow=TRUE)

> x

[,1] [,2] [,3]

[1,] 1 2 3 [2,] 4 5 6

行数列数要素

nrow

は

nr

、

ncol

は

nc

でもよい

byrow

は、

by

でもよい

TRUE

は、

T

でもよい

行

列

> dim(x)

[1] 2 3 > nrow(x)

[1] 2 > ncol(x) [1] 3

行列の大きさ

の確認

(16)

行列の操作

16

> x

[,1] [,2] [,3]

[1,] 1 2 3 [2,] 4 5 6

> x[1,2]

[1] 2

> x[1,]

[1] 1 2 3

> x[,2]

[1] 2 5

1

行目を

取り出す

2

列目を

取り出す

(1,2)

成分を取り出す

> y<-matrix(c(2,3,4,5,6,7),nrow=2, ncol=3, byrow=TRUE)

> y

[,1] [,2] [,3]

[1,] 2 3 4 [2,] 5 6 7

> x+y

[,1] [,2] [,3]

[1,] 3 5 7 [2,] 9 11 13

> x*y

[,1] [,2] [,3]

[1,] 2 6 12 [2,] 20 30 42

x

と

y

を足す

x

と

y

を要素毎に掛ける

「

x/y

」で要素毎に割れる行列の積を表す演算子は、

「

%*%

」

(17)

3. R における保存と読み込み

保存

(1) R

の作業ファイルを丸ごと保存する（

.Rdata

ファイル）

　　

[

メニュー

]

　⇒

[

ワークスペース

]

　⇒

[

ワークスペースをファイルに保存

]

(2)

コマンド履歴のみを保存する（

.txt

ファイル）

　　　

[

メニュー

]

　⇒

[

ファイル

]

　⇒

[

保存

]

(3)

変数を保存する（

.txt, .csv

ファイル等）

　

write.table(x, “output.txt”, row.names=FALSE, quote=FALSE, append=FALSE) “write.csv()”

や

”write()”

もある

17

次回、途中から作業を続けられる

次回、コピペで

同じ作業を実行出来る

計算結果の保存等

(18)

3 (3) 　変数の保存

> write.table(x, “output.txt”)

18

変数

x

を

output.txt

という名前のファイルに保存する。

R

が今参照しているフォルダに、

”output.txt”

というファイルが出来ている。

> getwd()

[1] "/Users/kayano"

R

の参照フォルダを確認する

> setwd(“/Users/kayano/workshop”) ^R

_変更する ^{の参照フォルダを}

"V1" "V2" "V3"

"1" 1 2 3

"2" 4 5 6

行名

(row.names

）はナシ、データはダブルコーテーション（

quote

）で囲まない、

既存ファイルへの追加の書き込み（

append

）はナシとする。

> write.table(x, “output.txt”, row.names=FALSE, quote=FALSE, append=FALSE)

(19)

ファイルの読み込み

19

> data<-read.table(“data1.txt”, header=T)

> data

Sex Height Weight 1 M 168 60 2 M 176 67 3  M 168 58

・

29 F 168 50

Sex Height Weight M 168 60

M 176 67 M 168 58 M 165 61

・

Excel

からコピペ（

Windows)

“sample_data.txt”

> data<-read.delim(“clipboard”, header=T)

列名をコピーした場合

Excel

からコピペ（

Mac)

> data<-read.delim(pipe(“pbpaste”), header=F)

^{列名をコピー} ^{しなかった場合}

Excel

で「データを選択」⇒「コピー」の後に、

以下を行う

(20)

4. 平均や分散などを求める

20

> a <- c(1,2,3,4,5)

> a

[1] 1 2 3 4 5

> mean(a) [1] 3

> var(a) [1] 2.5

> sd(a)

[1] 1.581139

>summary(a)

Min. 1st Qu. Median Mean 3rd Qu. Max.

1 2 3 3 4 5

最小値、

1^st quanNle

（下から

25%

）、中央値、平均値、

3^rd quanNle

（上から

25%

）、最大値最小値は

min(a)

、中央値は

median(a)

、最大値は

max(a)

でも求められる

不偏分散（

n-1

で割っている）

和は

sum

で求められる

> sum(a) [1] 15

(21)

読み込んだデータの平均や分散などを求める

21

> x <- data[,2]

> y <- data[,3]

> x

（

data[,2]

と同じ）

> y

（

data[,3]

と同じ）

x:

身長　

y:

体重　としておく

> mean(x) [1] 169.3448

> mean(y) [1] 61.89655

平均値

> var(x)

[1] 53.94828

> var (y)

[1] 87.88177

分散

> sd(x)

[1] 7.344949

> sd (y)

[1] 9.374528

標準偏差

data$Height

でも可

data$Weight

でも可

> cor(x,y)

[1] 0.6774206

> cor (x,y,

method=“spearman”) [1] 0.5897278

相関係数

(22)

5. 　作図：ヒストグラム

22

ヒストグラム

> hist(x)

これだけ！

Histogram of x

x

Frequency

150 160 170 180 190

0246810

身長のヒストグラム

> par(mar=c(5,5,1,1))

> hist(x,xlab="Height(cm)",main="",cex.axis=2,cex.lab=2)

Height(cm)

Frequency

150 160 170 180 190

0246810

par(mar=c(5,5,1,1)

：図の余白を下

5

、左

5

、上

1

、右

1

にする

xlab

：

x

軸の

label main

：タイトル

cex.axis

：軸目盛りの文字の大きさ

cex.lab

：軸ラベルの文字の大きさ

階級（データの区切り）は自動で設定してくれる

breaks=c(150,160,170,180,190

）

などとすれば階級を指定出来る

(23)

箱ひげ図（ boxplot ）

23

> boxplot(x)

これだけ！

155160165170175180185190

身長の箱ひげ図

> par(mar=c(1,5,1,1))

> boxplot(x,ylab="Height(cm)",cex.axis=2,cex.lab=2)

155160165170175180185190Height(cm)

> M<-x[1:24]

> F<-c(x[25:29],rep(NA,19))

> x<-cbind(M,F)

> par(mar=c(5,5,1,1))

> boxplot(x,

ylab="Height(cm)",cex.axis=2,cex.lab

=2)

M F

155165175185Height(cm)

Mと数を合わせるために、

NA（欠測値）を入れておく

男女別の

boxplot

M

と

F

を横に結合

(24)

散布図

24

> plot(x,y)

これだけ！

155 160 165 170 175 180 185 190

50607080

x

y

> plot(x,y,xlab=“Height(cm)”,

ylab=“Weight(kg)”,cex.axis=2,cex.lab=2)

155 165 175 185

50607080

Height(cm)

Weight(kg)

> plot(x,y,xlab=“Height(cm)”,

ylab=“Weight(kg)”,cex.axis=2,cex.lab=2, col=c(rep(4,24),rep(2,5)),cex=2)

男女で色分け

155 165 175 185

50607080

Height(cm)

Weight(kg)

赤：

col=2

、　青：

col=4

cex:

点の大きさ

pch

パラメータで◯を

×

（

pch=4)

などに変更可

h3p://cse.naro.aﬀrc.go.jp/takezawa/r-Nps/r/53.html

(25)

図の保存

図を選択して、

[

メニュー

]

⇒　

[

ファイル

]

⇒　

[

保存

]

pdf

ファイルで保存

eps

ファイルで保存

> plot(x,y)

> dev.copy2eps(ﬁle=“

ファイル名

.eps”) quartz

2

(26)

補足 1 ： R の記憶やショートカットキー、

関数の確認

R

の記憶

•  R

は、こちらが入力した数式や命令を全て覚えている。

　（

R

を閉じてしまうと忘れることもある）

• 

以前に計算した式を呼び出せる：　上矢印キー「

↑

」を押すショートカットキー（

Mac

の場合）

• 

行の最初に移動：　

Control-

「

a

」

• 

行の終わりに移動：　

Control-

「

e

」

関数（括弧を付けて使うコマンド）の確認

> ? “

関数名

”

　　　　　例：

> ?read.table

> help(“

関数名

”

）　　例：

> help(read.table)

> “

関数名

”

　　　　　　例：

> read.table

26

(27)

補足２：パッケージのインストール

•  R で追加パッケージ（拡張機能）をインストールするには以下のようにすればよい。

27

1. R

のメニュー

2. →

「パッケージ」

3. →

「パッケージのインストール」

4. →Tokyo (Japan

）など近いところを選択

5. →

パッケージ名を選択

警告メッセージが出たら対応する後は、勝手にインストールしてくれる

他の方法

> install.packages(“

パッケージ名

”,dependencies=TRUE)

後は、上記の

4

〜と同じ

(28)

補足 3 ： R におけるコメントと .r ファイル

R

におけるコメント

R

では、

”#”

以降はコメントとして扱われ、計算に影響しない

特に、右の

.r

ファイルを作って計算するときに、コメントを書いておくと便利。

（コメントが無いと忘れてしまう）

28

.r

ファイル

R

のコンソール画面に命令文を１つ１つ書かなくても、

まとめてテキストファイル（

.r

など）に書いておいて、コピペやショートカットキーを使って、まとめて実行出来る。

長い計算のときに便利。

> a<-1:5 # a

に

(1 2 3 4 5)

を入れる

# a

に

(1 2 3 4 5)

を入れる

> a<-1:5

data=read.table("data2.txt",header=T) x1<-data[1,]

x2<-data[2,]

x3<-data[3,]

x4<-data[4,]

#Odds raNo

x1[1]*x1[4]/x1[2]/x1[3]

x2[1]*x2[4]/x2[2]/x2[3]

x3[1]*x3[4]/x3[2]/x3[3]

x4[1]*x4[4]/x4[2]/x4[3]

r

ファイルの例

(29)

2-3 　配列を扱う

ベクトル：

1

次元行列：

2

次元

配列：

3

次元以上も可！

3

次元の例：

複数の（同じサイズの）行列をまとめて保存する

29

行列

1

行列

2

行列

3

2

行

2

列

配列

この場合の配列のサイズは、

2×2×3

> x1<-matrix(1:4,nc=2,by=T)

> z<-array(c(x1,x2,x3),dim=c(2,2,3))

> z , , 1

[,1] [,2]

[1,] 1 2 [2,] 3 4 , , 2

[,1] [,2]

[1,] 5 6 [2,] 7 8 , , 3

[,1] [,2]

[1,] 9 10 [2,] 11 12

z[,,1]

z[,,2]

z[,,3]

Advanced!

(30)

配列を詳しく見てみる

30

z= 1 2

3 4 5 6

7 8 9 10 11 12

z[,,1] z[,,2] z[,,3]

> z[1,1,1]

[1] 1

> z[1,2,1]

[1] 2

> z[2,1,1]

[1] 3

> z[2,2,1]

[1] 4

> z[1,1,2]

[1] 5

> z[1,2,2]

[1] 6

> z[2,1,2]

[1] 7

> z[2,2,2]

[1] 8

> z[1,1,3]

[1] 9

> z[1,2,3]

[1] 10

> z[2,1,3]

[1] 11

> z[2,2,3]

[1] 12

配列を使う場面：

Mantel-Haenzel

検定を行う場合

配列

=

層ごとの分割表をまとめたもの

層1 お菓子をよく食べる　層2 お菓子を食べない　 E not E 合計　　 E not E 合計　

歯磨き

をするしない

歯磨き

をするしない　 D 13 32 45 　 D 25 8 33 not D 7 14 21 　 not D 63 22 85 合計 20 46 66 　合計 88 30 118

D:

虫歯あり

z= 13 32

7 14 25 8 63 22

z[,,1] z[,,2]

Advanced!

(31)

統計ソフト R の応用推測統計学

1.   t 検定

2.   カイ二乗検定

3.   Fisher の正確確率検定

4.  マンテルヘンツェル検定

31

(32)

R で t 検定

> x<-data$Height[1:24]

> y<-data$Height[25:29]

> t.test(x,y)

Welch Two Sample t-test data: x and y

t = 2.8693, df = 6.112, p-value = 0.02787

alternaNve hypothesis: true diﬀerence in means is not equal to 0 95 percent conﬁdence interval:

1.34017 16.40983 sample esNmates:

mean of x mean of y

170.875 162.000 ₃₂

これだけ！

引数（パラメータ）

paired=T

（対応のある

t

検定）

var.equal=T

（等分散の

t

検定）

R

：　

t.test(x,y)

Excel: t.test(x,y,2,3)

(33)

R でカイ二乗検定

33

左の

case-control

研究の結果

について、暴露効果があるのかどうかを検定する

E

：暴露あり

E

：暴露なし合計

D 70 (a) 40 (b) 110 (m₁)

D 42 (c) 58 (d) 100 (m₀)

112(n₁) 98 (n₀) 210 (n)

> x<-matrix(c(70,40,42,58),nr=2,by=T)

> chisq.test(x,correct=FALSE) Pearson's Chi-squared test data: x

X-squared = 9.8523, df = 1, p-value = 0.001696

これだけ！

引数（パラメータ）：

correct=TRUE

（データ数が少ない

場合のイェーツの補正）

n(ad-bc)² n₁ n₀ m₁ m₀

χ²=

全てのセルが

5

以上なら、

correct=FALSE

にする

>x[1,1]*x[2,2]/(x[1,2]*x[2,1])

オッズ比！

>x[1,1]*sum(x[,2])/(x[1,2]*sum(x[,1]))

リスク比！

(34)

R でフィッシャーの正確確率検定

34

> ﬁsher.test(x)

Fisher's Exact Test for Count Data data: x

p-value = 0.002262

alternaNve hypothesis: true odds raNo is not equal to 1 95 percent conﬁdence interval:

odds raNo 2.406193

これだけ！

近似を使わない（

exact

な）、二因子の関連性を評価する検定

←

この値は条件付き最尤推定値。

　　サンプルオッズ比とは異なる（ヘルプより）

(35)

R で Mantel-Haenzel 検定

> data1<-matrix(c(13,32,7,14),nr=2,by=T)

> data2<-matrix(c(25,8,63,22),nr=2,by=T)

> data<- array(c(data1,data2),dim=c(2,2,2))

> mantelhaen.test(data, correct=FALSE)

Mantel-Haenszel chi-squared test without conNnuity correcNon data: data

Mantel-Haenszel X-squared = 0.008, df = 1, p-value = 0.9288

alternaNve hypothesis: true common odds raNo is not equal to 1 95 percent conﬁdence interval:

0.4745731 1.9737919 sample esNmates:

common odds raNo 0.967837

35

これだけ！！

基本的に、

“data”

は、

2 x 2 x k

の配列

“mantelhaen.test(data)”

にすると、

自動で

correct

するかどうか決めて

くれるよう

(36)

演習 1 ：身長・体重データの要約と検定

29

人の性別、身長、体重のデータを

”data1.txt”

から読み込み、以下を行って下さい。

1. 

男性の身長の平均値と（不偏）分散をそれぞれ求める。

2. 

男性の身長と体重の相関係数を求める。

3. 

男性と女性の身長の箱ひげ図を描き、

eps

形式で保存する。

4. 

男性と女性の身長について、

t

検定を行い、差があるかどうかを検定する。

M F 36

155165175185Height(cm)

(37)

演習 2: ﬁxed コホート研究における層化解析

641人のfixedコホート研究結果（上表,

“data2.txt”

）

について、暴露Eの疾病Dへの効果を検証して下さい。

(1)

各層とCrudeデータのそれぞれのリスク比を求め、

　　カイ二乗検定や正確確率検定を行う (2) Mantel-Haenzel検定を行う

37

層1 女性、年齢≦20 層2 女性、年齢>20 　 E not E 合計　 E not E 合計 D 4 30 34 D 5 7 12 not D 10 251 261 not D 18 61 79 合計 14 281 295 合計 23 68 91

層3 男性、年齢≦20 層4 男性、年齢>20 　 E not E 合計　 E not E 合計 D 23 29 52 D 19 5 24 not D 27 102 129 not D 36 14 50 合計 50 131 181 合計 55 19 74

Crude 　

　 E not E 合計 D 51 71 122 not D 91 428 519 合計 142 499 641

前回の演習

2

と同じデータを

R

で扱う

“data2.txt” ^a ^b ^c ^d

4 30 10 251 5 7 18 61 23 29 27 102 19 5 36 14

＊

crude

データは他の層の和

(38)

> data<-read.table(“data1.txt”,header=T)

38

[ 解答例 ] 演習 1: 身長・体重データの要約と検定

> xm<-data$Height[1:24]

> c(mean(xm),var(xm)) [1] 170.87500 44.80978

> ym<-data$Weight[1:24]

> cor(xm,ym) [1] 0.609144

> xf<-data$Height[25:29]

> xf<-c(xf,rep(NA,19))

> x<-cbind(xm,xf)

> par(mar=c(5,5,1,1))

> boxplot(x,

ylab="Height(cm)",cex.axis=2,cex.lab

=2) > dev.copy2eps

(ﬁle=“boxplot_height_M-F.eps”)

男性の身長

男性の体重

女性の身長

> t.test(xm,xf)

Welch Two Sample t-test data: xm and xf

t = 2.8693, df = 6.112, p-value = 0.02787

alternaNve hypothesis: true

diﬀerence in means is not equal to 0 95 percent conﬁdence interval:

mean of x mean of y 170.875 162.000

身長に有意差あり

(39)

[ 解答例 ] 演習 2:ﬁxed コホート研究における層化解析

> data<-read.table("data2.txt",header=T)

> data<-as.matrix(data)

39

> (RR1<-x1[1]*(x1[2]+x1[4])/(x1[2]*(x1[1]+x1[3]))) a

1 2.67619

> (RR2<-x2[1]*(x2[2]+x2[4])/(x2[2]*(x2[1]+x2[3]))) a

2 2.111801

> (RR3<-x3[1]*(x3[2]+x3[4])/(x3[2]*(x3[1]+x3[3]))) a

3 2.077931

> (RR4<-x4[1]*(x4[2]+x4[4])/(x4[2]*(x4[1]+x4[3]))) a

4 1.312727

> x=c(sum(data[,1]),sum(data[,2]), sum(data[,3]),sum(data[,4]))

> x[1]*(x[2]+x[4])/(x[2]*(x[1]+x[3])) [1] 2.524202

> x1<-data[1,]

> x2<-data[2,]

> x3<-data[3,]

> x4<-data[4,]

層

1

のリスク比

層

2

のリスク比

層

3

のリスク比

層

4

のリスク比

crude

のリスク比

←data

を行列形式にしておく

(40)

> chisq.test(matrix(x1,nr=2,by=T),correct=FALSE)

Pearson's Chi-squared test data: matrix(x1, nr = 2, by = T)

警告メッセージ

:

In chisq.test(matrix(x1, nr = 2, by = T), correct = FALSE) :

カイ自乗近似は不正確かもしれません

40

>ﬁsher.test(matrix(x1,nr=2,by=T)) Fisher's Exact Test for Count Data data: matrix(x1, nr = 2, by = T)

p-value = 0.06381

odds raNo 3.326825

層

1

(41)

41

> chisq.test(matrix(x2,nr=2,by=T),correct=FALSE) Pearson's Chi-squared test

data: matrix(x2, nr = 2, by = T)

警告メッセージ

:

In chisq.test(matrix(x2, nr = 2, by = T), correct = FALSE) :

カイ自乗近似は不正確かもしれません

> ﬁsher.test(matrix(x2,nr=2,by=T)) Fisher's Exact Test for Count Data data: matrix(x2, nr = 2, by = T)

p-value = 0.171

odds raNo 2.393584

層

2

(42)

42

層

3, 4

Crude > chisq.test(matrix(x,nr=2,by=T),correct=FALSE) Pearson's Chi-squared test

data: matrix(x, nr = 2, by = T)

X-squared = 33.7381, df = 1, p-value = 6.305e-09

(43)

43

> mantelhaen.test(array(c(x1,x2,x3,x4),dim=c(2,2,4)),correct=FALSE) Mantel-Haenszel chi-squared test without conNnuity correcNon data: array(c(x1, x2, x3, x4), dim = c(2, 2, 4))

Mantel-Haenszel X-squared = 14.3635, df = 1, p-value = 0.0001507 alternaNve hypothesis: true common odds raNo is not equal to 1 95 percent conﬁdence interval:

common odds raNo 2.544875

参考：　調整リスク比

> arr1<-0

> arr1<-arr1+x1[1]*(x1[2]+x1[4])/sum(x1)

> arr1<-arr1+x2[1]*(x2[2]+x2[4])/sum(x2)

> arr1<-arr1+x3[1]*(x3[2]+x3[4])/sum(x3)

> arr1<-arr1+x4[1]*(x4[2]+x4[4])/sum(x4)

> arr2<-0

> arr2<-arr2+x1[2]*(x1[1]+x1[3])/sum(x1)

> arr2<-arr2+x2[2]*(x2[1]+x2[3])/sum(x2)

> arr2<-arr2+x3[2]*(x3[1]+x3[3])/sum(x3)

> arr2<-arr2+x4[2]*(x4[1]+x4[3])/sum(x4)

> arr1/arr2 a

1 1.948444

分子

分母

(44)

反復をする : 　 1 から 10 までを足す

> a <- 0

> a <- a + 1

> a <- a + 2

> a <- a + 3

> a <- a + 4

> a <- a + 5

> a <- a + 6

> a <- a + 7

> a <- a + 8

> a <- a + 9

> a <- a + 10

> a [1] 55

> a <- 0

> for (i in 1:10){

> a <- a + i

> }

> a [1] 55

i

が

1

から

10

までのとき、

a <- a + i

を、ひたすら繰り返す！

地道にやる

for

文を使う

> 1+2+3+4+5+6+7+8+9+10 [1] 55

これを使えば、例えば、

検定を

100

万回行うのも難しくない

補足資料

for

文について：『

The R Tips

』（第３版）

13.2

節

繰り返し文

44

(45)

今日の内容：統計ソフト R の基礎と応用

導入

•  R

の紹介

R

の基礎

•  R

のメイン画面と電卓としての利用

•  R

における「変数」と「代入」

•  R

における保存と読み込み

• 

記述統計学（平均や分散などを求める、図を描く）

R

の応用

• 

推測統計学：検定をする

　　　　　　　　（

t

検定、カイ二乗検定、フィッシャーの正確確率検定、

　　　　　　　　　マンテルヘンツェル検定）

45

(46)

このセミナーについて

内容：　疫学と統計を復習し、交絡因子とその調整方法、

　　　　　ロジスティック回帰等を紹介する　　

目標：　交絡因子調整の検定やロジスティック回帰を理解し、

　　　　　

R

等で実行できるようになる！

ポイント：疾病の規定要因（リスク因子）を正しく同定する

日時（予定）：　毎月下旬月曜

or

火曜の午後

5

時から

1.5

時間程度

スケジュール（予定）：　全４回

　　第

1

回

(11/28)

：　疫学と統計の基礎

　　第

2

回

(12/20)

：　交絡因子とその調整方法

　　第

3

回

(1/24)

：　統計ソフト

R

の基礎と応用

　　第

4

回

(2/21?)

：　ロジスティック回帰（仮）＋

α

？

46

疫学統計セミナー