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時%M分") 　　[1] "2014-12-27 JST" %Y： 4桁の西暦 %m： 2桁の月 %d： 2桁の日 %H： 2桁の時 %M:  2桁の分 となっているので、それぞれの要素が文字列のどこに登場しているかを第2引数で指定しているわけですね。 「パーセントほにゃらら」とし

まずは今回のサンプル用として、身長っぽいデータを作ってみます。 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関数の戻り値のオブジェクトに各階級の度数が入っているんで

↓こちらで公開した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ファ

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 となる、というのを期待したのですが、その部分（lo

グラフに円を描くような関数が用意されていてもよさそうなものですが、どうもないみたいですので、自力で描く必要がありそうです。 簡単のため、半径１の円の例です。 まず思いつくのは、 　　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 ( -

例えば↓こんな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

以前のコマンド入出力を確認したいときは、スクロールバーでコンソール画面をスクロールさせればいいですよね。 でも、要素数が非常多いリストやベクトルを表示させたりしたあとだと、スクロールバーのつかむところがちっちゃくなっていて、行き過ぎたり戻り過ぎたりと、なかなか該当箇所が見つからなかったりしますよね。 で、いったんRのコンソール画面をリセットしたいってときは、 　　Ctrl + L を入力すれば、コンソール画面がクリアされます。 という分かり切った操作ですが、あえてスクリーンショットをつけてみました。 コンソールがごちゃごちゃしてきたイメージ 「Ctrl + L」で、↑が、↓になります すっきりとクリア ちなみにショートカットキーでなくても、メニューの「編集」→「コンソール画面を消去」としても同じことができます。「コンソール画面を消去」だと、コンソールのウィンドウ自体が閉じられてしまいそうですが、そんなことはないので大丈夫です。 （以下は蛇足） 私は結構長い間このコマンドを知らなくて、過去の出力結果が邪魔な時は、たくさんリターンを押して、表示結果に「切れ目」を入れて、該当箇所を見つけやすくするという、面倒なことをしていました。 で、放送大学の「データからの知識発見」のテレビ放送を見ているときに、講師の秋光先生が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軸に追記したい場合はsi

今回は、普段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 ) ) 多次元尺度構成法を使ってみました。データを（年齢、身長、バスト、ウエスト、ヒップ）の５次元だととらえて距離を算出し、なるべくその距離感が再現されるように、２次元に落として（５次元のままじゃ表現しづらいので）プロットされた図という感じでしょうか。 結果はあまりあてにできないと

↓こんなデータフレームがあって、 > 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

例えば、複数の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 { # 空のファイ

もともと存在する画像（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) # 現状パラメタの退避

放射能の測定結果の表し方にはお約束があって、検出下限値未満だったときは、検出下限値の数値の左側に不等号（小なり記号）をつけて表現します。 つまり 　　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])   } } &