東京農工大学

生物統計学演習１

　以下の問題の解答は，グラフや表をワードに貼り付けて，考察を書いて提出．R のスクリプトは 特に必要ない．提出は 10月23日（金）の昼頃までに，クラスルームを用いた課題提出か、
fs7263@go.tuat.ac.jp
まで、メールの添付ファイルで送付して下さい． 解答例は次回講義の開始時もしくは次週にアップする予定です．

問題１：12 人のきょうだい中の女児数のデータ

19 世紀末のドイツの病院のデータによると，同じ両親で 12 人きょうだいがいる 6115 家族の女児数 の数は以下のようであった．

   女児数       0    1    2    3    4    5    6
度数    7       45       181       478       829       1112       1343   
   女児数       7    8    9    10    11    12    合計
度数    1033       670       286       104       24       3   6115   

「Weldon のサイコロ実験」と同様な解析を「12人のきょうだい中の女児数のデータ」 に対して行え．
具体的には，Weldon サイコロ実験の解析結果の表のようなものをつくり，結果の解釈を行う．

ヒント：サイコロでは 5 か 6 の目が出る確率は p0 = 1/3 と考えたが， このデータでは男か女なので p0 = 1/2 である．

問題２：捕獲再捕獲法の推定精度

ここでは捕獲（m = 50），再捕獲（k = 50）合わせて100頭を捕獲している．捕獲総数100頭 を固定して，捕獲頭数と再捕獲頭数を変えたときの推定精度はどうなるか． また，予算の関係で捕獲総数を半分の50頭にした場合の推定精度はどうなるか．

注）推定精度は，推定値の平均 mv や標準偏差（標準誤差）sd を求め，そこから正規近似により，未知 母数 θ の 95 ％信頼区間（Wald 信頼区間）として，

mv - 1.96×sd < θ < mv + 1.96×sd
を出すのが一般的である．
　しかし，この例では，推定値が無限大に発散してしまうことがあるので推定値の 平均や標準誤差を求めることができない．このため，シミュレーションにより行った．
　推定精度に必要な項目として，推定値の発散割合，メディアン（中央値），90％もしくは95％信頼区間を求める．

ヒント：m = 50，k = 50 のときは， 

> m <- 50	# 標識をつけた数　
> n <- 450	# 標識をつけられていない数　
> k <- 50	# 再捕獲数
> x <- rhyper(1000, m, n, k)	# 超幾何分布に従う乱数1000抽出　
> estimate <- k * m / x			# 個体数推定値列ベクトル　
> table(estimate) 		#推定値の階級分け　
estimate
192.307692307692 208.333333333333 227.272727272727              250 
               1                2                5                7 
277.777777777778            312.5 357.142857142857 416.666666666667 
              26               59              110              187 
             500              625 833.333333333333             1250 
             193              182              137               58 
            2500              Inf 
              30                3 
> quantile(estimate, c(0.05, 0.1, 0.5, 0.9, 0.95)) # 分位点（パーセンタイル）　
       5%       10%       50%       90%       95% 
 312.5000  352.6786  500.0000  833.3333 1250.0000

となり，真値は 500頭なのに推定頭数が無限大（Inf）になったケースが 3回あった．無限大に 発散するケースは少ないほどよい．また，推定値分布の分位点から全体のケースのうちの 90％が 312.5 から 1250までの間に含まれていた．この区間の幅が短いほど推定精度がよいことになる．

たとえば，m = 30，k = 70 のときは， 

> m <- 30	# 標識をつけた数　
> n <- 470	# 標識をつけられていない数　
> k <- 70	# 再捕獲数　
> x <- rhyper(1000, m, n, k)	# 超幾何分布に従う乱数1000抽出　
> estimate <- k * m / x			# 個体数推定値列ベクトル　
> table(estimate) 		#推定値の階級分け　
estimate
190.909090909091              210 233.333333333333            262.5 
               1                4                7               30 
             300              350              420              525 
              71              129              168              208 
             700             1050             2100              Inf 
             187              137               48               10 
