統計実務の好著が出た：朝野煕彦『入門 多変量解析の実際』

統計実務の好著がちくま学芸文庫に入った：

朝野煕彦『入門 多変量解析の実際』（ちくま学芸文庫，2018年5月）

以前，同じくちくま学芸文庫に入っていた，林知己夫の『調査の科学』を紹介したことがある（参照）が，ちくま学芸文庫は統計学をはじめ，こういう数学の背景にある思想を伝える名著を出してくれるので，とても助かる。

◆　　　◆　　　◆

この本はもともとは講談社から2000年に発刊されていた。それから幾星霜，いまやデータサイエンス全盛期である。統計実務とでもいおうか，多変量解析のような分析手法を使いこなすことが文系理系関係なく多くの人に求められている。

分析自体はRのほか，SASやSPSSのような統計プログラム・パッケージが自動的に行ってくれる。問題は分析の目的や分析後の解釈の部分である。そういった「思想」に関わる部分をしっかり教えてくれる本が必要である。そういう意味で，『入門 多変量解析の実際』のような本は貴重。

これまで，Rをさんざん駆使して重回帰分析だのクラスター分析だのをやってきたわけだが，こういう本を読むと，自分のやってきたことの意味を再確認することができる。それだけではなく，今まで勘違いしてきたこと，知らなかったことも判明する。

◆　　　◆　　　◆

理工系書籍で知られる，朝倉書店からは朝野煕彦先生ほか監修の「シリーズ　マーケティング・エンジニアリング」が出ている。R使いのみならず，統計実務に関わっている人たちにはお勧め：

タグクラウドもしくはワードクラウドの練習（続）

Rを使ってワードクラウド（タグクラウド）を作る練習の第2弾。

今回は文章（テキストデータ）を使って，その中のキーワードを抽出し，タグクラウドで出現回数を表示してみる。

まず，今回，対象とするテキストは今まで翻訳に取り組んできた（参照）：

ウィンストン・チャーチル『マラカンド野戦軍物語』

の全文である。テキストデータはグーテンベルク・プロジェクトのページから手に入る：

"The Story of the Malakand Field Force: An Episode of Frontier War by Churchill"

テキストデータをダウンロードしたら，次は保存先に注意。

Rのワーキングディレクトリの下に"corpus"というディレクトリを作り，さらにその中に"target"を作っておき，この"target"ディレクトリの中に保存する。

そしてRを立ち上げて，以下の通り作業に移る：

＜パッケージの読み込み＞

まず，パッケージを読み込む。

# # load packages # library(ctv) library(tm) library(SnowballC) library(wordcloud)

＜テキストデータの読み込み＞

つぎに，分析対象とするテキストデータを読み込む。今回は"./corpus/target"の中に保存したので，そのディレクトリに到達できるようにパスを設定する。

そして，そのテキストデータをもとに「コーパス」を作る。コーパスは"docs"というオブジェクトに保存される：

# # Store all text files in the location(folder) "./corpus/target" # Identify text files and create a corpus # cname <- file.path(".", "corpus", "target") docs <- Corpus(DirSource(cname))

＜コーパスの整形＞

次に，コーパス内のデータを整形する。例えば，"/"とか"@"とかいった記号を半角スペースに変換したり，大文字を小文字にしたり，"a"とか"the"とか"I"とか，頻出する単語を削ったりする。

# # Data processing on the text files # # Replace symbols "/", "@", "-", with a blank space for (j in seq(docs)) { docs[[j]] <- gsub("/", " ", docs[[j]]) docs[[j]] <- gsub("@", " ", docs[[j]]) docs[[j]] <- gsub("-", " ", docs[[j]]) } # # Remove words like "a", "an", "the", "I", etc.. docs <- tm_map(docs, tolower) docs <- tm_map(docs, removeWords, stopwords("english")) docs <- tm_map(docs, removePunctuation) docs <- tm_map(docs, stripWhitespace) docs <- tm_map(docs, PlainTextDocument)

＜コーパスの構造化＞

コーパスを構造化して，wordcloud関数が読めるようにする。

# # Create structured data # dtm <- DocumentTermMatrix(docs)

＜ワードクラウドの作成＞

単語とそれぞれの登場回数をまとめたデータフレームを作り，それをwordcloud関数に読み込ませて，ワードクラウドを作成する。

