理学療法士がまとめたノート

統計学備忘録 since2016

Rを使って統計学を勉強するブログです

データフレームからの抽出

例)iris アヤメ

head(iris,5)
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1           5.1         3.5          1.4         0.2  setosa
2           4.9         3.0          1.4         0.2  setosa
3           4.7         3.2          1.3         0.2  setosa
4           4.6         3.1          1.5         0.2  setosa
5           5.0         3.6          1.4         0.2  setosa

Sepal.Length(へたの長さ)とSpecies(種)
2列のみを抽出したフレームにします

irisSL <- subset(iris,select = c(Sepal.Length, Species))
head(irisPW,5)

  Sepal.Length Species
1          5.1  setosa
2          4.9  setosa
3          4.7  setosa
4          4.6  setosa
5          5.0  setosa

つぎに列を種(3種類)にしてへたの長さを提示します

iun2 <- unstack(irisSL,formula = Petal.Width ~ Species)
head(iun2,10)

  setosa versicolor virginica
1    5.1        7.0       6.3
2    4.9        6.4       5.8
3    4.7        6.9       7.1
4    4.6        5.5       6.3
5    5.0        6.5       6.5

Sepal.Length(へたの長さ)とSpecies(種)の2列に戻します

x<-c(iun2$setosa,iun2$versicolor,iun2$virginica)
Species <- c(rep("setosa",50),rep("versicolor",50),rep("virginica",50))
xy <- data.frame(Sepal.Length=x, Species)

  Sepal.Length Species
1          5.1  setosa
2          4.9  setosa
3          4.7  setosa
4          4.6  setosa
5          5.0  setosa
6          5.4  setosa


irisのvirginicaだけを抜き出す方法

ise <- split(iris,iris[5])     #irisの5列目の項目で分割
head(ise$virginica)

    Sepal.Length Sepal.Width Petal.Length Petal.Width   Species
101          6.3         3.3          6.0         2.5 virginica
102          5.8         2.7          5.1         1.9 virginica
103          7.1         3.0          5.9         2.1 virginica
104          6.3         2.9          5.6         1.8 virginica
105          6.5         3.0          5.8         2.2 virginica
106          7.6         3.0          6.6         2.1 virginica



Sepal.Lengthの順に並べ換えます

Sepal.Lengthの昇順に並び替え
ise <- split(iris,iris[5])     #irisの5列目の項目で分割
isev <- ise$virginica
isevo <- isev[order(isev$Sepal.Length),]
head(isevo)

    Sepal.Length Sepal.Width Petal.Length Petal.Width   Species
107          4.9         2.5          4.5         1.7 virginica
122          5.6         2.8          4.9         2.0 virginica
114          5.7         2.5          5.0         2.0 virginica
102          5.8         2.7          5.1         1.9 virginica
115          5.8         2.8          5.1         2.4 virginica
143          5.8         2.7          5.1         1.9 virginica
````

二項分布の最尤推定法 

  母集団Aの成功率0.7をパラメータとして、この母集団のサンプルから考えていきます.成功率0.7、試行回数10回の二項分布の乱数を1000個生成して母集団Aとします.

y <- rbinom(1000,10,0.7)
head(y)   #最初の6データを確認してみます
[1] 6 7 6 4 9 8  #成功率0.7なので7に近い値になっています
hist(y, xlim=c(0,10))    #ヒストグラムを作成してみます

母集団Aの成功率は0.7なので、10回試行で7回成功したデータが最も多く生成されることになります.
f:id:yoshida931:20180110135858p:plain:w400

次に、試行回数10回で成功率0.3(黒)、0.5(赤)で、成功回数0~10回の確立密度のグラフを描いてみます

x <- rep(0:10)
prob <- dbinom(x,10,0.3)
plot(prob,type="b",xaxt="n")
axis(side=1,at=1:11,labels=x)
prob2 <- dbinom(x,10,0.5)
lines(prob2,type="b",xaxt="n",col="2")    #赤色

f:id:yoshida931:20180110140135p:plain:w400
  やはりどちらも母集団Aの分布とは異なります.本来であればパラメータは未知な数値です.どうのようにしてサンプルから推定するのでしょうか?.実際に10回中7回が成功したら, 成功確立は 7 ÷ 10 = 0.7となり, 真のパラメートも0.7になりそうですが…


