2014年12月27日土曜日

Rで文字列を日付時刻型に変換する

日付時刻と観測値が対応付いたデータってありますよね。例えば、原子力規制委員会が公開している放射線モニタリング情報だと下記のようなCSV形式で、測定時刻と測定値が対応付いています。


                V5,      V6
  2014/12/26 23:50,   0.033
  2014/12/26 23:40,   0.034
  2014/12/26 23:30,   0.033
  ................,   .....


こういうデータを見ると、V5を横軸、V6を縦軸にして時系列グラフを描きたくなりますよね。でも、V5が文字列型で読み込まれていると、そのままplotのx軸に指定してもうまくいきません。

日付や時刻を表現する型に変換してやる必要があります。

上記のような書式(日付がスラッシュ、時刻がコロン)ならば、POSIXlt関数に引数で渡すだけで、簡単に変換することができます。

  > as.POSIXlt("2014/12/26 23:50")

  [1] "2014-12-26 23:50:00 JST"

まとめて処理するなら、↓こんな感じでしょうか。

  d$日付時刻 <- as.POSIXlt(d$V5)


文字列が標準の書式に従っていれば、上記のように簡単です。

では、文字列が独自の書式で書かれている場合はどうでしょうか。例えば、日付時刻の情報が「2014年12月27日 00時00分」のような書式で入っているとします。このような場合は、strptimeを使うと便利です。

  > strptime("2014年12月27日 00時00分", "%Y年%m月%d日 %H時%S分")

  [1] "2014-12-27 JST"

%Y: 4桁の西暦
%m: 2桁の月
%d: 2桁の日
%H: 2桁の時
%M:  2桁の分

となっているので、それぞれの要素が文字列のどこに登場しているかを第2引数で指定しているわけですね。

「パーセントほにゃらら」として何が指定できるかは、ヘルプに書かれています。コマンドプロンプトに「?strptime」と打てば確認できます。

as.POSIXltやstrptimeをつかって日付時刻型に変換してしまえば、plotに指定することで、時系列のグラフを描けるようになります。■

[2015年1月2日追記]
実際のデータの時系列プロットで使ってみたサンプルです↓
Rで放射線モニタリング情報を時系列プロットする - Rプログラミングの小ネタ




2014年11月22日土曜日

Rで度数分布表を作る

まずは今回のサンプル用として、身長っぽいデータを作ってみます。

a <- round(rnorm(30, mean=170, sd=5), 1)
a
 
 [1] 173.0 168.6 168.5 164.2 167.3 170.6 162.7 168.0 170.8
[10] 175.7 166.0 166.5 162.4 172.6 170.1 170.2 164.2 167.1
[19] 163.8 163.4 168.1 168.7 171.1 164.6 166.3 177.0 170.0
[28] 173.3 170.7 169.0


このaに対して、度数分布表を作りたい。

質的データ(カテゴリカルデータ)なら、table関数を使って度数を集計できるのですが、量的データに対してtableを使うと、↓こんな感じになってしまいます。

table(a)
 
a
162.4 162.7 163.4 163.8 164.2 164.6   166 166.3 166.5 167.1 
    1     1     1     1     2     1     1     1     1     1 
167.3   168 168.1 168.5 168.6 168.7   169   170 170.1 170.2 
    1     1     1     1     1     1     1     1     1     1 
170.6 170.7 170.8 171.1 172.6   173 173.3 175.7   177 
    1     1     1     1     1     1     1     1     1


量的データだとうまくいかないですね。

Rには、お馴染みhist関数があるんで、ヒストグラムは簡単に描けます。

hist(a)




実は、このhist関数の戻り値のオブジェクトに各階級の度数が入っているんです。知ってました?

h <- hist(a)
h
 
$breaks
[1] 162 164 166 168 170 172 174 176 178
 
$counts
[1] 4 4 5 6 6 3 1 1
 
$density
[1] 0.06666667 0.06666667 0.08333333 0.10000000 0.10000000
[6] 0.05000000 0.01666667 0.01666667
 
$mids
[1] 163 165 167 169 171 173 175 177
 
$xname
[1] "a"
 
$equidist
[1] TRUE
 
attr(,"class")
[1] "histogram"


breaksが階級を区切る値で、countsが度数ですね。

つまり、155~160の範囲の度数が1、160~165の範囲の度数が7、・・・のような意味です。

細かく言うと、境目のところは

  155 < 階級に含まれる値 ≦ 160

のように扱われています。

これらの値を使えば、度数分布表が作れそうです。

前置きが長くなりましたが、↓こんな感じでどうでしょう?

h <- hist(a)
n <- length(h$counts) # 階級の数
class_names <- NULL # 階級の名前格納用
for(i in 1:n) {
  class_names[i] <- paste(h$breaks[i], "~", h$breaks[i+1])
}
frequency_table <- data.frame(class=class_names, frequency=h$counts)
frequency_table
 
       class frequency
1 162164         4
2 164166         4
3 166168         5
4 168170         6
5 170172         6
6 172174         3
7 174176         1
8 176178         1


データフレームの形で度数分布表ができました。

階級の取り方を指定したいときには、hist関数を呼ぶときにbreaksオプションで指定すればOKです。↓こんな感じ。

h <- hist(a, breaks=seq(155, 185, 10))


データフレームの形になってしまえば、あとは好きなように使えますよね。