# # Making the word cloud # m <- as.matrix(dtm) v <- sort(colSums(m), decreasing=TRUE) head(v, 14) words <- names(v) d <- data.frame(word=words, freq=v) wordcloud(d$word, d$freq, min.freq=100, colors=brewer.pal(8, "Dark2"))

さて，こういうのができたら完成：

"enemy"がトップ。そして軍人の階級やら兵科やらの単語が良く登場する。やはり戦記物らしい結果。

タグクラウドもしくはワードクラウドの練習

今頃知ったのだが，キーワードの出現回数を文字の大小で表現するグラフのことを「タグクラウド」もしくは「ワードクラウド」と呼ぶのだそうだ。こういうやつ：

多分，"R"でできるだろうと思って，調べてみると，"wordcloud"というパッケージを使えばできることがわかった。

そこで，ちょっと練習してみる。

◆　　　◆　　　◆

こういう例題を考えてみた：

alpha, beta, gamma, deltaという4つのキーワードがある。
ある文章中に，それぞれ40回，30回，20回，10回出現した。
これをワードクラウドで表現してみよう。

まず，R上で，alpha, beta, gamma, deltaの出現回数をまとめたデータフレームを作ってみる（参考："R-Source 39. データフレーム事始"）。

wo <- c("alpha", "beta", "gamma", "delta") fr <- c(40, 30, 20, 10) dd <- data.frame(word=wo, freq=fr)

すると，こんなデータフレーム"dd"ができる。

word freq 1 alpha 40 2 beta 30 3 gamma 20 4 delta 10

つぎに，"wordcloud"というパッケージを呼ぶ。これはRに初めから入っているわけではないので，あらかじめどこかのサイトからパッケージをインストールする。そしてライブラリーから読み込む。

install.packages("wordcloud") # ライブラリの読み込み library(wordcloud)

すでにデータフレーム"dd"ができているので，これを関数"wordcloud"に放り込めばいい。

wordcloud(dd$word, dd$freq)

すると，こんなワードクラウド（タグクラウド）が表示されるはず（並び方は実行する度にランダムに変わる）：

フォントの大きさが出現回数の多さを表している。

さて，このワードクラウド，黒一色で面白くない。そこで，色の指定をしてみる。まずは，グレースケール。

colors = c("grey80", "grey60", "grey40", "grey20")

出現回数が少ないほど，色が薄くなる設定の4段階のグレースケールである。先ほどの関数"wordcloud"に，色指定のパラメータを加えてみる：

wordcloud(dd$word, dd$freq, color=colors)

実行すると，こんな感じになる：

今度はヒートカラー「黄色→オレンジ→赤→濃赤」という指定をして実行してみる：

colors = c("yellow", "orange", "red", "darkred") wordcloud(dd$word, dd$freq, color=colors)

実行すると，こんな感じ：

ということで，wordcloudの最も簡単な例題終了。

近いうちに，自然言語処理をしてワードクラウド（タグクラウド）を作ってみる。

予習用資料：
"Build a word cloud using text mining tools of R"
「R言語で小説を読もう累計１位の「無職転生- 異世界行ったら本気だす - 」をタグクラウドしてみた。 」

Rでネパールのdistrict別人口密度の塗分け地図を作る

およそ1年ぶりにRの記事を書く。

今回はネパールの郡（district）別の人口密度の塗分け地図を作ってみた。

スクリプトは次の通り：

