群馬大学 | 医学部 | Top
医学情報処理演習
2010年度 | 2009年度 | 2008年度 | 2007年度 | 2006年度 | 2005年度 | 2004年度
キーワード
データ解析，データの水準，帰無仮説，検定，推定，信頼区間
授業の目標および期待される学習効果
- （目標）
- 統計ソフトウェアを使った実際のデータ解析を通して，医学統計学の講義で学んだ統計手法を，よりしっかりと身に付ける
- （学習効果）
- 実験データや調査データを適切に分析できる能力を身につけることが期待される
授業の概要
実験計画や調査デザインに適したデータ構造の設計を行い，表計算ソフトで入力したデータを統計処理ソフトで適切な手法を用いて解析し，その結果を解釈する，という一連の流れを，さまざまな種類のデータについて繰り返し実践する。講義と演習を情報処理演習室で行う。
授業内容のレベル
教養教育でコンピュータ操作の基礎は身に付いているはずであり，前期に「医学統計学」で統計学の基礎は理解しているはずなので，実際の統計処理を実践することを主眼におく。
履修資格
医学部医学科２年生（「医学統計学」履修を前提とする）及び生命医科学専攻修士１年生
関連授業科目
「医学統計学」（２年前期）
テキスト／参考書
中澤 港(2007)『Rによる保健医療データ解析演習』ピアソン・エデュケーション。草稿をオンラインで公開しているので，必ずしも買う必要はない。演習課題はweb上にその都度提示する。
- 中澤 港(2003)『Rによる統計解析の基礎』（ピアソン・エデュケーション）1,800円
- 間瀬 茂・神保 雅一・鎌倉 稔成・金藤 浩司(2004)『工学のための数学３ 工学のための データサイエンス入門 ― フリーな統計環境Rを用いたデータ解析 ― 』（数理工学社）2,300円
- 岡田昌史（編）(2004)『The R Book −データ解析環境Rの活用事例集−』（九天社）
- 舟尾暢男 (2005) 『The R tips −データ解析環境Rの基本技・グラフィックス活用集』（九天社）
- U.リゲス（著），石田基広（訳）『Rの基礎とプログラミング技法』（シュプリンガー）
- P.Dalgaard（著），岡田昌史（監訳）『Rによる医療統計学』（丸善）
また，以下のウェブサイトが参考になるであろう。
- Rによる統計処理（本学社会情報学部・青木繁伸教授による）
- RjpWiki（筑波大学の岡田さんが管理されている，日本におけるRユーザコミュニティの代表的なサイト）
授業の形式
学務課奥の情報処理演習室のコンピュータを使って，サンプルデータを分析してみる。
評価
毎回，課題レポートを提出することとし，その結果によって評価する。期末試験も実施する。なお，この演習は生命医科学専攻修士課程の講義を兼ねているが，修士の学生については期末試験は行わず，平常の課題レポート点だけで評価するので注意されたい。
メッセージ
統計は，実際にデータを扱って初めてわかってくるものです。慣れるまでは大変かもしれませんが，やっているうちにわかってくると思います。
オフィスアワー
基本的に平日昼はオフィスにいますが，できれば月曜から木曜の17:00〜17:30が好都合です。電子メール（nminato[atmark]med.gunma-u.ac.jp // [atmark]は半角の＠です）で事前に連絡ください。
授業の展開（講義日程表）
- 第1回：データ入力（10月3日1-3時限）：出席登録はこちら。作業フォルダの設定等について：課題はこちら。出てくる関数・文の主なものの一覧はこちら。課題解答例はこちら。
- 第2回：基本的な図示（10月17日1-3時限）：出席登録はこちら。課題はこちら。出てくる関数・文の主なものの一覧はこちら。課題解答例はこちら。
- 第3回：記述統計（10月24日1-3時限）【pdfファイル49ページ参照】ベクトルの一部の操作について。出席登録はこちら。課題はこちら。出てくる関数・文の主なものの一覧はこちら。課題解答例はこちら。
- 第4回：標本統計量と母数推定（10月31日1-3時限）：出席登録はこちら。課題はこちら。出てくる関数・文の主なものの一覧はこちら。課題解答例はこちら。
- 第5回：データの分布（11月7日1-3時限）：出席登録はこちら。課題はこちら。出てくる関数・文の主なものの一覧はこちら。課題解答例と主なコメント・質問への回答はこちら。
- 第6回：平均値の差の検定（11月14日1-3時限）：出席登録はこちら。課題はこちら。出てくる関数・文の主なものの一覧はこちら。課題解答例と主なコメント・質問への回答はこちら。
- 第7回：一元配置分散分析と多重比較（11月21日1-3時限）：出席登録はこちら。課題はこちら。出てくる関数・文の主なものの一覧はこちら。課題解答例と主なコメント・質問への回答はこちら。
- 第8回：相関と回帰（11月28日1-3時限）：出席登録はこちら。課題はこちら。出てくる関数・文の主なものの一覧はこちら。課題解答例と主なコメント・質問への回答はこちら。
- 第9回：計数データと比率の解析[1]（12月5日1-3時限）：出席登録はこちら。課題はこちら。出てくる関数・文の主なものの一覧はこちら。課題解答例と主なコメント・質問への回答はこちら。
- 第10回：計数データと比率の解析[2]（12月12日1-3時限）：出席登録はこちら。課題を除く演習内容はこちら。課題はこちら。出てくる関数・文の主なものの一覧はこちら。課題解答例と主なコメント・質問への回答はこちら。
- 第11回：クロス集計[1]（12月19日1-3時限）：出席登録はこちら。課題を除く演習内容はこちら。課題はこちら。出てくる関数・文の主なものの一覧はこちら。課題解答例と主なコメント・質問への回答はこちら。
- 第12回：クロス集計[2]（1月16日1-3時限）：出席登録はこちら。課題を除く演習内容はこちら。課題はこちら。出てくる関数・文の主なものの一覧はこちら。課題解答例と主なコメント・質問への回答はこちら。
- 第13回：量的データのノンパラメトリックな分析法（1月23日1-3時限）：出席登録はこちら。課題はこちら。出てくる関数・文の主なものの一覧はこちら。課題解答例と主なコメント・質問への回答はこちら。
- 第14回：モデルの当てはめと生存時間解析（1月30日1-3時限）：出席登録はこちら。課題を除く演習内容はこちら。今回，課題はありませんが，こちらから感想やコメントを送ってください。
期末試験（2月6日1-3時限）Aクラス | Bクラス（要，この試験のためのユーザ名とパスワード）
※学外からは見えません。
基本的な事項についての補足
Rというよりもコンピュータ操作上の作法に近い，基本的なところで躓かないために，以下注意事項を補足します。不明なことがあれば気軽に問い合わせてください。
- 作業フォルダの設定等について
- ●Vistaのせいだと思いますが，演習室ではRのファイルメニューの「ディレクトリの変更」がうまく動作しません。Vistaの場合，管理者権限がないために制限される機能が非常に多いのですが，暫定的な対処法を考えたので説明しておきます。
- ●簡単に言うと，マイドキュメントを作業フォルダとしてRを起動させるための専用のアイコンを作ってしまえばいいのです。手順は以下の通りです。
- スタートメニューの「すべてのプログラム」「R」「R-2.11.1」の上でマウスの右ボタンをクリックしたままデスクトップまでドラッグし，手を離したときに出てくるウィンドウで「コピー」を選ぶ
- コピーされた起動アイコンの上でマウスの右ボタンをクリックし，「プロパティ」を選ぶ→プロパティウィンドウ（＊）が表示される
- スタートメニューの「ドキュメント」の上でマウスの右ボタンをクリックし（もしかするとキーボードの[Shift]キーを押しながらクリックする必要があるかもしれません），「パスとしてコピー」を選ぶ
- 先に表示されたプロパティウィンドウ（＊）の，「作業フォルダ」のところにマウスカーソルをもっていき，右ボタンをクリックして「貼り付け」を選ぶ
- （参考）プロパティウィンドウ（＊）の「リンク先」のところで，"****\Rgui.exe"となっている後にスペースを空けてオプションをつけておくと，表示言語やウィンドウタイプの設定ができる。デフォルトではMDIといって大きなウィンドウの中に複数の小さなウィンドウが開く形で動作するが，--sdiオプションを付けておけば完全に独立したウィンドウが複数開く形で動作する。また，LANG=C LC_ALL=Cオプションを付けると表示言語が英語になる。
- [OK]ボタンをクリックしてウィンドウを閉じる
デスクトップに残っている「R」の文字のアイコンをダブルクリックしてRを起動すれば，作業フォルダが，Excelなどでファイルを保存できる場所と共通の，マイドキュメントになるはずです。
- 講義日程表第1回の項に戻る。
- 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人の人がいて，以下のデータがあるとする。
| PERSON | VALUE | SEX |
|---|
| 1 | 11 | M |
| 2 | 13 | M |
| 3 | 15 | F |
| 4 | 17 | F |
| 5 | 19 | M |
| 6 | 51 | F |
| 7 | 52 | M |
| 8 | 53 | F |
| 9 | 54 | F |
| 10 | 55 | M |
- ・あなたのブラウザがMozilla Firefoxならば上の表をマウスで選択してコピーし，dat<-read.delim("clipboard")でもdatというオブジェクトにデータフレームとしてデータを読み込むことができるが（注：ブラウザがInternet Explorerではうまく行かなかった。ただ，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)とすると一度にできる。
- 講義日程表第3回の項に戻る。
- ２群の平均値の差の検定について補足
- 独立２標本の平均値の差の検定は，古典的には，まずF検定をして２群の分散に差がないかを調べ，差がない場合は通常のt検定，有意差があったら（この場合の有意水準は，第二種の過誤を見ていることになるので，0.2程度にすることが多い）Welchの検定と段階を踏む方式でなされてきたので，テキストにもそのやり方を説明してある。しかし，シミュレーションしてみると，常にWelchの方法が最良の結果が得られることが明白である。
- 第11回に説明する参考事項
- ●表の操作について，簡単にまとめておく。下の生データファイル（http://phi.med.gunma-u.ac.jp/medstat/sample11.txt）があるとする。
| AGE | EXPOSURE | DISEASE |
|---|
| 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回の項に戻る。
リンクと引用について