東京農工大学
2020.10.9

生物統計学１

東京大学大学院農学生命科学研究科　大森宏

この講義の目的

　統計解析ソフトＲを用いて，統計解析の理論と実践を学ぶ

１．R の基本的使い方

基礎知識

  1. 作業ディレクトリ（フォルダ）
    作業ディレクトリに入力したいデータ等を置いておく。Rからの出力もここに保存される。
  2. 作業（ワーク）スペース
    R上で行った作業結果を保存しているファイル。Rを終了したときに保存するかを聞いてくる。作業ディレクトリに .RData で保存。
  3. Rの実行とエディタの利用
    Rを開いた後、コンソール画面に直接コードを打ち込むか、ウェブページなどに記載されているRコードをコンソールにコピペすれば実行し、結果が表示される。
    ↑ボタンを押すと、1回前のコードが実行される。
    しかし、エディタを立ち上げ、ここにコードを書き込み、Rに実行させた方がコードの編集ができるので便利である。 Rコンソール画面でコードと実行結果が表示される。
    通常、１コードを１行で記述する。１行にコードを連続して記述したいときは、間に";"を入れる。
    エディタでコード複数行をマウスで選択状態にして上記と同様な操作を行うと、複数コードが連続して実行される。
  4. ファイル（データ）入出力
    エクセルとの連携を考えると、csvファイルが使いやすい。
  5. グラフの保存
    グラフウィンドウにグラフが表示された後、
    ファイル -> 別名で保存 -> png

基本演算

1+1		# "#"以降はコメント文で、実行されない。
x <- 3		# "<-" は代入の意味
x
x + 5*x		# x を用いた演算
# 特別な演算子
5 ^ x  		# 5の3乗
5 %/% x 	# 商
5 %% x 		# 余り
# 値の大小の確認
x < 3  		# 小なり 
x <= 3  	# 小なりイコール
x == 3  	# 等しい
x != 3  	# 等しくない
# 種々の数学関数が定義されている。
abs(x)  	# 絶対値
sqrt(x)  	# 平方根
sin(x)  	# sine 正弦
atan(x)  	# arc tangent 逆正接
log(x)  	# log e 底eの対数
factorial(x) 	# x! xの階乗

ベクトルや行列の演算

1:10					# 数列生成
seq(0,1, by=0.05)		# 数列生成
rep(1, 5)				# 繰り返し
rep(1:4, 3)				# 繰り返し
rep(1:4, each=3)		# 繰り返し
x <- c(2,-1,3,1)		# ベクトル生成
y <- 1:4
x 
y 
2*x + y					# ベクトル演算
x * y					# 要素同士の積
x / y					# 要素同士の割り算
t(x) %*% y				# 内積、t()は転置
x %*% t(y)				# 行列生成
z <- c(-2,2,-5,-1)
a <- cbind(x,y,z); a 	# 列ベクトル－＞行列生成
a[2,3]					# 要素取り出し
a[1,]					# 行取り出し
a[,1:2]					# 2列取り出し
a[2:3,2:3]				# 部分行列取り出し
a[-(1:2),c(1,3)]		# 複雑な部分行列取り出し
aa <- t(a) %*% a; aa	# 行列の積	
diag(aa)				# 対角成分
solve(aa)				# 逆行列
b <- matrix(1:9, nrow=3); b	# 行列作成
a %*% b					# 行列の積

データ形式

グラフ表示

散布図はplot()で簡単に描ける。回帰直線も簡単にひける。 
x <- 1:10
y <- c(3,5,6,6,9,10,6,9,10,7) 
plot(x, y)			# x, y 散布図
reg <- lm(y ~ x)		# 回帰の計算
reg$coefficients		# 回帰係数
abline(reg, col="red")		# 回帰直線の表示
cor(x, y)			# 相関係数
title(main="散布図と回帰直線") 	# グラフタイトル
#
curve(sin(x), 0, 2*pi)	# sin カーブ
abline(h=0)	
title(main="Graph of sin(x)")		
# 
#グラフの重ね合わせ
curve(sin(x), 0, 2*pi)	# sin カーブ
abline(h=0)				# x軸
abline(v=0)				# y軸
curve(cos(x), 0, 2*pi, col="red", add=T)	# add=T で重ねる。
# 
# グラフの重ね合わせ２
x <- seq(0, 2*pi, by=0.01)
plot(x, sin(x), type="l")
lines(x, cos(x), col="red")
abline(h=0)	
abline(v=0)
# 
# 自作関数のグラフ
# 関数定義
f2 <- function(x){ 
x*log(x)-x 
}  
curve(f2(x),0,4)		# 
abline(h=0)
abline(v=0)