#おまじない rm(list=ls(all=TRUE)) # # library(maptools) par(family="sans") # x <- readShapePoly("c:\\rtemp\\NPL_adm\\NPL_adm3.shp") #シェープファイルの読み込み # ix <- length(x) # データのサイズ（行＝districtの数） zdata <- x$POP_2011/x$AREA # 人口密度の計算；POP_2011はdistrict別人口，AREAはdistrict面積 colnam <- rep("white", length=ix) # 塗分けのための変数 t <- c(50, 100, 500, 1000) # 塗分けのための閾値 color.set <- c("lemonchiffon", "yellow", "orange", "red", "darkred") # カラーパレット(heat color) unit <- "Persons/km^2" # zdataの単位 # # zdata（人口密度）に応じて，カラーパレットから色を選んで # 塗分けのための変数colnamに代入する for (ii in 1:ix) { if (zdata[ii] > t[4]) { colnam[ii] <- color.set[5] } else if (zdata[ii] > t[3]) { colnam[ii] <- color.set[4] } else if (zdata[ii] > t[2]) { colnam[ii] <- color.set[3] } else if (zdata[ii] > t[1]) { colnam[ii] <- color.set[2] } else { colnam[ii] <- color.set[1] } } # plot(x, col=colnam) # 人口密度で色分け # # 以下は凡例の描画用 oo <- par("usr") # プロット範囲（経度・緯度）をooに保存 ox <- oo[1] # 経度の値の小さい方を原点x座標に oy <- (oo[3] + oo[4])/2 # 緯度の値の小さい方と大きいほうの中間を原点y座標に oxl <- oo[2] - oo[1] # oyl <- oo[4] - oo[3] # par(new = T) # 現在の図に新たな図を重ね書きする for (ii in 1:5) { aa <- c(ox + oxl * 15 / 20, ox + oxl * 15 / 20, ox + oxl * 16 / 20, ox + oxl * 16 / 20) bb <- c(oy + oyl / 20 * ii, oy + oyl / 20 * (ii + 1), oy + oyl / 20 * (ii + 1), oy + oyl / 20 * ii) polygon(aa, bb, col=color.set[ii]) if (ii < 5) text(ox + oxl * 16.5 / 20, oy + oyl / 20 * (ii + 1), labels = t[ii], adj = 0) else text(ox + oxl * 15 / 20, oy + oyl / 20 * (ii + 1.5), labels = unit, adj = 0) }

詳細はここに記載。

出来上がりはこれ：

首都圏（カトマンドゥ郡，バクタプル郡，ラリトプル郡）の人口密度が非常に高いことがわかる。カトマンドゥなんか，1平方キロメートル当たり4500人以上である。

この地図では75郡だが，実は，2015年の憲法で行政区分がいろいろと変わり，現在ではナワルパラシ郡とルクム郡は東西に分割されて77郡となっている。

Rで雨温図(Hythergraph)を書く（その3）

前回，前々回に続き，統計解析フリーソフトRで雨温図(Hythergraph)を書くというシリーズの3回目。

今回は4都市分の雨温図を並べて書いてみる。

前回，気象データとして，ASEAN各国の首都の月別最高気温・最低気温・雨量をまとめたcsvファイルを作ってみたのでそれを再び使用する：

「aseanclimate.csv」をダウンロード

前回も述べたが，データの出典は，カンボジアやラオスに関してはそれぞれの国の統計年鑑，その他の国については世界気象機関(WMO)のデータベースである。

スクリプトは以下の通りであるが，これまでとの大きな違いは，こんな感じ：

• ASEAN諸国のうち，カンボジア，ラオス，タイ，ベトナムの首都の雨温図を書くことにした
• 関数par()を使って画面を4分割し，4都市分がまとめてみられるようにした
• for文を使って描画作業を繰り返すことにした

