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

統計学備忘録 since2016

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

二元配置分散分析 (1)

二元配置分散分析 投稿日2017.3.8更新日2017.8.17 2要因ともに対応なし、繰り返し数がすべて等しい場合 例) 要因1(水準A1、水準A2)、要因2(水準B1、水準B2、水準B3) (仮想データ使用) Rの関数で分散分析表を作成してみます 各セルにベクトル名をつけ…

一元配置分散分析

一元配置分散分析 投稿日2016.11.2 更新日2017.8.15 言葉の整理 要因:実験結果に影響を与える要素 因子:研究対象となる要因( 母平均に差をもたらすと考えられる ) 水準:要因を分ける条件( 要因に含まれる項目 ) 群:=水準 平方和:偏差の二乗の合計 …

標準偏回帰係数の求め方

標準偏回帰係数の求め方 投稿日2017.8.12 更新日2017.8.13 忘れないように載せておきます 例)予測変数が二つの場合の重回帰分析 母回帰方程式 Y = b0 + b1*X1i + b2*X2i + Ei 予測方程式 y = β0 + β1*x1i + β2*x2i +ei 標準偏回帰係数 sβ1、sβ2 sβ1 = ( yと…

トービット回帰直線

トービット回帰直線 投稿日2017.8.4 打ち切りデータの場合、つまり天井効果や床効果が生じているデータの場合には、トービット回帰直線で分析します. 参考書はもちろん豊田秀樹 (著, 編集);回帰分析入門 (Rで学ぶ最新データ解析) ,東京図書 ,2012 Rのサンプ…

残差分析

残差分析 Rのサンプル ChickWeight を使用y<-ChickWeight$weight[1:10]x<-ChickWeight$Time[1:10]summary(lm(y~x))plot(x,y,xlab = "生後日数",ylab = "体重")abline(30.327 ,7.030) #回帰直線の挿入 残差 ei をプロットすることで、回帰モデルからのズレを…

単回帰分析

投稿日2017.2.2更新日2017.8.1線形回帰分析説明変数と目的変数を直線関係で傾向を示す。説明変数と目的変数との関係を直線でモデル化する回帰分析。 非線形回帰分析非線形関係でモデル化する回帰分析。説明変数と目的変数を非直線的関係で傾向を示す。 以下…

感度、特異度、ROC曲線

感度と特異度 投稿日2017.7.31 感 度(陽性反応的中度)= b / ( a + b ) 疾患に罹患していて、検査も陽性 特異度(陰性反応的中度)= c / ( c + d ) 疾患に罹患してなくて、検査も陰性 偽陽性率 a / ( a + b ) 疾患に罹患していて検査が陰性偽陰性率 d / (…

回帰係数の区間推定

投稿日2017.2.3更新日2017.7.31 回帰係数の区間推定 標本回帰係数の不偏推定量から母回帰係数の信頼区間を求めてみます. 次式をRにペーストしてx<-read.table("clipboard",header=T) 次のデータをコピーしてRで実行します東京 福岡1019.4 1018.41005.7 1007…

ユールのQ

ユール ( Yule ) のQユールの関連係数、ユールの連関係数とも呼ばれています.関連が強いほど1または-1に近い値をとります。 Q= ( a*d - b*c ) / ( a*d + b*c ) = ( オッズ比 - 1 ) / ( オッズ比 + 1 )例)2つの質問 ( q1 , q2 ) に関する答え ( yes , no ) …

確率点と累積分布

確率点と累積分布 上側確率が0.025の場合 lower.tail logical; if TRUE (default), probabilities are P[X ≤ x] otherwise, P[X > x]. Rのhelpより 正規分布の確率点(quantile) qnorm( 確率 ) 上側確率が0.025になるZ値( p=0.975とおくことで上側確率は0.02…

適合度の検定 (ピアソンの 𝜒^2 適合度検定)

適合度の検定 (ピアソンの 𝜒^2 適合度検定)投稿日2016.11.4更新日2017.7.21理論上の確率分布から得られる期待度数 ( Expected frequency ) を求めることが前提となります.その期待度数を利用して観測値の度数 ( Observed frequency ) が適合するかどうか(…

二項分布からの最尤推定

投稿日2017.6.13更新日2017.7.20 二項分布から最尤推定に挑戦 二項分布のグラフ - 統計学備忘録 since2016 例)試行回数3回・確率1/2の二項分布 𝐵𝑖(3 , 0.5) 表が出る確率をp=1/2とすると、次のような分布が考えられます.確率pのときに表が〇回出る確率(同…

Rstudioでの data.frame 表示の不具合

データフレームを表示したときに、ラベルが見えなくなりました RStudio-1.0.143 から RStudio-1.0.153 へバージョンアップすることで不具合は解消しました.https://www.rstudio.com/products/rstudio/download/preview/ Rstudioのサポートページの掲示板で…

ファイ係数

連関係数 クロス集計表における2つの変数間の関連性の程度を表す指標として,連関係数が提案されています. 連関係数<|0.2| 連関はほとんどない|0.2|≦連関係数<|0.4| やや連関がある|0.4|≦連関係数<|0.7| 強い連関がある|0.7|≦連関係数 かな…

イェーツの補正

イェーツの補正 / イェーツの連続修正(Yate's continuity correction) 離散型分布を連続型分布に近似させて統計的検定を行う際に使用する修正です.2×2分割表のデータに対して行われ、より正確な検定が可能になります.「連続的なカイ二乗分布(自由度1のχ…

独立性の検定

投稿日2016.11.4更新日2017.7.14 この分割表において独立とは Ai ∩ Bj の各確率に対して全ての i j(どの i j でも)に対して P( Ai ∩ Bj )=P( Ai ) * P( Bj ) となるAとBが完全に独立な場合の例 A1 : A2 : A3 がどのBでも 1 : 3 : 6 B1 : B2 : B3 がどのA…

並べ替え検定

正確なp値 ( exact p value )特定の確率分布をもとに推定を行うのではなく、p値の計算において母集団の未知のパラメータやサンプリング誤差が入らないため計算上も正しいp値が得られる.並べ替え検定 ( permutation test )例)x群とy群を比較するx<-c ( 5…

符号検定

符号検定 ( sign test ) 一標本問題 n個のうちx個以下の符号が他の符号と異なる確率を求める検定です.正と負が同じ確率で出るときに、正(または負)の値の数がxとなる確率が B ( n , 0.5 ) の2項分布に従うことを利用します.例y<- c (1.5, 1.1, 1.2, -1.…

比率の信頼区間

投稿日2017.6.19更新日2017.7.9 比率の信頼区間例)有権者から2,400人を無作為抽出した結果、1,250人は支持していたことがわかった.有権者の支持率の95%信頼区間を求めよ.出典 日本統計学会 (編集); 日本統計学会公式認定 統計検定2級対応 統計学基礎, …

なぜt検定を繰り返してはいけないのか?

なぜt検定を繰り返してはいけないのか投稿日2017.2.24修正日2017.7.9 sample は「石村卓夫;分散分析のはなし,東京図書,1992,p137」 A<-c(12,14,16)B<-c(13,15,17)C<-c(15,17,19)D<-c(17,19,21) x<-c(rep("A",3),rep("B",3),rep("C",3),rep("D",3))y<-c(A,B,…

二項分布

投稿日2017.6.1更新日2017.7.6 二項分布 B ( n , p ) X=確率変数 nは、ベルヌイ試行の実行回数 Xは、ベルヌイ試行の成功回数 pは、成功確率 例)確率0.5で5回実施した場合の、それぞれの成功回数に対する確率 choose=二項係数 dbinom=二項確率choose ( 5…

F分布の確率密度関数

F分布の確率密度関数df(x, df1, df2, ncp, log = FALSE) 自由度1、自由度5のF分布の確率密度関数のグラフを描いてみます # par()を使用する# グラフィックパラメータ:グラフを構成する点や線の形状や色、大きさや余白などを決定する# par()は余白などの土台…

t検定のまとめ

投稿日2017.1.27 最終更新日2017.6.28 Rでは次の式を使用しますt.test ( x, alternative=c("two.sided","less","greater") ) mu=0, paired=FALSE, var.equal=FALSE, conf.level=0.95) 引数alternative:両側検定("two.sided")か片側検定("less","greater"…

ベクトルとデータフレームの関係

更新日2017.6.26ここのページは、ちょこちょこ更新する予定物忘れが激しいので何でも書いときますブログ名通りの使い方ですベクトルをフレームにしたら列になる!data.frame (列名1 = ベクトル1, 列名2 = ベクトル2, ... )まずはここを忘れないように… A<-c(…

Rstudioでのデータ保存、読み込み

よく使うので書き留めておきます 事前準備Rstudio 四つの画面から編成されていますので、勝手に名前をつけます.①書き込み画面、②Consol画面、③環境画面、④ファイル画面 まずRスクリプト上に書いているベクトルを保存してみます①にデータを書き込みますt3.4<…

並行箱ひげ図

エクセルでデータのみコピーt<-read.table("clipboard")x<-t$V1y<-t$V2boxplot(x,y,xaxt="n",main="タイトル") #x軸のラベルを消去name<-c("ラベル1","ラベル2") #x軸のラベルを定義axis(side=1,at=c(1,2),labels=name) #x軸のラベルを更新 このような小技が…

t分布からの信頼区間

自由度5のt分布より確率点:qt(0.025,5,lower.tail = F)累積分布:pt(2.570582,5,lower.tail = F)確率密度:dt(2.570582,5) 拡張期血圧を6回測定して95%信頼区間を求めるDBP<-c(86,92,88,94,89,88)vDBP<-var(DBP) t値 { mean(DBP)-X } ÷ {sqrt(vDBP/6)} 95%…

シミュレーションの準備

まずはforの使い方(なかなか慣れません) 変数xに1から5を足す0+1+2+3+4+5=15をfor関数で実行x<-0for(i in 1:5) x<-x+i ベクトルに繰り返し数を付けるx<-c()i<-c(1,2,3,4,5)x<-c(x,i)for関数で実行x<-c()for(i in 1:5) x<-c(x,i) 関数として定義するとmyfu…

正規分布から解く

前回の掲載がやや分かりにくい内容でしたので、r-de-r様からの助言を基に修正しました(2017.6.11)平均体重がN(70,25^2)に従う集団があります.その集団から10人ランダムに選択します.その合計が800㎏を超えてしまう確率は? 求める確率はP(合計≧800)= P{Z=(…

データの要約

irisx<-iris$Sepal.Length[1:50]y<-iris$Sepal.Length[51:100]xy<-data.frame(x,y)summary(xy)boxplot(x,y)t<-c(1,2)ave<-c(mean(x),mean(y)) #平均値points(t,ave,pch=16) #平均値(・)を挿入 x y Min. :4.300 Min. :4.900 1st Qu.:4.800 1st Qu.:5.600 Me…

x軸操作いろいろ

覚えきれないことは素直にかくことにしましたサンプル標準正規分布N(0,1)x<-seq(-3,3,0.01)y<-dnorm(x)plot(x,y,type="l",ylab ="確率密度",xlab="") 目盛り操作でN(3,1)に変更する方法plot(x,y,type="l",ylab ="確率密度",xlab="",xaxt="n") #xaxt="n"で軸…

確率変数の同時分布

確率変数X、Yの同時分布 出典)柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010,p60 二次元確率変数の期待値X=1の周辺分布P(X=1,Y) = 0.1 + 0.2 = 0.3確率変数Xの期待値E(X) = 1*(0.1+0.2) + 2*(0.4+0.3) = 1.7確率変数Yの期待値E(Y…

独立

「二つの事象A,Bが独立である」ということ一方の事象が起こるかどうかが、他方の事象の起こる確率に影響しないこと条件付き確率P(A│B)=P(A∩B)/(P(B)) より事象Bが事象Aの起こる確率に影響しないのでP(A│B)=P(A) ということを意味していますつまり P(A∩B)=P(A…

共分散と相関係数

(修正日2017.6.5)何度やっても忘れてしまう(悩)ので、大切なことは、真面目にコツコツ書くことにしましたたとえばエクセルからコピーしたデータをベクトルにする方法などx<-read.table("clipboard")x<-x$V1これすらも覚えられないので重症(´;ω;`)共分散デ…

X連鎖劣性遺伝のキャリアの診断

臨床の課題を統計学で読み解く練習前提)高校の生物学の復習.デュシェンヌ型筋ジストロフィー(Duchenne muscular dystrophy: DMD)はX連鎖劣性遺伝に該当します.母親が保有者、父親が正常である場合に、生まれてくる男児がDMDである確率は0.5となります.ち…

ポアソン分布の最尤推定

忘れないうちに書いときますポアソン分布とは…二項分布において、n(→∞)が大きく、p(→0)が小さい場合の確率分布. 平均1,3,5,7のポアソン分布y<-0:15prob<-dpois(y,lambda = 1) #ラムダ=1plot(y,prob,type = "b",lty=2,xlim=c(0,15),ylim=c(0,0.4),xlab = …

異なる周期の三角関数の合成いろいろ

フーリエ変換の予備知識として色々な三角関数の合成を可視化しておきますあくまで準備段階ですので…弧度法を度数法に変換しますx<-pi/180y<-0:360x<-x*y sin(x)、sin(2x) の合成plot(sin(x)+sin(2*x),type="l",ylim = c(-2,2),ann = F)lines(sin(x),type="l"…

cos(nx)とsin(nx)の合成

グラフでイメージするa*cos(x)とb*sin(x)の合成ベクトル(振幅)は√(a^2+b^2)sun(nx)、cos(nx)を合成すると位相が変化するが、周期は変化なし 振幅 √(2^2+2^2)=2.828427 #sqrt(8)1周期plot(2*cos(y*x)+2*sin(y*x),type="l",ylim = c(-3,3),ann = F)lines(2*c…

三角関数のグラフ

角度の概念を復習します 弧度法ラジアン(弧度)とは、「弧の長さ=半径」となる角度を「1」とする角度の単位のこと。単位記号は「rad」R言語ではこの「弧度法」を使用しています. 度数法円周を360等分したものを1とする角度の単位のこと。単位記号は「°」 つ…

RでBland Altman Plotを描く

データ国際医療福祉大学 下井研究室信頼性の検討方法:相対信頼性と絶対信頼性例1(http://shimoi.iuhw.ac.jp/reliability_3_systematicbias.html)を使用させていただきましたx125.1 23.3 33.2 20.9 22.8 34.9 26.5 26.1 29.4 29.8 29.2 26.5 14.3 14.6 17.…

1標本t検定のスクリプト

只今練習中~2標本xとyについてsumxy <- function(x,y,sxy) { xs<-summary(x) ys<-summary(y) xt<-t.test(x) yt<-t.test(y) print(xs);print(ys);print(xt);print(yt)}sumxy (x,y,sxy)より x<-c( );y<-c( )にデータを入れると、各標本の要約・t値・p値・95%C…

繰り返し for

いつも忘れるので書いておこう#x=3に2を3回足す式#3+2+2+2=9 myf<-function(){ x<-3 for(i in 1:3){ x<-x+2} return(x)}myf() #Xに2を3回加える式#x+2+2+2#xに5を代入すると5+2+2+2=11 myf<-function(x){ a<-x for(i in 1:3){ a<-a+2} return(a)}myf(5) #Xに…

相関係数の区間推定

サンプル:日本統計学会 (編集); 日本統計学会公式認定 統計検定2級対応 統計学基礎, 東京図書, 2012,p125 身長のデータ大学生x<-c(172,167,184,175,176,175,170,180,170,179,167,175,174,162,165,163,170,169,165,175)父親y<-c(165,165,178,176,150,171,17…

等分散性の検定

準備としてカイ二乗分布とF分布が必要になります カイ二乗分布確率変数Z1,Z2,…Znが互いに独立で標準正規分布にしたがうときW=Z1^2+Z2^2+…+Zn^2の従う分布を自由度nのカイ二乗分布とよび、χ^2(n)で表します.期待値E[W]=n, 分散V[E]=2n カイ二乗分布の性質① 二つの確…

一元配置分散分析 (対応なし) F値の算出方法

更新日2017.3.13更新日2017.6.22 一元配置分散分析とは、一つの因子に基づく“3つ以上の母平均の差の検定”です 因子Factor:実験結果に影響を与えると考えられる要因水準Lebel:因子をいくつかに分ける基準 (水準=群) 例題)治療因子に基づく、患者20名の…

正規分布のグラフ

グラフ初歩の初歩 手探り状態でやってます 赤文字のみRで実行 正規分布#1正規分布の乱数を100個生成x<-rnorm(100)#ヒストグラム,freq=Fで確率密度,ylimでy軸範囲hist(x,freq=F,ylim=c(0,0.6))#枠の処理、lty線種は実践、btyで左と下のみbox(lty=1,bty = "l")…

二項分布のグラフ

こんな説明でどうでしょうか? コイン投げ(いかさまなし)表の出る確率確率1/23回コイン投げを実施してみた場合・・・表が出る回数(x)は0回、1回、2回、3回xの確率は𝐵𝑖(3 , 0.5) に従う.その確率分布を二項分布といいます. 村上 正康,安田 正実; 統計…

グラフに色をつける

library(RColorBrewer)#RColorBrewerパッケージのサンプルdisplay.brewer.all() ヒストグラムを塗ってみるどの色セットを使用するかを指定するcols <- brewer.pal(8,"Pastel1") # brewer.pal(何色、パレット名) y<-c(1,2,3,4,5,6,7)p<-c(2,3,4,5,4,3,2)q<-…

棒グラフの中央に散布図をプロットする方法 (pos.x)

y<-c(1,2,3,4,5,6,7)p<-c(2,3,4,5,4,3,2)q<-c(2,3,4,5,4,3,2)par(mar = c(5, 4, 1, 4)) #余白 底辺、左、上、右の順 pos.x <- barplot(q,ylim=c(0,6))points(pos.x, p)

Rことはじめ

1、Rのインストール方法と使い方 2、以下のデータの母平均の95%信頼区間を求めてみます 73.0 75.8 74.0 72.3 76.9 82.5 69.1 80.9 74.9 71.6 73.5 84.6 81.3 82.578.1 71.3 76.3 74.9 74.8 72.5 78.5 71.5 70.0 81.1 76.6 77.9 75.5 70.881.9 75.1 85.4 8…