> quantile(estimate, c(0.05, 0.1, 0.5, 0.9, 0.95)) # 分位点（パーセンタイル）　
  5%  10%  50%  90%  95% 
 300  300  525 1050 2100
となり，推定値が無限大に発散したケースは 10回で先ほどより増えた．また， 推定値の 90％区間は，300から2100で，これも先ほどより幅が長くなった．これより， この組み合わせはよくないことがわかる．

問題３：工場労働者の事故数の分布

　工場労働者一人当たりの遭遇事故数の分布データ（Greenwood et al. 1920）が以下のようであった．

上のデータをポアソン分布と負の二項分布にあてはめ，どちらの分布のあてはまりがよいか 考察せよ．ただし，'5 以上'は 5 とみなすとする．（正確ではない．）

ヒント：ポアソン分布へのあてはめは，上の死亡記事件数の例と同じである．
　負の二項分布へのあてはめは，データの平均と分散が負の二項分布の平均と分散と等しいと考えることにより， 負の二項分布のパラメータの推定ができる．これをモーメント法による分布パラメータの推定という．
　なお，負の二項分布の確率密度は dnbinom(x, n, p) で与えられる．

　すなわち，データの平均が m，分散が v であったとする．負の二項分布の平均は n(1 - p)/p， 分散は n(1 - p)/p2 であるので，これらを等しいとおいた連立方程式，

m = n(1 - p)/p, v = n(1 - p)/p2

を解けば負の二項分布のパラメータ n，p の値が m，v で表せる．

問題４：ダイズ地上部乾物重データの正規分布へのあてはめ

　長野県中信試験場（現在は，長野県野菜・花卉試験場）で，畦間 75 cm，株間 15 cm で栽培された ダイズ品種エンレイの個体ごとの地上部 乾物重データの平均、標準偏差、25％分位点、メディアン、75％分位点を求めよ．
　ボックスプロットとデータヒストグラムにあてはめた正規分布を上書したグラフを書き、正規分布への当てはまりの程度を考察せよ。．

　夏作物では周辺部分は，隣接個体が少ないことから競合が少なく，また受光や風通し がよくなったりして圃場内部の個体より生育条件が有利になることがある．これを周辺効果 (border effect) という．
　周辺効果の影響を避けるため，圃場実験では周囲の１〜２畦や４〜５株のデータは解析から外すことが多い． ここで取り上げたデータは，畦１が西側，１行目が北側境界であった．
　周辺効果の有無を確認してみよう．これは、データの行平均、列平均を求め境界に近い部位が他の部位より大きいかを比較検討すればよい．
　特定した周辺部分のデータを除いたときの、 平均、標準偏差、25％分位点、メディアン、75％分位点を求めよ． また、このデータを正規分布にあてはめ，あてはまりの変化を考察せよ．
　また，この栽培試験では 2 粒撒きであったが，NA のデータ欠損は何らかの個体の異常により、発芽しなかったか生育途中で枯死したものである．従って，異常に発育の悪い個体は、その個体自体に何らかの問題が あったため発育不良になったものと考えられる．このため，正常に生育した個体の平均値 を知るためには極端に生育の悪い個体を取り除く必要がある．極端に生育悪い個体（たとえば20以下）を取り除くと，平均、標準偏差、25％分位点、メディアン、75％分位点 はどうなるか．さらに、正規分布へのあてはめはどうなるか．

ヒント：

# 周辺部分の削除として、たとえば、 1 - 4 行目と 1 - 2 列のデータ以外を使うときは
daizu2 <- as.vector(daizu.mat[-(1:4), -(1:2)])
# とすればよい．
# また，異常に生育の悪い個体を除く（普通のを残す）には，
daizu3 <- daizu2[daizu2>20]
# のようにすればよい．
Copyright (C) 2011, Hiroshi Omori. 最終更新日：2020年9月28日