#おまじない rm(list=ls(all=TRUE)) # # csvデータの読み込み d <- read.csv("c:\\rtemp\\aseanclimate.csv") attach(d) names(d) # # データフレームを雨量，月平均最高気温，月平均最低気温の3つに分割 dtmax <- subset(d, Item == "tempmax") dtmin <- subset(d, Item == "tempmin") drain <- subset(d, Item == "rainfall") # # 棒グラフの横軸タイトル＝ラベル（月名）の定義 xlabel <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec") # # 対象国群の指定 countryname <- c("Cambodia", "Lao PDR", "Thailand", "Vietnam") # # 画面を2×2に分割し，4つのグラフを表示できるようにする # また，軸表示用の余白を定義する： # ベクトルの数値は余白を表し，単位は"行"。底辺，左側，上側，右側の順で定義 # フォントのサイズはpsで指定 par(mfrow=c(2,2), mar=c(4, 4, 4, 4), ps = 9) # ##### # グラフの作成 # 対象国群の数だけ繰り返す ##### for (i in 1:length(countryname)) { # dtmp <- subset(dtmax, Country == countryname[i]) cityname <- data.matrix(dtmp$City) # # 雨量，月平均最高気温，月平均最低気温のデータを代入 tmax <- data.matrix(subset(dtmax, Country == countryname[i], xlabel)) tmin <- data.matrix(subset(dtmin, Country == countryname[i], xlabel)) rain <- data.matrix(subset(drain, Country == countryname[i], xlabel)) # # # 棒グラフを表示 # 戻り値は、各バーのx座標位置 barx <- barplot(rain, xlim=c(0, 14), ylim=c(0, 800), yaxp = c(0, 800, 8), names.arg=xlabel, las=3) # # グラフの上書きを指定し， # ラベルや座標軸を非表示のまま月平均最高気温をプロット par(new=T) plot(barx, tmax, xlim=c(0, 14), ylim=c(0, 40), type="o", pch=1, axes=F, xlab="", ylab="") # # グラフの上書きを指定し， # ラベルや座標軸を非表示のまま月平均最低気温をプロット # 最後に折れ線グラフのy軸を表示 par(new=T) plot(barx, tmin, xlim=c(0, 14), ylim=c(0, 40), type="o", pch=19, axes=F, xlab="", ylab="") gtitle <- paste(cityname, ", ", countryname[i]) mtext(gtitle, side=1, line=2) mtext("Mean Total Rainfall [mm]",side=2, line=2) mtext("Mean Daily Max and Min Temp. [C]",side=4, line=2) axis(side=4, yaxp = c(0, 40, 8)) box() # ##### # グラフの作成 # 対象国群の数だけ繰り返す ##### }

こんな図ができるはず。

今回のスクリプトを応用すれば，一画面に何都市分の雨温図でも表示できる。もちろん図が小さくなるので見えにくくなるが。

Rで雨温図(Hythergraph)を書く（その2）

前回からの続き。

今回は気象データを読み込んで雨温図(Hythergraph)を書くことにする。

まず，気象データであるが，ASEAN諸国の首都の月別最高気温，最低気温，雨量をまとめたcsvファイルを作ってみた：

「aseanclimate.csv」をダウンロード

データの出典として，カンボジアやラオスに関してはそれぞれの国の統計年鑑，その他の国については世界気象機関(WMO)のデータベースを用いた。

スクリプトは次の通りで，前回との違いは，csvファイルの読み込みを行うところである。

あと，ブルネイではなく，カンボジアのプノンペンの雨温図を書くことにする。

ちなみに，小生の場合，csvファイルはcドライブ直下に作った"rtemp"というフォルダーに置いている："c:\\rtemp\\aseanclimate.csv"

rm(list=ls(all=TRUE)) # # # csvデータの読み込み d <- read.csv("c:\\rtemp\\aseanclimate.csv") attach(d) names(d) # # データフレームを雨量，月平均最高気温，月平均最低気温の3つに分割 dtmax <- subset(d, Item == "tempmax") dtmin <- subset(d, Item == "tempmin") drain <- subset(d, Item == "rainfall") # # 棒グラフの横軸タイトル＝ラベル（月名）の定義 xlabel <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec") # # 対象国の指定 countryname <- "Cambodia" dtmp <- subset(dtmax, Country == countryname) cityname <- data.matrix(dtmp$City) # # 雨量，月平均最高気温，月平均最低気温のデータを代入 tmax <- data.matrix(subset(dtmax, Country == countryname, xlabel)) tmin <- data.matrix(subset(dtmin, Country == countryname, xlabel)) rain <- data.matrix(subset(drain, Country == countryname, xlabel)) # # # 軸表示用の余白の定義 # ベクトルの数値は余白を表し，単位は"行"。底辺，左側，上側，右側の順で定義 # フォントのサイズはpsで指定 par(mar=c(4, 4, 4, 4), ps = 9) # # 棒グラフを表示 # 戻り値は、各バーのx座標位置 barx <- barplot(rain, xlim=c(0, 14), ylim=c(0, 800), yaxp = c(0, 800, 8), names.arg=xlabel, las=3) # # グラフの上書きを指定し， # ラベルや座標軸を非表示のまま月平均最高気温をプロット par(new=T) plot(barx, tmax, xlim=c(0, 14), ylim=c(0, 40), type="o", pch=1, axes=F, xlab="", ylab="") # # グラフの上書きを指定し， # ラベルや座標軸を非表示のまま月平均最低気温をプロット # 最後に折れ線グラフのy軸を表示 par(new=T) plot(barx, tmin, xlim=c(0, 14), ylim=c(0, 40), type="o", pch=19, axes=F, xlab="", ylab="") gtitle <- paste(cityname, ", ", countryname) mtext(gtitle, side=1, line=2) mtext("Mean Total Rainfall [mm]",side=2, line=2) mtext("Mean Daily Max and Min Temp. [C]",side=4, line=2) axis(side=4, yaxp = c(0, 40, 8)) box()

