統計分析ソフトRのミニソース集

かゆいところに手が届くようなRのミニソースを作成

ggplotでF分布と棄却域を描く自作関数

RでF分布と棄却域を描画したいとき,片所強先生が公開されている関数がたいへん便利である。

引数に"norm", "t", "chisq", "F"のいずれかを指定するだけで簡単に分布曲線が描かれます。描こうとする分布によっては、x軸の値が適当でないこともあるので、その場合は引数myself = TRUEとすることによって分布を描くための関数自体が返されます。*1

分布のプロット(正規分布、t分布、F分布、χ2分布)

片所先生の関数を使えば,例えば,df(2,100)のF分布と有意水準を5%とした場合の棄却域は,次のように書けば良い。

#※予め関数の定義を読ませておく
d <- my.dist("F", myself=TRUE)
#自由度を聞かれるので2と100を入力する
plot(d, 0, 10) #第2,3引数は,x軸の描画範囲
my.paint(d, qf(1-0.05, 2, 100), 10, col="grey") #第2,3引数は,塗りつぶしに対応するx軸の範囲

f:id:mini-tool:20151220020637p:plain

これで十分であるが,場合によっては,もう少し丁寧に描画したいことがある。
今回は,ggplotを使った関数を書いた。
以下に関数の定義コードを示す。

F.dist <- function(df1,df2,alpha=0.05,fmax=20){
my.F <- data.frame(X=x<-seq(0,fmax,len=10000), Y=df(x,df1,df2)) #F分布の確率密度関数データフレーム
p <- ggplot()
p <- p + geom_path(data=my.F, aes(x=X, y=Y), size=1.0)
p <- p + labs(x = "F value", y = "Probability density")

#以下,棄却域
upper <- qf(1-alpha, df1, df2) #臨界値
ymin <- par("usr")[3]
ymax <- par("usr")[4]
pos1 <- ((ymin+ymax)/100)*20*0.70 #矢印↓のy座標上限
pos2 <- ((ymin+ymax)/100)*20 #テキストのy座標

p <- p + geom_ribbon(data=subset(my.F, X>=upper), aes(x=X, ymin=0, ymax=Y), fill="gray")
p <- p + annotate("text",size=4.5,label=format(upper,digits=2,nsmall=2,scientific=F),x=upper,y=pos2, fill="white")
p <- p + annotate("segment",size=1, x=upper, xend=upper, y=pos1, yend=0.0,arrow=arrow())
p <- p + scale_x_continuous(expand = c(0, 0))
p <- p + scale_y_continuous(expand = c(0, 0))
p <- p + annotate("text", x = Inf, y = Inf, label = paste("F(",df1,",",df2,")\nα=",alpha,sep=""),
            hjust=1.1, vjust=1.1, col="black", cex=6,
            fontface = "bold", alpha = 0.8)

#以下,書式の調整
p <- p + theme(axis.title.x = element_text(vjust = 0), #vjustでラベル位置調整
		axis.title.y = element_text(vjust = 2, angle = 90), #vjustでラベル位置調整
		panel.background = element_rect(fill = "white", colour = NA),
		plot.background = element_rect(colour = "white"),
		axis.line = element_line(colour = "#000000", size = 1, linetype = "solid", lineend = "round"),
		axis.text = element_text(size = rel(1)),
		axis.text.x = element_text(vjust = 1,colour="#000000"),
		axis.text.y = element_text(hjust = 1,colour="#000000"),
		axis.ticks = element_line(size = 1,colour="#000000"),
		axis.ticks.length = unit(0.2, "cm"),
		axis.ticks.margin = unit(0.2, "cm"),
		panel.grid = element_blank(),
		text = element_text(face = "plain", colour = "black", size = 12, hjust = 0.5, vjust = 0.5, angle = 0, lineheight = 0.9))

plot(p)
return(paste("F分布(",df1,",",df2,")における,α=",alpha,"の棄却域はF>=",format(upper,digits=2,nsmall=2,scientific=F),"です。",sep=""))
}

この関数を使って,df(2,100)のF分布と有意水準を5%とした場合の棄却域を描くには,次のようにすれば良い。

#※予め関数の定義を読ませておく
library(ggplot2)
library(grid)
F.dist(2,100,fmax=10) #fmaxはx軸の描画上限値(省略時20)

f:id:mini-tool:20151220021505p:plain

有意水準を変えたい場合は,引数にalphaを追加すれば良い。

F.dist(2, 100, fmax=10, alpha=0.1)

f:id:mini-tool:20151220022335p:plain

一行で描画できるため,次のようなことを理解する際に良い。

  • 自由度を変えた場合に,F分布の形状がどう変化するか?
  • 有意水準を変えた場合に,棄却域はどう変化するか?

次回は,F分布に形状が近い,χ2乗分布の関数を作成する予定である。

*1:2015年12月19日閲覧