2014年11月3日月曜日

Rを使ってシェープファイルをKMLファイルに変換する(土砂災害危険箇所マップ)

↓こちらで公開したKMLファイルは、シェープファイルをRで変換処理したものなのですが、その手順です。

神奈川県の土砂災害危険箇所ハザードマップ(Google Earth閲覧用KML形式): 主張

元にしたシェープは↓こちらで入手しました。

国土数値情報ダウンロードサービス

<災害・防災>のカテゴリに、「土砂災害危険箇所」というのがあります。シェープファイルは都道府県別に分かれて公開されています。

データ形式のタブ(?)で、「JPGIS2.1」のところを選択すると、JPGIS2.1形式とともになぜかShapeファイルも置いてあるという、微妙なトラップがあるので要注意です。

このShapeファイルをRでKMLに変換するんですが、スクリプトはごく簡単です↓

# 必要なパッケージの読み込み
library(maptools)
library(plotKML)
 
# 座標系の情報を持つオブジェクトを作る
# longlatは緯度経度指定、WGS84は世界測地系であることを表す
pj <- CRS("+proj=longlat +datum=WGS84")
 
# シェープファイルの読み込み
hazard.area <-readShapePoly("xxx.shp", proj4string=pj)
 
# もし、表示して確認したいときはplotすればOK
# plot(hazard.area)
 
# KMLファイルとして出力する
kml(hazard.area, file = "xxx.kml")


あとは、GoogleEarthで開けば閲覧できます。

↓こんな感じですね。
国交省が公開している「土砂災害危険箇所」シェープをRでKMLに変換後、Google Earthで表示したもの



閲覧時の色の変更などは、冒頭のリンク先を参照してみてください。■

[2014年11月23日追記]東京都、埼玉県、千葉県のKMLファイルも追加しました。あと、Google Earthの設定で色を変更しなくて済むよう、KMLファイル内のタグで色を赤に変えときました↓
土砂災害危険箇所のKMLファイル(東京、神奈川、埼玉、千葉): 主張




2014年10月5日日曜日

Rで数式を微分(多項式関数、三角関数、指数関数、対数関数)

Rで関数を微分することができます。数値的にではなく解析的にできます。

やり方はexpression関数で数式をセットしておいて、Dという関数を呼び出せばOK(プログラミング言語の「関数」と、数学の「関数」がごっちゃになる・・・)

高校で習うようなやつをいろいろ試してみました。

> # 多項式の関数の微分
> f <- expression( x^4 + 2*x^3 + 3*x^2 + 4*x + 5 )
> D(f, "x")

4 * x^3 + 2 * (3 * x^2) + 3 * (2 * x) + 4

できてる、みたいですね。数式を簡単にするところまではやってくれないので、ごちゃごちゃしているけど。

Dに「どの文字を変数とするか」("x"のところ)を指定するので、変数以外の文字を含んでいてもOKです。

> # 微分する変数以外の文字を含んでいてもOK
> f <- expression( a*x^n )
> D(f, "x")

a * (x^(n - 1) * n)

↓三角関数もいけます。

> # 三角関数の微分
> f1 <- expression( sin(x) )
> f2 <- expression( cos(x) )
> f3 <- expression( tan(x) )
> D(f1, "x")

cos(x)
> D(f2, "x")
-sin(x)
> D(f3, "x")
1/cos(x)^2

↓指数関数も。

> # 指数関数の微分
> f1 <- expression( a^x )
> f2 <- expression( exp(0)^x )
> D(f1, "x")

a^x * log(a)
> D(f2, "x")
exp(0)^x * log(exp(0))

↑exp(0)ってのは自然対数の底のeのことなので(e^0のことね)、exp(0)^x を微分したら exp(0)^x となる、というのを期待したのですが、その部分(log(exp(0)) = 1)は計算してくれませんでした。

> # 対数関数の微分
> f1 <- expression( log(x) )
> D(f1, "x")

1/x

↑自然対数はいけたのですが、常用対数はなぜか無理でした↓

> # 常用対数はできない
> f2 <- expression( log10(x) )
> D(f2, "x")

 以下にエラー D(f2, "x") :  関数 'log10' は導関数の表中にありません

以下は、計算時のテクニックの類ですね。これもやってくれます(文字列的に淡々と、という感じですが)。

> # 合成関数の微分
> f <- expression( sin(x^2) )
> D(f, "x")

cos(x^2) * (2 * x)

