群馬大学 | 医学部 | サイトトップ | 医学情報処理演習

医学情報処理演習:2011年度第4回課題解答例

課題

テキスト表4-1に示す成人男性100万人の母集団についての身長データは,次のコードで入力できる。Xに100万人の身長データが付値される。

HT <- 144:191 # Set heights into variable HT in cm.
NUM <- c(2,3,6,17,19,56,125,219,463,915,1609,2649,4550,7214,11005,16081,
 22098,29903,39048,48312,57703,66639,73332,78051,79829,77866,73767,66321,
 57993,48410,39081,29967,22055,15810,10875,7309,4596,2726,1519,939,462,
 224,128,50,31,14,5,4) # Set the number of individuals with each height.
X <- rep(HT,NUM) # Make a new variable X as heights for a million people.

次に,ランダムな標本抽出を行うために,まずRNGkind("Mersenne-Twister")によって乱数生成アルゴリズムをメルセンヌツイスター法に指定し(注:乱数については,このページを参照されたい),set.seed(54321)によって乱数を初期化する。このことによって,解答者全員がまったく同じサンプリングを行うことができる。

RNGkind("Mersenne-Twister") # Specify how to make pseudo-random numbers (random number generator; RNG)
set.seed(54321) # Initialize the RNG

なお,ここまでのRコードは,c04enter.Rにも示すので,コピー&ペーストできる。

ここからが問題だが,10人,50人,10000人の非復元抽出を2回ずつ行い,それぞれの平均と不偏標準偏差を,母集団の平均と標準偏差とともに求め,それぞれヒストグラムを描画するコードは,次のようになる。空白に適切な文字列または数字を入れよ。また,結果から何がいえるか考察せよ。

meansdhist <- function(Y, ...) { # ヒストグラムを描き平均と不偏標準偏差を計算する関数を定義
  BoxA(Y, right=FALSE, BoxB=c(140, 200), ...) # 横軸目盛140以上200未満でヒストグラム描画
  RES <- c(mean(Y), sd(Y)) # RESという名前の数値ベクトルに平均と不偏標準偏差を付値
  BoxC(RES) <- c("mean","sd") # RESの2つの要素に"mean"と"sd"という名前を付ける
  return(RES) # ベクトルRESを返す
}

S10M <- sample(X,10,rep=FALSE) # Randomly sampling 10 from the population X without replacement
S10N <- sample(X,10,rep=FALSE) # Randomly sampling 10 from the population X without replacement
S500M <- sample(X,500,rep=FALSE) # Randomly sampling 500 from the population X without replacement
S500N <- sample(X,500,rep=FALSE) # Randomly sampling 500 from the population X without replacement
S10000M <- sample(X,10000,rep=FALSE) # Randomly sampling 10000 from the population X without replacement
S10000N <- sample(X,10000,rep=FALSE) # Randomly sampling 10000 from the population X without replacement
windows(width=600,height=768) # 横600ピクセル,縦768ピクセルのグラフィックウィンドウを開く
layout(1:7) # グラフィックウィンドウを縦に7等分する
meansdhist(S10M, main="A histogram of randomly sampled 10 values of heights, 1st trial")
meansdhist(S10N, main="A histogram of randomly sampled 10 values of heights, 2nd trial")
meansdhist(S500M, main="A histogram of randomly sampled 500 values of heights, 1st trial")
meansdhist(S500N, main="A histogram of randomly sampled 500 values of heights, 2nd trial")
meansdhist(S10000M, main="A histogram of randomly sampled 10000 values of heights, 1st trial")
meansdhist(S10000N, main="A histogram of randomly sampled 10000 values of heights, 2nd trial")
hist(X, right=FALSE, xlim=c(140, 200), main="The histogram of heights for the population X")
print(EX <- mean(X)) # Mean of X, assigned to EX
sqrt(sum((X-EX)^BoxD)/length(X)) # Standard deviation for the population X

(The program listed above is to enter the heights for one million people (the whole population) as a variable X, to initialize random number generator, and to do several kinds of random sampling with viewing the nature of those samples by calculating means and sds and by drawing histograms. Please fill the boxes from BoxA to BoxD and discuss the results.)

解答例

解答となるプログラムは,it2011-04.Rとしてダウンロードできる。

項目解答例
BoxAhist
BoxBxlim
BoxCnames
BoxD2
結果からいえること(What the result means)サンプルサイズが大きくなるにつれて,標本の分布は母集団の分布に近づくようにみえる。2回の標本抽出で得られるサンプル間の違いは,サンプルサイズが大きくなるにつれて,あまり目立たなくなる。

主なコメントへの回答

プログラムのコメントとグラフに表示されるメッセージが間違っている→ごめんなさい,その通りです。直しそびれました。解答例のコードでは修正済みです。

エラー回避の方法をもっと早く説明してほしかった→実はここに書いた通り,layoutをする前にサイズを指定してグラフィックウィンドウを開いておけば,最初からエラーは出ないので,そういうコードにしておくべきでした。申し訳ありません。解答例のコードではこの点も修正済みです。

もっと詳しい関数の引数一覧が欲しい→東工大の間瀬先生が訳してCRANで公開されている,R基本統計関数マニュアル(リンク先は筑波大ミラーのpdf,404ページあります)が参考になると思います。間瀬先生のサイトにRの基本統計関数のマニュアルについてという概要説明があります。

コンピュータの起動がこんなに遅いのに遅刻扱いはひどい→遅いといっても,実測で5分はかからなかったので,少なくとも5分過ぎまでは遅刻にしません,という扱いでは不十分でしょうか? 英語が延長したなどの理由で遅くなる場合は全体に遅くなるでしょうから,その場合は考慮します。

適切なサンプルサイズはどのように決めればよいのでしょうか→いい質問です。基本的には,目的に応じて計算する必要があります(例外はあります)。例えば,ある指標値Xの母平均値を,決まった95%信頼区間で推定したいとき,母集団での標準偏差が先行研究などからわかっていれば,サンプルサイズがいくつ以上あればいいかは計算できます。2群の平均値に統計学的に有意な差があるかどうかを検定したいとき,有意水準と検出力と母集団での標準偏差と,どれくらいの差があれば臨床的に意味があると考えていいかという4つの値が決まっていれば計算できます。これら2つの計算式は異なります。Rでもpower.t.test()などの関数がありますが,例えばヴァンダービルト大学のWilliam D. DupontとWalton D. Plummerが開発し公開しているPS (Power and Sample Size Calculations)という専用ソフトを使ってもいいと思います。


リンクと引用について