こんな図ができるはず。

今回，csvファイルの読み込みを行ったのは，複数都市の雨温図を並べて書くという作業を行うための伏線である。

次回（その3）では4都市分の雨温図を書くことにする。

Rで雨温図(Hythergraph)を書く（その1）

久々に統計解析フリーソフトRの話をする。

雨温図 (Hythergraph)というのは，ある地域の降水量と気温をあらわすグラフのことである。月別の気温を折れ線グラフで，雨量を棒グラフで表し，両者を重ねて表現する。

これをRで書こうとすると，ちょっと手間取る。

折れ線グラフだけ，または棒グラフだけならplot()とかbarplot()といった関数を使ってあっという間に描けるのだが，両者を重ねるのは意外と難しい。

"Symfowareについての考察blog"や定番"統計解析フリーソフト R の備忘録頁"を参考にしながらスクリプトを書いてみたので紹介する。

最近ASEAN諸国の研究をしているので，ASEAN加盟国の一つ，ブルネイの首都，バンダルスリブガワンの雨温図を書いてみる。

雨温図に必要なデータは世界気象機関(WMO)のデータベースから入手：

"World Weather Information Service - Bandar Seri Begawan"

肝心のスクリプトは次の通りである：

# 軸表示用の余白の定義 # ベクトルの数値は余白を表し，単位は"行"。底辺，左側，上側，右側の順で定義 # フォントのサイズはpsで指定 par(mar=c(4, 4, 4, 4), ps = 9) # # 棒グラフの横軸タイトル＝ラベル（月名）およびデータ（雨量） xlabel <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec") rain <- c(292.6, 158.9, 118.7, 189.4, 234.9, 210.1, 225.9, 226.6, 264.4, 312.3, 339.9, 339.6) # # 棒グラフを表示 # 戻り値は、各バーのx座標位置 barx <- barplot(rain, xlim=c(0, 14), ylim=c(0, 800), yaxp = c(0, 800, 8), names.arg=xlabel, las=3) # # 月平均最高気温と月平均最低気温のデータを代入 tmax <- c(30.4, 30.7, 31.9, 32.5, 32.6, 32.5, 32.3, 32.4, 32.0, 31.6, 31.4, 31.0) tmin <- c(23.3, 23.3, 23.5, 23.7, 23.7, 23.4, 23, 23.1, 23.1, 23.2, 23.2, 23.2) # # グラフの上書きを指定し， # ラベルや座標軸を非表示のまま月平均最高気温をプロット par(new=T) plot(barx, tmax, xlim=c(0, 14), ylim=c(0, 40), type="o", pch=1, axes=F, xlab="", ylab="") # # グラフの上書きを指定し， # ラベルや座標軸を非表示のまま月平均最低気温をプロット # 最後に折れ線グラフのy軸を表示 par(new=T) plot(barx, tmin, xlim=c(0, 14), ylim=c(0, 40), type="o", pch=19, axes=F, xlab="", ylab="") mtext("Bandar Seri Begawan, Brunei", side=1, line=2) mtext("Mean Total Rainfall [mm]",side=2, line=2) mtext("Mean Daily Max and Min Temp. [C]",side=4, line=2) axis(side=4, yaxp = c(0, 40, 8)) box()

こんな図ができるはず。

スクリプトを読むとおわかりだと思うが，棒グラフと折れ線グラフのx座標がずれないようにするために，barxというベクトルを作って，棒グラフの各バー（棒）の中心の位置を記録し，折れ線グラフのx座標として利用している。

このあたりの知恵は"Symfowareについての考察blog"を参照。

今回は雨量，最高気温，最低気温のデータを全部ベタ打ちで入力したが，次回（その2）では，csvデータを取り込んでグラフを書く，というやや進化したものをお見せする予定。

Rでサンキ―・ダイアグラム的なものを書く

