Estimation of diameter-height curves

D-H.R

# 『たのしいベイズモデリング2』第5章のRコード
# CmdStanR版
# R 4.3.1 Mac, Stan 2.23.2, CmdStanR 0.5.3 で動作確認

library(readxl)
library(dplyr)
library(cmdstanr)
options(mc.cores = parallel::detectCores())

# データファイルのURL
data_url <- "http://www.ffpri.affrc.go.jp/pubs/bulletin/435/documents/tables-s1-s3.xlsx"

# データの取得と読み込み
destfile <- "tables-s1-s3.xlsx"
if (!file.exists(destfile)) { # データファイルがないならダウンロード
  download.file(data_url, destfile = destfile)
}
data <- read_xlsx(destfile, sheet = 1, skip = 2,
                  col_names = c("Jpn_name", "Sci_name", "No", 
                                "D", "H1", "H2", "H3", "H4"))

# データの表示
print(data)

# 集計表の作成
data %>%
  dplyr::rowwise() %>%
  dplyr::mutate(H = mean(c(H1, H2, H3, H4), na.rm = TRUE)) %>%
  dplyr::ungroup() %>%
  dplyr::group_by(Jpn_name) %>%
  dplyr::add_count() %>%
  dplyr::summarize(N = mean(n), mean_D = mean(D), mean_H = mean(H))

# Stan用に樹高データの処理
heights <- data[, c("H1", "H2", "H3", "H4")]
heights[is.na(heights)] <- 0  # NAを0に
nh <- sapply(1:nrow(heights), function(i) sum(heights[i, ] > 0))  # 測定回数

# 樹種をコード化
spc <- as.numeric(factor(data$Jpn_name))

# Stanに渡すデータ
stan_data <- list(N = nrow(data),
                  D = data$D,
                  Nh = nh,
                  MaxNh = max(nh),
                  H = heights[, 1:max(nh)],
                  Ns = max(spc),
                  S = spc)

# 初期値
init_fun <- function() {
  list(D_bar = runif(nrow(data), 10, 50),
       log_H_max_bar = rnorm(1, 4, 2),
       log_a_bar = rnorm(1, 0, 1),
       log_h_bar = rnorm(1, 0, 1),
       e_a = rnorm(max(spc), 0, 1),
       e_h = rnorm(max(spc), 0, 1),
       e_H_max = rnorm(max(spc), 0, 1),
       e = rnorm(nrow(data), 0, 1),
       sigma = runif(6, 0, 1))
}

# 擬似乱数系列の指定
set.seed(123)

# あてはめ
model <- cmdstan_model("D-H.stan")
fit <- model$sample(data = stan_data, init = init_fun,
                    seed = 12, chains = 4,
                    iter_sampling = 8000, iter_warmup = 8000,
                    refresh = 800,
                    adapt_delta = 0.9, max_treedepth = 15)

# 事後分布の要約の表示
fit$print(c("a", "h", "H_max"))
fit$print(c("log_H_max_bar", "log_a_bar", "log_h_bar"))
fit$print("sigma")
fit$draws("sigma") |>
  bayesplot::mcmc_trace()

# グラフ表示
library(ggplot2)

# 樹高を計算する関数の定義
get_height <- function(D, a, h, H_max, H_D) {
  a * D^h * H_max / (a * D^h + H_max) + H_D
}

# 樹種のリスト
spp <- levels(factor(data$Jpn_name))

# 取り出す樹種を指定する
s <- c("カナメモチ", "コナラ")
n_sp <- length(s)

# 指定した樹種の樹高測定値の平均を計算する
data2 <- dplyr::filter(data, Jpn_name %in% s) %>%
  dplyr::rowwise() %>%
  dplyr::mutate(H_mean = mean(c(H1, H2, H3, H4), na.rm = TRUE))

# 指定した樹種の位置を取り出す
pos <- match(s, levels(factor(data$Jpn_name)))

# MCMCサンプルの取出
log_H_max_bar = fit$draws("log_H_max_bar", format = "draws_matrix") |>
  c()
log_a_bar = fit$draws("log_a_bar", format = "draws_matrix") |>
  c()
log_h_bar = fit$draws("log_h_bar", format = "draws_matrix") |>
  c()
e_H_max = fit$draws("e_H_max", format = "draws_matrix")[, pos] |>
  c()
e_a = fit$draws("e_a", format = "draws_matrix")[, pos] |>
  c()
e_h = fit$draws("e_h", format = "draws_matrix")[, pos] |>
  c()

# 樹種ごとのパラメータのMCMCサンプル
H_max <- exp(replicate(n_sp, log_H_max_bar) + e_H_max)
a <- exp(replicate(n_sp, log_a_bar) + e_a)
h <- exp(replicate(n_sp, log_h_bar) + e_h)