２．データの図示と基本統計量の計算

要約統計量

　大きさ n のデータ、 x1、x2、‥、xn が得られたとき、 その標本平均（average, sample mean）x- と 標本分散（sample variance）s2 はそれぞれ、
mean
で与えられ、標本標準偏差（sample standard deviation : SD）は 分散の平方根で与えられる。
　平均は、データの最もありそうな値であり、標準偏差は、データのちらばりの程度を表している。

　平均とほぼ同様な値としてメディアン(median)（中央値、50％点(percentile)、1/2分位点(1/2 quantile)） がある。データのちらばりの程度は、 25％点（1/4分位点、第一四分位数(1st quartile)）と 75％点（3/4分位点、第三四分位数(3rd quartile)）の幅である 四分位範囲(interquartile range : IQR)でも測ることができる。

データの図示と要約統計量の取得

英語得点データを用いて，データの基本統計量の計算演習を行う． 
#英語の得点	
eigo <- c(	
36,70,56,68,76,60,50,63,62,42,64,60,50,68,71,67,	
50,65,67,57,72,64,61,66,46,80,46,51,59,32,55,65,	
65,52,57,64,23,57,53,54,38,71,57,69,77,61,51,64,	
63,43,65,61,51,69,72,68,53,66,68,58,73,65,62,67,	
47,81,47,52,59,33,56,66,67,52,58,65,24,58,54,55)	
length(eigo)	     	 #データ数　
summary(eigo)			# データの要約統計量
mean(eigo)	     	 	#標本平均　
var(eigo)	     	 	#標本分散　
sd(eigo)	     	 	#標本標準偏差　
boxplot(eigo, main="英語得点の箱ヒゲ図")	     	 #箱ひげ図
boxplot.stats(eigo)	     	 #箱ひげ図用統計量　
hist(eigo, breaks=seq(0, 100, by=5), xlab="英語得点", ylab="頻度", main="")
title(main = "英語得点のヒストグラム")	     	 #グラフタイトル
stem(eigo, scale=2)	     	 #幹葉表示

ダイズ地上部乾物重データ

　ダイズの植物体に圃場環境の不均一性がどのように現れるかの実験が、1971年に、長野県塩尻市 にある中信試験場（現在は，長野県野菜・花卉試験場）で実施された。一枚の圃場（南北約 8m、東西約 18m）で、畝は南北に作り、畝間は約75cm、株間は約15cmの間隔に作られた植穴に2粒ずつ播種し、間引きと補植で1本仕立てにした。各個体は、根際で切断され、地上部を乾燥した後、地上部全重(g)、粒重(g)、主茎長(cm)、分枝数(本)の4形質が測定された。
　今回使用するデータは、その北西部分にあたる12畝、25株の地上部全重データをそれらの位置に対応して並べたものである。すなわち、データ行列の左側と上側は道路に面しているが、右側と下側は他のダイズ個体に面していた。