最尤原理
  最尤原理は、「現在起きている事象の起きる確率」が最大値となるとき、母数の推定値として「最も尤もらしい」とする考え方です.

  例)成功率0.7の母集団から10回試行した場合の成功回数のサンプルを5つ用意して、「現在起きている事象」とします.

#母集団から抽出したサンプルを仮定します
set.seed(5)
rb <- rbinom(5,10,0.7)
[1] 8 6 5 8 9                    #現在起きている事象=成功回数

  試行回数は10回なので、それぞれの事象の起きる確率Y_iは、「0.8, 0.6, 0.5, 0.8, 0.9」 となります.
  Y_1=0.8,  Y_2=0.6,  Y_3=0.5,  Y_4=0.8,  Y_5=0.9

母集団の確率密度関数 f\theta(y)=_nC_y{\theta}^y{(1-\theta)}^{n-y}

#母集団から抽出したサンプルの確率密度
#これは成功率を固定して、成功回数が変化しているので確率です
choose(10,8)*0.7^8*0.3^2
choose(10,6)*0.7^6*0.3^4
choose(10,5)*0.7^5*0.3^5
choose(10,8)*0.7^8*0.3^2
choose(10,9)*0.7^9*0.3^1

#Rの関数を使用
dbinom(rb,10,0.7)

[1] 0.2334744 0.2001209 0.1029193 0.2334744 0.1210608


母集団から抽出したサンプルの同時確率密度
サンプル(測定値)「8, 6, 5, 8, 9 」は、互いに独立しています.したがって、「8, 6, 5, 8, 9」が揃う確率は,

f\theta(y_1,y_2,...,y_n)= \prod_{i=0}^n f\theta (y_i)

という同時確率になり、最尤推定値を求めるときに必要になります.

#Rを使用して、母集団から抽出したサンプルの同時確率密度を求めてみます
prod(dbinom(rb,10,0.7))
[1] 0.0001359164
#母集団の成功率0.7が前提だったので算出可能となります


ここからが本題です

  成功回数Y_ia_1, a_2, … ,a_iが観察されたとき、その同時確率密度を母数(成功率)\thetaの関数とみて、尤度関数といいます. 尤度関数は、成功回数(y1=8,  y2=6 , y3=5 , y4=8 , y5=9 )を固定して、成功率\thetaが変化する関数ということになります.つまり、データを定数として母数を変数とした、「観測されたデータが得られる確率」を「未知のパラメータの関数」とみた関数です.

尤度関数 L(θ)=f\theta(a_1,a_2,...,a_n)= \prod_{i=0}^n f\theta (a_i)

# 測定値「8,6,5,8,9」, 母数θは未定
L(θ) = choose(10,8)*θ^8*0.3^2 ×
        choose(10,6)*θ^6*0.3^4 ×
        choose(10,5)*θ^5*0.3^5 ×
        choose(10,8)*θ^8*0.3^2 ×
        choose(10,9)*θ^9*0.3^1



最尤推定
  上述の例題の\thetaは0.7だったのですが、母集団の母数が分かっていることはほとんどありません.本来であれば不明です.したがって、成功率\thetaが不明の場合は、尤度関数から\thetaを求めます.

尤度関数 L(θ)= \prod_{i=0}^n f\theta (a_i)

  \thetaを尤度関数から推定する場合に、尤度関数を最大にする\thetaa_1, a_2, … ,a_iが観測されたときの最尤推定値といいます. したがって、尤度関数を微分して L'(θ)=0 となるθを求めます.
  とはいっても、尤度関数から直接計算するのは、なんとなく大変そうなので、対数をとって考えていきます.

二項分布の対数尤度関数
上述の例題では試行回数n=10回、成功回数a_i   (i = 1,2,3,4,5    )

f\theta(ai) =_{10}C_a{_i}{\theta}^{ai}{(1-\theta)}^{10-ai}

=log[\frac{n!}{x!(n−x)!}θ^{ai}(1−θ)^{n−x}]

=log(n!)−log(x!)−log(n−x)!+ailogθ+(10-ai)log(1−θ)

対数尤度関数を微分
=log(n!)−log(x!)−log(n−x)!をθで微分すると0

したがって f' \theta(ai) = \sum(\frac{ai}{\theta}-\frac{10-ai}{1−θ})

ai=(8,6,5,8,9)を挿入

