Rプログラミングの小ネタ

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

表題は正確に言うと、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 +datu...

例えば、セルの1つに「2014年7月22日」と入力してあったとしましょう。 これをXLConnectパッケージを使ってRで読み込むと、 > library(XLConnect) > wb <- loadWorkbook("input.xlsx") > ws <- readWorksheet(wb, sheet=1, header=F) > ws    Col1 1 41842 のように、「2014年7月22日」と入力してあったセルの値が、「41842」となっています。 これはシリアル値といって、1900年1月1日を「1」、1月2日を「2」、・・・として日付を数値化したもので、エクセルは内部的にはこのように日付データを保持しています。というのを知っている人は多いと思います。 じゃあ、この値を日付としてRから使いたいときにどうやるのか？ というのが今回のエントリの目的です。 Rには日付のシリアル表現のための関数があって、整数値と起源日付を指定することによって、シリアル値を日付（Date）型に変換することができます。 ということは、↓これでいけるはずだが・・・ > as.Date(41842, origin="1900-01-01") # 2014年7月 22 日のつもり [1] "2014-07- 24 " あれ？ 2日ほどずれている。 このずれの 原因の1つ目 は以下のコードを実行することで、すぐに理解できます。 > as.Date( 0 , origin="1900-01-01") [1] "1900-01-01" > as.Date( 1 , origin="1900-01-01") [1] "1900-01-02" > as.Date( 2 , origin="1900-01-01") [1] "1900-01-03" エクセルでは1900年1月1日を「 1 」としているのに対し、Rのas.Date関数はoriginオプションで指定された日付（つ...

Rでsom関数を使って自己組織化マップを作った後、各ユニットの座標に文字列をplotした場合、同じユニットに配置された文字列の重なりが気になることってありますよね。 下記は、放送大学の 「データからの知識発見」という授業のテキスト からの引用です。 重なり対処を全く行っていない表示例 入門者向けの講義なので、サンプルコードのシンプルさを優先したのか、文字の重なりについての対処はありません。用いられていた入力のデータ数（動物の数）が14、同じ場所に配置されてしまったデータが最大で2個なので、なんとか読み取れます。ライオンとトラが重なると「ライトオン」っぽくなるというのはトリビアですね。 ただ、データ数が多くなって、3つも4つも重なってくると、読み取れなくなって実害が出てくると思います。 「 データマイニング入門 」（わたせせいぞうの表紙のシリーズのやつ）では、そこに一工夫されていて、$visualで得られる座標をそのまま使うのではなく、乱数を加えることによって、文字が全く同じ位置で重なることを避けています↓ 乱数による重なり対処を行った表示例 でも、そこはやはり乱数なので、うまく分離されているところもあれば、これじゃ読み取れないでしょってレベルのところもあります。 ↓以下は データマイニング入門 に載っていたサンプルをそのまま動かしたものです。 乱数による重なり対処を行った表示例（私の環境で実行させたもの） library ( som ) 動物データ 1 <- read.csv ( "animal1.csv" , header=T ) 標準化動物データ <- normalize ( 動物データ 1 [ , 2 : 14 ] , byrow=F ) 動物 SOM <- som ( 標準化動物データ , xdim= 10 , ydim= 10 , topol= "rect" ) 乱数 <- cbind ( rnorm ( nrow ( 動物データ 1 ) , 0 , 0.15 ) , rnorm ( nrow ( 動物データ 1 ) , 0 , 0.15 ) ) 動物マップ <- ...

ジオコーディングというのは、 　入力： 住所の文字列 　出力： 緯度経度の数値 という処理のことです。YahooやGoogleがAPIを公開しています。 YahooのAPIの仕様は↓下記のページにあります。 YOLP(地図):Yahoo!ジオコーダAPI - Yahoo!デベロッパーネットワーク このページの説明にしたがって、住所文字列をエンコードしてURLに埋め込んで、ブラウザのアドレス欄に入れれば結果が返ってきます。 でも、こういうのはプログラムから使いたいですよね。で、私が多少なりとも使いなれているRでのやり方を書きとめておきます。 普通に使うときは、一連の処理を関数の形にして、それを呼び出すような書き方をする方がすっきりするのですが、今回は使い方の紹介がメインなので、あえてベタ書きしています。 あと、アプリケーションIDは各人で取得する必要があります。YahooのIDを持っていれば（持っていない場合は作れば）、無料で取得できます。 library ( RCurl ) library ( RJSONIO )   # appidのところは自分で取得したものを入れてください root <- "http://geo.search.olp.yahooapis.jp/OpenLocalPlatform/V1/geoCoder?appid=xxxxxxxxxx&output=json" addr <- "東京都墨田区押上1-1-13" url <- paste ( root , "&query=" , addr , sep= "" ) url.u <- URLencode ( iconv ( url , "" , "UTF-8" ) ) # UTF-8に変換   # プロキシへの対処（使っていない場合はgetURL(url.u)とそのまま呼んでOK） curl <- getCurlHandle ( ) curlSetOpt ( .opts= list ( proxy = "ZZZZZ:8080" ) , curl=...

まずはShapeファイルを入手するところから。 総務省統計局がe-Statというサイトを運営していて、そこに全国各地域のShapeファイルが公開されています。 地図で見る統計（統計ＧＩＳ) (1) 「平成２２年国勢調査（小地域） 2010/10/01」を選択 (2) 「男女別人口総数及び世帯総数」にチェック (3) 「統計表各種データダウンロードへ」ボタンをクリック (1) 「都道府県」を選択 (2) 「市区町村」を選択 (3) 「検索」ボタンをクリック (4) 境界データの欄に表示されている     「世界測地系緯度経度・Shape形式」のデータをクリック ダウンロードしたzipファイルを展開すると、 　xxx.shp 　xxx.dbf 　xxx.shx 　xxx.prj という、拡張子だけが違うファイルが4つできるんですが、これがshapeデータの1セットになります。 それぞれのファイルの中身は、 シェープファイルとは？ | 用語集とGISの使い方 | 株式会社パスコ によると、 　shp：    図形の座標が保存 　dbf：    属性の情報が保存 　shx：    shpの図形とdbfの属性の対応関係が保存 だそうです。あと、「.prj」は測地系の情報が入っています。（テキストとして開くとそれっぽい文言が入っている） では単純にshapeファイルをRで表示させてみましょう。 境界情報だけをプロット # GIS関係のパッケージ library ( maptools )   # 座標系の情報を持つオブジェクトを作る # longlatは緯度経度指定、WGS84は世界測地系であることを表す pj <- CRS ( "+proj=longlat +datum=WGS84" )   # シェープファイルの読み込み x <- readShapePoly ( "h22ka14137.shp" , proj4string=pj )   # 表示 plot ( x ) あら簡単。実質4行で表示できまし...