2012年度 バイオスタティスティクス基礎論 テキスト1
Introduction to Biostatistics 2012, Text 1
西田洋巳Hiromi Nishida
講義・実習日:2012年4月11日（水）
[1-1] 推薦書Recommendation books
R による医療統計学
Introductory Statistics with R, Springer-Verlag
Peter Dalgaard著、岡田昌史監訳、丸善
統計学:R を用いた入門書
Statistics: An Introduction using R, John Wiley & Sons
Michael J. Crawley著、野間口謙太郎・菊池泰樹訳、共立出版
[1-2] Rのダウンロード、起動
ダウンロードはCRANより行う
Rを起動させる
ディレクトリの変更が可能
終了するときはq()
Rでは大文字と小文字を区別している
Rではスペースは無視される
[1-3] ベクタVector
各要素を連結Concatenate
数字の場合
c(1, 2, 3)
文字の場合
c("A", "B", "C")
なお、a, b, c, ..はlettersに、A, B, C, ..はLETTERSに収納されているので
LETTERS[1:3]
もc("A", "B", "C")に同じである
1:3は1, 2, 3を示し、[1:3]は1番目から3番目までを抽出することを意味している
アルファベットの中でTはTruth、FはFalseを示す
例えば
c(T, T, F)
要素が文字の場合には"T"のようにする必要がある
数値の連続Sequence
seq(1, 10)
は1から10まで1間隔で示す
seq(1, 10, 2)
は1から10まで2間隔で示す
1間隔の場合には
1:10
とすることができる
繰り返しReplicate
rep(c(1, 2, 3), 3)
はc(1, 2, 3)を3回繰り返すことである
rep(c(1, 2, 3), 1:3)
はc(1, 2, 3)の1番目を1回、2番目を2回、3番目を3回繰り返すことである
rep(1:3, c(3, 2, 1))
は1:3の1番目を3回、2番目を2回、3番目を1回繰り返すことである
[1-4] 算術演算Arithmetic operations
a <- c(1, 2, 3)
はaにc(1, 2, 3)を収納することである
a = c(1, 2, 3)
でも収納することができる
a+a
は各要素の同じ位置における和を示す
a-a
は各要素の同じ位置における差を示す
a*a
は各要素の同じ位置における積を示す
a/a
は各要素の同じ位置における商を示す
c(a, a)
はaとaの連結を示し、和ではない
a+a; a-a
;を使うと区切ることができる
sqrt(a) #square root
は平方根を意味し、#を用いると切断できる
log(a)
は自然対数である
log2(a)
は底が2の対数である
a > 2
はaの中で2より大きいものを判別する
a >= 2
はaの中で2以上のものを判別する
a == 2
は2と同じものを判別する
a != 2
は2ではないものを判別する
aから特定条件を満たすものを抽出するときは
a[a > 2]
のように[ ]を使う（LETTERSを参照）
[1-5] 平均、分散Mean, Variance
a <- c(1, 20, 300, 4000)
a[1]-a[4]
はaの1番目の要素から4番目の要素を引くことを意味している
abs(a[1]-a[4])
は絶対値Absolute valueを示す
sum(a)
はaの要素の総和Sum of elementsを示す
length(a)
はaのサンプルサイズSample sizeを示す
sum(a)/length(a)
は総和/サンプルサイズであり、平均を意味している
mean(a)
も平均を示す
median(a)
は中央値を示す
要素数が偶数の場合には上位半分の底と下位半分の頂の平均を示す
aの場合には(20+300)/2 = 160
max(a)
は最大値を示す
min(a)
は最小値を示す
range(a)
はaの最小値と最大値を示す
sum((a-mean(a))^2)/(length(a)-1)
は平方和/自由度Sum of squares/degrees of freedomであり、不偏分散を意味している
var(a)
も不偏分散Varianceを示す
sqrt(var(a))
は不偏分散の平方根であり、標準偏差を意味している
sd(a)
も標準偏差Standard deviationを示す
[1-6] 正規分布Normal distribution
f(x) = (1 / sqrt(2π(variance)))exp(-(x – mean)2/2(variance))
平均が0、分散が1の場合を標準(基準)正規分布Standard normal distributionと呼ぶ
f(x) = (1 / sqrt(2π))exp(-x2/2)
seikibunpu <- function(x) 1/sqrt(2*pi)*exp(-x^2/2)
function()を使うことのより関数を定義することができる
a <- seq(-5, 5, 0.1)
b <- seikibunpu(a)
plot(a, b)
なお、本関数はdnormとして入っている
par(new=T)
は次の図を重ねて示す場合Overwritingに使う
plot(a, dnorm(a), type="l", col="red")
プロットのタイプは他にも"b", "c", "s", "S", "h"などがある
色については
colors()
を参照すると良い
a <- rnorm(1000)
は標準正規分布に従う乱数1000個を発生させaに収納することを意味している
rnorm(1000, mean=0, sd=1)としてもよい
mean(a); sd(a)
aの平均と標準偏差
hist(a)
はaの分布をヒストグラムHistogramで表す
hist(a, breaks=seq(-4, 4, 0.2), main="rnorm", xlab="")
はヒストグラムを-4から4の範囲で0.2間隔表記、タイトルはrnorm、X軸ラベルはなしで表す
Separate at 0.2 intervals between -4 and 4; main, title of the histogram; x-axis, no label
hist(a, prob=T)
はY軸を頻度から確率密度への変更を示すChange the y-axis to the density from the frequency
density(a)
は密度評価Density estimationを意味している
lines(density(a))
は近似曲線Fitted curveを描く
lines(density(a, bw=0.5), col="red")
lines(density(a, bw=0.1), col="blue")
などと条件を変えて見ることができる
quantile(a)
はクォンタイル、四分位を示す
summary(a)
はaの分布の様子の要約を示す
quantile(a, prob=c(0.025, 0.975))
はaの分布における上位2.5%および下位2.5%の値を示す
なお、正規分布における理論値Theoretical valueはqnorm(c(0.025, 0.975))により得られる
qqnorm(a)
はQ-Qプロットを示し、正規分布にaが従っていると直線に乗る
分布に対しての最初の文字については
d, density
p, probability
q, quantile
r, random
をそれぞれ意味している
[1-7] ブートストラップBootstrapping
Bradley Efron (1938-)により提案された
例えば、下記aの標本の大きさは5であることから、5回の復元無作為抽出Random sampling with replacementにより一つの仮想データを作成する。この作業を繰り返す。
a <- c(2, 3, 6, 18, 108)
これがオリジナルデータとして
d <- c() #またはd <- numeric()
dをここで定義しておく
for(i in 1:1000) d[i] <- mean(sample(a, 5, replace=T))
は1000の仮想データを生成し、個々の平均を求め、dに収納することを意味しているGenerate 1000 virtual data
hist(d)
によりdのヒストグラムを作成する
mean(d)
これは仮想データの平均値1000個の平均値を示す
この値を実際のデータの平均を比較するCompare it with mean(a)
ブートストラップ解析は系統樹における分岐の評価に使われているIn phylogenetics, the bootstrapping has been used frequently.
Felsenstein J (1985) Confidence-limits on phylogenies – an approach using the bootstrap. Evolution vol. 39, pp. 783-791
Efron B et al. (1996) Bootstrap confidence levels for phylogenetic trees. Proc Natl Acad Sci USA vol. 93, pp. 13429-13434
[1-8] 中心極限定理Central limit theorem
サンプル平均値の分布はサンプリングを繰り返すことにより正規分布に近づく
When we obtained a mean from a sample, repeated it many times, the distribution of those means tend to follow the normal distribution.
母集団が一様分布でサンプルサイズ1、サンプリング回数1000の場合
Population, uniform distribution
Sample size, 1
Sampling number, 1000 times
a1 <- c()
for(i in 1:1000) a1[i] <- mean(runif(1))
1000のデータをa1に収納
母集団が一様分布でサンプルサイズ2、サンプリング回数1000の場合
Population, uniform distribution
Sample size, 2
Sampling number, 1000 times
a2 <- c()
for(i in 1:1000) a2[i] <- mean(runif(2))
a2に収納
母集団が一様分布でサンプルサイズ10、サンプリング回数1000の場合
Population, uniform distribution
Sample size, 10
Sampling number, 1000 times
a10 <- c()
for(i in 1:1000) a10[i] <- mean(runif(10))
a10に収納
a1, a2, a10の分布を比較するCompare the distributions of a1, a2, and a10
par(mfrow = c(3, 1))
は3つの図を3行1列で表記することを意味しているDisplay three figures (3 rows, 1 column)
hist(a1)
hist(a2)
hist(a10)
Q-Qプロットで比較すると
par(mfrow = c(1, 1))
qqnorm(a1, ylim=c(0, 1))
par(new=T)
qqnorm(a2, ylim=c(0, 1), col="blue")
par(new=T)
qqnorm(a10, ylim=c(0, 1), col="red")
直線に乗っている方が正規分布に近いことを示す
[1-9] 集合、並び替えなどUnion, Intersect, Sort, etc.
a <- runif(100, 1, 100)
は1から100の範囲で一様分布に従う乱数を100個生成し、aに収納する
100 random numbers from a uniform distribution (in the range from 1 to 100)
b <- runif(100, 1, 100)
aと同様の内容であるが、中身は違う
union(a, b)
はaとbの結合集合を意味している
intersect(a, b)
は共通集合を意味している
a <- trunc(a); b <- trunc(b)
は要素の整数部分Integers of elementsだけを要素とすることを意味している
union(a, b); intersect(a, b)
を先ほどの結果と比較する
特定の要素を抽出するには
a[30]; a[101]
はaの30番目、aの101番目の抽出を意味しており、ここでは101番目は存在していない
NAは欠損Not Availableを意味しており、NaN非数Not a Numberと間違わないように
1/0; -1/0
は無限大Infinityとなる
a[a>20]; a[a>20 & a<=30]
は20より大きな要素の抽出、20より大きく30以下の要素の抽出を意味している
a[1] <- 30
はaの1番目を30に換えることを意味しているChange the first element to 30
rev(a)
は反転Reverseを示す
order(a)
はaの中で小さい値から何番目に存在しているかを示す、デフォルトDefaultはdecreasing=F
sort(a)
は要素の値が小さいものから並び換える、a[order(a)]に同じ
[1-10] 行列Matrix
matrix(c(1, 2, 3, 4, 5, 6), 2, 3)
は要素1, 2, 3, 4, 5, 6を2行3列に行列させる
matrix(c(1, 2, 3, 4, 5, 6), 2)
も同じこと
matrix(c(1, 2, 3, 4, 5, 6), nrow=2)
も同じこと
matrix(c(1, 2, 3, 4, 5, 6), ncol=3)
も同じこと
matrix(c(1, 2, 3, 4, 5, 6), 2, byrow=T)
は要素を並べる際に行を埋めて次の行へということを示す
a <- c(1, 2, 3, 4, 5)
b <- c(2, 4, 6, 8, 10)
ab <- data.frame(a, b)
はベクトルaとbをデータフレーム化（aとbは列となる）
ab[ ,1]
はデータフレームabの1列目を抽出する
ab[ ,2]
は2列目を抽出する
ab[1, ]
は1行目を抽出する
ab[3, ]
は3行目を抽出する
a <- matrix(1:6, 2, byrow=T)
rownames(a) <- c("A", "B")
は1行目をA、2行目をBと名づける
colnames(a) <- letters[1:3]
は1列目、2列目、3列目をa, b, cと名づける
t(a)
は転置関数Transposition functionを意味する
cbind(A=1:3, B=4:6)
は列で結合する
rbind(A=1:3, B=4:6)
は行で結合する
a <- matrix(1:10, ncol=2)
a[ ,1]
は1列目を抽出する
a[ ,1]>3
は1列目の要素で3より大きいものを示す
a[a[ ,1]>3, ]
は1列目の要素で3より大きいものを持つ行を抽出する
a[a[ ,1]>3, ][ ,2]
は1列目の要素で3より大きいものを持つ行を抽出し、その2列目の要素を抽出する
a[a[,1]>3, 2]に同じ
apply(a, 1, mean)
は各行の平均Mean of rowを示す
apply(a, 2, median)
は各列の中央値Median of columnを示す
apply(a, c(1,2), sqrt)
はすべての要素の平方根を示す、sqrt(a)でもできる