rworldmap:Mapping Worldwide Suicide Rates

R Journal Volume 3/1, June 2011 で位置データを可視化する “rworldmap” という R パッケージが紹介されていた。

Andy South : rworldmap: A New R package for Mapping Global Data. The R Journal Vol. 3/1, June 2011. ISSN 2073-4859

rworldmap は R プログラムが本職でない人でもお手軽に使えるよう設計されており、同解説記事では rworldmap のビジョンが次のように宣言されている。

The vision for rworldmap is to produce a package to facilitate the visualisation and mapping of global data. Because the focus is on global data, the package can be more specialised than existing packages, making world mapping easier, partly because it doesn’t have to deal with detailed local maps.

R Journal の解説記事では “The Happy Planet Index(HPI)” というポジティブなインデックスが例で取り上げられている。
一方で数か月前に “Dark Contrasts: The Paradox of High Rates of Suicide in Happy Places” (Mary C. Dalya, , Andrew J. Oswaldb, Daniel Wilsona, and Stephen Wuc. Journal of Economic Behavior & Organization. 2011)というネガティブな報告が発表され、一部で話題になっていた。(The New York Times “Happiest Places Post Highest Suicide Rates”  by TARA PARKER-POPE. 2011/04/22)
対比のために、後者の論文で引用されている WHO による国別自殺率(Suicide rates per 100,000 by country, year and sex)を題材に rworldmap を使ってみた。


まずはパッケージのインストール

> install.packages('rworldmap')
> library(rworldmap)
Loading required package: sp
Loading required package: maptools
Loading required package: foreign
Loading required package: lattice

次にWHOサイトから取得したデータを読み込む

> dat <-read.csv('suicide.csv', header=T, as.is=T)
> head(dat)
              Country Year Males Females
1             ALBANIA   03   4.7     3.3
2 ANTIGUA AND BARBUDA   95   0.0     0.0
3           ARGENTINA   05  12.7     3.4
4             ARMENIA   06   3.9     1.0
5           AUSTRALIA   04  16.7     4.4
6             AUSTRIA   07  23.8     7.4

この Males/Females の rate を rworldmap で Country 別にマップするのが最終ゴール

次に読み込んだデータでいきなりマップしてみる
データフレームの国名(dat$Country)を rworldmap の国名と JOIN させる。(joinCode = ‘NAME’)
国の指定は、NAME 以外にも ISO3(ISO 3166-1 alpha-3 で定義されているアルファベット3文字コード。 e.g. JPN for Japan)、UN(ISO 3166-1 で定義されている数値コード。e.g. 392 for Japan)が指定可能

sPDF <- joinCountryData2Map(dat,
          joinCode='NAME',
          nameJoinColumn='Country',
          verbose='TRUE')
0 codes from your data successfully matched countries in the map
105 codes from your data failed to match with a country code in the map
       failedCodes                        failedCountries
  [1,] "ALBANIA"                          "ALBANIA"
  [2,] "ANTIGUA AND BARBUDA"              "ANTIGUA AND BARBUDA"

246 codes from the map weren't represented in your data
       codesMissingFromUserData
  [1,] "Antigua and Barbuda"
  [2,] "Algeria"
  [3,] "Azerbaijan"
  [4,] "Albania"
  [5,] "Armenia"

       countriesMissingFromUserData
  [1,] "Antigua and Barbuda"
  [2,] "Algeria"
  [3,] "Azerbaijan"
  [4,] "Albania"
  [5,] "Armenia"
  [6,] "Angola"

国名の修正1
ログには「Country カラムで指定された国名にマッチする国名が見つかりません」というエラーが大量に出ている。期待される国名をよく見ると、1文字目が大文字で、2文字目以降が小文字でないといけないのに、WHOデータはすべて大文字。
仕方なしに、上記仕様に国名変換する。
Python の title がぴったりなのだけど、あいにく R にはそのような関数がない。
グーグルで検索すると ggplot2 に firstUpper という近い関数があると分かったので、これを流用して mytitle 関数を作る。

> firstUpper # ggplot2 function
function (s)
{
    paste(toupper(substring(s, 1, 1)), substring(s, 2), sep = "")
}
<environment: namespace:ggplot2>
mytitle <- function(s) #
{
  paste(toupper(substring(s, 1, 1)), tolower(substring(s, 2)), sep = '')
}
> mytitle('JAPAN')
[1] "Japan"

この関数を使って国名を正規化して再トライ

> dat$Con <- title(dat$Country)
> head(dat)
              Country Year Males Females                 Con
1             ALBANIA   03   4.7     3.3             Albania
2 ANTIGUA AND BARBUDA   95   0.0     0.0 Antigua and barbuda
3           ARGENTINA   05  12.7     3.4           Argentina
4             ARMENIA   06   3.9     1.0             Armenia
5           AUSTRALIA   04  16.7     4.4           Australia
6             AUSTRIA   07  23.8     7.4             Austria
sPDF <- joinCountryData2Map(dat,
            joinCode='NAME',
            nameJoinColumn='Con',
            verbose='TRUE')
