須通り_統計_Rで分割表を100%積み上げ棒グラフ

ホーム | 統計 Top | 何が何でもRで分割表を100%積み上げ棒グラフで表示する

Rのbarplot()関数で100%積み上げ棒グラフを表示できないので

分割表データの入力から検定、グラフ作成までを一貫してRでやろうとすると少し面倒である。基本的には以下のようなマトリックスの形でデータを作成・保持するのだが、

> example <- matrix( c(67,128, 155,180), nrow=2, byrow=T) ##テストデータ作成
> example ##表示してみう
      [,1] [,2]
[1,]    67  128
[2.]   155  180

これを100%積み上げ棒グラフで表示したいとき、単純に「棒グラフだったらbarplot()関数でいいんじゃろ」と考えて

> barplot(example, beside=F) ##棒グラフを表示
## なおbeside=F（論理値FALSE）のとき、複数系列が積み上げグラフとして処理される。
## TRUEならば横に並べられ、下の例だと「黒→灰色→黒→灰色」の順に4本のバーが描かれる。

高さが揃わない棒グラフの例
などとやっても各棒の高さがバラバラで、そんなん考慮しとらんよ……ということになる。これだと100%ではない、ただの積み上げ棒グラフである。各列の度数が違うのだから当たり前だ。

【自作関数】分割表形式の行列データを割合に変換する【しょぼい】

そこで、分割表を構成する行列要素を、生の数字から割合に変換する関数を作ってみた。ひょっとして既に同等の機能が用意されているのかもしれないが、ことRに関しては調べるより自作したほうが早いので気にしない。

↓↓以下の行をRコンソールにコピペしてリターンを押すことで、関数がインストールされる（"x2ratio"というオブジェクト名称を既に使っている場合は、適宜変更）

##function "x2ratio" (C) 2012 SUDO Masaaki
##分割表形式の行列の各要素を、横(デフォルト)もしくは縦方向の合計が1になるような割合値に変換する
##（たぶん）行列要素が全て実数の場合のみ有効
x2ratio <- function(x, byrow=FALSE) {
  if ( !is.matrix(x) ) return(NA)  #引数xが行列以外ならば実行しない
  if ( !is.logical(byrow) ) return(NA) #引数byrowが論理値以外ならば実行しない
      if (byrow == FALSE)  { return(     t(  t(x) / apply(x , 2, sum)  )  )
      }
      else return(    x / apply(x , 1, sum)  )
}
##ここまで

##使用例
> example ##サンプルデータ
      [,1] [,2]
[1,]    67  128
[2.]   155  180

> x2ratio(example) ##早速実行してみう
          [,1]      [,2]
[1,] 0.3018018 0.4155844
[2.] 0.6981982 0.5844156

> barplot ( x2ratio(example), beside=F ) ##実行結果をネストして棒グラフを表示

高さを揃えた積み上げ棒グラフの例
これで100%積み上げ棒グラフを表示でき、大勝利である。

なお本関数のオプションとして、縦に（列について：デフォルト）集計するか、横に（行について）集計するかを論理値で選べる。

> x2ratio(example, byrow=F) ##縦に集計（省略時はこちら）
          [,1]      [,2]
[1,] 0.3018018 0.4155844
[2.] 0.6981982 0.5844156
  
> x2ratio(example, byrow=T) ##横に集計
          [,1]      [,2]
[1,] 0.3435897 0.6564103
[2.] 0.4626866 0.5373134

barplot()が返す座標について

さて、上記の方法で何とかグラフを表示できたとする。これにエラーバーを追加したい場合、今のところbarplot()関数自体にはそのようなオプションが無いため、低水準作画関数を用いてシコシコ描画することになる。

##2×4分割表の例
> example24 <- matrix( c(67,128,3,22, 155,180,9,49), nrow=2, byrow=T) 
> example24
     [,1] [,2] [,3] [,4]
[1,]   67  128    3   22
[2,]  155  180    9   49

> barplot(x2ratio(example24), beside=F )

ここまでの処理により、積み上げ棒グラフが4本横に並んだものが描かれるはずだ。
4×2の積み上げ棒グラフ

Y座標について

さて、このグラフにおける座標値を考えよう。y座標は簡単で、グラフの目盛りの値をそのまま読めば良い。

