群馬大学 | 医学部 | Top

医学情報処理演習2009年度

2008年度

2007年度

2006年度

2005年度

2004年度

キーワード

データ解析，データの水準，帰無仮説，検定，推定，信頼区間

授業の目標および期待される学習効果

（目標）
統計ソフトウェアを使った実際のデータ解析を通して，医学統計学の講義で学んだ統計手法を，よりしっかりと身に付ける
（学習効果）
実験データや調査データを適切に分析できる能力を身につけることが期待される

授業の概要

実験計画や調査デザインに適したデータ構造の設計を行い，表計算ソフトで入力したデータを統計処理ソフトで適切な手法を用いて解析し，その結果を解釈する，という一連の流れを，さまざまな種類のデータについて繰り返し実践する。講義と演習を情報処理演習室で行う。

授業内容のレベル

教養教育でコンピュータ操作の基礎は身に付いているはずであり，前期に「医学統計学」で統計学の基礎は理解しているはずなので，実際の統計処理を実践することを主眼におく。

履修資格

医学部医学科２年生（「医学統計学」履修を前提とする）

関連授業科目

「医学統計学」（２年前期）

テキスト／参考書

中澤　港(2007)『Rによる保健医療データ解析演習』ピアソン・エデュケーション。草稿をオンラインで公開しているので，必ずしも買う必要はない。演習課題はweb上にその都度提示する。

また，以下のウェブサイトが参考になるであろう。

授業の形式

学務課前の情報処理演習室のコンピュータを使って，サンプルデータを分析してみる。

評価

毎回，課題レポートを提出することとし，その結果によって評価する。期末試験も実施する。

メッセージ

統計は，実際にデータを扱って初めてわかってくるものです。慣れるまでは大変かもしれませんが，やっているうちにわかってくると思います。

オフィスアワー

基本的に平日昼はオフィスにいますが，できれば火曜から金曜の17:00〜17:30が好都合です。電子メール（nminato@med.gunma-u.ac.jp）で事前に連絡ください。

授業の展開（講義日程表）

  1. 第1回：データ入力（10月5日1-3時限）：課題はこちら。解答例はこちら。
  2. 第2回：基本的な図示（10月19日1-3時限）：課題はこちら。解答例はこちら。
    【参考】ESRIジャパンの日本地図データサイト | 群馬県市町村別65歳以上人口割合データ | 塗り分け地図を画面に表示するためのRコード（日本地図データを作業ディレクトリに展開し，群馬県データを作業ディレクトリにコピーした状態で実行する） | 精子数データ | 精子数と年次の散布図を描くRコード（精子数データを作業ディレクトリにコピーして実行する）
  3. 第3回：記述統計量（10月26日1-3時限）：テキスト第3章（37ページ），pdfでは49ページ。実行上の参考事項。課題はこちら。解答例はこちら。
  4. 第4回：標本統計量と母数推定（11月2日1-3時限）：課題はこちら。解答例はこちら。
  5. 第5回：データの分布と検定の概念（11月9日1-3時限）：課題はこちら。解答例はこちら。
  6. 第6回：２群の平均値の差の検定（11月16日1-3時限）：課題はこちら。解答例はこちら。
  7. 第7回：一元配置分散分析と多重比較（11月30日1-3時限）：Dunnettの多重比較の実行方法をまとめたpdfファイルはこちら。課題はこちら。解答例はこちら。
  8. 第8回：相関と回帰（12月7日1-3時限）：課題はこちら。解答例はこちら。
  9. 第9回：計数データと比率の解析(1)（12月14日1-3時限）：課題はこちら。解答例はこちら。
  10. 第10回：計数データと比率の解析(2)（12月21日1-3時限）：3群以上の比率の差の解析と，スクリーニングにおけるROC分析を行う（テキストでいうと，9.8，9.9と10.4を扱う）。課題はこちら。解答例はこちら。
  11. 第11回：クロス集計（1月18日1-3時限）：表の操作についてはこちらを参考に。課題はこちら。解答例はこちら。
  12. 第12回：量的データのノンパラメトリックな分析法（1月25日1-3時限）：課題はこちら。解答例はこちら。
  13. 第13回：モデルの当てはめと生存時間解析（2月1日1-3時限）：課題はこちら。解答例はこちら。