エネルギーの流れを理解するためには，サンキー・ダイアグラム（Sankey diagram）という図表を用いると良い。サンキー・ダイアグラムというのはフロー図の一種で，工程間のエネルギー・物資・資金等の流量を表現するのに用いられる。

小生が関わっている分野で言うと，環境工学やエネルギー関連の研究で良く用いられる。

サンキー・ダイアグラムはわかりやすい図表だが，作るのは大変である。

何か良い手はないかと思ったら，Rにちょうどいいパッケージがあった。

"riverplot"

というパッケージがそれである。実は"sankey"というそのものズバリのパッケージもあるのだが，riverplotの方が，説明がわかりやすく，出力される図もきれいだったので，こちらを使ってみる。

◆　　　◆　　　◆

準備

まず，何としてもパッケージriverplotを入手する。

CRANにあるはずなので，探してみること。

◆　　　◆　　　◆

riverplotの練習 (1)

A, B, Cという3つのノードがあり，AからCに5，BからCに10が流れるとする。

これをriverplotでフロー図（サンキー・ダイアグラム）にしてみよう。

まず，なにも考えずに，以下のコマンドを打ち込んでいく。

library(riverplot) nodes <- data.frame( ID = c("A", "B", "C"), x= c( 1, 1, 2 ), col= c( "yellow", NA, NA ), labels= c( "Node A", "Node B", "Node C" ), stringsAsFactors= FALSE ) edges <- data.frame( N1 = c("A", "B"), N2 = c("C", "C"), Value = c(5, 10)) r <- makeRiver( nodes, edges ) plot( r )

最後に"Enter"キーを押して，次のような図が出てきたらOK。

◆　　　◆　　　◆

nodesの説明

さて，コマンドの解説を始める。

どんなフロー図を書くのか，という定義はmakeRiver()という関数を用いて行う。

makeRiver()には最低限，nodesとedgesという2つの変数が必要である。

nodesというのはノード，つまり結束点を定義したデータフレームである。

今回はA, B, Cという3つのノードがある。それぞれのID（識別コード），横軸(x)座標，色(col)，名前(labels)を次のように決めよう：

• ノードA：IDを"A"，横軸(x)座標を1とする。色(col)を黄色とし，名前(labels)を"Node A"とする。
• ノードB：IDを"B"，横軸(x)座標を1とする。色(col)をデフォルトの灰色とし，名前(labels)を"Node B"とする。
• ノードC：IDを"C"，横軸(x)座標を2とする。色(col)をデフォルトの灰色とし，名前(labels)を"Node C"とする。

ノードごとにIDやx座標を定義するのではなく，ID，x座標といった属性ごとに定義を行うことにする。

すると，

ID = c("A", "B", "C")
x= c( 1, 1, 2 )
col= c( "yellow", NA, NA )
labels= c( "Node A", "Node B", "Node C" )

となる。色(col)に関しては，NAと書いておけば，デフォルトの灰色が指定される。

これら4つの属性をnodesに収納するため，data.frame関数を使用する：

nodes <- data.frame( ID = c("A", "B", "C"),
x= c( 1, 1, 2 ),
col= c( "yellow", NA, NA ),
labels= c( "Node A", "Node B", "Node C" ),
stringsAsFactors= FALSE )

最後に

stringsAsFactors= FALSE

というのがあるが，今回はこれを「まじない」とでも言っておく。

これで，ノードに関する定義が終わった。

edgesの説明

つぎに，edgesという変数の説明。

edgesというのはノード間の流れを定義したデータフレームである。

今回の例では，AからCに5，BからCに10の何か（水かもしれないし，金かもしれない）が流れるという設定である。

edgesというデータフレームでは，流れの出発点をN1と呼び，流れの終着点をN2と呼ぶ。N1からN2への流量をValueと呼ぶ。

すると，

• AからCに流量5の流れがあることを：N1 = "A", N2 = "C", Value = 5
• BからCに流量10の流れがあることを：N1 = "B", N2 = "C", Value = 10

と表すことになる。

流れごとにN1, N2, Valueを定義するのではなく，N1, N2, Valueといった属性ごとに定義を行うことにする。

すると，

N1 = c("A", "B")
N2 = c("C", "C")
Value = c(5, 10)

となる。これら3つの属性をedgesに収納するため，data.frame関数を使用する：

