2018年10月8日（火）

Rmarkdownは便利だけど数式相互参照は発展途上

またもや日記を。こんどは細大漏らさず記録するというより、気が向いたことを記録しよう。3日坊主でおわるのもなんだか悔しいし。

Rmarkdownを使って論文を書きはじめた。いい感じ。

数式（\(Y_{it} \sim N(h_{it}, \sigma^{2}_{yi})\)）を相互参照しようとして手間取る。markdown記法で数式を書くと、相互参照がうまく使えない。Rstudio上でプレビューされるはずの数式も表示されなくなる。図表については相互参照がスムースに作動するから、ちょっと意外。

$$
Y_{it} \sim N(h_{it}, \sigma^{2}_{yi})
$$

LaTeX記法（）を使って書けば相互参照も使えるようになるかもとおもったが、なぜだかできなかった。本文中に??と表示されておわり。

\begin{equation}
Y_{it} \sim N(h_{it}, \sigma^{2}_{yi}) \label{eq:hoge}
\end{equation}

どうやらBookdownが数式を相互参照することに対処していないみたい。これは記法を考えたくらいでは解決しないみたい。

結局、提出（提出はなぜかWord）する前にWord側で数式番号を振り直すことに。

Rmarkdownだけで論文をすべて完成させるにはまだまだ。ただ、MSWordで書くよりは断然楽できれい。図表を扱う作業は楽。書式も整えてくれる。LaTeXで書いてもいいけれど、あとからMSWordに変換するのがとっても面倒。その点、Rmarkdownは融通がきいて便利。ElsevierあたりはすべてRmarkdownで完結するんだろうか。

ところで、いつもおもうのだけれども、経営とか労働について論文を書く先生方は、どうしてあんなに長い論文を書けるのだろう。自然科学みたいにせいぜい15ページくらいでもいいようにおもう。

とど

  • HRMタイランド：論文執筆。次調査準備。そろそろ機材調達しないと。
  • プラスチックごみ関係：データ整理。とくにログデータ。
  • 統計学習会関係：スライド作成。相手が想定よりもRが苦手かもという一大事。
  • その他もろもろ：いまのところ特になし。

2018年10月9日（水）

執筆

本日はまぁまぁ順調。書きたいことを\(\frac{1}{4}\)半分くらいは書いたろうか。考えきれてない箇所は当然書けない。以前書いた内容もあるけれど、忘れているところ多し。猛省。それでも、相互参照がうまく機能することやbibtexがうまく動作することも確認できた。まぁまぁの一日だったとしよう。寝る前に論文でも1,2本読むことができれば上等。ほんとうはもう少し読みたいけれど、講義やらなんやらいろいろあってなんともはや。人的資源を水産資源に代表される再生産可能な資源であると認識してもらう説明がうまくできるかがポイントな気がする。

その講義、ことしははじめて日本人学生が2名も混在している。見る限り、英語についてきてない。

さて、\(\LaTeX\)用IDEみたいに補完機能があるとうれしい。いまのところは実装されてないみたい。将来的にRmarkdownを使って論文執筆する方々が増えれば、Rstudioに補完が実装されることは考えられる。なにしろRコードについては十分すぎる補完があり実績十分。補完に加えて、できれば日本語エディタとしてもう少しRstudioによる日本語表示（特に入力時）が改善されるとうれしい。

formattable関数を使って美麗な作図を自動化するため、MCMCを回し直す。前回パラメータ名に適当な記号をつけたせいで、いま執筆している論文とパラメータ名に使う記号が整合しない。記号を変更するだけだから結果はおなじはずで、時間がかかるだけではある。昨晩からぶん回しているけれど、まだ終わらない…結果を取り出して記号だけを置き換えるという手もあるが、時間もありそうなので回し直すことに。

12月に出かけるRおじさんの旅について手続きはほぼ順調。懸念はビザ取得料支弁に関する細々としたことと以下。

  • R初心者が多すぎはしないか：R経験者と明記してもらったけれど、CTUから到来したメールを読む限りR初心者が混ざっている気がする。RとRstudioのインストールやらCUI全般について一応準備をしつつ、内容を少々圧縮しないといけないかも。頭がいたい。
  • 教材を完成させるには時間が少々不足：水産学部対象Rおじさん終了後、E-4対象Rおじさんがある。こちらはRやRstudioにはじめて触れる人々向け。E-4向けにアレンジできるあたりが前回よりも助かる。部分的に参加してもなんとかやりたいことができるような内容を考えないといけないところはなんともはや。

