医学情報処理演習

キーワード

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

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

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

授業の概要

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

授業内容のレベル

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

履修資格

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

テキスト／参考書

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

中澤　港(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@med.gunma-u.ac.jp）で事前に連絡ください。

授業の展開

第1回：データ入力（10月6日1-3時限）：課題：p01-2008.xlsを読み込み，データ構造について記述せよ。（解答と解説の例）
第2回：基本的な図示（10月20日1-3時限）：課題はこちら。解答例はこちら。
ESRIジャパンの日本地図データサイト | 群馬県市町村別65歳以上人口割合データ | 塗り分け地図を画面に表示するためのRコード（日本地図データを作業ディレクトリに展開し，群馬県データを作業ディレクトリにコピーした状態で実行する） | 精子数データ | 精子数と年次の散布図を描くRコード（精子数データを作業ディレクトリにコピーして実行する）
第3回：記述統計量（10月27日1-3時限）：課題はこちら。解答例はこちら。
第4回：標本統計量と母数推定（11月10日1-3時限）：課題はこちら。解答例はこちら。
第5回：データの分布と検定の概念（11月17日1-3時限）：課題はこちら。解答例はこちら。
第6回：２群の平均値の差の検定（12月1日1-3時限）：課題はこちら。解答例はこちら
第7回：一元配置分散分析と多重比較（12月8日1-3時限）：課題はこちら。解答例はこちら
第8回：相関と回帰（12月15日1-3時限）：課題はこちら。解答例はこちら
第9回：計数データと比率の解析（12月22日1-3時限）：課題はこちら。解答例はこちら
第10回：クロス集計（1月5日1-3時限）：課題はこちら。解答例はこちら
第11回：量的データのノンパラメトリックな分析法（1月19日1-3時限）：課題はこちら。解答例はこちら
第12回：一般化線形モデル（1月26日1-3時限）：課題はこちら。解答例はこちら
第13回：生存時間解析（2月2日1-3時限）：課題はこちら。解答例はこちら

試験（2月9日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人の人がいて，以下のデータがあるとする。

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ではうまく行かなかった），以下のコードを実行しても同じことである。

	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の方法が最良の結果が得られることが明白である。

リンクと引用について