# 予測につかう直径
D <- seq(0.1, 65, length = 100)

# MCMCサンプルの値を使って樹種ごとに直径-樹高関係を計算する
samp_h <- lapply(1:n_sp, function(i) {
  sapply(D, function(x) {
    get_height(x, a[, i], h[, i], H_max[, i], 1.3)
  })
})

# 中央値と95%確信区間を計算する
H_median <- lapply(1:n_sp, function(i) {
  apply(samp_h[[i]], 2, median)
})
H_lower <- lapply(1:n_sp, function(i) {
  apply(samp_h[[i]], 2, quantile, probs = 0.025)
})
H_upper <- lapply(1:n_sp, function(i) {
  apply(samp_h[[i]], 2, quantile, probs = 0.975)
})

gg_data <- data.frame(Jpn_name = rep(s, each = length(D)),
                      D = rep(D, n_sp),
                      H_median = unlist(H_median),
                      H_lower = unlist(H_lower),
                      H_upper = unlist(H_upper))

# グラフ表示
ggplot(gg_data) +
  geom_ribbon(aes(x = D, ymin = H_lower, ymax = H_upper),
              colour = "grey", alpha = 0.5) +
  geom_line(aes(x = D, y = H_median),
            colour = "black") +
  geom_point(data = data2, aes(x = D, y = H_mean),
             size = 2) +
  labs(x = "直径(cm)", y = "高さ(m)") +
  xlim(0, max(D)) + ylim(0, 30) +
  facet_wrap(~Jpn_name) +
  theme_bw(base_family = "IPAexGothic", base_size = 10)

# グラフの保存
ggsave("fig_d-h.pdf", device = cairo_pdf,
       width = 12, height = 6, units = "cm")

D-H.stan

/*
  Allometric equation by Ogawa et al. (1965)
  1/H = 1/(aD^h) + 1/H_max
  
  Modeled for Stan 2.32
*/
data {
  int<lower=0> N;              // 木の本数
  vector<lower=0>[N] D;        // 直径の測定値(cm)
  array[N] int<lower=0> Nh;    // 木ごとの樹高測定回数
  int<lower=0> MaxNh;          // 木ごとの樹高測定回数の最大値
  matrix<lower=0>[N, MaxNh] H; // 樹高の測定値(m)
  int<lower=0> Ns;             // 樹種の数
  array[N] int<lower=1> S;     // 樹種のインデックス
}
transformed data {
  matrix[N, MaxNh] H2 = H - 1.3; // 樹高-1.3(m)
}
parameters {
  vector<lower=0>[N] D_bar;   // 「真の」直径
  real log_H_max_bar;         // 全体的なH_maxの対数 (最大樹高-1.3m)
  real log_a_bar;             // 全体的なaの対数
  real log_h_bar;             // 全体的なhの対数
  vector[Ns] e_a;             // 種によるaのランダム効果
  vector[Ns] e_h;             // 種によるのhのランダム効果
  vector[Ns] e_H_max;         // 種によるH_maxのランダム効果
  vector[N] e;                // 全体の残差
  vector<lower=0>[6] sigma;   // 各標準偏差
}
transformed parameters {
  vector<lower=0>[Ns] a;      // それぞれの種のaパラメータ
  vector<lower=0>[Ns] h;      // それぞれの種のhパラメータ
  vector<lower=0>[Ns] H_max;  // それぞれの種H_maxパラメータ
  vector[N] log_H_bar;        // 「真の」高さ-1.3mの対数

  a = exp(log_a_bar + e_a);
  h = exp(log_h_bar + e_h);
  H_max = exp(log_H_max_bar + e_H_max);
  for (i in 1:N) {
    log_H_bar[i] = log(a[S[i]]) + h[S[i]] * log(D_bar[i])
                   + log(H_max[S[i]])
                   - log(a[S[i]] * pow(D_bar[i], h[S[i]]) + H_max[S[i]])
                   + e[i];
  }
}
model {
  // 直径の誤差
  D ~  normal(D_bar, sigma[1]);
  // 高さの誤差
  for	(i in 1:N)
    H2[i, 1:Nh[i]] ~ lognormal(log_H_bar[i], sigma[2]);

  e ~ normal(0, sigma[3]);
  e_a ~ normal(0, sigma[4]);
  e_h ~ normal(0, sigma[5]);
  e_H_max ~ normal(0, sigma[6]);
  
  log_H_max_bar ~ normal(4, 4);
  D_bar ~ normal(0, 100);
  log_a_bar ~ normal(0, 10);
  log_h_bar ~ normal(0, 10);
  sigma ~ normal(0, 10);
}