2012/08/26

文系のための「数の可視化」(3)

今回は、少し、長いかもしれないが、実際の説明は長くない。
長く見えるのは、今回はグラフ、データ、コマンドの類が多いため。

今回の話では、色々と細かいテクニックが出てくる。
落ち着いて、前から順に進めて行けば、なんとか理解できるハズ。

さて、棒グラフの話では、「数量」の可視化について述べた。
棒グラフは、変数間の値の大小を差を視覚的に把握することが目的であった。

今回は、もう一歩踏み込んで、「比率」の可視化というのを考えたい。
では、「比率」とは何か?ここでは、これが大きな問題である。

そもそも、比率というのが何故、重要かというと、
数値だけの比較では、全体と個々の状況を比較できないからである。

全体に対して、個々の状況がどうのようになっているか?
それを知ることは、数量の大小だけでは見えなかったことが解るかもしれない。

では、具体的にどのようにすれば良いか、という疑問が出てくる。

要するに、全部あわせて「1」となるような単位を考え、
その単位「1」をどのように配分するか、というのが比率である。

当然、全体が「1」となるような単位であるので、全体(総数)の逆数を取ればよい。
ある変数における対象の値に掛けると、全体を「1」とした単位で、
個々の値を揃えることができる。例えば、次のような状況を考える。



xという変数が 1〜n 個並んでいて、その不特定番目が「i」だったとする。
これはベクトルであり、実際には縦と横の方向は関係無いのだけれど、
後の話に合わせるために、ここでは縦ベクトルとして書いておく。

ここで、全体が1となるような逆数を考えるので次のようになる。



ここまでの基礎を積み上げてきた人は、容易に理解できるハズ。

全体を1とした単位で基準化した単位を p としたときの i 番目の値p(i)と置いている。
分数の分母の部分は、変数 x における総和を表しているに過ぎない。

ところで、似たような話を小学校で聞いたことがあった。
その時には、この単位をさらに「100倍」し、「%」という記号を付けたものであった。
そうそう、確か、「百分率」とか言う名前だった。これを思い出せば良い。

さて、説明だけでは物足りない。実際にやってみよう。
データは前回から引き続き、発電量のデータ。データについては前回の話を参照。

# Windows と Linux の人は次のコマンド 
X <- read.table("clipboard", sep=",") 

# Mac の人は次のコマンド 

X <- read.table(pipe("pbpaste") , sep=",")

利用するデータは、以下の通り。

pnum,max,hydro_pnum,hydro_max,thermal_pnum,thermal_max,atomic_pnum,atomic_max
Hokkaidou,66,7419,53,1234,11,4065,1,2070
Touhoku,228,17206,209,2423,13,11286,2,3274
Tokyo,192,64988,162,8981,25,38696,3,17308
Chubu,197,32828,183,5219,11,23969,1,3617
Hokuriku,138,8057,127,1904,6,4400,1,1746
Kansai,165,34877,149,8196,12,16907,3,9768
Chugoku,110,11986,97,2906,12,7801,1,1280
Shikoku,65,6963,58,1141,4,3797,1,2022
Kyushu,194,20330,139,3279,45,11577,2,5258
Okinawa,22,1919,0,0,21,1919,0,0

まずは、簡単なことから。最初に、各変数の合計を出したい。
確か、colSums()関数 を使うのだった。

ステップ3までの操作は、文系のための「データの観察」(4)とほぼ同じ。
colMeans()関数の部分が、colSums()関数になっているだけ。

# ステップ1:データの総数(n)を求める。
n <- nrow(X)
# ステップ2:各変数の列和を求める。転置が必要。
m <- t(colSums(X))
# ステップ3:行列の成分の割り算ができるように、列和を総数個複製する。
M <- matrix(rep(m, each=n), nrow=n)
# ステップ4:比率のデータにするために行列の成分の割り算を行う。
A <- X/M

何度も言っていることであるが、
いわゆる、行列の計算行列の成分の計算混同してはいけない。

さて、上手く行っているかを確認するのは簡単。
先程の説明から、行列Aの列和は「1」になって居るハズ。

> colSums(A)
        pnum          max   hydro_pnum    hydro_max thermal_pnum  thermal_max
           1            1            1            1            1            1
 atomic_pnum   atomic_max
           1            1

ふむ。その通りになっている。上手く行っている。
では、この行列Aを改めて観察してみよう。

> A
                pnum         max hydro_pnum  hydro_max thermal_pnum thermal_max