仮説山盛り論文は…

仮説が明示されこれが検証されるような論文は読みやすい。もっとも、経営学論文にありがちな仮説が10も20もあるような論文は読みにくい。あまりに読みにくいため、近頃は読まないことにしている。Hypotheses \(n (n \in \mathbb{N})\)みたいな論文は、なにを真に明らかにしたいかわからないではないか。たいていは既存研究を適当につないでつくった仮説で、理論的背景も因果推論もなにもあったものではない。

とど

  • HRMタイランド：昨日同様。論文執筆。次調査準備。そろそろ機材調達しないと。

2018年10月10日（木）

Rから官庁統計

きょうは昨日読んだ文献メモを入力する以外、ほとんど研究関連用務がこなせなかった。こういう日もあるだろうとはおもうが、なんともはや。

おもな用務はASEANのいろいろと講義準備。ASEAN社会指標を扱うときになにかいいパッケージでもないか探していたら、WDIというパッケージが見つかった。おなじみILOのデータベースのAPIをたたくILOStatや日本の官庁統計をさわるパッケージ（estatapi）みたいに、APIをたたいてデータベースから必要なデータをいただくパッケージみたい。E-4ワークショップでも紹介してみようか。

きょうはASEAN各国における総人口推移を表現する折れ線グラフをつくる。応用すればいろいろな指標について折れ線グラフが簡単につくれる。自動化してあらゆる変数について折れ線グラフをつくることさえ可能だろう。もっとも、その場合は1,000を超える指標について作図する相当な時間を覚悟しないといけないし、ひとつひとつの指標について勉強するという目が回る作業が待っている。

では、手順を。

  1. ライブラリを読み込む：使うライブラリはWDI。World Development Indicatorsを省略した表現であろう。作図にggplot2パッケージを使う都合、tidyverseもついでに読み込む。あとでRmarkdownを使ってレポートを書くことを考えて、本通りにヘッダをつけておく。今回はデータ読み込みと作図までのセクションとプロットのセクションを分ける。表示は図だけなので、こうすると余計な場所を表示しなくて済む。
# ---- readlibrary ----
# read library #####
library(WDI)
library(tidyverse)
  1. データを取得する：WDIパッケージを用いて世銀からデータを引っ張る。必要な変数は国名と指標名と開始・終了年。国名はisoによる2レターコード。BN（ブルネイ）、KH（カンボジア）、ID（インドネシア）、LA（ラオス）、MM（ミャンマー）、MY（マレーシア）、PH（フィリピン）、SG（シンガポール）、TH（タイ）、VN（ベトナム）があれば、当座はよろしい。指標名は略称。WDIsearch関数を使えば見つかる。そうでなければWDIのページで探すと見つかるはず。コード自体は先程のセクションのつづき。
# read data from the world bank
## total population
poptotal<- WDI(country = c("BN","KH","ID","LA","MM","MY","PH","SG","TH","VN"),  
                    indicator = "SP.POP.TOTL", 
                    start=1970, 
                    end=2015,
                    extra=TRUE, 
                    cache=NULL
)

  1. データを加工する：読み込んだデータは、けっこうそのまま使える。特に加工する必要はない。塗り分け地図を書くためか、各国首都緯度・経度まで取得できる。

  2. 作図する:おなじみggplot2関数を使う。作図速度は遅いが美麗に作ることができる。折れ線グラフ右端に国名を付した（参考ページ）。

## total population
fig.poptotal <- ggplot(poptotal, aes(x = year, y = SP.POP.TOTL,colour = country)) + 
  geom_line(size = 1, stat = "identity") + 
  geom_point(size = 2) +
  geom_text(
    data = subset(poptotal, year == max(year)),
    aes(label = country),
    nudge_x = 5
  )

あとはプロットすればいい。軸ラベルを変更し背景を削除した。