f' \theta(ai) =\frac{8}{\theta}-\frac{10-8}{1−θ}+\frac{6}{\theta}-\frac{10-6}{1−θ}+\frac{5}{\theta}-\frac{10-5}{1−θ}+\frac{8}{\theta}-\frac{10-8}{1−θ}+\frac{9}{\theta}-\frac{10-9}{1−θ}=0
  = \frac{36-50\theta}{\theta(1-\theta)}=0
\theta=0.72
パラメータは0.72と推定することができました.

#Rを使って確認してみます
LL <- function(p) {
  prod(dbinom(x,10,p))           # 同時密度分布
}
optimize(f=LL, c(0, 1), maximum=TRUE)      #optimizeでパラメータの最適化(この場合は最大値を求める)
$maximum
[1] 0.7200026

参考)合成関数の導関数
y = f(u)  ,  u = g(x) のとき  y = f( g(x) ) となる y = f( g(x) ) 導関数は、
\frac{dy}{dx}=\frac{dy}{du}*\frac{du}{dx}


参考web
norimune.net
www.yasuhisay.info

参考文献
柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010
東京大学教養学部統計学教室 (編集); 統計学入門 (基礎統計学), 東京大学出版会, 1991

ベイズの定理

ベイズの定理とは「観測値(データ)で条件付けられた、母数の分布を与える定理 」です

#Rのサンプル women$height を使用します
x <- women$height
#標本平均
me <- mean(women$height)
=65
#標本標準偏差
sd <- sqrt(sum((x - mean(x))^2) / length(x))
=4.32

  女性の伸長(インチ)の平均65、標準偏差4.32を正規分布の母数と見立ててしまえば、95%予測区間を導き出すことは簡単です.
  しかし、x~N(65,4.3)の分布が真の分布とは言えません.x~N(64,4.5)やx~N(66,4.6)などの母分散も考えられます.つまり、母数(平均や標準偏差)も分布することが考えられます.

確率密度関数 f(母数|観測値) の分布

例として母数を平均値, 標準偏差とします

f (θ|x) = f ( 平均値, 標準偏差 | x )

平均値, 標準偏差 とデータxは独立していないので条件付き分布となります

条件付き分布

f (x,θ) ≠ f(x)*f(θ)
# x が与えられた場合
  f (x,θ) = f(x)*f (θ|x)
# θ が与えられた場合
  f (x,θ) = f (x|θ)*f(θ)   


分布に関するベイズの定理

観測値(データ)で条件付けられた、母数の分布を与える定理

f (x,θ) = f(x)*f (θ|x)より
f (θ|x) = f (x,θ) ÷ f(x)
        = f (x|θ)*f(θ) ÷ f(x)

上記式より
f (θ|x) =事後分布
f(x|θ) =尤度
f(θ) =事前分布
f(x) =正規化定数



確率に関するベイズの定理

Aを得られた結果とします.Bをその原因とします.
例)5つの袋に赤白の玉が入り混ざっています.どれかの袋から玉が取られたとします.

取り出した玉=結果A
どの袋から取り出したか=原因B

 結果Aが得られた原因をBを推定します.袋が5つなのでB1, B2, B3, B4, B5 が原因となります. ベイズの定理では、結果Aであったときに原因が、B1またはB2またはB3またはB4またはB5 である確率を算出することになります.

ベイズの定理
P( Bi | A ) = \frac{P(A)P(A|Bi)}{P(A)} = \frac{P(A)P(A|Bi)}{ΣP(Bi)P(A|Bi)}

証明

A = A∩Ω =A∩(B1∪B2∪…∪Bi)
=(A∩B1)∪(A∩B2)∪...∪(A∩B1)
 P(A) = ΣP(Bi)P(A|Bi)

例)
袋B1または袋B2から一つの球を取り出します.
B1には白玉3個、赤玉2個
B2には白玉1個、赤玉3個
取り出した球が赤だった場合(事象A)、袋がB1である確率、B2である確率を求めよ.
f:id:yoshida931:20180117184454p:plain:w300
P(赤)=\frac{1}{2}*\frac{2}{5}+ \frac{1}{2}*\frac{3}{4}

# 袋がB1である確率   
  P(B1|A)= (1/2)*(2/5) / ((1/2)*(2/5) + (1/2)*(3/4)) = 0.3478261
# 袋がB2である確率
  P(B2|A)= (1/2)*(3/4) / ((1/2)*(2/5) + (1/2)*(3/4)) = 0.6521739