edges <- data.frame( N1 = c("A", "B"),
N2 = c("C", "C"),
Value = c(5, 10))

これで，ノード間の流れを定義できた。

nodesとedgesの定義が終わったので，makeRiver()関数を使って，全体の流れの情報をrという変数に代入する：

r <- makeRiver( nodes, edges )

これで，あとは

plot( r )

と打ち込めば，フロー図が出力される。

Rで宇部市の町丁字別の人口密度を塗り分ける（making a choropleth map）

久々にRの話をする。

今回はRを使った地図の塗り分け（階級区分図，コロプレス図，choropleth map）を行う話。

実はだいぶ前にも同じ話題を取り扱った（「Rでカンボジアの塗り分け地図を作る」2011年11月28日）が，今回は自作コード（青木繁伸先生作成のコードを改造したもの）を利用するのではなく，Rのパッケージ"foreign"，"maptools"，"sp"を利用する。

◆　　　◆　　　◆

準備

コロプレス図を作成するためのデータを，あらかじめ政府統計の総合窓口にて入手するところから始める。

政府統計の総合窓口／地図で見る統計（統計GIS)：
http://e-stat.go.jp/SG2/eStatGIS/page/download.html

• 「Step1:統計調査（集計）を選択」の中から「平成22年国勢調査（小地域）2010/10/01」を選択する。
• 「Step2:統計表を選択」の中から，人口密度が求められそうなデータを選ぶ。「男女別人口総数及び世帯総数」がいいだろう
• 「統計表各種データダウンロードへ」をクリック

ここまで進むとまた新たなページ「統計表各種データダウンロード」に移動する。

• 「Step3:地域選択」で都道府県と市区町村を選択する。今回は山口県，35202宇部市の組み合わせで選択
• 「検索」ボタンを押すと統計データと境界データが表示される
• 境界データには町丁字の境界データ（形状データ，ポリゴン）と国勢調査による町丁字別人口，世帯数などのデータがまとめて入っている。どんなデータが含まれているのかは「定義書」をちゃんと読むこと
• 今回は「日本測地系緯度経度・Shape形式」の「宇部市(620KB)」を選択（クリック）してzipファイルをダウンロードする

zipファイルをダウンロードしたら，解凍して適当なところに保存する。

小生の場合はRによる作業のディレクトリとして"c:/rtemp/"というフォルダを使っているので，"c:/rtemp/JapanGIS/"というフォルダを作って，解凍したものを保存している。

正確に言うと，"c:/rtemp/JapanGIS/A002005212010DDSJC35202/"というのが宇部市のshapeデータを含むファイル群が保存されているフォルダとなる。

◆　　　◆　　　◆

ここからがRでの作業になる。

あらかじめRのパッケージ"foreign"，"maptools"，"sp"を入手しておくこと。

そうしたら以下のコマンドを打ち込んでいく。

ls() rm(list=ls(all=TRUE)) library(foreign) library(maptools)

いよいよ宇部市のshapeデータの取り込みである。

maptemp <- readShapePoly("c:/rtemp/JapanGIS/A002005212010DDSJC35202/h22ka35202.shp") plot(maptemp)

するとこういう図が登場するはず（↓）

困ったことに海面まで含まれている。「定義書」を読むとこれを削るためのヒントが得られる。

今作ったデータセット"maptemp"には"HCODE"という陸地と水面とを区別する変数が入っているので，これを利用すれば海面を取り除くことができる。HCODE=8154が水面調査区であるから，

maptemp1 <- subset(maptemp, maptemp$HCODE != 8154) plot(maptemp1)

とやれば，海面が消えたことがわかるだろう。

次にデータセットmaptemp1に人口密度のデータを加える。「定義書」によれば，変数"AREA"が町丁字の面積[m^2]，変数"JINKO"が町丁字別人口を表しているから，新たに人口密度を表す変数"DENS"を作って，データセットmaptemp1に加えてやればよい：

maptemp1$DENS <- maptemp1$JINKO / maptemp1$AREA * 1000000

ここで，"* 1000000"としたのは人口密度の単位を[人/km^2]にするため。

さて，ついに町丁字別人口密度の塗り分けを行う。

spplot(maptemp1, "DENS", col.regions = heat.colors(100))