試験（2月8日1-3時限）（要，この試験のためのユーザ名とパスワード）
※学外からは見えません。

基本的な事項についての補足

Rというよりもコンピュータ操作上の作法に近い，基本的なところで躓いてしまうことがあるとわかったので，以下注意事項を補足します。不明なことがあれば気軽に問い合わせてください。

Rの文法の基本の補足
・何かの動作をする関数に，データオブジェクト（変数）を与えると，その関数を変数に適用した結果を返すのが基本動作。（例）合計を計算する関数はsum()。オブジェクトxの合計を計算させたければsum(x)とすればいいし，オブジェクトyの合計を計算させたければsum(y)とすればいい。
・関数は入れ子にできる。棒グラフをグラフィックウィンドウに表示させる関数はbarplot()，カテゴリ変数について度数分布を計算させる関数はtable()なので，カテゴリ変数xについて，度数分布の棒グラフである度数分布図をグラフィックウィンドウに表示させるには，barplot(table(x))とすればいい。
・関数の値は<-を使って別の変数に保存しておける（定義できる）。既に関数名などとして予約されている名前以外のすべての名前は変数名にできる。実は，既に使われている名前でも，大抵の場合再定義できてしまうが，それは避けた方が無難。<-の左辺に使いたい変数名を書き，右辺に中身を書けばいい。
・同じ関数でもオブジェクトの種類によって異なる動作をすることがある。例えば，plot(Y~X)と書いたとき（注：~はキーボードの右上の方にあるShiftキーを押しながら^キーを打つことで入力できる記号で，チルダと読む），XとYの両方が量的変数なら散布図が描かれるが，Xがカテゴリの場合は層別箱ひげ図が描かれる。意図した通りの操作をさせたい場合は，強制的に変数の型を変えることで実現できる。この例だと，Xが量的変数でも，plot(Y~as.factor(X))とすれば層別箱ひげ図になる。ただし，明示的に箱ひげ図を描く関数を使ってboxplot(Y~X)としてもいい。
Rのオブジェクトについての基本的な注意
・文字列オブジェクトは，"身長のヒストグラム"のように半角の二重引用符（Shiftキーを押しながら数字キーの2を押しながら入力する。英語ではダブルクォーテーションマークという）で括る。二重引用符なしだと，変数名として扱われてしまう。グラフの表題をつけるときなど，表題文字列は半角二重引用符で括らないと表示されない。
・関数にオプションを付け加えるときは半角コンマ(,)で区切って追加する。例えば，身長のヒストグラムの区間の区切りを「〜を超えて〜以下」ではなく，「〜以上〜未満」にしたいとき，オプションright=Fを付け加えるには，hist(dat$HT)となっていたところを，hist(dat$HT,right=F)とする。
・Rでは半角と全角は区別されるし，大文字と小文字も区別される。基本的に関数名や予約語はすべて半角である。hist(dat$HT,xlab="身長(cm)")と書くべきところを，ｈｉｓｔ（ｄａｔ＄ＨＴ，ｘｌａｂ＝”身長（ｃｍ）”）と書くと，全角のｈｉｓｔという関数などないのでエラーが起こる。コマンド中に１つでも全角が混ざっているとうまくいかない。例えば，hist(dat$HT，xlab="身長(cm)")では，HTの次のコンマが全角になっているため，何も表示されない。
・データフレームの中の変数を参照するには"$"（半角ドル記号）を用いる。例えば，datというデータフレームの中にWTという変数が含まれているとき，WTを参照するには，dat$WTとする。
ベクトルの一部の操作（部分集合の参照）
・ベクトルの一部を操作するには，[ ]を使ってマッチング条件を入れるか，subset()を使って部分集合を作って，そのオブジェクトを関数に渡すのが基本。
（例）10人の人がいて，以下のデータがあるとする。
PERSONVALUESEX
111M
213M
315F
417F
519M
651F
752M
853F
954F
1055M
・あなたのブラウザがMozilla Firefoxならば上の表をマウスで選択してコピーし，dat<-read.delim("clipboard")でもdatというオブジェクトにデータフレームとしてデータを読み込むことができるが（注：ブラウザがInternet Explorerではうまく行かなかった），以下のコードを実行しても同じことである。 
	PERSON <- 1:10
	VALUE <- c(11,13,15,17,19,51,52,53,54,55)
	SEX <- c("M","M","F","F","M","F","M","F","F","M")
	dat <- data.frame(PERSON,VALUE,SEX)