# ---- lineplot.pop ----
# plot
fig.poptotal + xlab("Year (1970-2015)") + ylab("Total population (Unit: person)") + theme(legend.position = "none")

実行結果。図をクリックすると大きく表示される。ggplotの型通りの折れ線グラフ。

  1. 備考：官庁統計全般に言えることだけれど、信頼に足りるデータか常に確認すること。どうやら2008年以降のデータしか信用ならないという風のうわさもあったりなかったり…。

きょうの格言

立てばテイラー　すわればガント　歩く姿はギルブレス。

立てば研究　すわれば作図　歩く姿も分析に

経営学はかくあるべし

とど

  • HRMタイランド：昨日同様。論文執筆。次調査準備。そろそろ機材調達しないと。

2018年10月11日（金）

お買い物は自分で

サミュエルソンと秘書にまつわる逸話によれば、秘書がどれほど信頼できなくても秘書ができることを任せたほうがいいらしい。が、そうではないときょうは悟った。自分でやるか信頼できるプロに都度外注したほうがいい。

ということで、10月末から11月上旬に出かけるタイ出張準備。もっとも手がかかる物品調達についてきょうは考えた。というより、これくらいしか仕事が手につかなかった。頭をつかう仕事は静謐であることが必要条件ですな。

まずは物品。

  • USBフラッシュメモリ128GB:TS128GJF810。大きなデータを収納する用。統計資料や写真撮影した統計資料を保管する用。16GBもあれば十分ではあるが、OS入れ替えやらなんやらなにかと使う。防水らしい。あらっぽい環境で仕事をすることも多いし使ってみることに。
  • USBフラッシュメモリ32GB:TS32GJF720S。予備用。手持ちUSBメモリがなぜか海外で使えず往生したことがあった。MLCパッケージという読み込み速度が速いみたい。
  • 非破壊型ブックスキャナ:ET-16 plus。高い。が、いろいろ悩んでこれにした。撮影時にトリミングや台形補正や分厚い本の膨らみ補正やらを自動的にしてくれるソフトが魅力。なにしろ統計資料は分厚い。カメラで撮影するには工夫が必要。本体は1.5kgあるらしい。ケースがついているみたいだが、渡航時は発泡スチロールかなにかで梱包にもうひと工夫したほうがよさそう。
  • 4ポートUSBハブ:Anker Ultra Slim 4-Port USB 3.0 Data。 Hub。手持ちが相手を選ぶことがわかった。これなら薄いし相手を選ぶこともないだろうと信じたい。

つづいて文書。

  • NSOあてレター:前回、データをもらうときにレターが必要と要請を受けた。今回ひそかに目をつけているデータをいただくことができるならば、やはりレターが必要だろう。あらかじめ用意しておけば、手間も省けるというもの。
  • 研究計画書:CMU向け。簡単なものでいいので。手元資料に。

これなんだろう

論文をいじっていておもしろい現象を発見。発見事実としておもしろい。無論、理論的にも示唆に富んでいる。どうやら、特定時期に年齢選択的に発生する現象があるみたい。その現象はしだいにおさまる/ひどくなることもわかった。これらは解析してはじめてわかった。統計学って素敵。統計に逃げてよかったよかった。

わが国事例でおもいつくのは、希望退職や採用抑制だろうか。経時的におさまる/ひどくなることもある程度理解される。でも、希望退職だけ、採用抑制だけ、とそれぞれだけを実行するだろうか。何らかの人員整理をするなら同時に実施することが合理的であろうし、別個の時期に別個の方策だけを実施する理由も考えないといけないだろう。

特定時期に年齢選択的に現象が発生する理由ってなんだろう?理由がわかったとして、機序はどうだろう。機序がわかったとして、これをどうモデルに組み込んでまた検証しよう。

さしあたって、手元にあるデータをもう少し詳しく分析する、ねらっているデータで新たな観点から検討するなどだろうか。

それにはまず今書いている論文を仕上げないと。10月はなかなか忙しい。

文献を読む

少しずつでも文献を読むと気が休まる。

手持ちも少なくなってきたので、方法論に関する論文をいくつか集めた。社会科学においても適用例はけっこうあった。おもに計量経済学だった。経営学にはあまりなかった。まぁいいか。