Hokkaidou 0.04793028 0.035914665 0.04502974 0.03497435      0.06875  0.03267238
Touhoku   0.16557734 0.083292589 0.17757009 0.06867330      0.08125  0.09071108
Tokyo     0.13943355 0.314600650 0.13763806 0.25454185      0.15625  0.31101859
Chubu     0.14306463 0.158917187 0.15548003 0.14791826      0.06875  0.19265052
Hokuriku  0.10021786 0.039003161 0.10790144 0.05396367      0.03750  0.03536494
Kansai    0.11982571 0.168836198 0.12659303 0.23229317      0.07500  0.13588979
Chugoku   0.07988381 0.058023072 0.08241291 0.08236261      0.07500  0.06270043
Shikoku   0.04720407 0.033707212 0.04927782 0.03233852      0.02500  0.03051834
Kyushu    0.14088598 0.098415572 0.11809686 0.09293427      0.28125  0.09304999
Okinawa   0.01597676 0.009289694 0.00000000 0.00000000      0.13125  0.01542394

          atomic_pnum atomic_max
Hokkaidou  0.06666667 0.04466694
Touhoku    0.13333333 0.07064713
Tokyo      0.20000000 0.37347604
Chubu      0.06666667 0.07804846
Hokuriku   0.06666667 0.03767559
Kansai     0.20000000 0.21077617
Chugoku    0.06666667 0.02762014
Shikoku    0.06666667 0.04363118
Kyushu     0.13333333 0.11345834
Okinawa    0.00000000 0.00000000

小数点を全て表示する意味はあまりない。小数点第三位くらいで止めておく。
それには、round()関数というものを使う。ようするに、丸め込み。

> round(A,3)
           pnum   max hydro_pnum hydro_max thermal_pnum thermal_max atomic_pnum
Hokkaidou 0.048 0.036      0.045     0.035        0.069       0.033       0.067
Touhoku   0.166 0.083      0.178     0.069        0.081       0.091       0.133
Tokyo     0.139 0.315      0.138     0.255        0.156       0.311       0.200
Chubu     0.143 0.159      0.155     0.148        0.069       0.193       0.067
Hokuriku  0.100 0.039      0.108     0.054        0.038       0.035       0.067
Kansai    0.120 0.169      0.127     0.232        0.075       0.136       0.200
Chugoku   0.080 0.058      0.082     0.082        0.075       0.063       0.067
Shikoku   0.047 0.034      0.049     0.032        0.025       0.031       0.067
Kyushu    0.141 0.098      0.118     0.093        0.281       0.093       0.133
Okinawa   0.016 0.009      0.000     0.000        0.131       0.015       0.000

          atomic_max
Hokkaidou      0.045
Touhoku        0.071
Tokyo          0.373
Chubu          0.078
Hokuriku       0.038
Kansai         0.211
Chugoku        0.028
Shikoku        0.044
Kyushu         0.113
Okinawa        0.000

これを見て、なるほど、そういうことか!と思う人はほぼ皆無であろう。
このように数値が並んでいても、解釈は困難である。
やはり、この状況に対応した可視化が必要である。

そこで登場するのが、「円グラフ」と呼ばれるものである。
安易に円グラフを使う人も少なくはないが、比率のイメージを持つことが重要。

とりあえず、都道府県別の発電所数を観察してみる。
Rでは、円グラフを描くためのpie()関数があるので、これを使う。
pie(X$pnum)

棒グラフの回にも書いたが、このような図は意味が無い。単なる自己満足。
このようなグラフを見せられても、誰も理解できないのである。

パッと見て、タイトルラベル名の追加が必要である。では、そのように。
pie(X$pnum, labels=rownames(X))
title(main="Number of electric plants in Japan")

とりあえずは、見れる図にはなった。しかし、まだ色々と不備がある。
現在の順番は、3時(北海道)を起点に、反時計周りで沖縄までが並んでいるだけ。

これは、与えられたデータの順番であって、値の大小は関係ない
出来れば、データの並べ替えを行い、より見易いデータにしたい。

Rでは、並び替えの関数が二つ存在している。sort()関数order()関数である。
sort()関数は、昇順でデータを並び替え、その結果の値をベクトルで返す
order()関数は、昇順でデータを並び替え、元の順番の番号をベクトルで返す

とにかく、言葉だけだと解りにくい。まずは、pnumの変数で実験する。

> sort(X$pnum)
 [1]  22  65  66 110 138 165 192 194 197 228
> order(X$pnum)
 [1] 10  8  1  7  5  6  3  9  4  2

この結果から解るように、sort()関数では、値が小さい方から大きい方に並んでいて、
order()関数では、元の行番号が、sort()関数に対応するように並んでいる。