82 codes from your data successfully matched countries in the map
23 codes from your data failed to match with a country code in the map
      failedCodes                        failedCountries
 [1,] "Antigua and barbuda"              "ANTIGUA AND BARBUDA"
 [2,] "Bosnia and herzegovina"           "BOSNIA AND HERZEGOVINA"
 [3,] "China (hong kong sar)"            "CHINA (Hong Kong SAR)"
 [4,] "Costa rica"                       "COSTA RICA"
 [5,] "Czech republic"                   "CZECH REPUBLIC"
 [6,] "Dominican republic"               "DOMINICAN REPUBLIC"

164 codes from the map weren't represented in your data
       codesMissingFromUserData
  [1,] "Antigua and Barbuda"
  [2,] "Algeria"
  [3,] "Angola"

       countriesMissingFromUserData
  [1,] "Antigua and Barbuda"
  [2,] "Algeria"
  [3,] "Angola"
  [4,] "American Samoa"
  [5,] "Bermuda"
  [6,] "Bangladesh"

国名の修正2

ルックアップに失敗した国名を確認すると、スペースを含んだ国名で、スペースの次の文字が大文字ではなく小文字になっているのが原因とわかる。また、”and”, “or'” は例外的にすべて小文字としなければならない。

上記仕様変更を反映したのが次の mytitle2

mytitle2 <- function(x) {
  tmp <- sapply(str_split(x, ' '),
                 function(val) ifelse (tolower(val) %in% c('and', 'of'), tolower(val), mytitle(val)))
  str_c(tmp, collapse=' ')
}
mytitle2('UNITED STATES OF AMERICA')
   United    States        of   America
 "United"  "States"      "of" "America"
dat$Con <- sapply(dat$Country, mytitle2)

sPDF <- joinCountryData2Map(dat,
          joinCode='NAME',
          nameJoinColumn='Con',
          verbose='TRUE')
96 codes from your data successfully matched countries in the map
8 codes from your data failed to match with a country code in the map
     failedCodes
[1,] "China (selected Rural & Urban Areas)"
[2,] "China (hong Kong Sar)"
[3,] "Iran"
[4,] "Republic of Korea"
[5,] "Russian Federation"
[6,] "Saint Vincent and The Grenadines"
[7,] "Tfyr Macedonia"
[8,] "United States of America"
     failedCountries
[1,] "China (selected rural & urban areas)"
[2,] "China (hong kong sar)"
[3,] "Iran"
[4,] "Republic of korea"
[5,] "Russian federation"
[6,] "Saint vincent and the grenadines"
[7,] "Tfyr macedonia"
[8,] "United states of america"
150 codes from the map weren't represen

国名の修正3

この段階でも、いくつかの国名がこけている。プログラムで対応することはあきらめ、getMap()[[‘NAME’]] で rworldmap での国名を確認し、元のCSV ファイルのデータをクレンジング。楽しくないけど仕方ない。

データ修正後’、再度ファイルを読み込む

dat <-read.csv('suicide.csv', header=T, as.is=T)
sPDF <- joinCountryData2Map(dat,
          joinCode='NAME',
          nameJoinColumn='Con',
verbose='TRUE')
mapDevice() #create world map shaped window
mapCountryData(sPDF, nameColumnToPlot='Males')
savePlot('male', type='png')

mapCountryData(sPDF, nameColumnToPlot='Females')
savePlot('female', type='png')

ようやくワールドマップの完成

こうしてようやく完成した男女別 Suicide rates per 100,000 の国別マップ画像が以下。


件の論文は西ヨーロッパ、および、アメリカを対象に happiness と suicide rate の相関を見ていて、その他のエリア(日本とか Japan とか Japon とか Giappone)はモデルあてはめの対象外。

Tagged with: , , , ,
Posted in life, R

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Archives
  • RT @__apf__: How to write a research paper: a guide for software engineers & practitioners. docs.google.com/presentation/d… /cc @inwyrd 6 days ago
  • RT @HayatoChiba: 昔、自然と対話しながら数学に打ち込んだら何かを悟れるのではと思いたち、専門書1つだけ持ってパワースポットで名高い奈良の山奥に1週間籠ったことがある。しかし泊まった民宿にドカベンが全巻揃っていたため、水島新司と対話しただけで1週間過ぎた。 それ… 3 weeks ago
  • RT @googlecloud: Ever wonder what underwater fiber optic internet cables look like? Look no further than this deep dive w/ @NatAndLo: https… 3 weeks ago
  • @ijin UTC+01:00 な時間帯で生活しています、、、 6 months ago
  • RT @mattcutts: Google's world-class Site Reliability Engineering team wrote a new book: amazon.com/Site-Reliabili… It's about managing produc… 9 months ago
%d bloggers like this: