Rのggplot2で福島県南相馬市の地図を描く～もうすぐ野馬追
ggplot2で南相馬市の地図を描きたい
こちらでは、南相馬市と、福島第一原子力発電所を地図上にプロットします。
パッケージは、
jpndistrict package | R Documentation
sfとggplot2、geosphereで作図を行います。
#install these packages if needed
#install.packages("jpndistrict")
#install.packages("sf")
#install.packages("geosphere")
#install.packages("tidyverse")
library(jpndistrict)
library(sf)
library(geosphere)
library(tidyverse)
都道府県・市区町村の情報を取得
地方公共団体コードは以下を参照し、はじめの五桁を使用する。
コードは6桁の数字からなり、上2桁はJIS X 0401に定められた都道府県コード、続く3桁と合わせた上5桁がJIS X 0402に定められた市区町村コード及び一部事務組合等コード、最後の1桁は誤記入・誤入力を防ぐための検査数字(チェックディジット)です。
#Fukushima: Prefectural code 7 FKS <- jpn_pref(pref_code = 7) #これでは宮城県が含まれないので、足し合わせたものを作成する FKS <- c(4,7) %>% map(~jpndistrict::jpn_pref(pref_code = .)) %>% reduce(rbind) %>% st_simplify(dTolerance = 0.001) #Minamisoma (JIS code 07212) Minamisoma <- jpn_cities(jis_code = "07212")
地図に掲載する点を作成する
緯度経度の取得は以下などを参考に。
##Fukushima Daiichi nuclear power plant
FDNPP <- data.frame(longitude=c(141.033611), latitude=c(37.421389))
FDNPP_sf <- st_as_sf(FDNPP, coords = c("longitude", "latitude"),
crs = 4326, agr = "constant")
##MMGH
MMGH <- data.frame(longitude=c(140.9851348), latitude=c(37.6375094))
MMGH_sf <- st_as_sf(MMGH, coords = c("longitude", "latitude"),
crs = 4326, agr = "constant")
位置を微調整
gglocatorなどもありますが、 手作業で位置をずらしてみました。
#deciding map range using Miharu Town location mapRange <- c(range(st_coordinates(Minamisoma)[,1]),range(st_coordinates(Minamisoma)[,2])) mapX <- 0.1 mapY <- 0.15 #Main figure FKS %>% ggplot() + geom_sf(color="white") + geom_sf(data=Minamisoma, color="black")+ geom_sf(data = FDNPP_sf, color = "black", pch=13, size=4)+ geom_sf(data = MMGH_sf, color = "black", pch=10, size=4)+ coord_sf(datum = NA, xlim = c(mapRange[c(1)]-4*mapX, mapRange[c(2)]+1.5*mapX), ylim = c(mapRange[c(3)]-1.5*mapY, mapRange[c(4)]+1*mapY)) + theme_bw()+theme(axis.title = element_blank(), axis.text = element_blank()) -> p_fks #p_fksで出力して、coord_sf内の引数を微調整
特定の場所を中心とする円を描く
福島第一原子力発電所からの距離に応じた円を付記する
#Distance from FDNPP
FDNPP_equidistant <-
map_dfr(c(2.5,5),
function(x){
destPoint(FDNPP, -180:180, x*10000) %>%
as.tibble() %>%
# どの等距離線るかを示すグループ番号
mutate(group = as.character(x))
})
# labels: distance from FDNPP
FDNPP_labels <-
map_df(c(2.5,5),
function(x){
destPoint(FDNPP, 260, x*10000) %>%
as.tibble() %>%
mutate(distance = paste0(x*10, " km"))
})
完成
地図として出力する。
p <- p_fks +
# plotting FDNPP 25km and 50km circles
geom_path(data = FDNPP_equidistant,
aes(x=lon, y = lat, group = group),
colour = "black", alpha=0.3, size = 1.5) +
# plotting labels: distance from FDNPP
geom_label(data = FDNPP_labels,
aes(x=lon, y = lat, label = distance), size = 5) +
annotate("text", x = as.numeric(FDNPP[1]), y = as.numeric(FDNPP[2])+0.04, label = "FDNPP", size = 5)+
annotate("text", x = as.numeric(MMGH[1])-0.06, y = as.numeric(MMGH[2])+0.04, label = "MMGH", size = 5)+
annotate("text", x = as.numeric(FDNPP[1])-0.10, y= as.numeric(FDNPP[2])+0.36, label = "Minamisoma", size = 5)
p+theme_bw(base_family = "HiraKakuProN-W3")+theme(axis.title = element_blank(), axis.text = element_blank())
#saveするときは
#ggsave(file = "Figure1_jp.png", plot = p)
もうすぐ野馬追
あたたかくなってきました。もうすぐ野馬追です。
今年はコロナウイルスの影響で無観客なのか。
Rで質問紙調査のサンプルサイズ計算を行う
今日は、質問紙調査でのサンプルサイズ計算をrで行います。
を用います。
早速、関数からみていきましょう。
install.packages("samplesize")
library(samplesize)
sample.size.prop(e=, P=, N=, level=)
Arguments
e positive number specifying the precision which is half width of confidence interval
P expected proportion of events with domain between values 0 and 1. Default is P=0.5.
N positive integer for population size. Default is N=Inf, which means that calculations are carried out without finite population correction.
level coverage probability for confidence intervals. Default is level=0.95.
引数が4つありますが、特に大切な数字は2つ。
e: precision（精度）
たけのこの里が好きな人が標本内に90%いるとしたときに、精度が5%だと、実際は、85~95%のたけのこの里ファンが居ることになります。 5%が一般的な数値であり、1~10%に収めるのがよいです。 小さくするとたくさんのサンプルサイズが必要になります。
level: confidence level（信頼水準）
標本が、母集団を反映することをどのレベルで信頼したいか、という指標。
95%が最も一般的です。デフォルトもlevel=0.95になっています。
大きくすると、たくさんのサンプルサイズが必要になります。
残り２つは、不明であればデフォルトで構いません。
P: 期待されるイベントの割合
実際に予想されるイベント頻度、prevalenceがわかっていればそれを入力します。 デフォルトではP=0.5になっています。
N: サーベイをおこなう全対象者数
わかれば、入力を。わからなければ、入力しないでOKです。 デフォルトはN=Infになっています。
いくつか試してみる
eを0.01~0.1で動かす
以下のようになりました。 許容誤差が小さいと、必要サンプルサイズがとても大きくなります。
eは、標準的には0.05、一歩譲って0.075、百歩譲って0.1くらいがよいでしょう。
install.packages("samplingbook")
library(samplingbook)
sample.size.prop(e=0.01, P=0.5, N=Inf, level=0.95)
sample.size.prop(e=0.05, P=0.5, N=Inf, level=0.95)
sample.size.prop(e=0.075, P=0.5, N=Inf, level=0.95)
e0<-seq(0.01,0.1, length=19)
res0<-rep(0,19)
for(i in 1:19){
res0[i]<-sample.size.prop(e=e0[i], P=0.5, N=Inf, level=0.95)$n
}
plot(cbind(e0,res0), xlab = "許容誤差", ylab = "サンプルサイズ")
WindowsでRstanをインストールした話~ Compilation ERROR の解決
Windows PCでrstanを使うべく、進めてみた、がうまくいかない。
compilation error
その解決の備忘録として記しておく。
stan公式
RStan Getting Started (Japanese) · stan-dev/rstan Wiki · GitHub
スペック
Windows 10 64 bit
R version 3.6.1 (2019-07-05)
Rstudio 1.2.5001
rstan 2.19.3
ネットに山積するエラー報告
rstanのWindowsへのインストールは大きな関門であることがわかった。
行ったこと
いくつかポイントがあって、以下を行ってみた。
rstanの再インストール
remove.packages("rstan")
StanHeadersのversion
rstanと、StanHeadersのversionを揃える。 discourse.mc-stan.org
Rtools
35を用いていたが、どうにもうまく行かない。
何度か入れ直した。
蛇退治
Anacondaをアンインストールした。
見事に厄災は消え去った。
なんだこちらを参照していたのか。
晴れてWindowsでもStanし放題である。
レセプトデータから医療の実態に迫る
レセプト
レセプト - Wikipedia は、医療報酬の明細書です。
まさにレシート、のことですね。
診療報酬の請求のためにつける病名があり、 実態を反映していないとの意見もありますが、 診療行為や病名と組み合わせて使えば、 診療実態を明らかにすることが可能です。
このページでは、レセプトデータを見るときに必要な情報を随時更新していきます。
その公開
NDBオープンデータ として、国のHPで公開もされています。
薬品や診療行為マスターデータ
診療報酬情報提供サービス で、薬品名、コード、薬価などの情報が公開されています。
病名
社会保険表章用疾病分類.
https://www.mhlw.go.jp/stf/seisakunitsuite/bunya/kenkou_iryou/iryouhoken/database/hokensippei.html
厚労省通知
International shortlist for hospital morbidity tabulation (ISHMT).
https://www.who.int/classifications/icd/implementation/morbidity/ishmt/en/ http://apps.who.int/classifications/apps/icd/implementation/hospitaldischarge.htm
レセプトの落とし穴
レセプトを利用する際に注意すべきポイントについて、以下の文献が参考になります。
ベイズ統計分析によるデータモデリング入門 [ベイズ統計入門の決定版]
COI開示：特になし。
統計分析の教科書、データ分析の技術書、どんどん充実してきています。
私の答えはNoでした。理由は、StanやBUGSといったベイズモデリングのためのコードを扱うハードルがあったから。
今回、馬場真哉さんによるベイズ統計分析によるデータモデリング入門を読んでとてもためになったので記しておきます。
という構成になっています。
1章は理論の解説。2章はRやStanでのプログラミングの基本事項が解説してあり、初学者にもフレンドリーな構造です。
実装だけ試したければ3-5章だけでOK。
今回、一番の収穫は、
brms (Bayesian Regression Models using Stan) でした。
何がすごいのか？
ほとんど普通のglm関数と同じようにベイズ統計モデリングが出来ます。
これほど簡単になっていたのかと驚きました。
時系列分析については時系列分析と状態空間モデルの基礎で勉強させてもらいました。こちらもおすすめ。
馬場さんのwebsiteはこちらです↓↓
Rで表記ゆれを癒やす：全角から半角へ
日本語データを扱うときに困ること
日本語データを分析するときにもっとも困ることは何でしょうか？
それは、表記ゆれです。
送り仮名、漢字かな混じり、全角半角混じり、大文字小文字混じり・・・色々ありますが、
今回は、特に、全角と半角が混ざった表記ゆれについて考えていきます。
日本語には、全角と半角の文字が存在します。歴史を紐解くとわかるもので、
ウィキペディア にも詳しいのですが、東アジアの印刷物は、もともと文字が正方形に収まるように作られていました。
コンピュータが出現し、文字が扱えるようになり、はじめは、すべての文字が同じ幅の文字として用いられていました。
その後、東アジアの文字は、1文字を表すのに2バイト以上を用いる「マルチバイト文字」によって利用が可能となりました。
これ自体は、アルファベットが、漢字などよりも細く印刷される伝統によく合っています。
しかし、解析を行う上で、、Ｄr、Dr.、など、全角アルファベットや全角数字でタイプされた日本語情報を扱う場合には、半角アルファベットや半角数字と同じと認識してくれません。
このような表記ゆれに対する対処は、以下のブログにとてもまとまっていて、参考にさせてもらいました。
Rで全角英数記号を半角にするのに{stringi}が役に立った話 - Note of Pediatric Surgery
Nippon package
まず、全角半角ごちゃごちゃのデータを準備します。
data <- c("ＢＯＹＯＮ", "ｂｏｙｏｎ", "ボヨーン")
library(Nippon) zen2han(data)
まずは、もともとpackageとして入っていたNipponを用いてみます。
[1] "BOYONboyonボヨーン"
zen2hanではうまくいかない・・・
すべて連続する文字列として処理されてしまい、データフレーム全体を調整するには適していない模様です。
三重大学谷村先生のスライドを見ていると・・・
library(Nippon) sanitizeZenkaku(data)
[1] "BOYON" "boyon" "ボヨ-ン"
sanitizeZenkakuでいけました！！
ただ、これだと
伸ばし棒までsanitizeされてハイフンになってしまいます。。
halfwidthen
Rcppで全角英数を半角英数に変換するパッケージをつくった - Technically, technophobic.
こちらを参考にして、githubからpackageをインストールしまして・・・
library(devtools)
install_github("yutannihilation/halfwidthr")
library(halfwidthr)
halfwidthen(data)
[1] "BOYON" "boyon" "ボヨーン"
これでうまくいきました！
いまのところこちらが使い勝手がよいので、
当面はhalfwidthenを使用させていただきます。
stringrには実装されないかな・・・（他力本願）
Rで可視化する：可視化の方法とは
手持ちのデータを可視化する
データの特徴を示すためには、グラフや図にすることが重要である。 いくつも可視化があるのだが、実際にそれぞれのグラフを作成するには、
- データの性質を知る
- 可視化を実装する
- (実装の仕方を調べるために可視化の名称を知る)
が必要である。 Rでの可視化の実装について、以下のサイト（英語）にまとまっている。
データの性質
以下を押さえておくと良い。
- 連続変数
- カテゴリー変数
- 連続変数&カテゴリー変数
- 地図情報
- ネットワーク
- 時系列
可視化の実装
もちろんgoogle や書籍で調べるのだが、可視化方法の名称がわからないと調べ辛く、新たな手法に取り組む場合ここで躓く場合も多い。
可視化の名称
ここに具体例とともにまとまっている。
