2017-07-13

sfパッケージによる空間データの操作

r sf

利用パッケージ

library(sf)
library(tidyverse)
library(viridis)
library(jpndistrict)

データの準備

fukuoka <- jpndistrict::spdf_jpn_pref(code=40) 
class(fukuoka)

[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

spクラスからsfクラスへの変換

fukuoka_sf <- st_as_sf(fukuoka)
class(fukuoka_sf)

[1] "sf"         "data.frame"

CRSをセットする。

fukuoka_sf <- 
  fukuoka_sf %>% 
  st_set_crs(4612)

出力（作図）する

通常のプロット

fukuoka_sf %>%
  st_geometry() %>%
  plot()

f:id:jerrarrdan:20170713223601p:plain

ggplot2を利用したプロット

fukuoka_sf %>%
  ggplot() +
  geom_sf()

f:id:jerrarrdan:20170713223629p:plain

dplyr関連の操作を行う

filter()

北九州市のみを選択する

fukuoka_sf %>%
  filter(city_name_=="北九州市")

Simple feature collection with 23 features and 5 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 129.9819 ymin: 33.00003 xmax: 131.1906 ymax: 34.24994
epsg (SRID):    4612
proj4string:    +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
First 20 features:
   pref_name city_name_ city_name    city_name_full city_code
1     福岡県   北九州市    門司区   北九州市 門司区     40101
2     福岡県   北九州市    門司区   北九州市 門司区     40101
3     福岡県   北九州市    門司区   北九州市 門司区     40101
4     福岡県   北九州市    若松区   北九州市 若松区     40103
5     福岡県   北九州市    若松区   北九州市 若松区     40103
6     福岡県   北九州市    若松区   北九州市 若松区     40103
7     福岡県   北九州市    若松区   北九州市 若松区     40103
8     福岡県   北九州市    若松区   北九州市 若松区     40103
9     福岡県   北九州市    若松区   北九州市 若松区     40103
10    福岡県   北九州市    戸畑区   北九州市 戸畑区     40105
11    福岡県   北九州市  小倉北区 北九州市 小倉北区     40106
12    福岡県   北九州市  小倉北区 北九州市 小倉北区     40106
13    福岡県   北九州市  小倉北区 北九州市 小倉北区     40106
14    福岡県   北九州市  小倉北区 北九州市 小倉北区     40106
15    福岡県   北九州市  小倉北区 北九州市 小倉北区     40106
16    福岡県   北九州市  小倉北区 北九州市 小倉北区     40106
17    福岡県   北九州市  小倉北区 北九州市 小倉北区     40106
18    福岡県   北九州市  小倉北区 北九州市 小倉北区     40106
19    福岡県   北九州市  小倉南区 北九州市 小倉南区     40107
20    福岡県   北九州市  小倉南区 北九州市 小倉南区     40107
                         geometry
1  POLYGON((131.003197471 33.8...
2  POLYGON((131.008333333 33.8...
3  POLYGON((130.977422581 33.8...
4  POLYGON((130.744441712 33.8...
5  POLYGON((130.808888418 33.9...
6  POLYGON((130.759679559 33.9...
7  POLYGON((130.71166083 34.00...
8  POLYGON((130.72547917 34.01...
9  POLYGON((130.72615917 34.02...
10 POLYGON((130.851154955 33.9...
11 POLYGON((130.91506677 33.89...
12 POLYGON((130.855854449 33.9...
13 POLYGON((130.846127782 33.9...
14 POLYGON((130.852372218 33.9...
15 POLYGON((130.832898612 33.9...
16 POLYGON((130.813251115 34.0...
17 POLYGON((130.81524166 34.00...
18 POLYGON((130.809503606 34.0...
19 POLYGON((130.977270078 33.8...
20 POLYGON((131.038985383 33.8...

select()

必要な列を選択する

fukuoka_sf %>%
  select(city_name) %>%
  head()

Simple feature collection with 6 features and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 130.6732 ymin: 33.84096 xmax: 131.0239 ymax: 33.96988
epsg (SRID):    4612
proj4string:    +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
  city_name                       geometry
1    門司区 POLYGON((131.003197471 33.8...
2    門司区 POLYGON((131.008333333 33.8...
3    門司区 POLYGON((130.977422581 33.8...
4    若松区 POLYGON((130.744441712 33.8...
5    若松区 POLYGON((130.808888418 33.9...
6    若松区 POLYGON((130.759679559 33.9...

rename()

列名を変更する

fukuoka_sf %>% 
  rename(Prefecture = pref_name) %>% head()

Simple feature collection with 6 features and 5 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 130.6732 ymin: 33.84096 xmax: 131.0239 ymax: 33.96988
epsg (SRID):    4612
proj4string:    +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
  Prefecture city_name_ city_name  city_name_full city_code
1     福岡県   北九州市    門司区 北九州市 門司区     40101
2     福岡県   北九州市    門司区 北九州市 門司区     40101
3     福岡県   北九州市    門司区 北九州市 門司区     40101
4     福岡県   北九州市    若松区 北九州市 若松区     40103
5     福岡県   北九州市    若松区 北九州市 若松区     40103
6     福岡県   北九州市    若松区 北九州市 若松区     40103
                        geometry
1 POLYGON((131.003197471 33.8...
2 POLYGON((131.008333333 33.8...
3 POLYGON((130.977422581 33.8...
4 POLYGON((130.744441712 33.8...
5 POLYGON((130.808888418 33.9...
6 POLYGON((130.759679559 33.9...

mutate()

新規列の作成する。新たに面積の列を追加する。

fukuoka_sf %>%
  mutate(area = st_area(fukuoka_sf)) %>% 
  head()

Simple feature collection with 6 features and 6 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 130.6732 ymin: 33.84096 xmax: 131.0239 ymax: 33.96988
epsg (SRID):    4612
proj4string:    +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
  pref_name city_name_ city_name  city_name_full city_code
1    福岡県   北九州市    門司区 北九州市 門司区     40101
2    福岡県   北九州市    門司区 北九州市 門司区     40101
3    福岡県   北九州市    門司区 北九州市 門司区     40101
4    福岡県   北九州市    若松区 北九州市 若松区     40103
5    福岡県   北九州市    若松区 北九州市 若松区     40103
6    福岡県   北九州市    若松区 北九州市 若松区     40103
              area                       geometry
1     2700.535 m^2 POLYGON((131.003197471 33.8...
2     2229.020 m^2 POLYGON((131.008333333 33.8...
3 73715391.181 m^2 POLYGON((130.977422581 33.8...
4 62644478.607 m^2 POLYGON((130.744441712 33.8...
5  8136861.212 m^2 POLYGON((130.808888418 33.9...
6     5817.998 m^2 POLYGON((130.759679559 33.9...

arrange()

並び変える

fukuoka_sf %>%
  mutate(area = st_area(fukuoka_sf)) %>% 
  arrange(desc(area)) %>%
  head()

Simple feature collection with 6 features and 6 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 130.039 ymin: 33.10404 xmax: 130.9911 ymax: 33.86936
epsg (SRID):    4612
proj4string:    +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
  pref_name city_name_ city_name    city_name_full city_code          area
1    福岡県       <NA>    八女市            八女市     40210 482531309 m^2
2    福岡県       <NA>    朝倉市            朝倉市     40228 246641750 m^2
3    福岡県       <NA>  久留米市          久留米市     40203 229978807 m^2
4    福岡県       <NA>    糸島市            糸島市     40230 215059938 m^2
5    福岡県       <NA>    飯塚市            飯塚市     40205 213853733 m^2
6    福岡県   北九州市  小倉南区 北九州市 小倉南区     40107 169355048 m^2
                        geometry
1 POLYGON((130.724037224 33.3...
2 POLYGON((130.791228885 33.4...
3 POLYGON((130.666793891 33.3...
4 POLYGON((130.291798612 33.5...
5 POLYGON((130.743817224 33.7...
6 POLYGON((130.977422581 33.8...

sample_n()

サンプリングする

fukuoka_sf %>%
  sample_n(size = 5, replace=T)

Simple feature collection with 5 features and 5 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 130.4423 ymin: 33.10404 xmax: 130.8875 ymax: 33.90448
epsg (SRID):    4612
proj4string:    +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
    pref_name city_name_ city_name city_name_full city_code
100    福岡県     朝倉郡    筑前町  朝倉郡 筑前町     40447
50     福岡県       <NA>    八女市         八女市     40210
65     福岡県       <NA>    宗像市         宗像市     40220
95     福岡県     遠賀郡    岡垣町  遠賀郡 岡垣町     40383
77     福岡県       <NA>    朝倉市         朝倉市     40228
                          geometry
100 POLYGON((130.662384721 33.5...
50  POLYGON((130.724037224 33.3...
65  POLYGON((130.443467471 33.9...
95  POLYGON((130.635642568 33.8...
77  POLYGON((130.791228885 33.4...

group_by(), summarise()

以下の二つの図を比較することで group_by(), summarise()を利用可能性を理解することができる。

# 左の図
fukuoka_sf %>%
  filter(city_name_=="福岡市") %>%
  #group_by(city_name_) %>%
  #summarise() %>%
  ggplot() +
  geom_sf()

# 右の図
fukuoka_sf %>%
  filter(city_name_=="福岡市") %>%
  group_by(city_name_) %>%
  summarise() %>%
  ggplot() +
  geom_sf()

f:id:jerrarrdan:20170713223704p:plain:w300 f:id:jerrarrdan:20170713223716p:plain:w300

join関数

平成28年経済センサス‐活動調査結果のデータを結合する

eco <- readxl::read_excel("28_kigyou_03_0.xlsx",skip=10,na = "-")
eco

# A tibble: 8 x 34
    地域 総数_企業数 総数_事業所数 総数_常用雇用者数 `300万円未満_企業数`
   <chr>       <dbl>         <dbl>             <dbl>                <dbl>
1 福岡市       20737         49288            642559                 1715
2 東  区        2728          5971             94885                  198
3 博多区        5664         19689            292437                  406
4 中央区        5317         12845            165969                  457
5 南  区        2719          4630             38766                  247
6 城南区         958          1341              9125                  112
7 早良区        1900          2685             22685                  155
8 西  区        1451          2127             18692                  140
# ... with 29 more variables: `300万円未満_事業所数` <dbl>,
#   `300万円未満_常用雇用者数` <dbl>, `300～500万円未満_企業数` <dbl>,
#   `300～500万円未満_事業所数` <dbl>,
#   `300～500万円未満_常用雇用者数` <dbl>,
#   `500～1,000万円未満_企業数` <dbl>,
#   `500～1,000万円未満_事業所数` <dbl>,
#   `500～1,000万円未満_常用雇用者数` <dbl>,
#   `1,000～3,000万円未満_企業数` <dbl>,
#   `1,000～3,000万円未満_事業所数` <dbl>,
#   `1,000～3,000万円未満_常用雇用者数` <dbl>,
#   `3,000～5,000万円未満_企業数` <dbl>,
#   `3,000～5,000万円未満_事業所数` <dbl>,
#   `3,000～5,000万円未満_常用雇用者数` <dbl>,
#   `5,000～1億円未満_企業数` <dbl>, `5,000～1億円未満_事業所数` <dbl>,
#   `5,000～1億円未満_常用雇用者数` <dbl>, `1～3億円未満_企業数` <dbl>,
#   `1～3億円未満_事業所数` <dbl>, `1～3億円未満_常用雇用者数` <dbl>,
#   `3～10億円未満_企業数` <dbl>, `3～10億円未満_事業所数` <dbl>,
#   `3～10億円未満_常用雇用者数` <dbl>, `10～50億円未満_企業数` <dbl>,
#   `10～50億円未満_事業所数` <dbl>, `10～50億円未満_常用雇用者数` <dbl>,
#   `50億円以上_企業数` <dbl>, `50億円以上_事業所数` <dbl>,
#   `50億円以上_常用雇用者数` <dbl>

fukuoka_sf %>% 
  left_join(eco, by =c("city_name"="地域")) %>%
  filter(city_name_=="福岡市") %>%
  ggplot() +
  geom_sf(aes(fill=総数_企業数)) +
  scale_fill_viridis(na.value="white")

f:id:jerrarrdan:20170713223917p:plain

2017-07-09

ggplot2でグループごとの図を出力する

r ggplot

ggplotでグループごとに図を出力したい場合がある。意外と簡単な方法でできるのでメモ。

stackoverflow.com

library(ggplot2)
library(dplyr)

res <-
  iris %>% 
  group_by(Species) %>%
  do(
    plots=
      ggplot(data=.) +
      geom_point(aes(Sepal.Length, Petal.Length)) +
      ggtitle(.$Species)
  )

res

 Species
 <fctr>  
 plots  
 <list>  
 1  setosa  <S3: gg>      
 2  versicolor  <S3: gg>      
 3  virginica   <S3: gg>

res[[2]]

f:id:jerrarrdan:20170709154837p:plain f:id:jerrarrdan:20170709154846p:plain f:id:jerrarrdan:20170709154856p:plain

2017-07-08

また株価で遊んでみた

r 株

ここで取得した株価で遊んでみた。 ooooooha.hatenablog.com

建設業関連の株価の終値の増加率の累積地を可視化してみた。全体的に株価が上昇しているのがわかる。塩漬けでも投資しておけば小遣い程度は稼げたんだろうな。これをみると景気拡大期間バブル超えというのも納得できる。実感としては増税等のせい？で減ってきている印象しかありません。

load(file = "output/建設2017-06-18.Rdata")
code_list <- read.table("codelist.csv", sep=",", header = T, check.names = F)


data_list %>%
  as_tibble() %>%
  arrange(code, date) %>%
  group_by(code) %>%
  mutate(prev_end_price = lag(end_price),
         diff_end_price = end_price / prev_end_price -1) %>%
  mutate(abs_rate = diff_end_price %>% abs() %>% max(na.rm=T)) %>%
  filter(abs_rate < 2) %>%
  replace_na(list(diff_end_price = 0)) %>% 
  group_by(code) %>%
  mutate(cum_diff_end_price = cumsum(diff_end_price)) %>%
  ungroup %>%
  mutate(code = plyr::mapvalues(x = code, 
                                from = as.numeric(code_list$コード),
                                to = code_list$銘柄名,
                                warn_missing = F)) %>% 
  ggplot() +
  geom_line(aes(x = date, y = cum_diff_end_price)) + 
  geom_abline(slope = 0, linetype = 2) +
  facet_wrap(~ code, scale = "fixed") +
  theme(axis.text = element_text(size = 7), strip.text = element_text(size = 7)) +
  scale_x_date(labels=scales::date_format("%m"))

f:id:jerrarrdan:20170708001000p:plain

2017-07-03

Rでインタラクティブな地図表示はmapviewで十分

r sf mapview

たまたまmapviewパッケージのヘルプ(example)をみていたら、 mapview()関数を + でつないでいく方法があることを知った。

もう単なる地図表示だけとりあえず表示してインタラクティブに確認したい場合はこれを使えば十分そう。

ただし、データ容量の大きい場合には注意がいりそうなコメントをどこかで見たので、その点は注意が必要そうである。

library(mapview)
library(kokudosuuchi)
library(rgdal)
library(sf)
library(ggplot2)
library(dplyr)
library(sp)

rgdal::setCPLConfigOption("SHAPE_ENCODING", "")

# データの取得

get_url <- kokudosuuchi::getKSJURL(identifier="A38", prefCode = 18)
## 医療圏
A38 <- kokudosuuchi::getKSJData(zip_url = get_url$zipFileUrl)
A38 <- A38[[1]]

## 医療機関
get_url <- kokudosuuchi::getKSJURL(identifier="P04", prefCode = 18)
P04 <- kokudosuuchi::getKSJData(zip_url = get_url$zipFileUrl)
P04 <- P04[[1]]

plot(A38, axes=T)
plot(P04, add=T, col="green")

f:id:jerrarrdan:20170703232749p:plain

A38 <- A38 %>% st_as_sf() %>% st_set_crs(4612) %>% st_transform(2448)
P04 <- P04 %>% st_as_sf() %>% st_set_crs(4612) %>% st_transform(2448)

m <- mapview(A38, zcol="A38a_004") + mapview(P04, zcol="診療科目1")

原因不明だが、このmapview()関数の順序を逆にするとうまく表示されなかった。

mapshot(x = m, file = "fukui_iryou.png")

f:id:jerrarrdan:20170703232000p:plain

ついでに、医療圏域別に診療科目の開設？状況について表にしてみた。人口別にみる必要があると覆うがかなり偏りのある分布になった。徒歩圏内に医療がない地域も多そうで、医療の問題は大変そう。

歯科が多いのはよく聞くが、内科や小児科ってかなり多いんですね。安心して子供を産める地域といえるのかもしれませんが、小児科の数は少子化社会に対して妥当な数なんでしょうかね。

診療科目1に複数の科目が入力されているため切り分けて集計する。

st_intersection(x = P04, y = A38) %>%
  as.data.frame() %>%
  select(診療科目1, A38a_004) %>%
  tidyr::separate(col = 診療科目1, into = paste0("診療科目", 1:23), sep = "　") %>%
  mutate(rn = row.names(.)) %>% 
  gather(key = 科目, value = prof, starts_with("診療科目")) %>%
  xtabs(~A38a_004+prof, data = .) %>%
  as.data.frame() %>%
  spread(A38a_004, Freq) %>%
  mutate(計 = `奥越` + `丹南` + `福井・坂井` + `嶺南`) %>%
  arrange(desc(`計`))

prof	奥越	丹南	福井・坂井	嶺南	計
内科	25	93	191	62	371
歯科	21	64	158	46	289
小児科	15	56	113	34	218
消化器科	13	35	83	21	152
外科	10	39	64	23	136
整形外科	9	33	74	19	135
リハビリテーション科	8	27	70	15	120
循環器科	8	27	72	11	118
皮膚科	7	25	45	17	94
呼吸器科	6	18	50	7	81
眼科	4	14	35	9	62
耳鼻咽喉科	4	11	31	11	57
放射線科	3	21	25	8	57
小児歯科	0	12	31	6	49
肛門科	5	12	26	5	48
リウマチ科	2	7	34	4	47
泌尿器科	2	13	22	7	44
アレルギー科	4	2	30	5	41
胃腸科	2	12	19	4	37
精神科	1	4	22	9	36
矯正歯科	0	7	25	2	34
脳神経外科	1	9	18	4	32
産婦人科	2	6	16	6	30
神経内科	2	5	19	4	30
歯科口腔外科	0	4	20	5	29
麻酔科	2	7	15	5	29
形成外科	1	6	15	2	24
婦人科	1	8	9	5	23
心療内科	0	3	14	3	20
神経科	0	3	8	5	16
気管食道科	1	0	11	1	13
心臓血管外科	0	3	7	3	13
産科	0	3	6	1	10
呼吸器外科	0	1	4	3	8
性病科	0	4	1	0	5
美容外科	0	1	2	0	3
小児外科	0	0	2	0	2
皮膚泌尿器科	0	1	0	0	1

歯科はコンビニより多いということをどこかのTV番組できいたので、多いらしいということは知っていたけど、それを上回って内科が多い。ちなみに福井県の場合、コンビニ数はサイト確認時で353店舗なので、内科の病院のみがコンビニ店舗数を上回っている。

2017-06-28

ポリゴン内にテキスト（ラベル）を表示する

r sf sp plot

ポリゴン内に各ポリゴンの属性を表示したい。何回か調べているので、メモとして残す。

# using packages
library(ggplot2)
library(kokudosuuchi)
library(sf)
library(sp)
options(stringsAsFactors = F)

# preparation
url <- kokudosuuchi::getKSJURL("N03", prefCode = c(27))
url <- url %>% arrange(desc(year)) %>% .[1, ]
N03 <- kokudosuuchi::getKSJData(url$zipFileUrl)
N03 <- N03[[1]]
proj4string(N03) <- CRS("+init=epsg:4612")

# sp class
plot(N03, axes=T)
text(getSpPPolygonsLabptSlots(N03), labels = as.character(N03$市区町村名), cex=0.4)

f:id:jerrarrdan:20170628234700p:plain

# sf class
N03_sf <- st_as_sf(N03)
N03_sf %>% 
  st_geometry() %>%
  plot()
text(st_coordinates(N03_sf %>% st_centroid), labels=N03$市区町村名, cex = 0.3)

f:id:jerrarrdan:20170628234712p:plain

# using ggplot2
N03_sf %>%
  mutate(lon = st_centroid(.) %>% st_coordinates() %>% .[, 1],
         lat = st_centroid(.) %>% st_coordinates() %>% .[, 2]) %>%
  ggplot() +
  geom_sf(aes(fill = 市区町村名)) +
   ggrepel::geom_text_repel(aes(x=lon, y=lat, label=市区町村名), size=2)

f:id:jerrarrdan:20170629004605p:plain

こちらはstackoverflowにあった回答例からのコピペ。 ggplot2パッケージでもっと簡単にテキストを作図したいが、わからず。

2017-06-18

ここ一年のTOPIX30の株価（終値）の推移

r 株価遊び

ここ一年のTOPIX30の株価（終値）を可視化してみた。

便利なパッケージがあるらしいが、うまく機能しなかった。 yahooではスクレイピングを禁止しているらしい、また、させないように仕様を頻繁に？変更しているらしいので、そのほかでできそうな記事を参考にやってみた。

これだけで各株のトレンドが一目瞭然なので便利そう。

library(tidyverse)
library(lubridate)
options(stringsAsFactors = F)

# 証券コードと名前の変換用
code_list <- 
  rio::import("http://www.jpx.co.jp/markets/statistics-equities/misc/tvdivq0000001vg2-att/data_j.xls")

topix30 <- code_list %>% filter(規模区分 == "TOPIX Core30")
code <- topix30$コード

stock <- NULL
for (i in seq_len(length(code))){
  url <- paste("http://k-db.com/stocks/", 
               code[i], 
               "-T?1d/2017?&download=csv", 
               sep = "")
  stock[[i]] <- read.table(url, sep = ",", skip = 2)
  colnames(stock[[i]]) <- 
    c("date", "start_price", "max_price", "min_price",
      "end_price","volume","trading_value")
  # dateをでdate要素に変更
  stock[[i]] <- stock[[i]] %>% mutate_at(.vars = "date", ymd)
  Sys.sleep(5) 
}
names(stock) <- code

data_list <- NULL
for (i in seq_len(length(stock))){
  data_list <- rbind(data_list, 
                     stock[[i]] %>% mutate(code = as.character(code[i])))
}

p <- 
  data_list %>%
  mutate(code = plyr::mapvalues(x = code, from = topix30$コード, to = topix30$銘柄名)) %>% 
  ggplot() +
  geom_line(aes(x = date, y = end_price)) +
  facet_wrap(~code, scales = "free_y") 

p

点は月曜日である。 JR系の株価がうなぎのぼり。

f:id:jerrarrdan:20170618165811p:plain

データの取得に関してはほとんどこのサイトのパクリです。

2017-06-06

自作関数の高速化

r bytecompile 高速化

efficient R programingを読んでいて今後使えそうなことに関するメモ。

3.7 では自作関数の高速化について述べられている。

自作関数にcompiler::cmpfun()を適用すればよいだけらしい。

library(compiler)

mean_r = function(x) {
  m = 0
  n = length(x)
  for(i in seq_len(n))
    m = m + x[i] / n
  m
}

cmp_mean_r = cmpfun(mean_r)

平均を求める関数についての比較結果はこのようになる。結構差が出ているので、複雑な処理をする際には有効に使えそうである。

x = rnorm(1000)
microbenchmark(times = 10, unit = "ms", # milliseconds
          mean_r(x), cmp_mean_r(x), mean(x))
#> Unit: milliseconds
#>           expr   min    lq  mean median    uq  max neval cld
#>      mean_r(x) 0.358 0.361 0.370  0.363 0.367 0.43    10   c
#>  cmp_mean_r(x) 0.050 0.051 0.052  0.051 0.051 0.07    10  b 
#>        mean(x) 0.005 0.005 0.008  0.007 0.008 0.03    10 a

2017/7/23追記 R言語逆引きハンドブックにバイトコンパイルによる高速化として、しっかり書いてあったので、詳細はそちらを確認。