データダウンロード
#ダイズデータ読み込み
daizu.df <- read.csv("daizu.csv")     # daizu.csv ファイル読み込み
dim(daizu.df)                # データの大きさ
head(daizu.df)               # 最初の一部の行のみを表示
#daizu.dfはdata.frameなので行列に変換する。
daizu.mat <- as.matrix(daizu.df)
# その後ベクトル化して1次元数値データにする。
daizu <- as.vector(daizu.mat)
length(daizu)
summary(daizu)
#データに欠測値、NA(Not Avairable) が含まれていると計算してくれないので、
#NA部分を除く(remove)という命令、na.rm=Tを付け加える必要がある。
mean(daizu)
mean(daizu, na.rm=T)
sd(daizu, na.rm=T)
#ダイズデータをグラフにする。
boxplot(daizu)
hist(daizu, breaks=seq(0, 200, by=5))
　次に、ダイズデータの行平均や列平均を求めてみる。これは、apply() 関数を用いる。mean の他に max, sum, sd, など色々な操作もできる。na.rm=T も忘れずに入れる。 
rowmean <- apply(daizu.df, 1, mean, na.rm=T); rowmean		# 行平均
colmean <- apply(daizu.df, 2, mean, na.rm=T); colmean		# 列平均
plot(rowmean, type="b")
plot(colmean, type="b")
　左側ほど値が大きい、理由は後ほど考える。

主要な確率分布

１．確率分布

確率変数

　確率変数(random variable) X とは、起こりうる事象に数値 x を割り当て、その数値 x が起こる確率を Pr[X = x ]と表す。

離散型確率分布

　非負の整数 x = 0, 1, 2, … に対し、ある数値 p0, p1, p2, … が対応し、
fig1
であるとき、X = {0, 1, 2, … }$ を離散型確率変数 (discrete random variable)、 P = {p0, p1, p2, … } を離散型確率分布 (discrete probability distribution) といい、各 pi を確率密度（probability density）という。 すなわち、離散型確率変数は、

　X　 　0 　1 　2 　…
　P　 　p0　 　p1　 　p2　 　…　

のような表で表せる。たとえば、コイン投げにおいて表を0、裏を1としたならば、その確率分布は、

　X　 　0 　1
　P　 　1/2　 　1/2　

となる。

連続型確率分布

　関数 f(x) が，
fig2
という性質を持つとき，これを連続型確率分布（continuous distribution）という． また，f(x) を確率密度関数（probability density function）とも呼ぶ．なお， 確率密度関数 f(x) が， パラメータ θ を持つとき，f(x；θ) と表記することもある．

累積分布関数

　Pr[A] を事象 A が生起する確率とする．連続型確率変数 X に対し，X が x 以下である確率，
fig3
で定義される関数 F(x) を累積分布関数（cumulative distribution function）という． この関数を用いると，確率変数 X が区間 (a, b) に落ちる確率は，
fig4
で表せる． なお，離散型確率変数でも積分を和に変えることにより，同様に階段型の累積分布関数が定義できる．

２．分布の代表値

平均

　平均（mean）μ は，分布の中心的な位置（location）座標を表す．確率変数 X に対しては， X の期待値（expectation） E[X ] とも表記し，
fig5
と定義される．離散確率変数では，積分を総和に変えることにより括弧内のように定義できる．

分散

　分散（variance） σ2 は，分布の拡がり（dispersion）の程度を表す． 確率変数 X に対しては，Var[X ] と表記する． 確率変数 X の関数 (X - μ)2 の期待値でもあり，
fig6
と定義される．
　また，σ を標準偏差（SD : Standard Deviation）といい，平均と同じ次元で分布の拡がりの大きさ を表す量である．

独立な確率変数の平均と分散

　２つの独立な確率変数 X，Y の平均と分散がそれぞれ，
E[X ] = μx，Var[X ] = σx2， E[Y ] = μy，Var[Y ] = σy2，
であるとき，
E[X + Y ] = E[X ] + E[Y ] = μx + μy， Var[X + Y ] = Var[X ] + Var[Y ] = σx2 + σy2
E[X - Y ] = E[X ] - E[Y ] = μx - μy， Var[X - Y ] = Var[X ] + Var[Y ] = σx2 + σy2
である．一般に，スカラー a，b に対して
E[aX + bY ] = aE[X ] + bE[Y ] = aμx + bμy， Var[aX + bY ] = a2Var[X ] + b2Var[Y ] = a2σx2 + b2σy2
である．

離散型分布(discrete distribution)

１．二項分布（Binomial distribution）

