• 検索結果がありません。

疫学統計セミナー

N/A
N/A
Protected

Academic year: 2021

シェア "疫学統計セミナー"

Copied!
46
0
0

読み込み中.... (全文を見る)

全文

(1)

疫学統計セミナー

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

担当: 茅野光範

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

メール

: kayano@、内線:5521

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

H29.1.24

1

セミナー資料: 

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

(2)

このセミナーについて

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

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

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

     

R

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

ポイント: 疾病の規定要因(リスク因子)を正しく同定する

日時(予定): 毎月下旬月曜

or

火曜の午後

5

時から

1.5

時間程度

スケジュール(予定): 全4回

  第

1

(11/28)

: 疫学と統計の基礎

  第

2

(12/20)

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

  第

3

(1/24)

: 統計ソフト

R

の基礎と応用

  第

4

(2/21?)

: ロジスティック回帰(仮)+

α

2

(3)

1 回目の目標と内容

目標:コホート研究(追跡、

dynamic/fixed

)と症例対照研究(

case/

control

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

検証(検定)する。 (ただし、交絡は無視する)

内容:

• 

はじめに

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

• 

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

follow-up

研究)、症例対照研究

• 

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

   罹患率(

incidence raNo

)、有病率(

prevalence

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

3

復習

(4)

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

• 

復習

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

   暴露効果の指標と推測

• 

交絡の調整方法

1

 研究方法で調整する (考え方の紹介)

   

(1)

因子範囲の制限     

(2)

マッチング

• 

交絡の調整方法

2

 解析方法で調整する

   

(1)

層化解析(

StraNficaNon

):マンテル・ヘンツェル検定    

(2)

回帰分析 

 第4回目にやります

4

(交絡の影響は無視)

復習

(5)

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

導入

•  R

の紹介

R

の基礎

•  R

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

•  R

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

•  R

における保存と読み込み

• 

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

R

の応用

• 

推測統計学:検定をする

        (

t

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

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

5

(6)

R の紹介

•  R とは?

•  R の関連ソフト

•  R の資料

6

(7)

R とは?

• 

フリーの統計解析ソフト

利点

• 

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

• 

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

• 

パッケージの追加も簡単(フリー!)

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

• 

グラフィックスが豊富

• 

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

欠点

• 

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

• 

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

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

』(第3版) 

舟尾 暢男 (編集)

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

Excel

R

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

R

による生物医学統計』

  

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

•  R Tips h3p://cse.naro.affrc.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

年5月、6月)  

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

最小値、

1st quanNle

(下から

25%

)、中央値、平均値、

3rd 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.affrc.go.jp/takezawa/r-Nps/r/53.html

(25)

図の保存

図を選択して、

[

メニュー

]

⇒ 

[

ファイル

]

⇒ 

[

保存

]

pdf

ファイルで保存

eps

ファイルで保存

> plot(x,y)

> dev.copy2eps(file=“

ファイル名

.eps”) quartz

2

(26)

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

関数の確認

R

の記憶

•  R

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

  (

R

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

• 

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

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

Mac

の場合)

• 

行の最初に移動: 

Control-

a

• 

行の終わりに移動: 

Control-

e

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

> ? “

関数名

     例:

> ?read.table

> help(“

関数名

)  例:

> help(read.table)

> “

関数名

      例:

> read.table

26

(27)

補足2:パッケージのインストール

•  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

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

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

.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)

> x2<-matrix(5:8,nc=2,by=T)

> x3<-matrix(9:12,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 difference in means is not equal to 0 95 percent confidence interval:

1.34017 16.40983 sample esNmates:

mean of x mean of y

170.875 162.000 32

これだけ!

引数(パラメータ)

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 (m1)

D 42 (c) 58 (d) 100 (m0)

112(n1) 98 (n0) 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)2 n1 n0 m1 m0

χ2 =

全てのセルが

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

> fisher.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 confidence interval:

1.336666 4.376680 sample esNmates:

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 confidence 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: fixed コホート研究における層化解析

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

(file=“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

difference in means is not equal to 0 95 percent confidence interval:

1.34017 16.40983 sample esNmates:

mean of x mean of y 170.875 162.000

身長に有意差あり

(39)

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

> 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)

X-squared = 4.1881, df = 1, p-value = 0.04071

警告メッセージ

:

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

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

40

>fisher.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

alternaNve hypothesis: true odds raNo is not equal to 1 95 percent confidence interval:

0.717021 12.465966 sample esNmates:

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)

X-squared = 1.9665, df = 1, p-value = 0.1608

警告メッセージ

:

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

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

> fisher.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

alternaNve hypothesis: true odds raNo is not equal to 1 95 percent confidence interval:

0.531883 10.038545 sample esNmates:

odds raNo 2.393584

2

(42)

42

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

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

X-squared = 10.0638, df = 1, p-value = 0.001512

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

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

X-squared = 0.4364, df = 1, p-value = 0.5088

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 confidence interval:

1.544096 4.194290 sample esNmates:

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

』(第3版)

13.2

繰り返し文

44

(45)

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

導入

•  R

の紹介

R

の基礎

•  R

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

•  R

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

•  R

における保存と読み込み

• 

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

R

の応用

• 

推測統計学:検定をする

        (

t

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

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

45

(46)

このセミナーについて

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

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

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

     

R

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

ポイント: 疾病の規定要因(リスク因子)を正しく同定する

日時(予定): 毎月下旬月曜

or

火曜の午後

5

時から

1.5

時間程度

スケジュール(予定): 全4回

  第

1

(11/28)

: 疫学と統計の基礎

  第

2

(12/20)

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

  第

3

(1/24)

: 統計ソフト

R

の基礎と応用

  第

4

(2/21?)

: ロジスティック回帰(仮)+

α

46

次回

参照

関連したドキュメント

学校に行けない子どもたちの学習をどう保障す

これはつまり十進法ではなく、一進法を用いて自然数を表記するということである。とは いえ数が大きくなると見にくくなるので、.. 0, 1,

と言っても、事例ごとに意味がかなり異なるのは、子どもの性格が異なることと同じである。その

注)○のあるものを使用すること。

とである。内乱が落ち着き,ひとつの国としての統合がすすんだアメリカ社会

 英語の関学の伝統を継承するのが「子どもと英 語」です。初等教育における英語教育に対応でき

Q7 

 本計画では、子どもの頃から食に関する正確な知識を提供することで、健全な食生活