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

統計学備忘録 since2016

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

二元配置分散分析 分散分析表の理解

投稿日2017.3.8更新日2017.6.26 繰り返し数がすべて等しい場合 (=5)水準A1 × 5回 = (A11 ,A12 ,A13 ,A14 ,A15 ) 水準A2 × 5回 = (A21 ,A22 ,A23 ,A24 ,A25 )つまり水準B (1,2,3) はA1、A2をそれぞれ5回繰り返したことになります (仮想データ使用) B1<-c(…

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

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

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

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

比率の信頼区間

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

並行箱ひげ図

エクセルでデータのみコピー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軸のラベルを更新 このような小技が…

二項分布からの最尤推定

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

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 P(X=1,Y) = 0.1 + 0.2 = 0.3 E(X) = 1*(0.1+0.3) + 2*(0.4+0.3) = 1.7E(Y) = 1*(0.1+0.4) + 2*(0.2+0.3) = 1.5E(XY)=1*1*0.1+1*2*0.2+2*1*0.4+2*2…

独立

「二つの事象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これすらも覚えられないので重症(´;ω;`)共分散デ…

二項分布

臨床の課題を統計学で読み解く練習柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010,p50 演習問題3.2の答えdbinom=確率密度( ) =(成功回数、サイズ、確率)1-dbinom(0,50,1/80)-dbinom(1,50,1/80)-dbinom(2,50,1/80)-dbinom(3,50,1…

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…

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

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,C,D)stripchart(y~x,vertical=T) 視覚的にはAD間、BD間には差があるよう…

等分散性の検定

準備としてカイ二乗分布と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 赤文字の部分だけRにペーストして実行 一元配置分散分析とは、一つの因子に基づく“3つ以上の母平均の差の検定”です因子Factor:実験結果に影響を与えると考えられる要因水準Lebel:因子をいくつかに分ける基準 (水準=群) …

正規分布のグラフ

グラフ初歩の初歩 手探り状態でやってます 赤文字のみ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")…

回帰係数の区間推定

赤文字をRで実行 サンプル)統計学入門(基礎統計学),東京大学出版会,1991,p258 次のデータをコピーしてください東京 福岡1019.4 1018.41005.7 1007.61002 1006.21006.7 1009.91005.1 1010.81010.1 1013.21016.7 1016.21011 1009.1999.5 1003.11006.9 1012.51…

単回帰分析

赤文字のみRで実行 線形回帰分析説明変数と目的変数を直線関係で傾向を示す。説明変数と目的変数との関係を直線でモデル化する回帰分析。 非線形回帰分析非線形関係でモデル化する回帰分析。説明変数と目的変数を非直線的関係で傾向を示す。 ここでは単回帰…

t検定のまとめ

修正日2017.6.15修正日2017.6.19 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")かm…

Rで簡単(適合度検定)

例1)次のデータはメンデルの法則 1:2:1 の確率に適合しているか?AA Aa aa n52 102 46 200帰無仮説:標本はメンデルの法則に適合している 観測度数<-c(52,102,46)確率<-c(1/4,1/2,1/4) chisq.test(観測度数, p=確率) Chi-squared test for given probabilit…

Rで簡単(独立性の検定)

赤文字のみをRで実行 次のクロス表を検定してみます問)男性と女性で、yesの解答者数とnoの解答者数の比率に違いはあるか? Yes N0男性 n1 n2女性 n3 n4 帰無仮説:男性と女性では解答者数の比率に違いがない例)仮想データ n1<-15n2<-11n3<-4n4<-16 観測度…

二項分布のグラフ

こんな説明でどうでしょうか? コイン投げ(いかさまなし)表の出る確率確率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…

級内相関係数 ICC

修正日 2017.3.2修正日 2017.6.21 ICC (Intraclass correlation coefficients)ICCは信頼性の指標の一つとして利用されています.ICCには3つの分類があり、6Caseに分類されています. 注)ICCの記載方法は統一されていない ICCは分散分析の結果を活用していま…

Rで簡単(t検定の検定力分析)

パッケージpwrを使用しますlibrary(pwr) pwr.t.test(n = , d = , sig.level = , power = , type = c("two.sample", "one.sample", "paired"))n=標本数d=効果量()sig.level=有意水準power=検定力type=t検定の形("two.sample=対応なし", "one.sample=1標本"…

Rで簡単(t検定:対応あり)

例えばn=14の前後比較・・・ 下のデータをコピーして 前 後 90 98 90 95 92 77 91 84 83 75 93 88 89 89 81 88 80 86 77 86 78 79 77 90 76 88 79 99 Rに貼り付け x<-read.table("clipboard",header=T) 対応のあるデータなので「差」を検定することになる …

時系列データの可視化01

下のデータをコピーして時系列グラフを作ってみます 心拍変動を三次元加速度を同時測定した結果です 被験者は私です RRI TEM X Y Z HF LFHF LF activity HR793 21.5 0.16 -1.09 -0.12 75.531 4.247 80.941 0.11 76800 21.5 0.16 -1.09 -0.12 75.436 4.078 80…

Rで簡単(二元配置分散分析:対応なし)

サンプルをコピペで2式を実行 方法は Rで簡単(一元配置分散分析:対応なし)を参照 検定は帰無仮説が大切です! 商品のの主効果:帰無仮説=カップ麺でもインスタント麺でも評価の母平均は等しい スープの主効果:帰無仮説=スープが違っても評価の母平均は…

Rで簡単(一元配置分散分析:対応あり)

方法は「Rで簡単(一元配置分散分析:対応なし)」と同じです 最下部の式が少し違うだけです サンプル 人 種類 得点野原 A 12野原 B 8野原 C 9風間 A 8風間 B 5風間 C 8桜田 A 12桜田 B 10桜田 C 10佐藤 A 6佐藤 B 2佐藤 C 3東溝 A 8東溝 B 5東溝 C 6 x<-rea…

Rで簡単(一元配置分散分析:対応なし)

サンプルをコピーして(エクセルからでもOK)、下の2式を実行するのみ サンプル 治療方法 点数A 15A 9A 18A 14A 18B 13B 8B 8B 12B 7C 10C 6C 11C 7C 12D 10D 7D 3D 5D 7 x<-read.table("clipboard",header=T) summary(aov(x$点数~x$治療方法)) 実際の方法を…

Rで簡単(重回帰分析)

Rのデータセットtrees を使用して説明します Girth=y , Height=x1 , Volume=x2 次のようにエクセルに準備します x<-read.table("clipboard",header=T) をコピーしてR Consoleにペーストします 分析したい一覧表をコピー(上記の表など) x<-read.table("clip…

Rで簡単(t検定:対応なし)

RをPCにインストール https://cran.ism.ac.jp/ 次のような操作を行うだけ簡単にt検定ができますので、是非挑戦してみてください。 難しいことはさておき、とにかくやってみることが大切です R Consoleに直接入力する方法を試しましょう! 例)結果1と結果2…

確率とパーセント点

正規分布の累積分布(確率) pnorm(1.959964,lower=F) #Z≧1.959964となる上側確率=0.025 pnorm(-1.959964) #Z≦-1.959964となる上側確率=0.025 正規分布の確率点(パーセント点) qnorm(0.025,lower=F) #上側確率が0.025となるZ値=1.959964 qnorm(0.025) #…

平均と分散

標本例x<-c(1,2,3,4) 標本平均(次の2式は同じ解になります)sum(x)/length(x)mean(x) 二乗平均sum(x^2)/length(x) 標本分散(3式とも同じ解になります)var(x)*(length(x)-1)/length(x)sum( (x-mean(x))^2)/length(x)sum(x^2)/length(x)-(mean(x))^2 標準誤差 …

適合度の検定

k.ピアソンの適合度基準χ2乗値=Σ(観測度数-期待度数)^2)/期待度数)帰無仮説:各階級の発生確率は母集団の確率に適合する 例)次のデータはメンデルの法則 1:2:1 の確率に適合しているか? AA Aa aa n観測度数 52 102 46 200期待度数 50 100 50 200 観測度数<…

独立性の検定

赤文字のみRで実行 =======================================次のクロス表を検定してみます問)男性と女性で、yesの解答者数とnoの解答者数の比率に違いはあるか? Yes N0男性 n1 n2女性 n3 n4============================================================…