・・・なんかうまくいっていませんね。色の違いが判らない。全面的に真っ赤。

これは極端に人口密度が高い場所が生じているため。思い切って，表示する人口密度の上限を10000人/km^2とする処理を加えてみる。

maptemp1$DENS <- ifelse(maptemp1$DENS > 10000, 10000, maptemp1$DENS) spplot(maptemp1, "DENS", col.regions = heat.colors(100))

はい。ましになりました。

このあと，塗り分けのカラーパレットをいろいろと変えてみる。

spplot(maptemp1, "DENS", col.regions = topo.colors(100)) spplot(maptemp1, "DENS", col.regions = terrain.colors(100)) spplot(maptemp1, "DENS", col.regions = grey(0:99/100))

"terrain.colors"を使ってみた結果は次のようになる（↓）：

小生としては"terrain.colors"の方が"heat.colors"よりも気に入っている。

というわけで，結構楽にコロプレス図が書けることがわかった。

林知己夫『調査の科学』を読む

数学の中で難しいのは確率・統計ではないかと思っている。微分・積分，代数・幾何など，難しいようだが，計算結果などを巡って意見が対立することは無い。

ところが確率・統計では計算結果を巡ってもめることが多々ある。計算結果の何を巡ってもめているかというと，解釈を巡ってもめているのである。

例えば，小生は統計言語"R"を使っているわけだが，これを使えば，様々な統計的検定が簡単にできる。しかし，その検定結果に対しては，ネイマン・ピアソン流かフィッシャー流かで，解釈の仕方ががらりと変わるわけである。

結局，確率・統計分野では，数学の他の分野に比べて常々勉強しなくてはいけないことが山のようにあると思う。

というわけで，時折勉強のために確率や統計の本を読んでいるのだが，最近，質的データを量的に扱うという統計的手法「数量化理論」を生み出した偉大な統計学者である林知己夫の『調査の科学』を読んでいる。

調査の科学 (ちくま学芸文庫) 林 知己夫 筑摩書房 2011-05-12 売り上げランキング : 484579 Amazonで詳しく見る by G-Tools

この本で凄いと思ったのは統計用語のユニバースと母集団(population)を厳密に分けていることである。

この辺が詳しく書いてあるのが本書の第2章である。以下，主要部分を引用する。

「定義がはっきりし，しかもその構成要素が具体的にとらえられるような調査対象集団を，統計学ではユニバース universe と呼んでいる。」（42ページ）

「母集団とはユニバースに確率的な概念を加えたものである」（43ページ）

「ユニバースと母集団が同じではないか，と思われるかもしれない。＜中略＞　しかし，厳密な考え方としては，ユニバースと母集団はまったく別のものであることを認識していただきたい。」（43ページ）

「ユニバースに，要素を取り出すという考え方を加えたのが，「母集団」である。」（44ページ）

「日本人という調査集団では，ユニバースであった住民登録簿から機会均等に，つまり公平なくじ引きで，むずかしくいうと同じ抽出確率でつぎつぎと人を抽出すると決めたのが母集団ということになる」（44ページ）

さて，こうした説明を聞いたのち，ウェブ上でユニバースと母集団について解説している文章を読むと，なんか違う。

例えばこういう説明：

「統計学上の population と universe の違いは、前者が、検定や推定の対象であるところの個体から観測される値（小学校6年生のお小遣い） の集まりであるのに対して、後者は個体（小学校6年生）の集まりという点です。」（コラム『統計備忘録』バックナンバー第34話「sample, population, universe」，統計WEB，2008年1月15日）

『調査の科学』の解説を書いた吉野諒三によれば，ユニバースと母集団の区別について厳密に取り扱っているのは，本書を除くと林文ほか『統計学の基本』（朝倉書店，1991）や杉山明子編『社会調査の基本』（朝倉書店，2011）ぐらいしか存在しないようである。

社会調査の基本 杉山明子編著 朝倉書店 2011-03-18 売り上げランキング : 288849 Amazonで詳しく見る by G-Tools

通常はユニバースと母集団の対応関係が固定しているために，両者を区別する意義があまりないようだが，アンケートの回収率の低下などの問題が生じている現在では，ユニバースと母集団について，厳密な定義にまでさかのぼって慎重に検討する必要がある，ということを吉野諒三は述べている。