Rの関数:MCMCglmm {MCMCglmm}

# ライブラリの読み込み
library(ggplot2)
library(gridExtra)

# 乱数シードの設定（再現性確保のため）
seed <- 20251208
set.seed(seed)

# パラメータ設定
n_regions <- 30
n_plots_per_region <- 5 # 1地域あたりの区画数
n_trees_per_plot <- 20 # 1区画あたりの樹木数

# 総サンプルサイズ
total_obs <- n_regions * n_plots_per_region * n_trees_per_plot

# IDの生成
region_id <- rep(1:n_regions, each = n_plots_per_region * n_trees_per_plot)
# 区画IDを1から通し番号で作ることで、階層構造（地域ごとのユニークさ）を表現
plot_id <- rep(1:(n_regions * n_plots_per_region), each = n_trees_per_plot)

# 説明変数の生成
temperature <- rnorm(total_obs, mean = 15, sd = 3)
fertilizer <- sample(c("なし", "あり"), total_obs, replace = TRUE)

# 真のパラメータ（固定効果）
intercept_true <- 10.0 # 切片
beta_temp <- 0.5 # 気温の効果
beta_fert <- 2.0 # 施肥の効果

# 真の分散成分（変量効果と残差）
sigma_region <- 3.0 # 地域間の標準偏差（分散=9.0）
sigma_plot <- 1.5 # 区画間の標準偏差（分散=2.25）
sigma_resid <- 1.0 # 残差標準偏差（分散=1.0）

# 変量効果の実現値を生成
region_eff_vals <- rnorm(n_regions, 0, sigma_region)
region_effect <- region_eff_vals[region_id]

plot_eff_vals <- rnorm(n_regions * n_plots_per_region, 0, sigma_plot)
plot_effect <- plot_eff_vals[plot_id]

residuals <- rnorm(total_obs, 0, sigma_resid)

# 応答変数の生成
fert_dummy <- ifelse(fertilizer == "あり", 1, 0)
y_val <- intercept_true +
  beta_temp * temperature +
  beta_fert * fert_dummy +
  region_effect +
  plot_effect +
  residuals

# データフレームの作成
forest_data <- data.frame(
  TreeHeight = y_val,
  Temperature = temperature,
  Fertilizer = factor(fertilizer, levels = c("なし", "あり")),
  Region = factor(region_id),
  Plot = factor(plot_id)
)

cat("--- サンプルデータの一部を確認 ---\n")
print(str(forest_data))

# 事前分布の設定（弱情報事前分布）
prior_setting <- list(
  R = list(V = 1, nu = 0.002),
  G = list(
    G1 = list(V = 1, nu = 0.002), # Region用
    G2 = list(V = 1, nu = 0.002) # Plot用
  )
)

model_fit <- MCMCglmm(
  fixed = TreeHeight ~ Temperature + Fertilizer,
  random = ~ Region + Plot,
  rcov = ~units,
  family = "gaussian",
  data = forest_data,
  prior = prior_setting,
  nitt = 30000,
  burnin = 10000,
  thin = 10,
  verbose = FALSE
)

cat("--- モデルの要約結果 ---\n")
summary(model_fit)