また、次のパラメータを与えることで、降順に並び替えることもできる。

> sort(X$pnum, decreasing=TRUE)
 [1] 228 197 194 192 165 138 110  66  65  22
> order(X$pnum, decreasing=TRUE)
 [1]  2  4  9  3  6  5  7  1  8 10

デフォルトの状態では、 「decreasing =FALSE」になっているが、
この部分を明示的に指定し、値をTRUEにすることで降順になる。

さて、今回の場合は、どちらを使うべきか?確か、
行列の列と行の指定は、X[行数,列数]であった。
また、指定する行数列数は、ベクトルで与えても良かった

したがって、以下のようにすることで、特定の変数で並び替えた行列を作れる。
とりあえず、新しい変数を準備して、そこに入れてみる。

>(O.pnum <- X[order(X$pnum),])
          pnum   max hydro_pnum hydro_max thermal_pnum thermal_max atomic_pnum
Okinawa     22  1919          0         0           21        1919           0
Shikoku     65  6963         58      1141            4        3797           1
Hokkaidou   66  7419         53      1234           11        4065           1
Chugoku    110 11986         97      2906           12        7801           1
Hokuriku   138  8057        127      1904            6        4400           1
Kansai     165 34877        149      8196           12       16907           3
Tokyo      192 64988        162      8981           25       38696           3
Kyushu     194 20330        139      3279           45       11577           2
Chubu      197 32828        183      5219           11       23969           1
Touhoku    228 17206        209      2423           13       11286           2

          atomic_max
Okinawa            0
Shikoku         2022
Hokkaidou       2070
Chugoku         1280
Hokuriku        1746
Kansai          9768
Tokyo          17308
Kyushu          5258
Chubu           3617
Touhoku         3274

コマンド全体に括弧で囲まれているのは何故だったか?結果を表示するためであった。
そういったことも、少しずつ覚えていこう。

ふむ。では、発電所数で並び替えた状態で、再度、円グラフを書いてみる。
pie(O.pnum$pnum, labels=rownames(O.pnum))
title(main="Number of electric plants in Japan")

なるほど、かなり、解かりやすい図になってきた。
しかし、これでも、まだ不足していることがある。

まず、円グラフは、上を起点として時計回りに数値が大きくなるように配し、
それぞれの値をパーセンテージで表示する。
また、比率であるので、総数を明示することも忘れてはならない。

わざわざ、値を表示する理由は、人の目は角度を測るのが苦手で、
例えば、この例の場合、Chubu、Kyushu、Tokyoの大きさの差が解りにくい。
また、色に関しても、大小が区別しやすいようにしたい。では、そのように。

O.pnum.percent <- round((O.pnum$pnum/sum(O.pnum$pnum))*100,1)
O.pnum.labels <- paste(rownames(O.pnum), O.pnum.percent, sep="\n")
O.pnum.labels <- paste(O.pnum.labels, "%", sep="")
pie(O.pnum$pnum, labels=O.pnum.labels, clockwise=TRUE, col=rev(rainbow(nrow(O.pnum))))
title(main="Number of electric plants in Japan (in percentage)")
text(x=-0.8, y=-1, paste("Total number:", sum(O.pnum$pnum)))

今回は、少々、複雑に見えるかもしれない。
まず、最初の三行は、ラベルを作成するための処理である。

一行目は、百分率になるように計算している。単に、比率に100を掛けているだけ。
二行目は、paste()関数を使って、行名(対象名)と百分率の値を継ぎ合わせている。
ここでは、sep="\n" というのがあるが、\n は 改行コードで区切ることを意味している。
そして、三行目では、さらに、「%」の記号を追加して、変数を上書きしている。

円グラフの描画では、新たに「clockwise」と「col」の二つのパラメータが加わり、
labels」には、先程のプロセスで作成したラベル(O.pnum.labels)が指定されている。

clockwiseは、日本語に訳すと「時計回り」これを「TRUE」に指定すると、
0時の方向を起点として、時計回りに円グラフを描画する。

colは、その値に「rainbow(nrow(O.pnum))」が指定されている。
これは、rainbow()関数という、虹色のセットを作成する関数であり、
作成する色数を指定することで、虹色の階調で色を作成することができる。
この例の場合、対象数の数(行数)だけ色を作成すれば良いので、
nrow()関数で、行数を指定している。

最後のtext()関数では、この円グラフを作成した際の総数を表示している。
描画テキストに表示する文は、paste()関数で作られている。
円グラフ描画領域は、中心を(0, 0) とし、右上が(1, 1)左下が(-1, -1) となる領域である。
したがって、ここでは、左下付近に、総数が表示されることになる。

