[WIP]種間差や測定誤差を考慮した多樹種のアロメトリーのパラメータ推定をRStanで階層モデル化（伊東 2015）の追試
Shinya Uryu
November 08, 2015
⚠️ そもそも収束していなかったり、 欠損値の扱い、もしくはデータの渡し方が異なる可能性があり、元の論文とは異なる結果になってしまっている… 💣
森林総研の伊東さんが書かれた論文がネット上に公開された。
伊東宏樹 (2015) 樹種間差および測定誤差を考慮した胸高直径ー樹高関係のベイズ推定. 森林総合研究所研究報告 14(2) 73–74.
親切なことにStanのコードだけでなく利用したデータセットまで資料として公開してくれているので、RStanによる階層ベイズモデルの勉強がてら追試してみることにした。
🔧 準備
Rパッケージの読み込みと並列化設定
library(zoo)
library(VIM)
library(tidyr)
library(dplyr)
library(rstan)
library(ggmcmc)rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())データの取得
補足電子資料をダウンロードし、R上で利用できるデータセットにする。
base_url <- "http://www.ffpri.affrc.go.jp/pubs/bulletin/435/documents/"
data_file <- "tables-s1-s3.xlsx"
stan_file <- "list-s1.txt"
download.file(url = paste0(base_url, data_file),
destfile = paste(getwd(), "ito2015stan", data_file, sep = "/"))
download.file(url = paste0(base_url, stan_file),
destfile = paste(getwd(), "ito2015stan", "ito2015.stan", sep = "/"))🔨 RStanによるMCMCサンプリング
まず、公開されている計測データを読み込んで全体の概要を把握する。
df_ito <- readxl::read_excel(path = "ito2015stan/tables-s1-s3.xlsx",
sheet = 1,
skip = 2,
col_names = c("jp.name", "sp.name", "N", "D",
"Nh1", "Nh2", "Nh3", "Nh4"))
glimpse(df_ito)## Observations: 232
## Variables: 8
## $ jp.name (chr) "コシアブラ", "コシアブラ", "コシアブラ", "コシアブラ", "コシアブラ", "コシアブラ", ...
## $ sp.name (chr) "Chengiopanax sciadophylloides", "Chengiopanax sciadop...
## $ N (dbl) 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1, ...
## $ D (dbl) 27.60, 30.94, 9.07, 17.06, 1.10, 20.53, 23.14, 8.79, 1...
## $ Nh1 (dbl) 13.26, 17.43, 11.13, 11.92, 3.37, 12.44, 13.81, 8.80, ...
## $ Nh2 (dbl) 13.99, 17.12, 11.12, 11.46, 3.48, 11.69, 13.77, 8.99, ...
## $ Nh3 (dbl) 13.43, 17.27, NA, 11.35, NA, 12.01, NA, NA, NA, NA, NA...
## $ Nh4 (dbl) NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...Nh1からNh4というのが、最大で４回記録した樹高の値であるらしい。いくつかの欠損値が確認できる。
データ処理
読み込んだ元データに対して、いくつかの処理を加える。
Stanのコードをみると、樹高の値は H[N, maxNh]として与えられている。Nは幹の本数でNhは観測回数を指している。これの行列をRStanのMCMC実行関数に渡せば良いはず…。そのための処理として以下を実行する。
測定回数の列を作成
測定誤差をパラメータに取り組むために必要な各幹の測定回数の変数を用意する。Nh1からNh4までの列で欠損している箇所を数えて、もし欠損している場合には計測回数（４）から引いた値となる。
df_ito %<>% rowwise() %>%
dplyr::mutate(Nh = countNA(c(Nh1, Nh2, Nh3, Nh4))) %>%
dplyr::mutate(Nh = ifelse(Nh == 0, 4, 4 - Nh)) %>%
ungroup()
# ungroupを忘れないこと欠損値の補完
測定誤差をなくすために同一の幹において樹高は最大で４回計測されている。一方で計測が２回あるいは３回しか行われていない幹もあり、データに欠損が見られる。
Stanでは、欠損値を含んだ値があると計算がうまくできないので、欠損値の対処が必要である。論文のほうに欠損データをどう扱ったかという話がないので、ここでは{zoo}パッケージのna.approx()によって欠損値の補完を試してみた。
よくわかっていないのだけど、欠損箇所を０で置換したり、観測した幹の最大値を代入したりすると、うまく結果が一致しなかった。これはなぜなのか。後で調べる…。
# 0で置換する -> 推定値がおかしいかも
# df_ito %<>% replace_na(list(Nh2 = 0, Nh3 = 0, Nh4 = 0))
# 最大値を適用する -> 推定値がおかしいかも
# df_ito %<>% rowwise() %>%
# dplyr::mutate(obs.hmax = max(Nh1, Nh2, Nh3, Nh3, na.rm = TRUE),
# Nh2 = ifelse(is.na(Nh2), obs.hmax, Nh2),
# Nh3 = ifelse(is.na(Nh3), obs.hmax, Nh3),
# Nh4 = ifelse(is.na(Nh4), obs.hmax, Nh4)) %>%
# ungroup()
df_ito %<>% dplyr::mutate_each(funs(na.approx), Nh2, Nh3) %>%
dplyr::mutate(Nh4 = Nh3 + 0.8)データの構造（再）とリストへの格納
処理を加えたデータは次のような形となった。これをRStanに渡すため、Stanコードのdataブロックで宣言されている変数についてリストオブジェクトに格納する。
glimpse(df_ito)## Observations: 232
## Variables: 9
## $ jp.name (chr) "コシアブラ", "コシアブラ", "コシアブラ", "コシアブラ", "コシアブラ", "コシアブラ", ...
## $ sp.name (chr) "Chengiopanax sciadophylloides", "Chengiopanax sciadop...
## $ N (dbl) 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1, ...
## $ D (dbl) 27.60, 30.94, 9.07, 17.06, 1.10, 20.53, 23.14, 8.79, 1...
## $ Nh1 (dbl) 13.26, 17.43, 11.13, 11.92, 3.37, 12.44, 13.81, 8.80, ...
## $ Nh2 (dbl) 13.99, 17.12, 11.12, 11.46, 3.48, 11.69, 13.77, 8.99, ...
## $ Nh3 (dbl) 13.43, 17.27, 14.31, 11.35, 11.68, 12.01, 11.67, 11.32...
## $ Nh4 (dbl) 14.23, 18.07, 15.11, 12.15, 12.48, 12.81, 12.47, 12.12...
## $ Nh (dbl) 3, 3, 2, 3, 2, 3, 2, 2, 2, 2, 2, 2, 3, 2, 3, 3, 3, 3, ...list_d <- list(N = nrow(df_ito),
D = df_ito$D,
Nh = df_ito$Nh,
maxNh = max(df_ito$Nh),
H = df_ito %$% c(Nh1, Nh2, Nh3, Nh4) %>% matrix(nrow = nrow(df_ito)),
Ns = df_ito$jp.name %>% unique() %>% length(),
S = df_ito$sp.name %>% as.factor() %>% as.numeric())MCMCサンプリング
乱数を固定して、論文と同じ試行回数、連鎖数でsampling()を実行する。
comp_stan_code <- stan_model(file = "ito2015stan/ito2015.stan",
model_name = "ito2015_dh_allometry")
# MCMCの試行回数などや連鎖の数は伊東(2015)に従う（乱数は固定）
fit_sample <- sampling(object = comp_stan_code,
data = list_d,
seed = 71,
iter = 35000,
warmup = 5000,
chain = 4,
sample_file = "ito2015stan/allometry_")以下、MCMCサンプリングの結果を収めたfit_sampleというオブジェクト（stanfitクラス）を対象にして解析を進めていく。
🔬 推定結果の確認
収束判定
操作を楽にするため与えられた事後分布の値をデータフレームに変換しておき、その後、パラメータのRhat、MCMCの連鎖軌道を確認する。
df_param <- ggs(fit_sample)
# 推定したパラメータをデータテーブルにする
df_param %>% dplyr::filter(grepl("^D_bar.3.", Parameter)) %>%
ggs_Rhat() + xlab("R_hat")
# 一部の幹のRhatを確認あれ？きちんと収束していない…。特に33番目の幹が大きく1.0から外れている。
というわけで修正を加えた元データを確認してみる。
df_ito[33, ]## Source: local data frame [1 x 8]
##
## jp.name sp.name N D Nh1 Nh2 Nh3 Nh4
## (chr) (chr) (dbl) (dbl) (dbl) (dbl) (dbl) (dbl)
## 1 サカキ Cleyera japonica 3 1.78 3.31 3.4 3.26 NAやはり欠損していたのが問題になっているのだろうか？
連鎖の軌道はというと
df_param %>% dplyr::filter(grepl("^D_bar.3.", Parameter)) %>%
ggs_traceplot()
うーむ…
アロメトリー係数の事後分布
推定したこれらのパラメータ、特にアロメトリー式で用いているa, b, H_maxの３つについて見ていく。階層ベイズ化によって、これらのパラメータは樹種ごとに推定されている。
fit_sample@model_pars## [1] "D_bar" "log_H_max_bar" "log_a_bar" "log_h_bar"
## [5] "e_a" "e_h" "e_H_max" "e"
## [9] "sigma" "a" "h" "H_max"
## [13] "log_H_bar" "lp__"各パラメータの事後分布から信用区間を算出したものをデータフレームに変換する。
df_ci <- ci(df_param)資料として、「各樹種のパラメータa, b, Hmaxの事後分布の平均および分位点」があるのでそれと見比べてみる。
# 種ごとに共通なパラメータのみ抽出する
df_ci_sp <- df_ci %>% dplyr::filter(grepl("^H_max|^a|^b", Parameter))
# dplyr::mutate(sp.id = gsub("\\.", "\\_", Parameter) %>% extract_numeric())
df_ci_sp %>% dplyr::filter(grepl("^H_max", Parameter)) %>%
dplyr::arrange(median) %>%
DT::datatable(options = list(pageLength = 15))推定されたH_maxで、もっとも小さいのは10番目の種であり、その値は12.75となった。もっとも大きいのは４番目の種であり、その値は30.195であった。これがどの種かというと
# 種番号の確認
df_ito$jp.name %>% unique() %>% {
print(.)
print(.[10])
print(.[4])
}## [1] "コシアブラ" "リョウブ" "サカキ" "スギ" "ヒサカキ"
## [6] "タカノツメ" "アオハダ" "ソヨゴ" "ネズミモチ" "ネジキ"
## [11] "カナメモチ" "アラカシ" "コナラ" "エゴノキ" "クロバイ"
## [1] "ネジキ"
## [1] "スギ"ネジキとスギ。Hmaxの傾向は論文の結果と一致した。うまくいっているのかいっていないのかよくわからない。どこかが間違っているのは確かなのだが…。勉強が足りない。