> barplot(x2ratio(example24), beside=F ) ##グラフを描画
> abline(h=0.5, col=rgb(1,0,0,0.8)) ##y=0.5の位置に水平線を赤色で追加する

積み上げ棒グラフに水平線を追加
このabline()は、Rの描画デバイス内に直線を追加描画するための関数であり、上記処理によってy=0.5の位置に水平な直線が追加された。

X座標について

これに対して、x座標の読み取りは少しやっかいである。

barplot(x2ratio(example24), beside=F )
abline(v=0, col=rgb(1,0,1,0.8)) ##h=0の位置に垂線を紫色で追加する
abline(v=1, col=rgb(0,0,1,0.8)) ##h=1の位置に垂線を青色で追加する

積み上げ棒グラフに垂線を追加
つまりx=0となる原点は、1本目のプロットの左端ではなく、y軸線上である。そしてy軸と1本目プロットの左端との間にマージンがあり、棒が増えるごとにマージンの幅が積み増しされてゆく。

これはどういうことなの......とbarplot()関数のヘルプファイルを読むと、"space"という設定項目があり、"the amount of space (as a fraction of the average bar width) left before each bar. May be given as a single number or one number per bar."と記されている。さらにこの値はbeside=Fの場合、デフォルトでは0.2になるとも書かれている。
つまり、1本目のプロットより左にマージン分の0.2が存在し、そこから系列ごとに棒の幅である1と、マージンの幅である0.2が追加されていくってことか。

barplot(x2ratio(example24), beside=F )
abline(v=0.7, col=rgb(1,0,1,0.8), lwd=2) ##h=(0.2+0.5)の位置に垂線を紫色で追加する
abline(v=1.2, col=rgb(0,0,1,0.8), lwd=2) ##h=(0.2+1.0)の位置に垂線を青色で追加する 
abline(v=1.4, col=rgb(1,0,0,0.8), lwd=2) ##h=(0.2+1.0+0.2)の位置に垂線を赤色で追加する 
text(0.7,0.1,"0.7") ##(x,y)=(0.7,0.1)にラベルを追加。位置は標準で中央揃えとなる
text(1.2,0.2,"1.2")
text(1.4,0.3,"1.4")

積み上げ棒グラフの座標
これで、x,y座標を自在に指定して図形ないしテキストを追加表示できる。科学の勝利なのであった。

barplot()関数が返すx座標の再利用

実際には、さらに簡単なx座標の取得方法がある。それは、barplot()関数が返す値の再利用である。

hogehoge <- barplot(x2ratio(example24), beside=F ) ##barplot()関数の結果を、適当な名前で新規オブジェクトに代入すると......
hogehoge
[1] 0.7 1.9 3.1 4.3

このベクトルhogehogeが表すものは、barplot()関数によって描画された棒グラフにおける、4本のプロットの各中心に他ならない。

100%積み上げ棒グラフに信頼区間のエラーバーを追加表示する

上記を踏まえて、barplot()関数で描いた積み上げ棒グラフに、二項信頼区間のエラーバーを追加してみよう。

> example24
     [,1] [,2] [,3] [,4]
[1,]   67  128    3   22
[2,]  155  180    9   49

> binom.test( c(example24[1,1],example24[2,1]) )

        Exact binomial test

data:  c(example24[1, 1], example24[2, 1]) 
number of successes = 67, number of trials = 222, p-value = 3.299e-09
alternative hypothesis: true probability of success is not equal to 0.5 
95 percent confidence interval:
 0.2421815 0.3668047 
sample estimates:
probability of success 
             0.3018018

たとえば左端の系列である(67, 155)について割合は67/(67+155) = 0.3018、二項分布に基づく95%信頼区間は0.242-0.367と計算される。
エラーバーとして、今回は両端の「ひげ」を付けない線を描く。Rで線分を描くための低水準作画関数である、segments()を用いる。←segment()ではないので注意

barplot(x2ratio(example24), beside=F )
segments(0.7, 0.242, 0.7, 0.367, lwd=2) ##信頼区間の線分を追加する

##なお、座標は左から順に（始点のx）（始点のy）（終点のx）（終点のy）
##lwd=2というのは、線幅を標準の2倍にするという意味のオプション

積み上げ棒グラフの座標