さて、以上が円グラフの基本的な作成方法となる。
では、このグラフから、何を読み取ることができるか、考えてみる。

この比率を見てみると、沖縄、四国、北海道の割合が低く、
逆に、割合が大きいのは、東北、中部、九州の順になっている。

直感的に、沖縄は人口が少なく、また、建設用地が少ないので数が少ないために、
発電所が多く必要としていないと考えられる。四国も同様であろう。

北海道は、土地が広いが、人口集中区が限られているために、
発電所の数が少ないのかもしれない。

割合の大きい方で見ると、東北、中部、九州が高い。
確かに、人口という視点では、低くなりなのだが。

ふむ。中々、解釈が難しそうである。発電所の数ということで考えれば、
なんとなく、日本全国に、均等に配されているようにも思える。

これだけ見ていても、よく解らない。今度は、最大出力で考えてみる。

O.max <- X[order(X$max),]
O.max.percent <- round((O.max$max/sum(O.max$max))*100,1)
O.max.labels <- paste(rownames(O.max), O.max.percent, sep="\n")
O.max.labels <- paste(O.max.labels, "%", sep="")
pie(O.max$max, labels=O.max.labels, clockwise=TRUE, col=rev(rainbow(nrow(O.max))))
title(main="Maximum output of electric plants in Japan (in percentage)")
footer <- paste("Total output:", sum(O.pnum$max))
text(x=-0.8, y=-1, paste(footer, "Kw"))

なるほど。今度は、もう少し解釈しやすそうだ。
順番は、東京、関西、中部、九州、東北・・・という順番となっている。
直感的に、人口が多そうな順番である。中部地方と九州地方高いのは、
電気の供給エリアに、大規模な工業地帯が含まれているからだろうか?

この結果から、発電所の状況を分析するには、
単に、発電所の数ではなく、総出力で見た方が良さそうである。
となると、今度は、書く地域ごとの発電所のタイプも見てみたくなる。
とりあえずは、東京圏だけで考えてみる。

# 発電所の種別で色を分ける
col.h <- rgb(50, 100, 255, maxColorValue = 255)  # 水力発電所のイメージ色
col.t <- rgb(255, 100, 50, maxColorValue = 255)  # 火力発電所のイメージ色
col.a <- rgb(200, 255, 100, maxColorValue = 255)  # 原子力発電所のイメージ色

# 今度は、総出力の合計で、パイの大きさを変えてみる
# とりあえず、最大値と最小値を控えておく
omax <- max(X[order(X$max),2])
omin <- min(X[order(X$max),2])

# order()関数の並び替えによって、東京圏のデータは10番目に入っている。
i = 10

#  パイのサイズを調整する。最大値サイズが1、最小サイズが0.5になるように調整
size = 0.5+(O.max[i,2]-omin)/(omax-omin)*0.5

# 折角なので、タイトルもオートマチックに与える
title <- paste("Rate of electric plants type in", rownames(O.max)[i])

#  水力、火力、原子力を、それぞれ、h、t、aという変数に格納する
h <- round((O.max[i,4]/sum(O.max[i,c(4,6,8)]))*100,1)
t <- round((O.max[i,6]/sum(O.max[i,c(4,6,8)]))*100,1)
a <- round((O.max[i,8]/sum(O.max[i,c(4,6,8)]))*100,1)

# ラベル名を調整する。
lbl <- paste(c("Hydro","Thermal","Atomic"), c(h,t,a), sep="\n")
lbl <- paste(lbl, "%", sep="")

# 円グラフを描く
pie(c(h,t,a), labels=lbl, clockwise=TRUE, col=c(col.h, col.t, col.a),radius=size, main=title)

この処理では、色々と、凝った工夫がされいる。
最後のpie()関数のパラメータを見てみると、radiusが加わっている。
このパラメータは、円グラフの半径を指定するパラメータであり、既定では「1」である。

これを、総出力の合計の最大値最小値を使って、0.5から1に収まるようにしている。
また、タイトルをtitle()関数を使わずに、mianパラメータで指定している。

今回の処理の特徴は、変数のi値だけを決めてやれば、
全て、自動的に書いてくれるようになっている。

例えば、上記の一連の処理の中で「i=1」とすれば、
グラフタイトルや、円グラフのサイズも含めて、自動的に沖縄のグラフができる。
(沖縄は、order()関数昇順1番目に来ている。)

さて、わざわざ、このように回りくどいことをした理由は、

複数のグラフを一気に作りたい から。

