Rでシェープファイルを図示する（Mac）

Rでシェープファイルを図示する方法はいろいろなところで紹介されていますが、
Mac環境のRでシェープファイルが図示できない問題があります。

以前GADM(https://gadm.org/)から入手したシェープファイルを図示しようとしたのですが、うまくいきませんでした。 身近な知り合いのMacユーザーの方も同じ症状でしたが、Windowsだと問題なく動くようです。

結局WindowsでうまくいってMacでうまくいかない原因はよくわからなかったのですが、 誰かの助けになればと思い解決策を記事に残します。

環境はMacOS 10.16, R version 4.0.2 です。

# 必要なパッケージの読み込み
library(raster)
library(sf)
library(ggplot2)
library(magrittr)
library(rmapshaper)

GADMから日本の白地図をダウンロードします。

map <- raster::getData("GADM", country="JP", level=0)

sfオブジェクトに変更します。

map <- st_as_sf(map)

これをそのまま

ggplot(map)+
geom_sf()

で図示しようとすると、少なくとも私の環境ではPCがフリーズします（一度根気よく20分くらい待ったら図示してくれました）。

で、調べても解決法が出てこず困っていたのですが、この前たまたま対処法がわかりました。 どうやらシェープファイルの情報量が多すぎるのが問題らしいので、{rmapshaper}というパッケージで情報量を落とします。

keep=1で元のままです。 以下ではベクターデータの頂点の数を0.01％まで削減しているのだと認識しています。

jpn <- rmapshaper::ms_simplify(map, keep=0.0001, keep_shapes=TRUE)

これを図示すると...

ggplot(jpn)+
  geom_sf()

頂点が減って単純化されてるのがわかります。 0.01％の情報量（？）にしてもギリ原型を保っています。

バイト数を確認すると以下のようになります。

object.size(map)     # オリジナル

## 9897520 bytes

object.size(jpn)  　 # 頂点削除後

## 8048 bytes

keep=0.01にすると、少なくとも日本地図のまま使う分には問題ないほどになります。

rmapshaper::ms_simplify(map, keep=0.01, keep_shapes=TRUE) %>% 
  ggplot() + geom_sf()

{sf}パッケージのst_simplify関数でも同様のことができるけど、アルゴリズムが違くて{rmapshaper}の方が良いということです。

参考：https://www.r-bloggers.com/2021/03/simplifying-geospatial-features-in-r-with-sf-and-rmapshaper/

ggplot2：別々のデータセットから描画した図に凡例を追加する

Rで別々のデータセット、あるいはワイドフォーマットからggplot2で作図する話をします。 （どうってことなくない？）と思う方も多いかもしれませんが、凡例の付け方にかなり悩んでしまったので、記事にします。

Rではロングフォーマットが求められる

ということを、初めてRを習った時には刷り込まれました。

特にggplot2での描画時はそうで、必要に応じてデータの連結やpivot_longerによるフォーマット変換などを行っていました。

しかしこれは結構めんどくさいです。

別にロングフォーマットである必要はなく、別々のデータでも、ワイドフォーマットでも、
不自由なく描画できるということに気付きました。

たとえば以下のような別々のデータがあります。

America <- data.frame(x=seq(-1, 1, length=10), y=0)
Japan <- data.frame(x=seq(-1, 1, length=10), y=seq(-1, 1, length=10))

別々のデータセットのまま、色分けして図示しようとすると、以下のような案が思いつきます。

library(ggplot2)

ggplot()+
  geom_line(data=America, aes(x=x, y=y), color="green", size=3)+
  geom_line(data=Japan, aes(x=x, y=y), color="blue", size=3)

緑がAmericaで青がJapanなのですが、この方法では凡例を追加することができません（多分）。

以下のようにすると解決します。

ggplot()+
  geom_line(data=America, aes(x=x, y=y, color="America"), size=3)+
  geom_line(data=Japan, aes(x=x, y=y, color="Japan"), size=3)+
  scale_color_manual(values = c("America" = "green",
                                "Japan" = "blue"))

凡例の順番の変更は、以下のように工夫すると簡単に行えます。

# factor型の名前を"1","2"にして順番を調節し、後からlabelを変える
ggplot()+
  geom_line(data=America, aes(x=x, y=y, color="2"), size=3)+
  geom_line(data=Japan, aes(x=x, y=y, color="1"), size=3)+
  scale_color_manual(labels = c("1" = "Japan",
                                "2" = "America"),
                     values = c("1" = "blue",
                                "2" = "green"))

別々のデータでも、ワイドフォーマットでもドンと来いになります。

これでもう図示のためだけにデータ整形で悩まなくていい！！…かもしれない

カブトムシの在不在と見つけやすさ

こちらはQiitaのStanアドベントカレンダー2020 第19日目用の記事になります。
近くの山にカブトムシ調査に行った状況を想定します。

知りたいことは、

・どのくらいカブトムシがいるのか
・カブトムシの見つけやすさはどのくらいなのか
・見つけやすい時間帯はあるのか

です。

調査デザインとしては、山を適当なサイズのプロットに区切り、それぞれのプロットでカブトムシの在不在状況を記録します。

夜の山は危険なので、カブトムシ取りは朝〜夕方の間に行う予定で、事前調査もその時間帯で行います。

また、カブトムシの存在量を数量化することも可能ですが、ここでは簡単のため在・不在情報として記録することとします。

ここで、この調査において最も重要な考え方について、いくつか記述します。

1.  カブトムシの検出は「プロットiにカブトムシが存在する確率psi(i)とカブトムシの検出率p(i)」の積から決定される

　つまり、「カブトムシが検出された」という事実は、疑いもなくそこにカブトムシが存在することを意味するが、「カブトムシが検出されなかった」ことは「カブトムシが存在しなかった」状態と「存在したが見つけられなかった」状態の2つのパターンが考えられる、ということです。

　この対応として、同じプロットで複数回の調査を繰り返します。一回の「不検出」からは「不在」と「在→不検出」の区別ができませんが、例えば調査を3回繰り返し「不検出・不検出・検出」となれば、「存在はしたが1回目と2回目は見つけられなかった」という事実が浮かび上がり、存在率と検出率の分離が可能になります。

2.  各プロット内での閉鎖個体群を想定する 

　これは1の考え方を可能にする仮定になります。閉鎖個体群とは、「個体の誕生・死亡・移出・移入」が発生しない、ということです。  

　現実にこの仮定が満たされることはないので、実際上は「閉鎖個体群として近似できるほどの期間の間に調査する」ことが必要となります。幸い、甥っ子とのカブトムシ取りと今回の事前調査は日中に行うため、夜行性のカブトムシは１日の昼間は閉鎖個体群として近似できるでしょう。

まずは、もっともシンプルな状況でシミュレーションしましょう（といってもずっとシミュレーションですが）。

これは調査時間は考慮せず、一定存在率・一定検出率を仮定します。

 擬似データを生成します。

library(rstan)
set.seed(123)

R <- 100    # 調査サイト数
T <- 5      # 時間的調査反復数
psi <- 0.6  # 占有率（存在率）
p <- 0.4    # 検出率

z <- rbinom(R, 1, psi)   # R地点の真の在不在情報

y <- array(NA, dim=c(R,T))  # 調査記録枠組み

for (i in 1:R){                       # 調査記録（検出/不検出）格納
  for (t in 1:T){
    y[i,t] <- rbinom(1, 1, z[i] * p)  # 検出成功率は z(在不在) * p(検出率)
  }
}

調査から得られるデータは以下のようになります。
行がそれぞれのサイト、列が反復です（6行目まで表示しています）。
100プロットそれぞれについて5反復の調査を実施し、検出ならば1、不検出ならば0となります。

さて、存在率をpsi[i]、検出率をp[i,t]、プロットiの真の存在状態をz[i] (二値)、調査結果の観測値をy[i,t] (二値)としてモデルを整理すると、以下のようになります。
ただし、ここではpsiとpは一定なので、添字は必要ありません。

z[i] ~ Bernoulli(psi[i])

y[i,t] | z[i] ~ Bernoulli(z[i] * p[i,t])

Stanでは離散パラメータは扱えないので、以下のようなコードとなります。

data{
  int R;                         // プロット数
  int T;                         // 反復数
  int<lower=0, upper=1> y[R,T];  // 観測データ
}

transformed data{
  int sum_y[R];  // プロットiの合計検出回数
  
  for (i in 1:R)
    sum_y[i] = sum(y[i]);
}

parameters{
  real<lower=0, upper=1> psi;  // 存在率
  real<lower=0, upper=1> p;    // 検出率
}

model{
  for (i in 1:R){
    if (sum_y[i] > 0)  // 一度でも検出された場合（存在）
      target += bernoulli_lpmf(1 | psi) 
              + bernoulli_lpmf(y[i,] | p);
    else               // 一度も検出されなかった（不在 or 存在するが不検出）
      target += log_sum_exp(bernoulli_lpmf(0 | psi), 
                            bernoulli_lpmf(1 | psi) 
                            + bernoulli_lpmf(y[i,] | p) 
                            );
  }
}

キックするコードは以下で、結果を図示します。

data <- list(R=nrow(y), T=ncol(y),  y=y)
fit <- stan("beetle0.stan", data=data, seed=123)

真値はpsi=0.6, p=0.4ですので、まぁいい感じと思われます。

先の例はシンプルなシミュレーションなので、もう少し複雑にします。
具体的には、

・カブトムシの見つけやすい時間帯を推定したい
・複数人に調査を手伝ったもらったので、調査者ごとの検出率差も考慮したい

ということにします。

プロット数100、反復5、調査者4人で手分け、として生成された擬似データは以下のようになります。
なお、データ生成コードや図示コードが長くなってしまったので、ページの最後に載せます。

time_y1列とPerson_y1列が、y1列に対応する観測時刻と観察者、という見方になります。
timeについては、午前7時からの経過時間を分単位で記録しています（0 < time < 600）

また、観測時間及び4人の観察者による真の検出率pは以下のようになっています

12時前後が最も見つけづらい、という風になっています。

基本的には、先のシンプルな例をアレンジするだけです。
Stanコードは以下になります。

data{
  int R;                             // プロット数
  int T;                             // 反復数
  matrix[R,T] Time;                  // 観測時刻
  int<lower=1, upper=4> Person[R,T]; // 観測者ID(4人)
  int<lower=0, upper=1> y[R,T];      // 観測データ
  int T_new;                         // 時刻と検出率図示用
  vector[T_new] Time_new;            // 同上
}

transformed data{
  int sum_y[R];              // 検出回数合計
  matrix[R,T] Time2;         // Timeの2乗
  vector[T_new] Time2_new;   // Time_newの2乗(図示用)
  
  for (i in 1:R)
    sum_y[i] = sum(y[i]);
    
  Time2 = Time .* Time;
  Time2_new = Time_new .* Time_new;
}

parameters{
  real<lower=0, upper=1> psi;  // 存在率
  real beta0;
  real beta1;
  real beta2;
  vector[4] eps;               // 観測者ごと検出率誤差
  real<lower=0> sigma;         // その標準偏差
}

transformed parameters{
  matrix[R,T] logit_p;         // ロジット軸での検出率p
  
  for(i in 1:R)
    for(t in 1:T)
      logit_p[i,t] = beta0 + beta1 * Time[i,t]
                   + beta2 * Time2[i,t] + eps[Person[i,t]]; 
}

model{
  sigma ~ normal(0, 1);   // 事前分布
  
  eps ~ normal(0, sigma); // 観測者ごと誤差
  
  for (i in 1:R){
    if (sum_y[i] > 0)
      target += bernoulli_lpmf(1 | psi) 
              + bernoulli_logit_lpmf(y[i,] | logit_p[i,]);
    else
      target += log_sum_exp(bernoulli_lpmf(0 | psi),
                            bernoulli_lpmf(1 | psi)
                            + bernoulli_logit_lpmf(y[i,] | logit_p[i,])
                            );
  }
}

generated quantities{
  vector[T_new] p_pred;   // 検出率と時間の効果
  
  p_pred = inv_logit(beta0 + beta1 * Time_new + beta2 * Time2_new);
}

調査者4人には、事前にカブトムシの見つけ方講座を行っていたものとし、調査者ごとの検出率誤差に関わるsigmaにNormal(0,1)と幅の小さめの事前分布を適用しています。 見つけやすい、あるいは見つけにくい時刻を探すために、Timeの2乗を用いて回帰しています。
また、観測時刻と検出率の関係を知りたいので、generated quantitiesブロックで計算します。

推定結果は以下のようになります。
観察しやすい、あるいは観察しづらい時間については、Timeの2乗にかかる係数beta2がほぼ正の値を取ることから、ある時間で最も観察しづらい時間が存在することが分かります。これも図示してみましょう。

真値は
存在率psi = 0.6, beta0 = -1, beta1 = -0.3, beta2 = 0.5, eps[1] = 0.6, eps[2] = 0.3, eps[3] = -0.3, eps[4] = -0.6
と設定しており、悪くはなさそうです。

時間と検出率の関係も良さそうです。
ちゃんと12時前後でもっとも見つけづらくなっています。
時刻と検出率の図示については、下に凸であることを強調するため、MCMCサンプル300個を図示する方法を取りました。

このモデルはいろいろな拡張が可能です。
今回は閉鎖個体群を想定していますが、現実的には開放個体群を想定したいところでしょう。
その場合、開放個体群を想定する調査期間と、閉鎖個体群を想定する調査期間を組み合わせるロバストデザインが用いられるようです。
これにより、個体の死亡率や移入率なども推定できる場合があります。

今回の記事に当たっては以下の書籍を参考にさせていただきました。

BUGSで学ぶ階層モデリング入門: 個体群のベイズ解析

今回紹介した内容とその拡張も含めて、生態学（特に個体数や確率論的個体群動態）での様々な事例を、主に状態空間モデルから解析した内容になっています。
Winbugsを用いた解析ですが、全ての例についてStanでの実装がGithub上に挙げられています。
何より、日本語訳が素晴らしいです。考え方が分かりやすく、的確な統計学的用語で解説してくれていると思いました。

不満があるとすれば、Stanでは離散パラメータが簡単に扱えないことから、Winbugsだとこんなに簡単なコードなのに、この複雑なStanコードは一体...となることが度々あったことでしょうか。
まぁそれは仕方ないですね...一度理解できればなんとかなります。
基本的にはStanの方が高速なので、その引き換えと思う他ありません。

では最後に、擬似データ生成、Stanのキック、図示に関わるRコードをバーっと載せときます。
もっと簡単にできそうな気がしますが、技術力不足ゆえ長くなっている気がします。

library(rstan)
library(tidyverse)
set.seed(123)

#### 擬似データ生成 ####
R <- 100     # 調査サイト数
T <- 5       # 時間的調査反復数
psi <- 0.6   # 占有率（存在率）
Person_logit_p <- c(0.6, 0.3, -0.3, -0.6) # 調査者による検出率誤差（ロジット軸）

time <- array(NA, dim=c(R, T))  # 観測時間
Time <- array(NA, dim=c(R, T))  # 標準化された観測時間
Person <- array(NA, dim=c(R,T)) # 調査者ID
p <- array(NA, dim=c(R,T))      # 検出率
y <- array(NA, dim=c(R,T))      # 観測値


for(i in 1:R){
  for(t in 1:T){
    time[i,t] <- sample(0:600, 1) # 午前7時からの経過時間（分）
    Person[i,t] <- sample(1:4, 1) # 調査者ID (1~4)
  }
}

for(i in 1:R){
  for(t in 1:T){
    Time[i,t] <- (time[i,t] - mean(time)) / sd(time)   # 観測時間を標準化
    p[i,t] <- plogis( -1 - 0.3 * Time[i,t] + 0.5 * Time[i,t] * Time[i,t] 
                      + Person_logit_p[Person[i,t]])   # 観測時間と観察者による検出率
  }
}

z <- rbinom(R, 1, psi) # 真の在不在情報

for(i in 1:R){         # 観測データ生成
  for(t in 1:T){
    y[i,t] <- rbinom(1, 1, z[i] * p[i,t]) # 観測成功率 = 在不在(1/0) * 検出率p
  }
}

##### 検出率と時刻、調査者の関係図示 ####
P <- array(NA, dim=c(3*R, 2))
for(i in 1:R){
  P[i,1] <- time[i,1]
  P[i,2] <- p[i,1]
  P[i+R, 1] <- time[i,2]
  P[i+R, 2] <- p[i, 2]
  P[i+2*R, 1] <- time[i,3]
  P[i+2*R, 2] <- p[i,3]
}
P <- as.data.frame(P)
colnames(P) <- c("time", "P")

ggplot(P, aes(x=time, y=P))+
  geom_point()+
  ylim(0,1)


#### 観測データをまとめる ####
d <- cbind(y, time, Person)   # 調査で得られるデータ
colnames(d) <- c("y1", "y2", "y3", "y4", "y5",
                 "time_y1", "time_y2", "time_y3", "time_y4", "time_y5",
                 "Person_y1", "Person_y2", "Person_y3", "Person_y4", "Person_y5")


#### 推定 ####
time_new <- seq(0, 600, 1)  # generated quantitiesブロック用
Time_new <- (time_new - mean(time)) / sd(time) # 標準化
T_new <- length(Time_new)   # generated quantitiesブロック用

data <- list(R=nrow(y), T=ncol(y), K=4, Time=Time, Person=Person, y=y,
             Time_new=Time_new, T_new=T_new)

fit <- stan("beetle.stan", data=data, seed=123, 
            control=list(adapt_delta=0.995))
fit

#### 検出率と時刻の関係を図示 ####
ms <- rstan::extract(fit)
p_pred <- ms$p_pred

# 時刻ごと中央値を算出
p50 <- apply(p_pred, 2, quantile, probs=0.5) 
P50 <- data.frame(Time=time_new, p50)

# MCMCサンプル300個を取り出す
p_p <- t(p_pred)
P_p <- as.data.frame(p_p[,1:300])
P_p$Time <- time_new
P_t <- pivot_longer(P_p, cols=-Time, names_to="MCMC", names_prefix="V", values_to="p")

# 図示
ggplot()+
  theme_classic()+
  geom_line(data=P50, aes(x=Time, y=p50), colour="blue", size=1.6)+
  geom_path(data=P_t, aes(x=Time, y=p, group=MCMC), alpha=0.3, size=0.2)+
  ylab("Detection rate")

Stanで「丸め」の実装

Stan reference guideにも記載がある「丸め」について検証します。
丸められたデータとは、本当は連続値なのに、四捨五入（あるいは切り捨て）などで離散値のようになっているデータを指すようです。
打ち切りと似ている気がしますが、レファレンスでは別項目として記載されています。

library(rstan)
set.seed(123)

Y <- rnorm(50, 8, 2) # 平均8, 標準偏差2の正規分布
Y <- round(Y, 0) # 四捨五入

こんなデータの平均と標準偏差を求めます。

Stan reference guideの1つ目のやり方

data{
  int N;
  int Y[N];
}
parameters{
  real mu;
  real<lower=0> sigma;
}
model{
  for(n in 1:N)
    target += log(
      Phi((Y[n] + 0.5 - mu) / sigma) - Phi((Y[n] - 0.5 - mu) / sigma)
      );
}

データY[n]の"真の値"は Y[n]-0.5 ~ Y[n]+0.5 の区間にあるので、その確率を求めて、対数尤度を計算しています。
正規分布の標準化（平均0標準偏差1に変換）を行いPhiで標準正規累積分布関数として実装します。

modelブロックは以下のように書くこともできます。

model{
  for(n in 1:N)
    target += log(normal_cdf(Y[n]+0.5, mu, sigma) - normal_cdf(Y[n]-0.5, mu, sigma));
}

行っていることは同じです。
(Y[n]+0.5 | mu, sigma)と書くとエラーが出ます。~cdfはコンマでつなげるんですね。

Stan reference guideの2つ目のやり方

整数に丸め込まれたデータYに誤差項Y_err(-0.5 ~ 0.5) をパラメータとして宣言し、潜在変数 z[n] につなげます。

data{
  int N;
  int Y[N];
}
parameters{
  real mu;
  real<lower=0> sigma;
  vector<lower=-0.5, upper=0.5>[N] Y_err;
}
transformed parameters{
  vector[N] z;
  for(n in 1:N)
    z[n] = Y[n] + Y_err[n];
}
model{
  z ~ normal(mu, sigma);
}

キックするRコードと、事後分布を確認します（ほぼ同じ推定をしてくれるので方法2の結果だけ載せます）。
また方法2では、 丸め誤差 Y_err も推定してくれるわけなので、分布の「端」...ここではYの最大となるY[k]と付随する誤差項Y_err[k]の事後分布も確認してみましょう。
Y[k]は真の平均値から離れているので、Y_err[k]の事後確率密度は -0.5周辺 → 0.5周辺 に向かうにつれて小さくなることが予想されます。

data <- list(N=length(Y), Y=Y)
fit <- stan("marume1.stan", data=data, seed=123)
fit2 <- stan("marume2.stan", data=data, seed=123)

max <- which(Y==max(Y)) # Yが最大値をとる行を調べる
stan_hist(fit2, pars=c("mu", "sigma", "z[16]", "Y_err[16]")) # Y[16]でYが最大値を取る

zは潜在変数です。
パラメータは良い感じに推定されているようです。
真の平均が8なのにY=12になるのは、丸められる前の潜在変数zが11.5に近い方が確率が高いよね...となっており、こちらも予想された結果となりました。