とど

  • HRMタイランド：論文執筆。次調査準備。機材調達。

2018年10月17日（水）

日時データを取り扱う

Mo OであつめたGPSデータをつかって、移動時間と距離を算出する。算出した移動時間と距離から、ゴミ拾いにかかる費用を推定する。

きょうの目標は、

  1. kmlデータを読み込む：GPSロガーが記録したデータは、そのままだとRで読めない。したがってなにも計算できない。gpx形式やkml形式に変換すればRが読める。ただし、読み込みにちょっとしたコツがいるみたい。
  2. 読み込んだデータから必要な箇所を抽出する：読み込んだkml形式データ内にいくつか階層（レイヤ）がある。必要なデータがある階層から必要なデータを取り出す。
  3. 移動距離を算出する：kml形式で記録されるデータには、移動距離も書いてる。ただ、100m単位と非常に大きい。記録された位置情報からm単位で移動距離を測定したい。
  4. 日付データを取り扱う：簡単に見える（実際Rなら簡単）。が、ちょっとしたコツがある。ヒントはロケール。
  5. 使える形式にデータを整理する：後々利用できるように1.から3.までを含むデータセットを作りたい。せっかくなのでtidyverseな形式に整える。

かなり多いなぁ。それでははじまりはじまりー。

kmlデータを読み込む

必要なライブラリを読み込む。rgdalパッケージはkmlファイルを読み込むために使う。kmlファイルを読み込むことができるライブラリは他にもあるが、shpファイルを読むときにも使えて汎用性が高いから採用。geosphereパッケージは緯度・経度を2点与えて移動距離と移動方向を算出するために使う（参考）。tidyverseはいつもの通り。とにかく便利でいろいろ速い。

kmlファイルがひとつだけなら直接読み込めばいい。しかし今回は9つもある。たいていはいくつもあることが想定される。できるだけ自動化したい。いくつもあるkmlファイルを単一フォルダにまとめておいて、ファイル名を取得させてリストをつくる。つくったリストを使ってあとからファイルを読み込む。

list.filesはbaseパッケージについている関数。pasteも同様。うまくパスが通るように、適当に文字列を足してあげる。“../”をみると、あややをおもいだすこの頃。元気にしているかしらん。

# パスのリストをつくる
## kmlファイルがおいてあるフォルダのパスを指定
path <- "../scattered_item/selected_kml_sep_2018"
## kmlファイルリストを取得する
list.files <- list.files(path)
## kmlファイルリストを後々参照できるよう編集
list.files <- paste("../scattered_item/selected_kml_sep_2018/",list.files, sep = "")

いよいよ必要なデータを取り出す。古典的ではあるが、forループを使って繰り返し処理する。ほかにも早くてうまい手があるとはおもわれる（実はわたしが知らないだけ）が、確実に動作する方法を選ぶ。空のデータフレームを作る方法は、FOR-RESTを参照させていただいた。ありがとうございます。空のデータフレームの各列名は、一度データを読ませて手動で取得した。

readOGR関数は、kmlファイルも読める（もちろん他形式もいろいろ読める）関数。kmlはいくつかレイヤーがあるが、必要なデータはwaypointsというレイヤーに入っていた。レイヤー名はデータを読み込んでogrDrivers関数を実行すると表示される。なかなか全自動というわけにはいかないところがなんともはや。

distGEO関数は距離を算出する関数。単位はメートル。10進法による緯度・経度を与える。幸い、GPSは10進法で値を返してくれた。kmlはGoogleマップに適した形式で、Googleマップは10進法で座標を表示する。あたりまえといえばあたりまえ。bearing関数は移動方向を返す。右方向に移動する場合は正、左方向へ移動する場合は負の値を返す。単位は度。ラジアンではない。座標系は指定できるが、GPSが返す座標系は同関数のデフォルトみたい。

最後に適切な名称をつけた変数に取り出したデータを代入する。変数名を変えなくてもいいけれど、ループ内であやし挙動が発生しているかもしれないときに役に立つことがある。