確かに、i の部分を書き換えるだけ良いとして、
それでも、10回は書きなおさないといけない。
実際に上記の方法で10個も作るのは大変
Excelで作るよりかはマシであるが。

そこで、「i」を指定も自動的にやってしまうことにする。
最終的に、以下のようなグラフが出来上がる。

とりあえず、以下のコード最後までコピーし、Rコンソール上で実行する。
今回は、一行ずつ実行できない。

# おまじない。
par(mfrow=c(2,5))

# 魔法の言葉「For文(ループ文)」
for (i in 1:10) {
    # サイズを調整する
    size = 0.5+(O.max[i,2]-omin)/(omax-omin)*0.5

    # i番目の対象名(行名)からタイトルを自動生成する
    title <- paste("Rate of electric plants type in", rownames(O.max)[i], sep="")
    title <- paste(title, O.max[i,2], sep="\n Total Maximum Outpt:")
    title <- paste(title, "Kw", sep="")

    #  水力、火力、原子力を、それぞれ、h、t、aという変数に格納する
    h <- round((O.max[i,4]/sum(O.max[i,c(4,6,8)]))*100,1)
    t <- round((O.max[i,6]/sum(O.max[i,c(4,6,8)]))*100,1)
    a <- round((O.max[i,8]/sum(O.max[i,c(4,6,8)]))*100,1)

    # ラベル名を調整する。
    lbl <- paste(c("Hydro","Thermal","Atomic"), c(h,t,a), sep="\n")
    lbl <- paste(lbl, "%", sep="")

    # 円グラフを描く
    pie(c(h,t,a), labels=lbl, clockwise=TRUE, col=c(col.h, col.t, col.a),radius=size, main=title)
}
dev.off()

なんと、上のコードを実行するだけで、10個のグラフが一瞬で出来てしまった。
ここで新しく出てきたのは、par()関数for文dev.off() の三つ。

さて、ここは踏ん張りどころ。ここからが、今回のブログの山場

ここで登場するpar()関数はそれほど難しくはない。
この関数によって、グラフ領域を色々と設定することができる。
余白設定なども含めて、かなり細かい設定が可能である。

今回設定したのは、mfrowというパラメータのみ。値は c(2,5) というベクトル。
これは、グラフを2行、5列で並べて描くという設定である。
dev.off() が実行されるまで、この設定が反映される。

少々、難しいかもしれないのが「for文」。ここでは、for(i in 1:10) となっている。
これは、「i」という変数が「1から10まで繰り返されるという意味である。

何が繰り返されるかと言うと、「」と「」で挟まれている間の処理。
この括弧で挟まれた範囲を、一つの処理のまとまりとし「10回」実行され、
この範囲に登場する変数「i」は、1サイクルする度に、一つずつ繰り上がっていく。

この処理は、一般的にループ文あるいはFor文と呼ばれるもので、
プログラミングの最も基本的な処理である。覚えておけば、色々と便利。

なお、この他に、if文という分岐処理があるが、
ループ文とif文を理解できれば、ほぼ全てのプログラミング処理が可能となる。

さて、今回は、少々長かったので、疲れたかもしれない。
実を言うと、今回は、書いている方も辛かった。正直、長かった。
このまま、今日は止めにしたいところではあるが、最後に結果の観察を。

このグラフを見てみると、いくつかのことが読み取れる。

まず、各地域とも、火力発電による発電量の比率が最も高いことが解る。
関西圏を除く、全地域で50%以上が火力発電となっている。
原子力依存という点では、四国、北海道、関西圏、九州の順に高い。

一方、原子力への依存が少ないのは、沖縄、中部、中国のようである。
中部地方の依存度が低いのは、地震との関係があるのだろうか?

中国地方は、特に、電気の消費が激しい瀬戸内では、
建設できる場所が無いのかもしれない。何故、だろうか?

水力発電に関して見てみると、中国、北陸、関西の順で比率が高い。
知らなかったが、意外に関西が高い。

前回の棒グラフのデータ観察から、一施設辺りの総出力は原子力が高いことが解ったが、
全体に対する、総出力の比率という観点からは、火力発電が主力のようである。

さらに、地域別の消費電力や発電量などのデータも合わせて観察すれば、
より多くのことが解るのだろう。まぁ、その話は、本ブログの目的ではないが。

何やら、データを見てみると、色々と想像が尽きない。
重要なことは、このような基本的なデータを眺めながら、様々な想像を巡らすこと。
すぐに、高級な分析やら、検定がどうだの、という人も居るが、
まずは、基本的なことから見ていくのが重要。

0 件のコメント:

コメントを投稿