参考文献
豊田秀樹; はじめての統計データ分析 ベイズ的<ポストp値時代>の統計学, 朝倉書店, 2016

このブログの参考・引用文献

投稿日 2016-10-31
最終更新日 2017-12-29

文献:著者五十音順 巨人の肩に乗らねば・・・

石井一夫 ;Rとグラフで実感する生命科学のための統計入門, 羊土社 ,2017
石田 基広 ;改訂3版 R言語逆引きハンドブック ,シーアンドアール研究所; 改訂3版,2016

石田 基広 ;Rで学ぶデータ・プログラミング入門 ―RStudioを活用する―,共立出版 ,2012
石村卓夫;すぐわかる多変量解析,東京図書,1997
石村卓夫;入門はじめての時系列分析,東京図書,2012
石村卓夫;分散分析のはなし,東京図書,1992
稲葉 由之;プレステップ統計学I 記述統計学, 弘文堂, 2012
稲葉 由之;プレステップ統計学II 推測統計学, 弘文堂, 2013
北川 源四郎; 時系列解析入門, 岩波書店, 2005
金 明哲; Rによるデータサイエンス-データ解析の基礎から最新手法まで, 森北出版,2007
向後 千春, 冨永 敦子;統計学がわかる,技術評論社,2007
向後 千春, 冨永 敦子;統計学がわかる【回帰分析・因子分析編】,技術評論社,2008
田中 孝文; Rによる時系列分析入門, シーエーピー出版, 2008
東京大学教養学部統計学教室 (編集); 統計学入門 (基礎統計学), 東京大学出版会, 1991
東京大学教養学部統計学教室 (編集); 自然科学の統計学, 東京大学出版会, 1992

対馬栄輝; よくわかる医療統計 -「なぜ?」にこたえる道しるべ, 東京図書, 2015
豊田秀樹 (著, 編集);回帰分析入門 (Rで学ぶ最新データ解析) ,東京図書 ,2012
豊田秀樹; 検定力分析入門, 東京図書,2009
豊田秀樹; はじめての統計データ分析 ベイズ的<ポストp値時代>の統計学, 朝倉書店, 2016
西内 啓; 統計学が最強の学問である, ダイヤモンド社, 2013
西内 啓; 統計学が最強の学問である[実践編]---データ分析のための思想と方法, ダイヤモンド社, 2014
日本統計学会 (編集); 日本統計学会公式認定 統計検定2級対応 統計学基礎, 東京図書, 2012
南風原 朝和; 心理統計学の基礎―統合的理解のために, 有斐閣 , 2002
南風原 朝和; 続・心理統計学の基礎--統合的理解を広げ深める, 有斐閣 , 2014
舟尾 暢男; The R Tips―データ解析環境Rの基本技・グラフィックス活用集, オーム社; 第2版,2009
村上 正康,安田 正実; 統計学演習, 培風館,1989
柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010
柳川 堯 ; 観察データの多変量解析―疫学データの因果分析,近代科学社 ,2016

山田 剛史, 杉澤 武俊, 村井 潤一郎; Rによるやさしい統計学, オーム社 , 2008

Hadley Wickham (著),石田 基広 (翻訳) ;R言語徹底解説, 共立出版,2016
R.A. フィッシャー(原著)、遠藤 健児、鍋谷 清治(訳); 実験計画法, 森北出版株式会社, 2013
R.A. フィッシャー(原著)、遠藤 健児、鍋谷 清治(訳); 研究者のための統計的方法, 森北出版株式会社, 2013
Peter Dalgaard (著), 岡田 昌史 (監修, 翻訳) ;Rによる医療統計学 原書2版,丸善,2007


ホームページ

おしゃべりな部屋 http://aoki2.si.gunma-u.ac.jp/
竹中明夫のページ  (http://takenaka-akio.org/index.html)

梶山喜一郎;コピペで学ぶ Rでテクニカルデータプレゼンテーション      (http://monge.tec.fukuoka-u.ac.jp/R_analysis/0r_analysis.html
間瀬 茂;R 基本統計関数マニュアル(https://cran.r-project.org/doc/contrib/manuals-jp/Mase-Rstatman.pdf)
続・わしのページ(R-Tips_舟尾 暢男 http://nfunao.web.fc2.com/
Quick-R (http://www.statmethods.net/)
R-Tips (http://cse.naro.affrc.go.jp/takezawa/r-tips/r2.html)
RjpWiki (http://www.okadajp.org/RWiki/)

一様分布

確立密度関数
f(x|a,b)=\frac{1}{b-a} ,  a≦x≦b

積分布関数(分布関数)
F(x|a,b)=\frac{x-a}{b-a}  ,  F(b|a,b)=1

平均
\frac{b+a}{2}

分散
\frac{(b-a)^2}{12}

#例)離散型一様分布(確立変数が整数の場合)
#1~10の整数をランダムで発生させたとき、2~5の数が出る確率は?
(5-1)/(10-1) - (2-1)/(10-1) 
= (5-2)/(10-1)
= 0.33333


サイコロを例にグラフを作成してみます

#注)出力される乱数には引数の最大値(7)は含みません(引数の最大値は7とします)
#試しに20個生成してみます・・・6以上の実数が含まれますが7はふくまれません
 runif(20,1,7)
 [1] 5.684374 3.100493 1.276716 5.788635 1.131957 5.964903 6.869709 4.580407
 [9] 6.930938 2.175627 1.485392 1.791500 2.468878 3.537859 6.206818 2.513363
[17] 5.314561 3.055166 5.022207 2.632470

#次に乱数を6個出力して確率密度関数をグラフにしてみます
runif(6,1,7) 
   3.996601 5.632944 6.618261 1.698278 1.169685 1.670989
dunif(6,1,7) #確率密度=1÷(7-1)
 = 0.1666667
plot(dunif(1:6,1,7)) #密度関数のグラフ(確立密度を6個なので1:6とします)

f:id:yoshida931:20171228103553p:plain:w300

正規分布

確立密度関数
f(x|μ,σ)=\frac{1}{\sqrt{2π}σ}exp[\frac{-1}{2σ^2}(x-μ)^2] 、 -\infty\leqq\ x\leqq\infty

積分布関数(分布関数)
F(μ+1.96σ|μ,σ)-F(μ-1.96σ|μ,σ)\simeq0.95

pnorm ( 1.96, lower.tail=TRUE ) - pnorm ( 1.96, lower.tail=FALSE )  
= pnorm ( 1.96 ) - pnorm ( -1.96 )  
= 0.9500042

yoshida931.hatenablog.com

クラスター分析

Rのサンプルattitudeを使用して、クラスタ分析(階層的方法)を勉強します.
attitudeは、管理者態度のデータです.無作為に選ばれた35名の雇用者よりアンケート.好意的な割合が数値化されています.
rating全般的評価、 complaints雇用者からの苦情処理、 privileges特別な特権の付与、 learning学習機会の付与、 raises業績による昇進機械、 critical批判的態度、 advance進歩

# xxはdata.frame

(result<-hclust(dist(xx),method="single"))
plot(result, main = "最近隣法")

(result<-hclust(dist(xx),method="complete"))
plot(result, main = "最遠隣法")

(result<-hclust(dist(xx),method="average"))
plot(result, main = "群平均法")

(result<-hclust(dist(xx),method="centroid"))
plot(result, main = "重心法")

(result<-hclust(dist(xx),method="median"))
plot(result, main = "メディアン法")

(result<-hclust(dist(xx),method="ward.D"))
plot(result, main = "ウォード法 ward.D")

(result<-hclust(dist(xx),method="ward.D2"))
plot(result, main = "ウォード法 ward.D2")

実際にウォード法を実行してみます  まずは変数クラスタを求めてみます

xx <- t(attitude)
(result<-hclust(dist(xx),method="ward.D"))
plot(result, main = "ウォード法 ward.D")

f:id:yoshida931:20171219114804p:plain:w400
進歩、批判的態度を除くと2グループに管理者の評価は分かれるようです 進歩、批判的態度、向上意欲をくみ取る能力、雇用者の全般的な評価といった分類になりそうです (これかは私の勝手な見解です)

つぎにサンプルクラスタをもとめてみます

xx <- attitude 
(result<-hclust(dist(xx),method="ward.D"))
plot(result, main = "ウォード法 ward.D")

f:id:yoshida931:20171219114910p:plain:w700
雇用者が数グループに分けられそうです.
こんなデータが病院経営戦略にも役立つのかも…?

引用・参考 間瀬 茂;R 基本統計関数マニュアル(https://cran.r-project.org/doc/contrib/manuals-jp/Mase-Rstatman.pdf)