データフレームからの抽出
例）ChickWeightからDietの1と3だけ抜き出す
subset(ChickWeight,subset = Diet==c(1,3))
例）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 ````
ベイズの定理
ベイズの定理とは「観測値（データ）で条件付けられた、母数の分布を与える定理 」です
#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(母数|観測値) の分布
例として母数を平均値, 標準偏差とします
平均値, 標準偏差 とデータ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)
上記式より
確率に関するベイズの定理
Aを得られた結果とします．Bをその原因とします．
例）5つの袋に赤白の玉が入り混ざっています．どれかの袋から玉が取られたとします．
取り出した玉＝結果A
どの袋から取り出したか＝原因B
結果Aが得られた原因Ｂを推定します．袋が5つなのでB1, B2, B3, B4, B5 が原因となります． ベイズの定理では、結果Aであったときに原因が、B1またはB2またはB3またはB4またはB5 である確率を算出することになります．
ベイズの定理
証明
A = A∩Ω
=A∩(B1∪B2∪…∪Bi)
=(A∩B1)∪(A∩B2)∪...∪(A∩B1)
∴
例）ベイズの逆確率
袋B1または袋B2から一つの球を取り出します．
B1には白玉３個、赤玉2個
B2には白玉1個、赤玉３個
取り出した球が赤だった場合(事象A)、袋がB1である確率、B2である確率を求めよ．
正規化
# 赤玉を引いた袋が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(B1|A) + P(B2|A) = 1
このブログの参考・引用文献
投稿日 2016-10-31
最終更新日 2017-12-29
文献：著者五十音順 巨人の肩に乗らねば･･･
石井一夫 ;Rとグラフで実感する生命科学のための統計入門, 羊土社 ,2017
石田 基広 ;改訂3版 R言語逆引きハンドブック ,シーアンドアール研究所; 改訂3版,2016
石田 基広 ;Rで学ぶデータ・プログラミング入門 ―RStudioを活用する―,共立出版 ,2012
石村卓夫;すぐわかる多変量解析,東京図書,1997
石村卓夫;入門はじめての時系列分析,東京図書,2012
石村卓夫;分散分析のはなし,東京図書,1992
稲葉 由之;プレステップ統計学Ｉ 記述統計学, 弘文堂, 2012
稲葉 由之;プレステップ統計学II 推測統計学, 弘文堂, 2013
北川 源四郎; 時系列解析入門, 岩波書店, 2005
金 明哲; Rによるデータサイエンス-データ解析の基礎から最新手法まで, 森北出版,2007
向後 千春, 冨永 敦子;統計学がわかる,技術評論社,2007
向後 千春, 冨永 敦子;統計学がわかる【回帰分析・因子分析編】,技術評論社,2008
田中 孝文; Rによる時系列分析入門, シーエーピー出版, 2008
東京大学教養学部統計学教室 (編集); 統計学入門 (基礎統計学), 東京大学出版会, 1991
東京大学教養学部統計学教室 (編集); 自然科学の統計学, 東京大学出版会, 1992
対馬栄輝; よくわかる医療統計 -「なぜ?」にこたえる道しるべ, 東京図書, 2015
豊田秀樹 (著, 編集);回帰分析入門 (Rで学ぶ最新データ解析) ,東京図書 ,2012
豊田秀樹; 検定力分析入門, 東京図書,2009
豊田秀樹; はじめての統計データ分析 ベイズ的＜ポストp値時代＞の統計学, 朝倉書店, 2016
西内 啓; 統計学が最強の学問である, ダイヤモンド社, 2013
西内 啓; 統計学が最強の学問である[実践編]---データ分析のための思想と方法, ダイヤモンド社, 2014
日本統計学会 (編集); 日本統計学会公式認定 統計検定２級対応 統計学基礎, 東京図書, 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/)
クラスター分析
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")
進歩、批判的態度を除くと2グループに管理者の評価は分かれるようです
進歩、批判的態度、向上意欲をくみ取る能力、雇用者の全般的な評価といった分類になりそうです （これかは私の勝手な見解です）
つぎにサンプルクラスタをもとめてみます
xx <- attitude (result<-hclust(dist(xx),method="ward.D")) plot(result, main = "ウォード法 ward.D")
雇用者が数グループに分けられそうです．
こんなデータが病院経営戦略にも役立つのかも…？
引用・参考 間瀬 茂；R 基本統計関数マニュアル(https://cran.r-project.org/doc/contrib/manuals-jp/Mase-Rstatman.pdf)
シンプソンのパラドクス
最近、統計学の学習が進んでおりません
ちょっとお仕事が忙しくて…
今日は少しだけ勉強しておきます．
観察研究には重要な概念です．
x <- c(110,70,90,120) x <- matrix(x,2,2) rownames(x)<-c("治療A","治療B") colnames(x)<-c("改善","変化なし") addmargins(x) 改善 変化なし Sum 治療A 110 90 200 治療B 70 120 190 Sum 180 210 390
これだけ見て、「治療Aが良い！」と結論付けるのは早いのです．
患者さんの症状別にみた次の表を見てみましょう！
y <- c(3,10,40,93,17,31,42,20,90,29,8,7) dat <- array(y, dim = c(2,2,3)) name <- list("治療"=c("A","B"),"変化"=c("改善","変化なし"),"患者"=c("重度","中等度","軽度")) dimnames(dat)<-name dat , , 患者 = 重度 変化 治療 改善 変化なし A 3 40 B 10 93 , , 患者 = 中等度 変化 治療 改善 変化なし A 17 42 B 31 20 , , 患者 = 軽度 変化 治療 改善 変化なし A 90 8 B 29 7 #患者の重症度によって分類してみたら･･･必ずしも治療Aが良いとは言えません． #症状が重度、または中等度の患者さんには治療Aはあまり効果がないようです． #症状が中等度の患者さんには治療Bが効果的のようです． #観察研究では交絡因子の影響をブロックすることが重要になります．
今日は、久しぶりにRを操作しました．
上記のようなクロス表が作れたら、色々勉強できますよ！
医療関係者には最適な医療を提供してもらいたいですね！