↑ x^2を1つの文字みたいにみなしてsin(t)という形にしておいて、全体を微分したものと中身を微分したものを掛け合わせるというやつ( cos(t) * t' )ですね。

> # 積の導関数
> f1 <- expression( x^2 * sin(x) )
> D(f1, "x")

2 * x * sin(x) + x^2 * cos(x)

↑ これもおなじみ、(fg)' = f'*g + f*g' というやつですね。

> # 商の導関数
> f1 <- expression( x^2 / sin(x) )
> D(f1, "x")

2 * x/sin(x) - x^2 * cos(x)/sin(x)^2

↑ これは出てきた数式の形を見ると、商の導関数の公式を使ったというよりは、x^2 * (sin(x))^(-1) と解釈したうえで、積の導関数の公式を使ったような処理っぽいです。

というわけで、一通りできた感じですが、この機能を何に使えるかいまいちピンときてません。数学の宿題を解くのには、間違いなく使えないでしょう。自分でやったほうが楽だもん。

数値的に処理すると誤差が心配な場合に、解析的に導関数を求めておいて、変数に値を代入して使うような場面があるのかもしれません。ないのかもしれません・・・。




2014年10月4日土曜日

Rで円を描く方法あれこれ

グラフに円を描くような関数が用意されていてもよさそうなものですが、どうもないみたいですので、自力で描く必要がありそうです。

簡単のため、半径1の円の例です。

まず思いつくのは、

  x^2 + y^2 = 1

をyについて解いて、

  y = ±√(1 - x^2)

という関数の形にして、curve関数を使って描く方法でしょうか↓

curve関数で上側と下側を別々に描く

curve( sqrt(1-x^2), xlim=c(-1,1), ylim=c(-1,1), asp=1)
curve(-sqrt(1-x^2), xlim=c(-1,1), ylim=c(-1,1), add=T)


媒介変数表示を使うやり方もありますね。

 x = cosθ
 y = sinθ

の、x と y の組を plotしていく感じです↓

媒介変数で座標を求めたあとにプロット

theta <- seq(-pi, pi, length=100)
plot(cos(theta), sin(theta), type="l", asp=1)


ちなみに type=l オプションを忘れると↓こんな感じになります。

点でプロットすると理屈が分かりやすい

plot(cos(theta), sin(theta), asp=1)


何度も使う場合は、中心座標や半径を引数として関数化しておくといいですね↓
こんなところに隠れミッキーが!(自分で描いたくせに)

plot.circle <- function(x, y, r){
  theta <- seq(-pi, pi, length=100)
  points(x + r*cos(theta), y + r*sin(theta), type="l", asp=1)
}
plot(c(-3,3), c(-3,3), type="n", asp=1) # 作図領域を用意
plot.circle(0, -1, 2)       # 顔?
plot.circle(-1.9, 1.6, 1.2) # 左耳?
plot.circle( 1.9, 1.6, 1.2) # 右耳?


plot関数でなく、polygon関数を使ってもいいですね↓

polygon関数で線を引いたもの

theta <- seq(-pi, pi, length=100)
plot(c(-1,1), c(-1,1), type="n", asp=1)
polygon( cos(theta), sin(theta))


polygon関数だと塗りつぶす色も指定できますね↓

polygon関数だと塗りつぶし指定もできる

polygon( cos(theta), sin(theta), col="gray")


そして、↓これはどうだ?

textで丸の文字「○」をプロット

plot(c(-1, 1), c(-1, 1), type="n", asp=1)
text(0, 0, "○", cex=65)


text関数で「○」という文字を置いてみました。線が太いし、任意の半径にできないし、これはダメですね・・・

次は絶対におすすめしない、すごくまどろっこしい方法

contour関数の等高線で描かれた円

n <- 100 # メッシュ分割数
x <- seq( -1, 1, length = n )
y <- seq( -1, 1, length = n )
z <- matrix(0, n, n)
for (i in 1:n) {
  for (j in 1:n) {
    z[i, j] <- x[i] ^ 2 + y[j] ^ 2 # ここに陰関数を書く
  }
}
contour(x, y, z, drawlabels = F, levels=1, asp=1)


Rのcontour(等高線)を使って陰関数のグラフを描く - Rプログラミングの小ネタで紹介した方法なんですが、格子状の各点について x^2 + y^2 の値を算出したあとに、値が1のところに等高線を引くという複雑な方法です。


任意の陰関数に適用できるところが便利なんですが、円みたいな簡単な方程式に使うメリットは皆無ですね。

最後の2つは蛇足でした・・・




2014年9月19日金曜日

colClassesで特定の列の読み込み型だけを指定、他の列はRの判断にまかせる(read.csv、read.table)

例えば↓こんなCSVデータがあったとして、

ファイル名: "data.csv"


id,name,age
001,taro,41
002,hanako,40
003,ichiro,11
004,jiro,2

↓ごく普通に読み込むと・・・

> data <- read.csv("data.csv")
> data
  id   name age
1  1 ichiro  41
2  2 hanako  40
3  3   taro  11
4  4   jiro   2


idの先頭につけていた00がなくなってしまうんですよね。

str関数で型を見てみると、

> str(data)
'data.frame':   4 obs. of  3 variables:
 $ id  : int  1 2 3 4
 $ name: Factor w/ 4 levels "hanako","ichiro",..: 2 1 4 3
 $ age : int  41 40 11 2


idはint型として読み込まれています。idのような用途だと先頭のゼロは残しておきたいってときがありますよね。

そんなときは、colClassesで読み込み時の型を指定してやればいいです。

> data <- read.csv("data.csv", colClasses="character")
> data
   id   name age
1 001 ichiro  41
2 002 hanako  40
3 003   taro  11
4 004   jiro   2


文字列として読み込まれ先頭のゼロも残りました。ただ、上記のように一律「文字列」と指定してしまうと、ageまで文字列になってしまい、↓数値的な処理ができなくなってしまいます。

> mean(data$age)
[1] NA
 警告メッセージ: 
In mean.default(data$age) :
   引数は数値でも論理値でもありません。NA 値を返します 


で、それに対してはすべての列に対して読み込み時の型を指定してやることで対処できます。

> data <- read.csv( "data.csv",
                    colClasses=c("character", "character", "integer"))
> str(data)
'data.frame':   4 obs. of  3 variables:
 $ id  : chr  "001" "002" "003" "004"
 $ name: chr  "ichiro" "hanako" "taro" "jiro"
 $ age : int  41 40 11 2


ageは整数にしておいたので、平均も計算できます。

> mean(data$age)
[1] 23.5


が、このcolClassesにすべての型を明に指定するというのは、列が10個も20個もある場合はかなり面倒くさいです。

やっと本題に来ました。

やりたいこととしては、「idの列は文字列で読み込みたいけど、その他の列はRの判断に任せるよ」ってな感じですね。

そんなときは、Rのデフォルトで読み込んでほしい列に対応するベクトルの要素を、NAにしておけばOKです。

> data <- read.csv("data.csv", colClasses=c("character", NA, NA))
> str(data)
'data.frame':   4 obs. of  3 variables:
 $ id  : chr  "001" "002" "003" "004"
 $ name: Factor w/ 4 levels "hanako","ichiro",..: 2 1 4 3
 $ age : int  41 40 11 2


↑1列目は強制的に文字列で読み込みましたが、他の列の判断はRに任せました。

ちなみに、例はread.csvで書いていますが、read.tableでも同様の方法が使えます。

列がたくさんあってNAを書くのが面倒な場合は、NAのベクトルを作っておいて、指定したい箇所だけを上書きすれば楽かもしれません。

CC      <- rep(NA, 20) # NAだけのベクトル
CC[1]   <- "character" # 1つ目だけ"character"に書き換える
manycol <- read.csv("manycol.csv", colClasses=CC)


たくさん列のあるデータを読み込むんだけど、使うのはいくつかの列だけ、なんてことも結構ありますからね。




2014年9月18日木曜日

Rのコンソール画面をクリアするショートカットキー

以前のコマンド入出力を確認したいときは、スクロールバーでコンソール画面をスクロールさせればいいですよね。

でも、要素数が非常多いリストやベクトルを表示させたりしたあとだと、スクロールバーのつかむところがちっちゃくなっていて、行き過ぎたり戻り過ぎたりと、なかなか該当箇所が見つからなかったりしますよね。

で、いったんRのコンソール画面をリセットしたいってときは、

  Ctrl + L

を入力すれば、コンソール画面がクリアされます。

という分かり切った操作ですが、あえてスクリーンショットをつけてみました。

コンソールがごちゃごちゃしてきたイメージ

「Ctrl + L」で、↑が、↓になります

すっきりとクリア

ちなみにショートカットキーでなくても、メニューの「編集」→「コンソール画面を消去」としても同じことができます。「コンソール画面を消去」だと、コンソールのウィンドウ自体が閉じられてしまいそうですが、そんなことはないので大丈夫です。

(以下は蛇足)

私は結構長い間このコマンドを知らなくて、過去の出力結果が邪魔な時は、たくさんリターンを押して、表示結果に「切れ目」を入れて、該当箇所を見つけやすくするという、面倒なことをしていました。

で、放送大学の「データからの知識発見」のテレビ放送を見ているときに、講師の秋光先生がRでの処理を実演したりするんですが、その時に「画面がいっぱいになってきたので、いったん消去して・・・」みたいな感じで、コンソールをクリアしていたんですね。


で、「なんだ、あの操作は?」って感じで、この便利なショートカットキーを知りました。上級者の書いたスクリプトだけでなく、操作を見るのも勉強になるなあと思いました。■


 データからの知識発見 (放送大学教材)
データからの知識発見 (放送大学教材)




2014年9月16日火曜日

Rのグラフで軸の目盛りの刻み幅を変更する方法

Rでplotなどを使ってグラフを描くとき、x軸やy軸の目盛りは勝手に調整してくれて、大抵の場合はそれで問題ないのですが、たまにちょっと変えたい時があります。そのたびに必死で検索して調べているような気がするので、ここに書き留めておきます。

(最初の方の例はcurve関数を使ってますが、plot関数でも同様の方法でいけます)

もちろん、特に指定しなくても、目盛りと目盛りのラベル(例えば、x軸の -4, -2, 0, 2, 4)は書いてくれます↓
指定しなくても目盛りは表示される
# 目盛りの幅に関しては指定しない
curve(dnorm, xlim=c(-4,4))


上記では、2刻みになっていますが、例えば、1刻みにしたいなあという場合↓
目盛りが1刻みになるように指定
# -4から4までを8区間に分ける(1刻み)
curve(dnorm, xlim=c(-4, 4), xaxp=c(-4, 4, 8))


さらに刻みを細かくして0.5にしてみましたが、目盛りのラベルは適宜間引かれてしまうようです↓(作図領域を横に伸ばせばラベルが表示されるようになります)
目盛りが0.5刻みになるように指定
# 16区間に分けると0.5刻みになるが
# すべての目盛りに数値ラベルが書かれるわけではない
curve(dnorm, xlim=c(-4, 4), xaxp=c(-4, 4, 16))


いくつに分割するかちゃんと考えないと、↓こんなことになることも・・・
割り切れないところに目盛りがきてしまう例
# 8の幅を9で割るとか・・・
curve(dnorm, xlim=c(-4, 4), xaxp=c(-4, 4, 9))

次は違うやり方。

最初は軸の目盛りを書かないようにしておいて、あとから、axis関数をつかって目盛りを書き足すという方法です↓
axis関数で目盛りを後から追記したもの(見た目おんなじだけど・・・)
curve(dnorm, xlim=c(-4, 4), xaxt="n") # 最初は目盛りを書かない
axis(side=1, at=-4:4) # atに目盛りをベクトルで指定する

side=1というのがx軸の意味で、y軸に追記したい場合はside=2と指定すればOK。

xaxpでの指定は等間隔でしたが、axisのatでの指定は目盛り位置を列挙するような形ですので、任意の間隔で目盛りをつけることが可能です。

こんないいかげんな刻み幅も↓
こんな目盛りのグラフはイヤだ
# axis関数なら等間隔でない目盛りもいけちゃう
curve(dnorm, xlim=c(-4, 4), xaxt="n")
axis(side=1, at=c(-4, -2, 1, 2.5, 3.2))




上記はcurve関数の例でしたが、plot関数でも同様です↓
関数を定義したあと、plotで描いたもの(xaxp使用)
f <- function(x){dnorm(x)}
plot(f, xlim=c(-4, 4), xaxp=c(-4, 4, 8))

関数を定義したあと、plotで描いたもの(axis使用)
f <- function(x){dnorm(x)}
plot(f, xlim=c(-4, 4), xaxt="n")
axis(side=1, at=-4:4)

デフォルトの目盛り刻みで問題ない場合がほとんどでしょうから、滅多に使わない機能かなあと思います。滅多に使わないがゆえに、やり方を忘れてしまうというジレンマ。




2014年9月10日水曜日

RでAKBの年齢、身長、スリーサイズを分析する

今回は、普段Rで解析を行っているような人にはなんてことない話です。

データの題材を変えると興味が俄然湧いてくるような人もいるんじゃないの、って感じのテーマです。

データは↓こちらのものを使わせてもらいました。

アイドルプロフィール(スリーサイズ、カップ情報) - AKB48


ページ自体はHTMLのテーブルですが、大抵のブラウザではコピー&ペーストすれば、タブ区切りやカンマ区切りのテキストに簡単に変えられるんじゃないでしょうか。そんな感じでタブ区切りのデータファイルを作りました↓

"AKB.txt"

名前        年齢   身長     B     W     H
岩田華怜      15    159    82    62    85
菊地あやか    20    160    74    58    82
佐藤すみれ    20    166    78    57    84
篠田麻里子    27    168    87    57    85
高橋みなみ    22    148.5  74    56.5  81
・・・


サイトからデータを取得したのが2013年なので、情報が古くなっていると思われます。まず年齢はそうだし、あと激太りした人とか(いるのか?)。

平面上に各メンバーをマッピングしてみるなんてのが、なんとなく面白そうですね。

RでAKBの多次元尺度構成法
d1 <- read.table("AKB.txt", header=T, row.names=1)
d2 <- dist(d1)
d3 <- cmdscale(d2)
plot(d3, type="n")
text(d3, rownames(d3))


多次元尺度構成法を使ってみました。データを(年齢、身長、バスト、ウエスト、ヒップ)の5次元だととらえて距離を算出し、なるべくその距離感が再現されるように、2次元に落として(5次元のままじゃ表現しづらいので)プロットされた図という感じでしょうか。

結果はあまりあてにできないと思います。距離の算出方法もいろいろな手法があったり、正規化する/しないでも結果は変わるし。描画された結果を見て、適切な処理を選ぶといいと思いますが、なにせ私はAKBのメンバーは5人くらいしか知らないので、結果の妥当性を検証できません。AKBに詳しい人は、ぽっちゃりゾーン、のっぽゾーン、お局ゾーン(いちおう年齢という要素も入っているので)などを見つけられるかもしれません。

単純なplotは、x座標、y座標を指定すると散布図が描かれるような使い方ですが、このような多次元のデータの場合はいきなりplotを使うだけで、それぞれのデータの組み合わせの散布図をマトリックスで描いてくれます。痒いところに手が届く。

plot関数一発で複数の散布図が描ける
plot(d1)


ざっとみたところ、ウエストとヒップの相関が強そうです。場所が近いからね。

数値でも見てみましょう。

> cor(d1)
            年齢      身長         B           W          H
年齢  1.00000000 0.1496649 0.4605498 -0.02294316 0.01799732
身長  0.14966489 1.0000000 0.2171821  0.17997791 0.44760344
B     0.46054975 0.2171821 1.0000000  0.47559306 0.50405497
W    -0.02294316 0.1799779 0.4755931  1.00000000 0.67685625
H     0.01799732 0.4476034 0.5040550  0.67685625 1.00000000


予想通り、WとHの相関係数が約0.68で一番高いですね。

これだけ取り出してplotしてみましょう。

ウエストとヒップの散布図
plot(d1$W, d1$H)


見えるぞ、見えるぞ、右肩上がりの直線が。

○だと味気ないので、名前でプロットしてみました↓

ウエストとヒップの散布図(氏名でプロット)
plot(d1$W, d1$H, type="n")
rndm <- rnorm(nrow(d1))
text(d1$W, d1$H + rndm, rownames(d1), cex=0.7)


プロット位置が重ならないようにy軸方向にランダム要素を入れてみました。

もし、少しでも重なるのは嫌だという場合は、

Rの自己組織化マップ(SOM)で文字が重ならないようにする方法 - Rプログラミングの小ネタ


の方法を参考にしてみてください。

というわけで、練習でRを使うにしても、架空のデータより実際のデータの方がやる気が増したりするのではないでしょうか。




2014年9月4日木曜日

Rで比較対象がNAを含む場合の対処(データフレームの条件抽出あれこれ)

↓こんなデータフレームがあって、

> d
     name     club
1    taro baseball
2    jiro football
3  hanako     <NA>
4  ichiro baseball
5 yoshiko   tennis

clubがbaseballの行だけ抽出したい場合、↓こんなコードが最初に思い浮かびました。

d[d$club == "baseball", ]

でも実は、これではうまくいかなくて、↓こんな結果になってしまいます。

> d[d$club == "baseball", ]
     name     club
1    taro baseball
NA   <NA>     <NA>
4  ichiro baseball

"=="による比較の結果、TRUE/FALSEのベクトルが返ってきて、それでデータフレームを引くと、TRUEのところだけ残るという目論見なのですが、ベクトルにNAが混ざっており、その行はNAとして返ってきてしまうので都合が悪いです。

「NAではない」という条件を付け加えると、一応うまくいきます↓

> d[ !is.na(d$club) & d$club == "baseball", ]
    name     club
1   taro baseball
4 ichiro baseball

でもなんか、is.na関数とか使うのもなんだかなあと思っていたら、which関数が使えることに気づきました↓

> d[which(d$club=="baseball"), ]
    name     club
1   taro baseball
4 ichiro baseball

さらにいろいろ検索していると、いかにもそれ用な感じのsubset関数なるものがあることも知りました↓

> subset(d, d$club=="baseball")
    name     club
1   taro baseball
4 ichiro baseball

"=="や"which"でインデックス用のベクトルを作ってからデータフレームを引くというのは常套手段ではあるものの、慣れていない人にはピンと来にくい気もします。(絶対に覚えなきゃいけないやり方だとは思いますが・・・)

その点、subset関数はパッと見 何をやっているか分かりやすい印象がありますね。




2014年8月29日金曜日

Rで空のCSVファイルを読み込むとエラーが出てしまう場合の対処方法

例えば、複数のCSVファイルがあって、それをまとめて処理したいとしましょう。

"001.csv"は、

  taro,60
  hanako,100
  ichiro,90
  ...


"002.csv"は

  jiro,50
  keiko,70
  ...


みたいな感じでたくさんのCSVファイルがあって、それぞれを集計したいとか、全部をマージしたいとか、そんな場合。

以下のようなコードで処理しようとするんじゃないでしょうか、

files <- dir(pattern="\\.csv")
 
for(f in files){
  data <- read.csv(f, header=F)
 
  # このあとデータを集計するような処理
}


普通は上記でうまくいくんですが、たまに変なファイルが混じってたりする。

例えば、空のCSVファイルとか。

すると、

 以下にエラー read.table(file = file, header = header, sep = sep, quote = quote,  :
   入力中には利用可能な行がありません


というエラーメッセージが出て、処理が中断してしまいます。

人名と点数のデータみたいなものだと空のファイルなんてのは考えにくいですが、「複数の画像から顔検出して、検出位置の座標をCSVで出力」みたいな場合だと、不検出の場合は空のファイルができあがっている、なんてことはよくあります。

そういう場合の対処方法としては、try関数を使うと便利です。

for(f in files){
  # read.csvでエラーが出たかどうかはresに設定される
  res <- try ( data <- read.csv(f, header=F), silent=T )
 
  if( class(res) != "try-error" ) {
    # エラーが出なかったので、分析処理を続ける
    # ...
  } else {
    # 空のファイルがあったときの処理
    print("空のファイルがあったよ")
  }
}

try関数で囲んでおくと、エラーが出ても止まりません。エラーが出たかどうかは、resに設定されていますので、if文で処理を分けてやればOKですね。

silent=Fに変えると、発生したエラーメッセージを表示させる(処理は継続される)こともできます。




2014年8月28日木曜日

Rでjpeg画像を読み込み、加工する

もともと存在する画像(jpegファイルやpngファイル)を読み込んで、線や図形を書き足したり、文字を書き足したり、そんな感じの加工をRを使って行う方法です。

発想としては、plot関数で用意した作画領域に、rasterImage関数で画像を貼り付けしまえば、そのあとは使い慣れた(?)低水準作図関数でいろいろ落書きができるでしょ、という感じです。

Rでlenaさんの画像を読み込んで、rectで加工したもの
library(jpeg)
lena <- readJPEG("lena.jpg") # あの女性の画像
plot( 0,0, xlim=c(0,512), ylim=c(512,0), type="n", asp=1 )
rasterImage(lena, 0, 512, 512, 0)
rect( 220,200, 360,390, border="cyan" )

plot関数のylimオプションが上下逆になっているのがポイントです。jpeg上での座標の扱いは左上が原点になっているからなんですね。

rasterImage関数には、xleft, ybottom, xright, ytopの順で指定することになっています。つまり、貼り付けたい場所の左下(xleft, ybottom)と右上(xright, ytop)を指定しろということですね。対応する座標は、(0,512)と(512,0)になりますので、

rasterImage(lena, 0, 512, 512, 0)

で、正しく貼り付けられるわけです。

あとは低水準作図関数を駆使して、rect関数で矩形を描くなり、abline関数で線を引くなり、text関数で文字を書き込むなりすればいいんじゃないでしょうか。

で、お気づきだとは思いますが、x軸・y軸とかラベルとかが邪魔ですよね。これについては「余白なし」の設定になるように、plotの前にパラメタ指定しておけばOKです。軸やラベルは作画領域の外に追いやられ、見えなくなります。

余白なしにして、軸やラベルを表示させないようにしたもの

old.par <- par(no.readonly = TRUE) # 現状パラメタの退避
par( plt=c(0,1,0,1) ) # 余白なしの設定

plot( 0,0, xlim=c(0,512), ylim=c(512,0), type="n", asp=1 )
rasterImage(lena, 0, 512, 512, 0)
rect( 220,200, 360,390, border="cyan" )

par(old.par) # パラメタを元に戻す

パラメタの退避・復元は必須ではないんですが、やっておかないと、その後のplotでは軸やラベルが表示されなくなってしまいます。Rで他の作業をするときに困りますよね。(Rを起動しなおせば直りますが)

あとは、できあがった画像をファイルとして出力しましょう。

jpeg("out.jpg", width=512, height=512) # jpegファイルのオープン

old.par <- par(no.readonly = TRUE) # 現状パラメタの退避
par( plt=c(0,1,0,1) ) # 余白なしの設定

plot( 0,0, xlim=c(0,512), ylim=c(512,0), type="n", asp=1 )
rasterImage(lena, 0, 512, 512, 0)
rect( 220,200, 360,390, border="cyan" )

par(old.par) # パラメタを元に戻す

dev.off() # jpegファイルのクローズ

ファイルをオープンしているところのjpeg関数を、png関数に変更すればpngファイルに書き出すこともできます。

できあがったファイルを見てもらえば分かりますが、外周に白枠が追加されてしまっています。

これを消す方法は・・・分かりませんでした。誰か教えてください。




2014年7月29日火曜日

Rのifelse文が行う先行評価のせいで警告メッセージが出るときの対処法

放射能の測定結果の表し方にはお約束があって、検出下限値未満だったときは、検出下限値の数値の左側に不等号(小なり記号)をつけて表現します。

つまり

  8.6

と書いてあれば検出値が8.6だったということですが、

  <8.6

と書いてあれば、検出下限値が8.6という環境(測定器機の精度、試料の量、計測時間 etc.)で測定したが、値が小さいせいで検出されなかった、ということになります。

これを知らない人がやらかしてしまったのが、↓こちらの誤報だったりします。

新潟県:「女性セブン」誌の放射能関連記事に誤った内容が掲載されたため、訂正と謝罪を求める抗議文を送付しました。

「出典は厚生労働省のホームページであり、不等号(<)が付されている値が検出下限値であることを理解せず転記した」 by小学館

とのこと。記者一人が勘違いしたとしても、チェック機構はなかったんでしょうかね。不等号を見て、「これ何だろう?」って思わなかったんでしょうか。

で、話を戻しまして・・・

このような表記法となっておりますので、データはこんな感じになっていたりします。

> data
  食品 結果
1    A <7.8
2    B  8.9
3    C <6.7
4    D  9.0


下限値未満をどう扱うかは分析の方針次第ですが、仮に「下限値未満はゼロとして扱う」としたかったとしましょう。その場合、不等号付きの結果を加工してやる必要があります。

また、data$結果はcharacter型かfactor型(因子型)になっているので、このままでは平均を出したり、ヒストグラムを書いたりはできません。

「結果」をnumeric型に変換した「結果2」という列を追加してみましょう。

for ( i in 1:nrow(data) ) {
  if ( substr(data$結果[i], 1, 1) == "<" ) {
    # 不等号付きのときは、ゼロにする
    data$結果2[i] <- 0
  } else {
    # 不等号付きでないときは、そのまま数値にする
    data$結果2[i] <- as.numeric(data$結果[i])
  }
}

>data

  食品 結果 結果2
1    A <7.8   0.0
2    B  8.9   8.9
3    C <6.7   0.0
4    D  9.0   9.0


結果2はnumeric型になっていますので、mean(data$結果2) とか、hist(data$結果2)とかすることも可能です。

が、このようにfor文とif文で処理する人はあまりいないと思います。Rにはifelseという便利な記法があり、しかもifelseは「ベクトル化された演算」(詳しくは「アート・オブ・Rプログラミング」参照)であるため、for文で記述するよりも高速に実行できるというメリットもあります。

先ほどのfor文によるループをifelseで書き直すと以下のようになります。

data$結果2 <- ifelse( substr(data$結果, 1, 1) == "<", 0, as.numeric(data$結果) )

ifelseは第1引数を評価して、真ならば第2引数を、偽ならば第3引数を返すという動きをします。引数としてベクトルを渡せば、すべての要素について処理を繰り返してくれます。

実行結果として生成される、データフレームの列は、先ほどのfor文の場合と同じになります。

が、しかし、上記の書き方だと、↓こんなwarningが出てしまうんですよね。

 警告メッセージ:
In ifelse(substr(data$結果, 1, 1) == "<", 0, as.numeric(data$結果)) :
   強制変換により NA が生成されました


結果の中にNAなどないし、強制変換しなきゃいけないような処理もしていないはずだし・・・、なぜでしょうか?

どうやらこれは、ifelseが「第1引数の真偽にかかわらず、第2引数・第3引数を評価する」という動きをしているのことに起因しているようです。分岐はせずにすべて評価しておいて、不要な結果は捨てるという動きをしているのでしょう。

つまり「<7.8」に対しても、as.numericしてしまい結果がNAとなる、でも第1引数の条件式は偽になるので、その結果は使われない。なので、最終結果としては問題なしです。

でも、警告が出るのが気持ち悪い。なんとかならないものか・・・

で、結論としては、下記のように書けばOKということが分かりました。(前置きが長くなりました・・・)

data$結果2 <- as.numeric( ifelse( substr(data$結果, 1, 1) == "<", 0, data$結果 ) )

最後にas.numericすればいいだけですね。0に対してもas.numericすることになってしまいますが、これについては警告は出ないので、こちらの方がベターかなと思います。

【重要な補足】
「結果」の列がcharacter型ではなく、factor型になっている場合があります。いやむしろデフォルトでCSVファイルを読み込むとそうなるので、そちらのケースが多いかもしれません。

その場合は・・・

data$結果2 <- as.numeric( ifelse( substr(data$結果, 1, 1) == "<", 0, as.character(data$結果) ) )

としてやる必要があります。

ファクタをそのまま as.numericすると文字列の表す数字ではなく、水準の値が返ってきてしまうんですよね。

水準というのは、"男"、"女"、"その他" というファクタがあったときに、内部的にはそれぞれに 1、2、3が割り当たっていたりするんですが、その数字のことです。

ファクタの水準値を使われてしまうと、数値が明らかにおかしくなるので、変換後のデータフレームの内容をちゃんとチェックしてから、次の処理を行いましょう。




2014年7月23日水曜日

Rでshapeを出力するときはPDFファイルにすると便利

表題は正確に言うと、shapeファイルを読み込んで、加工やいろいろな可視化を施した後、画像(的なもの)として出力するときは、PDFにすると便利ですよ、という感じの趣旨です。

shapeファイルを読み込んで、plotする方法は↓以前に書きました。

e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ

ただ、shapeファイルってものによっては、ものすごく大量のポリゴンを含んでいますよね。地名などを重畳した場合、文字サイズが大きいと互いに重なってしまう、文字サイズが小さいと読めなくなってしまう、なんて問題が起こります。

前述のエントリでは、高解像度のpngに出力するという対処法を書いたんですけど、高解像度とはいえラスタである以上、拡大すれば文字はギザギザになるし、解像度を大きくするのにも限度がある。

使っている環境(メモリとか? CPUのビット数とか?)によるんでしょうけど、あまり大きくすると、

> png("output.png", width=10000, height=10000)
 以下にエラー png("output.png", width = 10000, height = 10000) :
   デバイス png() を開始できませんでした
 追加情報:  警告メッセージ:
1: In png("output.png", width = 10000, height = 10000) :
   ビットマップを割り当てられません
2: In png("output.png", width = 10000, height = 10000) :
  opening device failed


みたいな感じにエラーになってしまいます。

で、結論としてはPDFに出力するといい感じじゃん、ってことになりました。

コードは↓こんな感じ。

# GIS関係のパッケージ
library(maptools)

# 座標系の情報を持つオブジェクトを作る
# longlatは緯度経度指定、WGS84は世界測地系であることを表す

pj <- CRS("+proj=longlat +datum=WGS84")

# シェープファイルの読み込み
x  <- readShapePoly("h22ka14137.shp", proj4string=pj)

# PDFファイルをオープン
pdf("output.pdf", family="Japan1GothicBBB")

# プロット(出力先はPDFファイルとなる)
plot(x, lwd=0.1)
text(coordinates(x), as.character(x@data$MOJI), cex=0.1)

# ファイルを閉じる
dev.off()


ちなみにPDFファイルに出力するときは、フォントを指定(family="Japan1GothicBBB"の箇所)しないと、文字化けが起きるというワナがあります。

出力されたPDFファイルは↓こんな感じ。

川崎市麻生区の町丁目境界と町丁目名称(PDFファイル)

↓全体として見ると、「字ちっちゃっ!」って感じですが。

麻生区全体のPDFファイル


こまごまとしたところでも文字が重ならないようにするにはこれくらいのサイズにする必要がありました。text関数のcexオプションで調整できます。

線もデフォルトだと拡大したときに太く感じられたので、細めにしました。plot関数のlwdオプションで調整できます。

境界も地名も画像ではなく、ベクタの情報として埋め込まれているので、拡大しても粗くなりません↓

PDFファイルなので、拡大しても線や文字は滑らかなまま


「上麻生〇丁目」に混ざって、「上麻生」というのが点在するのは、たぶんこれらの地域に対して住居表示が実施されていないからですね。こんなにちっちゃく残すくらいなら、全部住居表示にしちゃえばよかったのに、って気もします。何か事情があったのかもしれませんが。■

[2014年9月10日追記] 
PDFファイルのサイズが大きくなりすぎた場合は、↓こんな方法もお試しください。
 SkyPDFを使ってPDFファイルのサイズを圧縮する: 主張