# 必要なデータを取り出す
## 取り出すデータのいれものをこしらえる
datname.null <- NULL
n <- 7
dat.null <- data.frame(matrix(rep(NA, n), nrow=1))[numeric(0), ]
colnames(dat.null) <- c("lonlat.coords.x1","lonlat.coords.x2","lonlat.coords.x3","distance","direction","period","quadrat")

## データ抽出
##  古典的なforループ。遅いけれど確実に動作するし、なにをしているか読める。
for(i in 1:length(list.files)){
  # リストにあるkmlファイル読み込み
  # rgdal関数が大活躍　
  target.file <- rgdal::readOGR(list.files[i], layer = "waypoints")
  # データのいれものをこしらえる
  distance.null <- NULL   # 距離
  direction.null <- NULL  # 方角　
  points <- length(target.file@coords[,1])-1
    # 距離と方角算出
    for(j in 1:points){
      # 前の点からの距離を算出
      distance <- geosphere::distGeo(c(target.file@coords[j,1], target.file@coords[j,2]), c(target.file@coords[j+1,1], target.file@coords[j+1,2]))
      # 前の点からの方角を算出
      direction <- geosphere::bearing(c(target.file@coords[j,1], target.file@coords[j,2]), c(target.file@coords[j+1,1], target.file@coords[j+1,2]))
      # 算出した距離と方角を各々代入。
      distance.null <- c(distance.null, distance)     # 距離
      direction.null <- c(direction.null, direction)  # 方角
    }
}
# assign the data frame into a new data frame
df.gps.01 <- dat.null
head(df.gps.01)
datname.gps.01 <- datname.null

なんかいろいろメッセージが出るけれど気にしない～気にしない～気にしない～。

日付データを取り扱う

GPSから緯度経度以外にも日時を表現するデータを得られる。これを使わない手はないだろう。地図上に点を描いたり移動距離折れ線グラフを描くだけなら、さほど面倒ではない。しかし時刻が重要であることは海象に影響される現象においてはよくある出来事である。

Rには日付データを取り扱う関数がいくつか準備される。今回はstrptime関数を使う。楽勝楽勝と考えていたが、意外や意外、NAばかりが返る結果に。調べてみると、ロケールが英語でないと動かないらしい。ヘルプにもそう書いてあった。ヘルプはよく読みましょう。

stringrライブラリのstr_sub関数は必要な箇所の文字列を取り出す関数。そういえばチートシート（英語）が出てたなぁ。チートシートとっても便利。問題は保管場所。

ロケールさえ英語になってしまえばあら不思議。あっという間でしたとさ。

# 日付データ
## ロケールを英語にする。
## 英語にしないと動かない。
Sys.setlocale("LC_TIME", "English")
## POSIXct形式に変換。
df.gps.01$period <- as.character(stringr::str_sub(df.gps.01$period, start = 5, end = 24)) %>% strptime(format = "%b %d %H:%M:%S %Y", tz = "Asia/Ho_Chi_Minh") %>% as.POSIXct()

tidyverseにて取り扱うには、strptimeが返すPOSIXlt形式よりもPOSIXct形式がいい。ゆえにデータ形式を変換。

tidyverseっぽく、 %>% なんて使ってみましたよ。なれると便利。だけど従来の代入記号<-とか=もまだまだ必要。結局両方覚えないといけないのね。

使える形式にデータを変換する

## as_tibble and rename
df.gps.01 <- as_tibble(df.gps.01) %>% rename(lon = lonlat.coords.x1, lat = lonlat.coords.x2, hoge = lonlat.coords.x3, distance = distance, direction = direction, period = period, quadrat = quadrat)
## reset the locale as default
Sys.setlocale("LC_TIME", "Japanese")
## [1] "Japanese_Japan.932"

今後の予定

  • 移動速度を計算してデータセットに加える。
  • いろいろ作図。
  • 昔考えたモデルを使って分析。

教材をつくりはじめる

ぼちぼちRワークショップ用教材を作成開始。『みんなのR』14章を研究室で勉強するがてら、14章説明資料をつくった。

交通事故に巻き込まれるの記（つづき）

一通りの手続きは終了。あとは療養に専念するのみ。

とど

  • LFSいろいろ。タイ調査用機材は発注した。