　独立な n 回のベルヌイ試行（成功確率 p，失敗確率 q = 1 - p）を行ったときの成功回数 X の分布．
　その確率密度は，
fig13
である．X ～ B(n，p) と書くこともある．
平均：　E[X ] = np，分散： Var[X ] = np(1 - p) = npq
これより，成功確率 p の推定値 X/n の平均と分散は，
E[X/n] = E[X ]/n = np/n = p， Var[X/n] = Var[X ]/n2 = np(1 - p)/n2 = p(1 - p)/n

２．超幾何分布（Hypergeometric distribution)

　m 個の白石と n 個の黒石が入った袋から k 個の石を無作為に取り出したとき， 白石の個数 X の従う確率密度は，
fig15
で与えられる．
平均：E[X ] = kp，分散：Var[X ] = kp(1 - p)(m + n - k)，p = m/(m + n)

３．ポアソン分布（Poisson distribution）

　正のパラメータ λ と，0 以上の整数 X に対し，確率密度が
fig16
となる分布．
平均：E[X ] = λ，分散：Var[X ] = λ

４．負の二項分布（negative binomial distribution）

　成功確率が p で，正のパラメータ n をもち，0 以上の整数 X に対し，確率密度が
fig17
となる分布．Γ 関数で確率密度を定義すれば，n は必ずしも整数である必要はない．
　成功確率 p のベルヌイ試行において，n 回の成功が起こるまでの失敗の回数 x の分布である。

平均：E[X ] = n(1 - p)/p， 分散：Var[X ] = n(1 - p)/p2

連続型分布 (continuous distribution)

１．一様分布（uniform distribution)

　2 つのパラメータ a，b（a＜b）をもつ確率密度関数が
dunif
で表される分布．

平均： E[X ] = (a + b)/2， 分散： Var[X ] = (b - a)2/12

２．正規分布（normal distribution)

　平均 μ，分散 σ2 の２つのパラメータをもつ確率密度関数が
dnorm
で表される分布で，N(μ，σ2）と表記する．μ は位置パラメータ（location parameter）で， スケールパラメータ σ を標準偏差（standard deviation）という．

３．χ2（カイ２乗）分布（chi-squared distribution)

　正の自由度パラメータ n をもつ確率密度関数が
dchi
で表される分布．ガンマ（Γ）分布で，シェープパラメータを a = n/2， スケールパラメータを s = 2， とおいた分布．
平均：E[X ] = as = n/2×2 = n， 分散：Var[X ] = as2 = n/2×4 = 2n．

４．F 分布（F distribution)

　正の２つの自由度パラメータ m，n をもつ確率密度関数が
dfdis
で表される分布．

５．t 分布（t distribution)

　正の自由度パラメータ n をもつ確率密度関数が
dtdis
で表される分布．標準正規分布より裾が重く（x が 0 より離れてもなかなか確率密度が 0 に近づかない）， 自由度が小さいほど裾が重くなる． 自由度 n = 1 のときはコーシー分布になり，n = ∞ のときは標準正規分布となる．

参考文献（古い順）

  1. Introduction to the Theory of Statistics, Mood, A. M., Graubill, F. A. & Boes, D. C., 1974, McGRAW-HILL
  2. 「実験」生産環境生物学，東京大学大学院農学生命科学研究科生産・環境生物学専攻編，1999，朝倉書店
  3. 工学のためのデータサイエンス入門－フリーな統計環境Rを用いたデータ解析－，間瀬茂ら，2004， 数理工学社
  4. 実践生物統計学－分子から生態まで－（第 1 章，第 2 章）， 東京大学生物測定学研究室編（大森宏ら）， 2004，朝倉書店
  5. The R Tips データ解析環境 R の基本技・グラフィックス活用集，船尾暢男，2005，九天社
  6. R で学ぶデータマインニング I －データ解析の視点から－，熊谷悦生・船尾暢男，2007，九天社
  7. R で学ぶデータマインニング II －シミュレーションの視点から－，熊谷悦生・船尾暢男，2007，九天社
Copyright (C) 2008, Hiroshi Omori. 最終更新日：2018年10月13日