Rを使ったデータ分析（実践編）１
データのダウンロード
厚生労働省のホームページの「(3)過去の食中毒事件一覧」から以下の３つのデータをダウンロードします。
- 令和元年（2019年）食中毒発生事例
- 平成30年（2018年）食中毒発生事例
- 平成29年（2017年）食中毒発生事例
ダウンロードする際、それぞれ「FP2019」、「FP2018」、「FP2017」という名前を付けて「CSV(コンマ区切り)(*.csv)」形式で保存します。
RStudioを立ち上げます。
まず必要となるパッケージを読み込む。
library(tidyverse) library(dplyr)
ggplotとjpndistrictで地図を書く①
まず必要となるパッケージをR studioに読み込みます。
library(tidyverse) #ggplot2などでグラフを作成する、データをきれいにするために使用 library(sf) #sf(simple features)のデータを扱うのに使用 library(jpndistrict) #元となるデータ
先ほど読み込んだパッケージ"jpndistrict"内に"jpn_cities"のデータセットが含まれています。これは区市町村レベルのデータセットで、その13番目が東京都です。
そこで13番目のデータを抽出します。ただし、島しょ地域が含まれていると、地図が見にくくなるため、今回は島しょ地域を除きます（54番目以降が島しょ地域）。
pref13 <- jpn_cities(13)[1:53, ] ggplot(pref13) + geom_sf()
次にmutate関数を使って、先ほど作成した"pref13"に新たに"HealthCenter"という名前の変数（列）を作成し、そこに管轄となる保健所名を入れていきます。作成したものを新たに"tokyo"とします。
tokyo <- pref13 %>% mutate(HealthCenter = case_when( city == "千代田区" ~ "千代田区", city == "中央区" ~ "中央区", city == "港区" ~ "港区", city == "新宿区" ~ "新宿区", city == "文京区" ~ "文京区", city == "台東区" ~ "台東区", city == "墨田区" ~ "墨田区", city == "江東区" ~ "江東区", city == "品川区" ~ "品川区", city == "目黒区" ~ "目黒区", city == "大田区" ~ "大田区", city == "世田谷区" ~ "世田谷区", city == "渋谷区" ~ "渋谷区", city == "中野区" ~ "中野区", city == "杉並区" ~ "杉並区", city == "豊島区" ~ "豊島区", city == "北区" ~ "北区", city == "荒川区" ~ "荒川区", city == "板橋区" ~ "板橋区", city == "練馬区" ~ "練馬区", city == "足立区" ~ "足立区", city == "葛飾区" ~ "葛飾区", city == "江戸川区" ~ "江戸川区", city == "八王子市" ~ "八王子市", city == "町田市" ~ "町田市", city == "立川市" ~ "多摩立川", city == "武蔵野市" ~ "多摩府中", city == "三鷹市" ~ "多摩府中", city == "青梅市" ~ "西多摩", city == "府中市" ~ "多摩府中", city == "昭島市" ~ "多摩立川", city == "調布市" ~ "多摩府中", city == "小金井市" ~ "多摩府中", city == "小平市" ~ "多摩小平", city == "日野市" ~ "南多摩", city == "東村山市" ~ "多摩小平", city == "国分寺市" ~ "多摩立川", city == "国立市" ~ "多摩立川", city == "福生市" ~ "西多摩", city == "狛江市" ~ "多摩府中", city == "東大和市" ~ "多摩立川", city == "清瀬市" ~ "多摩小平", city == "東久留米市" ~ "多摩小平", city == "武蔵村山市" ~ "多摩立川", city == "多摩市" ~ "南多摩", city == "稲城市" ~ "南多摩保健所", city == "羽村市" ~ "西多摩", city == "あきる野市" ~ "西多摩", city == "西東京市" ~ "多摩小平", city == "西多摩郡 瑞穂町" ~ "西多摩", city == "西多摩郡 日の出町" ~ "西多摩", city == "西多摩郡 檜原村" ~ "西多摩", city == "西多摩郡 奥多摩町" ~ "西多摩")) head(tokyo) city_code city geometry HealthCenter <chr> <chr> <S3: sfc_GEOMETRY> <chr> 13101 千代田区 <S3: sfc_GEOMETRY> 千代田区 13102 中央区 <S3: sfc_GEOMETRY> 中央区 13103 港区 <S3: sfc_GEOMETRY> 港区 13104 新宿区 <S3: sfc_GEOMETRY> 新宿区 13105 文京区 <S3: sfc_GEOMETRY> 文京区 13106 台東区 <S3: sfc_GEOMETRY> 台東区 6 rows
2.1 保健所ごとに色分けを行う
ggplot(tokyo) + geom_sf(aes(fill = HealthCenter), show.legend = FALSE)
2.2 保健所の凡例を追加します。
ggplot(tokyo) + geom_sf(aes(fill = HealthCenter), show.legend = FALSE) + geom_sf_label(aes(label = HealthCenter), size = 3, label.padding = unit(0.7, "mm"))
保健所数が多すぎてみずらいです。
3.1 新たに点をプロットしたいときのデータの作成方法
myhome <- tibble::tribble( ~city, ~lat, ~lon, "My Home", 139.65, 35.62, )
3.2 島しょ地域を除いた役場を抽出
df_office13 <- jpn_admins(13) %>% slice(1:172, 199)
3.3 点をプロットするのに２パターンある。一つ目は直接xyを入力する方法.２つめは事前に作っておいたデータを使う方法
ggplot() + geom_sf(data = tokyo, mapping = aes(fill = city_code), show.legend = FALSE) + geom_sf(data = df_office13, show.legend = FALSE, size = 0.1, alpha = 0.5) + geom_point(mapping = aes(x = 139.691, y = 35.689), color = "red") + geom_point(data = myhome, mapping = aes(x = lat, y = lon), color = "white", shape = 2)
maptoolsで日本の地図を描く
readShapePoly関数を使うため、maptoolsをlibraryで読み込みます。
library(maptools)
3種類の日本地図を描くため、GADM(https://gadm.org/index.html)から日本のデータ(shapefile)をダウンロード・展開し、使用しているdirectoryにデータを保存しておきます。
fortify (ggplot2の関数)でreadShapePolyで読み込んだshapefileを一般的なデータフレームに変換します。
Japan_shp0 <- readShapePoly("gadm36_JPN_shp/gadm36_JPN_0.shp") Japan_map0 <- fortify(Japan_shp0) ggplot(Japan_map0, aes(x = long, y =lat, group = group)) + geom_path()
参考までにJapan_map0の内容を見てみます。
> str(Japan_map0) 'data.frame': 562995 obs. of 7 variables: $ long : num 138 138 138 138 138 ... $ lat : num 34.6 34.6 34.6 34.6 34.6 ... $ order: int 1 2 3 4 5 6 7 8 9 10 ... $ hole : logi FALSE FALSE FALSE FALSE FALSE FALSE ... $ piece: Factor w/ 1398 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ... $ id : chr "0" "0" "0" "0" ... $ group: Factor w/ 1398 levels "0.1","0.2","0.3",..: 1 1 1 1 1 1 1 1 1 1 ... >
Japan_shp1 <- readShapePoly("gadm36_JPN_shp/gadm36_JPN_1.shp") Japan_map1 <- fortify(Japan_shp1) ggplot(Japan_map1, aes(x = long, y =lat, group = group)) + geom_path()
Japan_shp2 <- readShapePoly("gadm36_JPN_shp/gadm36_JPN_2.shp") Japan_map2 <- fortify(Japan_shp2) ggplot(Japan_map2, aes(x = long, y =lat, group = group)) + geom_path()
Rで地図を描いてみよう② 色分けする
USArrestsのデータを使って、アメリカの州ごとに殺人数で色分けします。
library(maps) library(plyr)
USArrestsの先頭列にstateという新しい列を作り、そこに小文字にした州名を入れ、crimesとして保存します。同様にstateのデータをstates_mapに保存します。その後、states_mapとcrimesを結合します。
crimes <- data.frame(state = tolower(rownames(USArrests)), USArrests) states_map <- map_data("state") crime_map <- merge(states_map, crimes, by.x = "region", by.y = "state") head(crime_map) region long lat group order subregion Murder Assault UrbanPop Rape 1 alabama -87.46201 30.38968 1 1 <NA> 13.2 236 58 21.2 2 alabama -87.48493 30.37249 1 2 <NA> 13.2 236 58 21.2 3 alabama -87.52503 30.37249 1 3 <NA> 13.2 236 58 21.2 4 alabama -87.53076 30.33239 1 4 <NA> 13.2 236 58 21.2 5 alabama -87.57087 30.32665 1 5 <NA> 13.2 236 58 21.2 6 alabama -87.58806 30.32665 1 6 <NA> 13.2 236 58 21.2
crime_mapをgroupとorderで昇順に並べ替えます。
crime_map <- arrange(crime_map, group, order)
Murderは人口10万人あたりの殺人数です。それを基にアメリカの州ごとに色付けします。
ggplot(crime_map, aes(x = long, y = lat, group = group, fill = Murder)) + geom_polygon() + coord_map("polyconic")
（参考：Rグラフィックス クックブック、Winston Chang著）
Rで地図を描いてみよう
今回はRを使って地図を描いてみます。
# ggplot2、maps、mapdataをRに読み込み library(ggplot2) library(maps) library(mapdata)
次にアメリカ合衆国のデータをstate.mapに保存します。map_dataはggplot2に含まれている関数の一つで、データフレームをggplot2でプロットしやすいフォーマットに変換します。stateはmapsに含まれているデータでアメリカの地図を描く基本情報を含んでいます。
state.map <- map_data("state")
ggplotを使ってstate.mapのデータからアメリカの地図を描きます。
# 投影法：直交座標、色付き ggplot(state.map, aes(x = long, y = lat, group = group)) + geom_polygon(fill = "green", color = "black")
# 投影法：直交座標、色なし ggplot(state.map, aes(x = long, y = lat, group = group)) + geom_path()
# 投影法：メルカトル図法、色なし ggplot(state.map, aes(x = long, y = lat, group = group)) + geom_path() + coord_map("mercator")
次に日本地図を描いてみます。
# 日本の地図 Japan <- map_data("world", region = "Japan") ggplot(Japan, aes(x = long, y = lat, group = group)) + geom_path()
# 日本の地図、色付き japan <- map_data("japan") ggplot() + borders("japan", fill = "lightgreen", size = 0.3) + theme_bw() + labs( y = "経度", x = "緯度")
データの部分集合化（subset関数）
今回はRにデフォルトで内蔵されているmtcarsのデータを使って、データの部分集合化を学びます。mtcarsは1973~1974年の車のモデルについてそのデザインとガソリンの消費量に対する性能を評価したデータです。まずmtcarsのデータがどのようなものか見てみましょう。結果から32のデータ（行）と11の変数（列）を持つデータフレームであることが分かります
> class(mtcars) [1] "data.frame" > dim(mtcars) [1] 32 11 > summary(mtcars) mpg cyl disp hp drat wt qsec Min. :10.40 Min. :4.000 Min. : 71.1 Min. : 52.0 Min. :2.760 Min. :1.513 Min. :14.50 1st Qu.:15.43 1st Qu.:4.000 1st Qu.:120.8 1st Qu.: 96.5 1st Qu.:3.080 1st Qu.:2.581 1st Qu.:16.89 Median :19.20 Median :6.000 Median :196.3 Median :123.0 Median :3.695 Median :3.325 Median :17.71 Mean :20.09 Mean :6.188 Mean :230.7 Mean :146.7 Mean :3.597 Mean :3.217 Mean :17.85 3rd Qu.:22.80 3rd Qu.:8.000 3rd Qu.:326.0 3rd Qu.:180.0 3rd Qu.:3.920 3rd Qu.:3.610 3rd Qu.:18.90 Max. :33.90 Max. :8.000 Max. :472.0 Max. :335.0 Max. :4.930 Max. :5.424 Max. :22.90 vs am gear carb Min. :0.0000 Min. :0.0000 Min. :3.000 Min. :1.000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:3.000 1st Qu.:2.000 Median :0.0000 Median :0.0000 Median :4.000 Median :2.000 Mean :0.4375 Mean :0.4062 Mean :3.688 Mean :2.812 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:4.000 3rd Qu.:4.000 Max. :1.0000 Max. :1.0000 Max. :5.000 Max. :8.000
①$と[ ]を使って部分集合化する。
以下の例では[ ]を使ってcylが６のデータについて、mpgの平均を計算しています。cylはシリンダーの数、mpgはマイル/ガロン(≒3.8リットル)です。
> mean(mtcars$mpg[mtcars$cyl == 6]) [1] 19.74286
また下のようにgear（ギアの数）が５のものだけを選んで部分集合化し、gear5として保存することもできます。忘れてはいけない点として最後にカンマ（,）が必要になります。
> gear5 <- mtcars[mtcars$gear == 5, ] > gear5 mpg cyl disp hp drat wt qsec vs am gear carb Porsche 914-2 26.0 4 120.3 91 4.43 2.140 16.7 0 1 5 2 Lotus Europa 30.4 4 95.1 113 3.77 1.513 16.9 1 1 5 2 Ford Pantera L 15.8 8 351.0 264 4.22 3.170 14.5 0 1 5 4 Ferrari Dino 19.7 6 145.0 175 3.62 2.770 15.5 0 1 5 6 Maserati Bora 15.0 8 301.0 335 3.54 3.570 14.6 0 1 5 8
②subset関数を使う
下の例ではsubset関数を使って、cylが４かつam（０：オートマ、１：マニュアル）が１のものだけを部分集合化しています
> cyl4 <- subset(mtcars, subset = cyl == 4 & am == 1) > cyl4 mpg cyl disp hp drat wt qsec vs am gear carb Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1 Fiat 128 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4 1 Honda Civic 30.4 4 75.7 52 4.93 1.615 18.52 1 1 4 2 Toyota Corolla 33.9 4 71.1 65 4.22 1.835 19.90 1 1 4 1 Fiat X1-9 27.3 4 79.0 66 4.08 1.935 18.90 1 1 4 1 Porsche 914-2 26.0 4 120.3 91 4.43 2.140 16.70 0 1 5 2 Lotus Europa 30.4 4 95.1 113 3.77 1.513 16.90 1 1 5 2 Volvo 142E 21.4 4 121.0 109 4.11 2.780 18.60 1 1 4 2
また、下の例は上記と同じ部分集合化ですが、列としてmpg、cyl、wt、amだけを選択しています。
> cyl4s <- subset(mtcars, select = c(mpg, cyl, wt, am), subset = cyl == 4 & am == 1) > cyl4s mpg cyl wt am Datsun 710 22.8 4 2.320 1 Fiat 128 32.4 4 2.200 1 Honda Civic 30.4 4 1.615 1 Toyota Corolla 33.9 4 1.835 1 Fiat X1-9 27.3 4 1.935 1 Porsche 914-2 26.0 4 2.140 1 Lotus Europa 30.4 4 1.513 1 Volvo 142E 21.4 4 2.780 1
上記と同じことを＄と[ ]を使って行うことができます。
> syl4r <- mtcars[mtcars$cyl == 4 & mtcars$am == 1, c("mpg", "cyl", "wt", "am")] > syl4r mpg cyl wt am Datsun 710 22.8 4 2.320 1 Fiat 128 32.4 4 2.200 1 Honda Civic 30.4 4 1.615 1 Toyota Corolla 33.9 4 1.835 1 Fiat X1-9 27.3 4 1.935 1 Porsche 914-2 26.0 4 2.140 1 Lotus Europa 30.4 4 1.513 1 Volvo 142E 21.4 4 2.780 1
③部分一致に対する条件文を作る
以下の例は、mpgが20以上かつgearが5又は4の条件に一致するものには１、一致しないものには０を割り当てる条件文です。mtcarsに新たにspcという変数（列）を追加しています。
> mtcars$spc <- ifelse(mtcars$mpg >= 20 & mtcars$gear == 5 | mtcars$gear == 4, 1, 0)
vapplyとtapplyの使い方について
この記事はRstudioのSwirlの「11: vapplyとtapply」を参照しています。
まずFlagsというデータセットをRに読み込みます（Swirlを使用している場合、自動的に読み込まれます。）。
参考までにflagsのsummaryの結果を示します。
> summary(flags) name landmass zone area Afghanistan : 1 Min. :1.000 Min. :1.000 Min. : 0.0 Albania : 1 1st Qu.:3.000 1st Qu.:1.000 1st Qu.: 9.0 Algeria : 1 Median :4.000 Median :2.000 Median : 111.0 American-Samoa: 1 Mean :3.572 Mean :2.211 Mean : 700.0 Andorra : 1 3rd Qu.:5.000 3rd Qu.:4.000 3rd Qu.: 471.2 Angola : 1 Max. :6.000 Max. :4.000 Max. :22402.0 (Other) :188 population language religion bars Min. : 0.00 Min. : 1.00 Min. :0.000 Min. :0.0000 1st Qu.: 0.00 1st Qu.: 2.00 1st Qu.:1.000 1st Qu.:0.0000 Median : 4.00 Median : 6.00 Median :1.000 Median :0.0000 Mean : 23.27 Mean : 5.34 Mean :2.191 Mean :0.4536 3rd Qu.: 14.00 3rd Qu.: 9.00 3rd Qu.:4.000 3rd Qu.:0.0000 Max. :1008.00 Max. :10.00 Max. :7.000 Max. :5.0000 stripes colours red green Min. : 0.000 Min. :1.000 Min. :0.0000 Min. :0.0000 1st Qu.: 0.000 1st Qu.:3.000 1st Qu.:1.0000 1st Qu.:0.0000 Median : 0.000 Median :3.000 Median :1.0000 Median :0.0000 Mean : 1.552 Mean :3.464 Mean :0.7887 Mean :0.4691 3rd Qu.: 3.000 3rd Qu.:4.000 3rd Qu.:1.0000 3rd Qu.:1.0000 Max. :14.000 Max. :8.000 Max. :1.0000 Max. :1.0000 blue gold white black Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:1.0000 1st Qu.:0.000 Median :1.0000 Median :0.0000 Median :1.0000 Median :0.000 Mean :0.5103 Mean :0.4691 Mean :0.7526 Mean :0.268 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.000 Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.000 orange mainhue circles crosses Min. :0.000 red :71 Min. :0.0000 Min. :0.0000 1st Qu.:0.000 blue :40 1st Qu.:0.0000 1st Qu.:0.0000 Median :0.000 green :31 Median :0.0000 Median :0.0000 Mean :0.134 white :22 Mean :0.1701 Mean :0.1495 3rd Qu.:0.000 gold :19 3rd Qu.:0.0000 3rd Qu.:0.0000 Max. :1.000 black : 5 Max. :4.0000 Max. :2.0000 (Other): 6 saltires quarters sunstars crescent Min. :0.00000 Min. :0.0000 Min. : 0.000 Min. :0.0000 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.: 0.000 1st Qu.:0.0000 Median :0.00000 Median :0.0000 Median : 0.000 Median :0.0000 Mean :0.09278 Mean :0.1495 Mean : 1.387 Mean :0.0567 3rd Qu.:0.00000 3rd Qu.:0.0000 3rd Qu.: 1.000 3rd Qu.:0.0000 Max. :1.00000 Max. :4.0000 Max. :50.000 Max. :1.0000 triangle icon animate text Min. :0.0000 Min. :0.0000 Min. :0.000 Min. :0.00000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.000 1st Qu.:0.00000 Median :0.0000 Median :0.0000 Median :0.000 Median :0.00000 Mean :0.1392 Mean :0.2526 Mean :0.201 Mean :0.08247 3rd Qu.:0.0000 3rd Qu.:0.7500 3rd Qu.:0.000 3rd Qu.:0.00000 Max. :1.0000 Max. :1.0000 Max. :1.000 Max. :1.00000 topleft botright black :12 red :69 blue :43 blue :47 gold : 6 green :40 green :32 white :17 orange: 4 black : 9 red :56 gold : 9 white :41 (Other): 3
vapply
前回の記事で以下のコードを実行すると、結果としてリストを返すことが分かりました
（それぞれのリストの要素の長さが異なるため、簡潔化（simplify）できずlapplyと同じ結果となります。）
> sapply(flags, unique)
ここで、もしunique関数の機能を忘れていて、ダブりを除いた結果を数値（numeric）で返すと認識していたらどうなるでしょうか。その場合コードの途中で間違いが起こってしまいます。そのような間違いを防ぐために、vapplyの出番となります。sapplyが結果のフォーマットをRがある程度予測するするのに対し、vapplyでは結果のフォーマットを自分で指定することができます。それにより、予期しない結果を得ることを防ぐことができます。
以下のコードは結果はnumericで長さ１を返すよう指定しています。コードを実行するとどうなるでしょうか。
> vapply(flags, unique, numeric(1))
結果はこのようになり、エラーが発生します。
Error in vapply(flags, unique, numeric(1)) : values must be length 1, but FUN(X[[1]]) result is length 194
例えば前回の記事で行った以下のコードを実行するとsapplyはキャラクターベクターを結果として返します。
> sapply(flags, class)
この結果をもとに、vapplyで以下のコードを実行するとsapplyと同じ結果を返します。これはvapplyにclass関数をflagsのそれぞれの列に実行し、結果はcharacterで長さ１と指定しているからです。
> vapply(flags, class, character(1)) name landmass zone area population language religion "factor" "integer" "integer" "integer" "integer" "integer" "integer" bars stripes colours red green blue gold "integer" "integer" "integer" "integer" "integer" "integer" "integer" white black orange mainhue circles crosses saltires "integer" "integer" "integer" "factor" "integer" "integer" "integer" quarters sunstars crescent triangle icon animate text "integer" "integer" "integer" "integer" "integer" "integer" "integer" topleft botright "factor" "factor"
このようにvapplyを使用することで予期せぬエラーを防ぐことができますが、コードの動作をその都度確認している場合にはsapplyで十分です。
tapply
データを処理する際、データをある変数を基にグループで分け、その分割されたグループに対し、ある計算を行いたいことがあります。tapplyはまさにそのようなことを行います。何を言っているのかわかりずらいと思いますので、例を見ていきましょう。
flagsデータのlandmassは1から6の整数の値を取り、それにより世界の異なる地域を表します。
> table(flags$landmass) 1 2 3 4 5 6 31 17 35 52 39 20
同様にflagsデータのanimateは、もしその国旗が何かキャラクター（鷲、木、人の手など）を含んでいる場合は1を、含んでいない場合は0となります。
> table(flags$animate) 0 1 155 39
この結果より、39の国旗が何かしらのキャラクターを含んでいることが分かります。もしこの平均値を取ると、全体の国旗に対するキャラクターを含んでいる国旗の割合を計算できます。下の結果から20%の国旗で何かしらのキャラクターを含んでいることが分かります。
> 39/194*100 [1] 20.10309
以上を基に、以下のコードを実行してみましょう。
> tapply(flags$animate, flags$landmass, mean) 1 2 3 4 5 6 0.4193548 0.1764706 0.1142857 0.1346154 0.1538462 0.3000000
これはどういう意味かというと、それぞれのlandmassグループ１から６に対し、animateの平均値（mean）を求めています。landmass＝１は北アメリカ大陸を表し、42%の国旗がキャラクターを含んでいます。
また以下の例では、国旗に赤色を含んでいる又はいない国の人口(100万)の概要を示しています。
> tapply(flags$population,flags$red, summary) $`0` Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00 0.00 3.00 27.63 9.00 684.00 $`1` Min. 1st Qu. Median Mean 3rd Qu. Max. 0.0 0.0 4.0 22.1 15.0 1008.0
国旗に赤色を含んでいない国の人口の中央値は3.00であることが分かります。