このとき，(1)１番目から５番目の人だけのVALUEの合計を求めたいとすると，以下のどのコードでも可能である。 
	sum(dat$VALUE[1:5])
	sum(dat$VALUE[dat$PERSON<6])
	sum(subset(dat,PERSON<6)$VALUE)
[ ]の中は，条件に合う行でTRUEとなり合わない行でFALSEとなれば，TRUEの行だけに対してその後の処理がなされるので，値自身も参照できる。(2)VALUEが50より大きい人だけのVALUEの合計を求めたい場合は，以下のどちらでもいい。 
	sum(dat$VALUE[dat$VALUE>50])
	sum(subset(dat,VALUE>50)$VALUE)
(3)SEXが"M"の人だけのVALUEの合計は，以下のどちらでもいい。等しいかどうかを調べるには，==と，"="（半角等号）を2つ続けなくてはいけないことに注意。
	sum(dat$VALUE[dat$SEX=="M"])
	sum(subset(dat,SEX=="M")$VALUE)
なお，SEXが"M"の人と"F"の人を別々にVALUEの合計を出すには，tapply(dat$VALUE,dat$SEX,sum)とすると一度にできる。
２群の平均値の差の検定について補足
独立２標本の平均値の差の検定は，古典的には，まずF検定をして２群の分散に差がないかを調べ，差がない場合は通常のt検定，有意差があったら（この場合の有意水準は，第二種の過誤を見ていることになるので，0.2程度にすることが多い）Welchの検定と段階を踏む方式でなされてきたので，テキストにもそのやり方を説明してある。しかし，シミュレーションしてみると，常にWelchの方法が最良の結果が得られることが明白である。
第3回に説明する参考事項
◎Rの実行環境は，設定しないでデフォルトでセットアップすると，MDI，即ち，大きなウィンドウの中に複数のウィンドウが開き，コンソールもグラフィックデバイスもスクリプトエディタも同じ大きなウィンドウのメニューで操作することになる。また，Rcmdrライブラリを使う場合，コマンダーのウィンドウが隠れてしまって表示できなくなることがある。そのため，SDI，即ちコンソールやスクリプトエディタのそれぞれのウィンドウが完全に独立した形にする方が操作しやすい。SDIへの変更は，R GUIを起動してから，メニューの「編集」の「GUIプリファレンス」を選んで開くウィンドウ（Rgui設定エディター）で，一番上の"Single or multiple windows"でSDIの左側にチェックを入れてSave（ファイル名やフォルダはデフォルトのままでいいはず）してからOKをクリックし，いったんRを終了して，再びRを起動するとSDIに変わっているはずである。うまくいかない場合は，起動時に読みこまれるプロファイルがSaveで保存したプロファイルと異なるのが原因なので，拡張子なしのRconsoleというファイルの保存場所を変えてみるとよい。
◎今回の課題を実行するには，この2項目上に書かれている「ベクトルの一部の操作」が役に立つはずなので，適宜参照されたい。
講義日程表第3回の項に戻る。
第11回に説明する参考事項
●表の操作について，簡単にまとめておく。下の生データファイル（http://phi.med.gunma-u.ac.jp/medstat/sample11.txt）があるとする。
AGEEXPOSUREDISEASE
69"YES""YES"
54"YES""NO"
76"YES""YES"
44"YES""NO"
50"YES""YES"
70"YES""YES"
40"YES""YES"
54"YES""YES"
50"YES""YES"
47"YES""NO"
58"YES""YES"
53"YES""YES"
43"YES""YES"
51"YES""YES"
56"YES""YES"
57"YES""YES"
59"YES""NO"
45"YES""YES"
49"YES""YES"
47"YES""YES"
76"NO""NO"
62"NO""YES"
40"NO""NO"
79"NO""YES"
43"NO""YES"
61"NO""NO"
55"NO""NO"
64"NO""YES"
40"NO""YES"
46"NO""NO"
44"NO""NO"
73"NO""NO"
69"NO""NO"
60"NO""NO"
78"NO""YES"
75"NO""NO"
60"NO""NO"
40"NO""YES"
67"NO""YES"
62"NO""NO"
●変数の意味は，AGEが年齢（整数），EXPOSUREがリスク要因への曝露の有無（"YES"が曝露あり，"NO"が曝露なし），DISEASEが疾病の有無（"YES"が疾病あり，"NO"が疾病なし）である。このデータファイルでは，"YES"と"NO"が文字列なので二重引用符（ダブルクォーテーション）で括ってあるが，テキストファイルからread.delim()関数でデータを読み込むときは，二重引用符なしでも，数値として解釈できない場合は自動的に要因型として読み込まれる。
◆読み込み（このファイルを，datというデータフレームに読み込む）：dat <- read.delim("http://phi.med.gunma-u.ac.jp/medstat/sample11.txt")
◆EXPOSUREの集計（やるまでもないが）：table(dat$EXPOSURE)と打てば，以下が表示される。
 NO YES 
 20  20
◆この結果を度数分布表ベクトルとしてオブジェクトEXCに付値するには，EXC <- table(dat$EXPOSURE)とする。
◆DISEASEの集計：table(dat$DISEASE)
 NO YES 
 16  24
◆曝露ありの人のDISEASEの集計：table(dat$DISEASE[dat$EXPOSURE=="YES"])
 NO YES 
  4  16
◆この結果を度数分布表ベクトルとしてオブジェクトEXDに付値するには，EXD <- table(dat$DISEASE[dat$EXPOSURE=="YES"])とする。
◆同様に曝露なしの人のDISEASEの集計結果をオブジェクトNEDに付値するには，NED <- table(dat$DISEASE[dat$EXPOSURE=="NO"])とする。
◆この２つのオブジェクトを行方向に結合すれば，曝露の有無と疾病の有無のクロス集計結果が得られる。rbind(NED,EXD)とする。
    NO YES
NED 12   8
EXD  4  16
◆しかし，実はいきなりtable(dat$EXPOSURE,dat$DISEASE)でクロス集計できる。
      NO YES
  NO  12   8
  YES  4  16
◆もう少し格好良くクロス集計する関数としてはxtabsがある。この場合，xtabs(~EXPOSURE+DISEASE,data=dat)とすればよい。
        DISEASE
EXPOSURE NO YES
     NO  12   8
     YES  4  16
◆rbind，table，xtabsのどれを使っても，クロス集計表オブジェクトとして使える。例えば次のようにする。
Table <- xtabs(~EXPOSURE+DISEASE,data=dat)
◇表の数字が最初からわかっていれば，Table <- matrix(c(12,4,8,16),2,2)でよい。ラベルをつけるには，rownames(Table) <- c("NO","YES"); colnames(Table) <- c("NO","YES")とするか，まとめて一気にdimnames(Table) <- list(c("非曝露","曝露"),c("健康","病気"))としてもよい。
☆クロス集計表オブジェクトに対してカイ二乗検定はchisq.test(Table)
☆クロス集計表オブジェクトに対してFisherの正確確率検定はfisher.test(Table)。この例は２×２なのでオッズ比も最尤推定される。曝露と疾病の間には有意水準5％で統計学的に有意な関連があり，オッズ比は5.7（95％信頼区間1.22〜32.8）である。
●3次元の表は少し難しいが，次にまとめる。
◆年齢を60歳以上か60歳未満かの2群に区分した変数ACを作って，同じデータフレームdatに含める。次のどちらかを実行する。
dat$AC <- cut(dat$AGE,c(min(dat$AGE),60,max(dat$AGE)+1),right=FALSE)
dat$AC <- ifelse(dat$AGE>=60,"[60,80)","[40,60)")
◆ACで元データを2群に分け（subsetを使うか[ ]を使う），それぞれについて曝露と疾病のクロス集計表を作り，クロス集計表オブジェクトとして，それぞれYTABとETABに付値する。。
YTAB <- xtabs(~EXPOSURE+DISEASE,data=subset(dat,AC=="[40,60)"))
ETAB <- xtabs(~EXPOSURE+DISEASE,data=subset(dat,AC=="[60,80)"))
☆60歳未満，60歳以上で別々にカイ二乗検定するには，chisq.test(YTAB); chisq.test(ETAB)とする。しかし，期待度数5未満のセルがあるので警告が出る。
☆60歳未満，60歳以上で別々にFisherの正確確率検定するには，fisher.test(YTAB); chisq.test(ETAB)。有意水準5％でどちらも有意でなくなる。しかしサンプルサイズが小さいために検出力が足りないだけかもしれない。
◆そこでクロス表の併合を行い，3次元のクロス表を作る。
Table3 <- array(c(YTAB,ETAB),dim=c(2,2,2))とすると，3次元のクロス表がTable3にできる（ラベルが全部消えてしまうが）。Table3を表示させる（ただTable3と打てばよい）と次のように見える。
, , 1

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

, , 2

     [,1] [,2]
[1,]    8    5
[2,]    0    3
◆実は，同じ表は，いきなりTable3 <- xtabs(~EXPOSURE+DISEASE+AC,data=dat)でできるし，こうした方がラベルが残っていい。
◇Table3 <- table(dat$EXPOSURE,dat$DISEASE,dat$AC)でも同じこと。
◇表の数字が最初からわかっていれば，Table3 <- array(c(4,4,3,13,8,0,5,3),dim=c(2,2,2))とする。名前を付けるには，dimnames(Table3) <- list(c("非曝露","曝露"),c("健康","病気"),c("[40,60)","[60,80)"))のようにする。
◆3次元の表から年齢層別に２×２クロス集計表を取りだすことは，YTAB <- Table3[,,1]とかETAB <- Table3[,,2]で可能である。
☆3次元の表について，「各層で曝露と疾病に関連はない」という帰無仮説（ただし３次の交互作用はないとする。つまり，対立仮説が「どの層でも共通した関連がある」となる）を検定する方法として，コクラン＝マンテル＝ヘンツェルの検定があり，mantelhaen.test(Table3)で可能。次の結果が得られ，帰無仮説が有意水準5％で棄却されるので，年齢で層別しても，この曝露と疾病の間には統計学的に有意な共通の関連があったといえる。また共通オッズ比は7.3 [1.29, 41.6]であり，年齢で層別した場合に，どの年齢層でも共通して非曝露群に比べて曝露群での疾病オッズが7.3倍と見ることができる。
        Mantel-Haenszel chi-squared test with continuity correction

data:  Table3 
Mantel-Haenszel X-squared = 3.9511, df = 1, p-value = 0.04684
alternative hypothesis: true common odds ratio is not equal to 1 
95 percent confidence interval:
  1.293228 41.584133 
sample estimates:
common odds ratio 
         7.333333
☆ちなみにマンテル＝ヘンツェルの検定の前提とした「３次の交互作用はない」自体を帰無仮説とする検定をWoolfの検定というが，vcdライブラリに入っていて，require(vcd); woolf_test(Table3)のようにすれば検定できる。
◎講義日程表第11回の項に戻る。

リンクと引用について