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

統計学備忘録(R言語のメモ)

since2016 ときどきTEXのメモ

正規分布のグラフ

グラフ初歩の初歩

手探り状態でやってます

赤文字のみ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")
#確率密度のグラフ,add=Tで重ねる、x軸範囲-3~3
curve(dnorm(x),add=T,-3,3)


これらをオリジナルの関数として記憶させます
mynorm<-function(n,g)
{
x<-rnorm(n)
hist(x,freq=F,ylim=c(0,0.6),xlab=g)
box(lty=1,bty = "l")
curve(dnorm(x),add=T,-3,3)
}


mynorm(100,"n=100")

 

f:id:yoshida931:20170212153358p:plain

 

#四つのグラフを2行2列のまとめてみます
par(mfrow = c(2,2))
mynorm(30,"n=30")
mynorm(100,"n=100")
mynorm(300,"n=300")
mynorm(1000,"n=1000")

f:id:yoshida931:20170212153327p:plain

グラフに色をつける

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<-c(2,3,4,5,4,3,2)
par(mfrow = c(3,3),mar = c(5, 4, 1, 4)) #余白 底辺、左、上、右の順 
pos.x <- barplot(q,ylim=c(0,6),col=cols[1])
pos.x <- barplot(q,ylim=c(0,6),col=cols[2])
pos.x <- barplot(q,ylim=c(0,6),col=cols[3])
pos.x <- barplot(q,ylim=c(0,6),col=cols[4])
pos.x <- barplot(q,ylim=c(0,6),col=cols[5])
pos.x <- barplot(q,ylim=c(0,6),col=cols[6])
pos.x <- barplot(q,ylim=c(0,6),col=cols[7])
pos.x <- barplot(q,ylim=c(0,6),col=cols[8])
pos.x <- barplot(q,ylim=c(0,6),col=cols[9])

 

f:id:yoshida931:20170114232836p:plain

 

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標本", "paired=対応あり")

n, d, power, and sig.level のいずれかを空白(NULL)にします
strict=T   #厳密な検定が可能になる
 (効果量(effect size)のはなし - 六本木で働くデータサイエンティストのブログ



必要な標本数を求めてみます
pwr.t.test(n =NULL, d =0.2 , sig.level =0.05 , power = 0.8, type ="two.sample")

Two-sample t test power calculation

n = 393.4057
d = 0.2
sig.level = 0.05
power = 0.8
alternative = two.sided

 

効果量0.2、有意水準0.05、検定力0.8の場合の標本の大きさは393必要になります
もちろんフリーソフトG*powerを使用した結果と同じになります

f:id:yoshida931:20170121182510p:plain

 


効果量をもとめてみます
pwr.t.test(n =150, d =NULL , sig.level =0.05 , power = 0.8, type ="two.sample")

Two-sample t test power calculation

n = 150
d = 0.3245459
sig.level = 0.05
power = 0.8
alternative = two.sided

 

標本数150、有意水準0.05、検定力0.8の場合の効果量は0.32と低くなります
フリーソフトG*powerを使用した結果

f:id:yoshida931:20170121183731p:plain


Rの予備知識
for(i in c(2,4,6)) でiに2,4,6を順に挿入する命令
print(・・・) 引数が数値ベクトルの場合にはベクトルの中身を表示する
print(c(1,2,3))
合わせて
for(i in c(0.26,0.24,0.22)){print(x<-c(i+1))}
round 四捨五入を導入して・・・
for(i in c(0.26,0.24,0.22)){print(round(x<-c(i+1),1))}

上記を利用して
有意水準が.05,.01,.001)のときの対応なしt検定のpowerを一気に算出してみます
pwr.t.test(n =NULL,d =0.2,sig.level=0.05,power=0.8,type ="two.sample")
for(i in c(.05,.01,.001)){print(round(pwr.t.test(n =394,d =0.2,sig.level=i,
,type ="two.sample")$power,2))}

[1] 0.8
[1] 0.59
[1] 0.31


対応のなし(2群の標本数が異なる場合)
pwr.t2n.test(n1 = , n2= , d = , sig.level =, power = )
例)
pwr.t2n.test(n1 =120, n2=NULL , d =0.5 , sig.level =0.05, power =0.8 )

t test power calculation

n1 = 120
n2 = 43.21666
d = 0.5
sig.level = 0.05
power = 0.8
alternative = two.sided

 

n1=120、効果量0.5、有意水準0.05、検定力0.8、両側検定の場合n1=43.2必要