【成長データ解析の基礎と再考 (2)】誰でも簡単、成長曲線推定：シミュレーションデータでやってみるLMS法

今回は、LMS法を使って実際に成長曲線を推定してみます。

　子どもの身長や体重の評価では、「同じ年齢の子どもたちの中で、その子がどの位置にいるか」を示すパーセンタイル（centile）が広く用いられます。たとえば、「50パーセンタイル」であれば同年齢の子どもの中でちょうど中央に位置することを意味し、「97パーセンタイル」であれば上位3%に相当することを表します。

　この年齢に伴うパーセンタイルの変化を線として描いたものがパーセンタイル曲線（centile curve）であり、複数のパーセンタイル値（例：3rd・10th・25th・50th・75th・90th・97th）に対応する曲線をまとめて描いたものが成長曲線（growth curve）と呼ばれます。

　このような成長曲線を統計的に推定するための代表的な手法が、Box-Cox-Cole-Green（BCCG）分布に基づくLMS法です。LMS法では、年齢ごとの身長・体重などの分布を、L（ゆがみ：Box-Cox変換のべき乗パラメタ）・M（中央値）・S（変動係数）の3つのパラメタで記述します。これら3つを年齢の滑らかな関数として推定することで、任意の年齢におけるパーセンタイルを算出できるようになります。

　本記事では、Rを使って以下の内容を体験します。

LMSパラメタを用いたデータ生成
スムージングによる連続モデルの構築
成長曲線の再構成と可視化
GAMLSSパッケージを用いたLMSパラメタの推定

1. データ生成の考え方
2. BCCG分布による成長曲線推定
3. まとめ：分布モデルを仮定するパラメトリック推定は正しい？

1. データ生成の考え方

1.1 LMSモデルとは

　LMS法では、年齢 $t$ ごとの身長・体重などの分布を、次の3つのパラメータで表現します。

$\displaystyle{L(t)}$ ：歪度パラメータ（Box-Cox変換のべき乗）——分布の左右の歪みを補正します
$\displaystyle{M(t)}$ ：中央値（50パーセンタイル）——その年齢での典型的な値を表します
$\displaystyle{S(t)}$ ：変動係数（標準偏差 ÷ 平均の近似値）——分布の広がりを表します

そして、観測値 $\displaystyle{Y}$ に対して次の変換を施すと、

$\displaystyle{ Z = \frac{(Y/M)^L - 1}{L \cdot S} }$

この $\displaystyle{ Z }$ が標準正規分布 $\displaystyle{ N(0,1) }$ に従うようにモデル化します（ $\displaystyle{ L \neq 0 }$ のとき。 $\displaystyle{ L = 0 }$ のときは $\displaystyle{ Z = \log(Y/M)/S }$ ）。

　つまりLMS法の本質は、「年齢ごとに異なる歪んだ分布を、Box-Cox変換によって正規分布に直す」という操作です。 $\displaystyle{ L }$ 、 $\displaystyle{ M }$ 、 $\displaystyle{ S }$ の3つを年齢の滑らかな関数として推定することで、任意の年齢における任意のパーセンタイル値を計算できるようになります。

逆変換を使えば、標準正規分位点 $z$ から元のスケールの観測値に戻すことも可能です。

$\displaystyle{ Y = M \cdot (1 + L \cdot S \cdot z)^{1/L} }$

この逆変換が、後述のデータ生成やパーセンタイル曲線の再構成で中心的な役割を果たします。

　以下が具体的なRスクリプトです。興味があれば、このコードをコピーして、Rで実行してみてください。ただし、正しく実行するには、

# 出力先ディレクトリ
DIR <- r"(... ここにフォルダパスを指定 ...)"

の部分に、出力ファイルの保存先となるフォルダのパスを指定する必要があります。

# =========================================================
# LMS パラメータに基づく成長データのモデルシミュレーション
# ---------------------------------------------------------
# このスクリプトは以下の処理を順に行う：
#
#   1) 性別・測定項目ごとに年齢依存型の LMS テーブルを定義する
#   2) L・M・S を年齢の連続関数としてスムージングする
#   3) スムージング済み LMS モデルから模擬観測データを生成する
#   4) モデルに基づくパーセンタイル曲線を再構成する
#   5) 模擬データと「真の」パーセンタイル曲線をエクスポートする
#
# ※ここで出力される「真の」パーセンタイル曲線は、元の離散的な LMS テーブル値
#   そのものではない。シミュレーションに使用したスムージング済み LMS 関数から
#   導出されるパーセンタイル曲線である。
# =========================================================


# ---------------------------------------------------------
# 1. グローバル設定
# ---------------------------------------------------------

# 出力先ディレクトリ
DIR <- r"(... ここにファイルパスを指定 ...)"

# 元の LMS 値が与えられている年齢点（ノット）
# ここでは 5 歳から 17 歳の整数年齢
age_knots <- 5:17

# 連続的なモデルベースパーセンタイル曲線を評価する年齢グリッド
# 0.1 年刻みで 5 歳から 17 歳まで出力する
age_grid <- seq(5, 17, by = 0.1)

# 推定対象のパーセンタイル（確率値）
probs <- c(0.03, 0.10, 0.25, 0.50, 0.75, 0.90, 0.97)

# パーセンタイルに対応する列名ラベル
prob_names <- c("p3", "p10", "p25", "p50", "p75", "p90", "p97")

# 性別・測定項目ごとのシミュレーションサンプルサイズ
n_sim <- 10000

# smooth.spline() に使用するスムージングパラメータ
# L・M・S に同じ値を使用（簡略化のため）
# spar が小さいほど曲線は元データに近くなり、大きいほど滑らかになる
spar_L <- 0.35
spar_M <- 0.35
spar_S <- 0.35

# 再現性確保のための乱数シード
seed_sim <- 123


# ---------------------------------------------------------
# 2. 入力 LMS テーブル
# ---------------------------------------------------------
# 各テーブルには年齢ごとの LMS パラメータが含まれる：
#
#   L : Box-Cox べき変換パラメータ（分布の歪みを補正する）
#   M : 中央値（50パーセンタイル）
#   S : 変動係数（分布の広がりを表す）
#
# これらのテーブルは入力データとして扱われる。スクリプト内では
# スムージングを施し、連続的な年齢依存 LMS 関数を定義する。
#
# LMS モデルの基本：
#   測定値 Y が LMS パラメータ (L, M, S) に従う場合、
#   正規化変換 Z = [(Y/M)^L - 1] / (L*S) は標準正規分布に従う（L ≠ 0 のとき）。
#   L = 0 のときは Z = log(Y/M) / S となる。

# ---- 男児：身長 ----
male_h_lms <- data.frame(
  age = age_knots,
  L = c(
     0.39189524,  0.07567997, -0.09136530,  0.23251045, -0.03336290,
    -0.60786430, -0.90515887,  0.78411177,  2.65704267,  2.37421076,
     0.64531977,  0.61711217,  0.98746908
  ),
  M = c(
    110.2860, 116.3997, 122.4224, 128.1819, 133.6935,
    139.1927, 145.6293, 153.6837, 161.2200, 166.1085,
    168.4837, 169.7087, 170.4979
  ),
  S = c(
    0.04298495, 0.04263699, 0.04222310, 0.04303110, 0.04308454,
    0.04491462, 0.04976466, 0.05243575, 0.04605124, 0.03887508,
    0.03502673, 0.03464069, 0.03467773
  )
)

# ---- 女児：身長 ----
female_h_lms <- data.frame(
  age = age_knots,
  L = c(
     0.39696980,  0.08285670,  0.06559417, -0.34334827, -0.58367403,
     0.41793428,  2.29001050,  2.20165372,  1.08364665,  0.74360945,
     0.68836865, -0.28622074,  0.53999018
  ),
  M = c(
    109.3860, 115.4997, 121.4318, 127.3382, 133.5641,
    140.6137, 147.4969, 152.3599, 154.8270, 156.2191,
    156.8107, 157.2234, 157.6595
  ),
  S = c(
    0.04333879, 0.04296948, 0.04290787, 0.04474366, 0.04709193,
    0.04892776, 0.04485439, 0.03825326, 0.03494484, 0.03435277,
    0.03436939, 0.03423403, 0.03439457
  )
)

# ---- 男児：体重 ----
male_w_lms <- data.frame(
  age = age_knots,
  L = c(
    -1.0183186, -1.1903911, -1.3547857, -1.2552963, -1.1764403,
    -1.0417090, -0.7758495, -0.5903966, -0.5258445, -0.6764802,
    -0.9351152, -0.9819073, -0.9584327
  ),
  M = c(
    18.63317, 20.86672, 23.49576, 26.49450, 29.91390,
    33.61381, 38.07091, 43.59493, 48.90715, 53.54038,
    57.40594, 58.93501, 60.67956
  ),
  S = c(
    0.1284525, 0.1435896, 0.1546339, 0.1718183, 0.1889550,
    0.2012785, 0.2116673, 0.2115569, 0.1944057, 0.1760894,
    0.1697971, 0.1564848, 0.1566104
  )
)

# ---- 女児：体重 ----
female_w_lms <- data.frame(
  age = age_knots,
  L = c(
    -0.9164836, -1.1261708, -1.1299641, -1.1134451, -1.0054783,
    -0.7095540, -0.3881036, -0.3576415, -0.6101772, -0.7935914,
    -0.9153502, -1.0291287, -0.9786760
  ),
  M = c(
    18.34246, 20.49235, 22.98892, 25.94783, 29.38421,
    33.67507, 38.78154, 43.47399, 46.61968, 48.81906,
    50.14582, 50.88365, 51.57595
  ),
  S = c(
    0.1302038, 0.1463061, 0.1547570, 0.1716172, 0.1859327,
    0.1988751, 0.1957782, 0.1777411, 0.1559749, 0.1458770,
    0.1483175, 0.1417884, 0.1447708
  )
)


# ---------------------------------------------------------
# 3. ユーティリティ関数
# ---------------------------------------------------------

# --- 関数: check_lms_table ---
# LMS テーブルの妥当性を検証する。
#
# 以下の条件を確認する：
#   - "age", "L", "M", "S" の列が存在すること
#   - NA・NaN・Inf が含まれないこと
#   - M（中央値）と S（変動係数）が正の値であること
#
# 引数:
#   dat  : 検証対象の data.frame
#   name : エラーメッセージ用のテーブル名（省略時は変数名）
check_lms_table <- function(dat, name = deparse(substitute(dat))) {
  req <- c("age", "L", "M", "S")

  if (!all(req %in% names(dat))) {
    stop(sprintf(
      "%s には次の列が必要です: %s",
      name, paste(req, collapse = ", ")
    ))
  }

  if (any(!is.finite(as.matrix(dat[, req])))) {
    stop(sprintf("%s に NA, NaN, または Inf が含まれています。", name))
  }

  if (any(dat$M <= 0) || any(dat$S <= 0)) {
    stop(sprintf("%s の M および S はすべて正の値でなければなりません。", name))
  }

  invisible(TRUE)
}


# --- 関数: qlms ---
# LMS モデルにおける分位点関数（逆 CDF）。
#
# 確率 p と LMS パラメータ (L, M, S) に対して、
# 対応するパーセンタイル値を返す。
#
# 変換式（LMS の標準的な定式化）：
#   z = Φ^{-1}(p)  （標準正規分布の分位点）
#   L ≠ 0 のとき： y = M * (1 + L*S*z)^(1/L)
#   L = 0 のとき： y = M * exp(S*z)
#
# ※(1 + L*S*z) > 0 が保証されない場合、対応する y は NA となる。
#
# 引数:
#   p   : 確率値（スカラーまたはベクトル）
#   L, M, S : LMS パラメータ（すべてスカラー、またはすべてベクトル）
#
# L, M, S がベクトルの場合、p はスカラーでなければならない（年齢軸方向の計算）。
# p がベクトルの場合、L, M, S はスカラーでなければならない（複数パーセンタイルの計算）。
qlms <- function(p, L, M, S) {
  z <- qnorm(p)  # 標準正規分位点に変換

  # --- L がスカラーの場合（複数の確率値 p に対してパーセンタイルを計算）
  if (length(L) == 1L) {
    if (abs(L) < 1e-8) {
      # L ≈ 0 のとき：対数正規分布に相当する変換
      return(M * exp(S * z))
    } else {
      # L ≠ 0 の通常ケース：Box-Cox 逆変換
      val <- 1 + L * S * z
      out <- rep(NA_real_, length(z))
      ok  <- val > 0  # 変換が有効な要素のみ計算
      out[ok] <- M * val[ok]^(1 / L)
      return(out)
    }
  }

  # --- L がベクトルの場合（年齢ごとに異なる LMS で 1 つのパーセンタイルを計算）
  if (length(p) != 1L) {
    stop("L, M, S がベクトルの場合、p はスカラーでなければなりません。")
  }

  out <- rep(NA_real_, length(L))

  # L ≈ 0 の年齢点：対数正規近似
  idx0 <- abs(L) < 1e-8
  out[idx0] <- M[idx0] * exp(S[idx0] * z)

  # L ≠ 0 の年齢点：Box-Cox 逆変換
  idx1 <- !idx0
  val  <- 1 + L[idx1] * S[idx1] * z
  ok   <- val > 0
  tmp  <- rep(NA_real_, sum(idx1))
  tmp[ok] <- M[idx1][ok] * val[ok]^(1 / L[idx1][ok])
  out[idx1] <- tmp

  out
}


# --- 関数: rlms ---
# LMS モデルに従う乱数を生成する。
#
# 単一の (L, M, S) パラメータセットから n 個の観測値を生成する。
#
# 生成メカニズム：
#   z ~ N(0, 1) を生成し、LMS 変換の逆変換を適用する。
#   L ≠ 0 のとき： y = M * (1 + L*S*z)^(1/L)
#
# ※(1 + L*S*z) ≤ 0 となる z の値は LMS モデルの定義域外なので、
#   そのような z は条件を満たすまで再サンプリングされる。
#   通常は極めて稀（分布の裾が非常に厚い場合のみ発生）。
#
# 引数:
#   n    : 生成するサンプル数
#   L, M, S : 単一の LMS パラメータ（スカラー）
rlms <- function(n, L, M, S) {
  z <- rnorm(n)  # 標準正規乱数を生成

  if (abs(L) < 1e-8) {
    # L ≈ 0 のとき：対数正規乱数
    return(M * exp(S * z))
  }

  val <- 1 + L * S * z

  # 定義域外の z を再サンプリングする
  while (any(val <= 0)) {
    bad      <- which(val <= 0)
    z[bad]   <- rnorm(length(bad))
    val      <- 1 + L * S * z
  }

  M * val^(1 / L)
}


# --- 関数: smooth_lms_curve ---
# 離散的な LMS テーブルを年齢の連続関数としてスムージングする。
#
# スムージング方針：
#   L : 元のスケールのままスプライン回帰
#   M : 対数スケールでスプライン回帰 → exp 変換（正値性を保証）
#   S : 対数スケールでスプライン回帰 → exp 変換（正値性を保証）
#
# M と S を対数スケールでスムージングすることで、
# 予測値が常に正になることが保証される。
#
# 引数:
#   age, L, M, S : 離散 LMS テーブルの各ベクトル
#   spar_L, spar_M, spar_S : smooth.spline() のスムージングパラメータ（0〜1）
smooth_lms_curve <- function(age, L, M, S,
                             spar_L = 0.5,
                             spar_M = 0.5,
                             spar_S = 0.5) {
  fit_L    <- smooth.spline(age, L,      spar = spar_L)
  fit_logM <- smooth.spline(age, log(M), spar = spar_M)  # 対数スケールでフィット
  fit_logS <- smooth.spline(age, log(S), spar = spar_S)  # 対数スケールでフィット

  list(
    fit_L    = fit_L,
    fit_logM = fit_logM,
    fit_logS = fit_logS
  )
}


# --- 関数: predict_lms_curve ---
# スムージング済み LMS モデルから任意の年齢における L, M, S を予測する。
#
# smooth_lms_curve() の出力オブジェクトを受け取り、
# 指定した年齢グリッド xout での予測値を data.frame として返す。
# M と S は対数スケールで予測した後に exp 変換して元のスケールに戻す。
#
# 引数:
#   fit_obj : smooth_lms_curve() の返り値（リスト）
#   xout    : 予測する年齢のベクトル
predict_lms_curve <- function(fit_obj, xout) {
  Lhat <- as.numeric(predict(fit_obj$fit_L,    x = xout)$y)
  Mhat <- exp(as.numeric(predict(fit_obj$fit_logM, x = xout)$y))  # exp で元スケールへ
  Shat <- exp(as.numeric(predict(fit_obj$fit_logS, x = xout)$y))  # exp で元スケールへ

  data.frame(
    age = as.numeric(xout),
    L   = Lhat,
    M   = Mhat,
    S   = Shat
  )
}


# --- 関数: build_model_from_lms_table ---
# 離散 LMS テーブルから連続 LMS モデルを構築する。
#
# 処理の流れ：
#   1. 入力テーブルの妥当性を検証する（check_lms_table）
#   2. L・M・S を年齢の連続関数としてスムージングする（smooth_lms_curve）
#
# 返り値はリストで、以下の要素を含む：
#   raw    : 元の離散 LMS テーブル
#   smooth : smooth_lms_curve() のフィットオブジェクト
#
# 引数:
#   lms_table : "age", "L", "M", "S" 列を持つ data.frame
#   spar_L, spar_M, spar_S : スムージングパラメータ
build_model_from_lms_table <- function(lms_table,
                                       spar_L = 0.5,
                                       spar_M = 0.5,
                                       spar_S = 0.5) {
  check_lms_table(lms_table)  # 入力検証

  fit_sm <- smooth_lms_curve(
    age    = lms_table$age,
    L      = lms_table$L,
    M      = lms_table$M,
    S      = lms_table$S,
    spar_L = spar_L,
    spar_M = spar_M,
    spar_S = spar_S
  )

  list(
    raw    = lms_table,
    smooth = fit_sm
  )
}


# --- 関数: simulate_from_lms_model ---
# スムージング済み LMS モデルから模擬観測データを生成する。
#
# アルゴリズム：
#   1. 年齢を [age_min, age_max] の一様分布からサンプリングする
#   2. サンプリングされた各年齢に対して、スムージング済み L, M, S を予測する
#   3. 各年齢の LMS パラメータから 1 つの観測値を乱数生成する（rlms）
#
# 返り値: sex・item・age・value および対応する L, M, S を含む data.frame
#
# 引数:
#   n        : 生成するサンプル数
#   model    : build_model_from_lms_table() の返り値
#   age_min, age_max : 年齢のサンプリング範囲
#   sex      : 性別ラベル（文字列）
#   item     : 測定項目ラベル（文字列）
#   seed     : 乱数シード（再現性のため）
simulate_from_lms_model <- function(n, model,
                                    age_min = 5,
                                    age_max = 17,
                                    sex,
                                    item,
                                    seed = NULL) {
  if (!is.null(seed)) {
    set.seed(seed)
  }

  # 年齢を一様分布からサンプリング
  age_sim <- runif(n, min = age_min, max = age_max)

  # サンプリングされた年齢での L, M, S を予測
  prm <- predict_lms_curve(model$smooth, age_sim)

  # 各サンプル点で LMS 乱数を 1 つずつ生成
  y <- mapply(
    FUN = function(L, M, S) rlms(1, L, M, S),
    L   = prm$L,
    M   = prm$M,
    S   = prm$S
  )

  data.frame(
    sex   = sex,
    item  = item,
    age   = age_sim,
    value = as.numeric(y),
    L     = prm$L,
    M     = prm$M,
    S     = prm$S
  )
}


# --- 関数: reconstruct_centiles ---
# スムージング済み LMS モデルからパーセンタイル曲線を再構成する。
#
# 処理：
#   1. age_grid における L, M, S を予測する
#   2. probs に対応するパーセンタイル値を qlms で計算する
#
# 返り値はリストで、以下の要素を含む：
#   age  : 評価年齢グリッド
#   qmat : パーセンタイル行列（行 = 年齢点、列 = パーセンタイル）
#   prm  : 年齢グリッド上のスムージング済み L, M, S
#
# 引数:
#   model     : build_model_from_lms_table() の返り値
#   age_grid  : パーセンタイルを評価する年齢のベクトル
#   probs     : 対象パーセンタイルの確率値ベクトル
#   prob_names: 各パーセンタイルの列名ベクトル
reconstruct_centiles <- function(model, age_grid, probs, prob_names) {
  prm <- predict_lms_curve(model$smooth, age_grid)

  # probs の各確率値についてパーセンタイル曲線を計算し、列として並べる
  qmat <- sapply(probs, function(p) {
    qlms(p, prm$L, prm$M, prm$S)
  })

  qmat <- as.matrix(qmat)
  colnames(qmat) <- prob_names

  list(
    age  = age_grid,
    qmat = qmat,
    prm  = prm
  )
}


# --- 関数: make_true_centile_df ---
# モデルに基づくパーセンタイル曲線を data.frame 形式に変換する。
#
# シミュレーション研究における「真の」パーセンタイル曲線の
# 主な出力形式として使用される。
# sex・item・age・L, M, S・各パーセンタイル値を列として持つ。
#
# 引数:
#   model     : build_model_from_lms_table() の返り値
#   age_grid  : パーセンタイルを評価する年齢のベクトル
#   probs     : 対象パーセンタイルの確率値ベクトル
#   prob_names: 各パーセンタイルの列名ベクトル
#   sex       : 性別ラベル（文字列）
#   item      : 測定項目ラベル（文字列）
make_true_centile_df <- function(model, age_grid, probs, prob_names, sex, item) {
  rec <- reconstruct_centiles(model, age_grid, probs, prob_names)

  data.frame(
    sex  = sex,
    item = item,
    age  = rec$prm$age,
    L    = rec$prm$L,
    M    = rec$prm$M,
    S    = rec$prm$S,
    rec$qmat,
    row.names = NULL
  )
}


# ---------------------------------------------------------
# 4. プロット設定
# ---------------------------------------------------------
# 色の割り当て：下位・上位パーセンタイルを視覚的に分離し、
# 中央値（50%）はより強調する配色を採用している。

col_p97 <- "#EF9F27"  # 97パーセンタイル：オレンジ系
col_p90 <- "#D85A30"  # 90パーセンタイル：赤橙系
col_p75 <- "#BA7517"  # 75パーセンタイル：黄茶系
col_p50 <- "#1A1A1A"  # 50パーセンタイル（中央値）：ほぼ黒
col_p25 <- "#378ADD"  # 25パーセンタイル：青系
col_p10 <- "#185FA5"  # 10パーセンタイル：濃い青
col_p3  <- "#85B7EB"  # 3パーセンタイル：淡い青

# パーセンタイル曲線の色と線幅（p3 から p97 の順）
centile_cols <- c(col_p3, col_p10, col_p25, col_p50, col_p75, col_p90, col_p97)
centile_lwds <- c(1.5, 2.0, 2.0, 3.0, 2.0, 2.0, 1.5)  # 中央値を最も太く


# --- 関数: plot_panel ---
# 模擬観測データとモデルベースのパーセンタイル曲線を 1 パネル描画する。
#
# 描画の流れ：
#   1. 模擬観測データを散布図として描画する（半透明グレー点）
#   2. 横・縦グリッド線を追加する
#   3. 各パーセンタイル曲線を重ね描きする
#   4. パーセンタイルの凡例を表示する
#
# 引数:
#   rec    : reconstruct_centiles() の返り値
#   sim_df : simulate_from_lms_model() の返り値（模擬データ）
#   main   : プロットのタイトル
#   ylab   : y 軸ラベル
#   ylim   : y 軸の範囲（2 要素ベクトル）
#   ygrid  : 横グリッド線を引く y 値のベクトル（NULL で非表示）
plot_panel <- function(rec, sim_df, main, ylab, ylim, ygrid = NULL) {
  # 散布図（模擬観測値）
  plot(
    sim_df$age, sim_df$value,
    pch  = 16, cex = 0.7, col = gray(0.7, 0.3),  # 半透明の小点
    xlim = c(5, 17), ylim = ylim,
    xlab = "年齢（歳）", ylab = ylab,
    main = main
  )

  # グリッド線の追加
  if (!is.null(ygrid)) {
    abline(h = ygrid, col = "gray92", lty = 1)       # 横グリッド（薄い灰色）
  }
  abline(v = 5:17, col = "gray96", lty = 1)           # 縦グリッド（各年齢）

  # パーセンタイル曲線の描画
  for (j in seq_len(ncol(rec$qmat))) {
    lines(rec$age, rec$qmat[, j], col = centile_cols[j], lwd = centile_lwds[j])
  }

  # 凡例の追加
  legend(
    "topleft",
    legend = c("3%", "10%", "25%", "50%", "75%", "90%", "97%"),
    col    = centile_cols,
    lty    = 1,
    lwd    = centile_lwds,
    bty    = "o",
    bg     = "white",
    cex    = 0.8
  )
}


# ---------------------------------------------------------
# 5. スムージング済み LMS モデルの構築
# ---------------------------------------------------------
# 各性別・測定項目の組み合わせに対して、
# 離散 LMS テーブルから連続 LMS モデルを作成する。

model_male_h   <- build_model_from_lms_table(
  male_h_lms,   spar_L = spar_L, spar_M = spar_M, spar_S = spar_S
)
model_female_h <- build_model_from_lms_table(
  female_h_lms, spar_L = spar_L, spar_M = spar_M, spar_S = spar_S
)
model_male_w   <- build_model_from_lms_table(
  male_w_lms,   spar_L = spar_L, spar_M = spar_M, spar_S = spar_S
)
model_female_w <- build_model_from_lms_table(
  female_w_lms, spar_L = spar_L, spar_M = spar_M, spar_S = spar_S
)


# ---------------------------------------------------------
# 6. 模擬データの生成
# ---------------------------------------------------------
# 乱数シードを設定してから各グループの模擬データを生成する。
# set.seed() はシミュレーション全体の先頭に一度だけ設定する。

set.seed(seed_sim)

sim_male_h   <- simulate_from_lms_model(
  n = n_sim, model = model_male_h,   sex = "Male",   item = "Height"
)
sim_female_h <- simulate_from_lms_model(
  n = n_sim, model = model_female_h, sex = "Female", item = "Height"
)
sim_male_w   <- simulate_from_lms_model(
  n = n_sim, model = model_male_w,   sex = "Male",   item = "Weight"
)
sim_female_w <- simulate_from_lms_model(
  n = n_sim, model = model_female_w, sex = "Female", item = "Weight"
)

# 全グループの模擬データを結合する
sim_all <- rbind(sim_male_h, sim_female_h, sim_male_w, sim_female_w)


# ---------------------------------------------------------
# 7. モデルベースのパーセンタイル曲線の再構成
# ---------------------------------------------------------
# 各グループのスムージング済み LMS モデルから、
# 指定した年齢グリッドと確率値に対応するパーセンタイル曲線を計算する。

rec_male_h   <- reconstruct_centiles(model_male_h,   age_grid, probs, prob_names)
rec_female_h <- reconstruct_centiles(model_female_h, age_grid, probs, prob_names)
rec_male_w   <- reconstruct_centiles(model_male_w,   age_grid, probs, prob_names)
rec_female_w <- reconstruct_centiles(model_female_w, age_grid, probs, prob_names)


# ---------------------------------------------------------
# 8. 「真の」パーセンタイル曲線の data.frame への変換とエクスポート準備
# ---------------------------------------------------------
# ここでの「真の」曲線とは、スムージング済み LMS モデルから導出された
# 理論的なパーセンタイル曲線を指す（元の離散テーブル値ではない）。

true_centiles_male_h <- make_true_centile_df(
  model      = model_male_h,
  age_grid   = age_grid,
  probs      = probs,
  prob_names = prob_names,
  sex        = "Male",
  item       = "Height"
)

true_centiles_female_h <- make_true_centile_df(
  model      = model_female_h,
  age_grid   = age_grid,
  probs      = probs,
  prob_names = prob_names,
  sex        = "Female",
  item       = "Height"
)

true_centiles_male_w <- make_true_centile_df(
  model      = model_male_w,
  age_grid   = age_grid,
  probs      = probs,
  prob_names = prob_names,
  sex        = "Male",
  item       = "Weight"
)

true_centiles_female_w <- make_true_centile_df(
  model      = model_female_w,
  age_grid   = age_grid,
  probs      = probs,
  prob_names = prob_names,
  sex        = "Female",
  item       = "Weight"
)

# 全グループの「真の」パーセンタイル曲線を結合する
true_centiles_all <- rbind(
  true_centiles_male_h,
  true_centiles_female_h,
  true_centiles_male_w,
  true_centiles_female_w
)


# ---------------------------------------------------------
# 9. 模擬データとパーセンタイル曲線のプロット
# ---------------------------------------------------------
# 2行 × 2列のパネル配置で 4 グループ（男女 × 身長・体重）を描画する。

op <- par(no.readonly = TRUE)  # 現在のグラフィックスパラメータを保存
par(
  mfrow = c(2, 2),             # 2行2列のパネル配置
  mar   = c(4, 4.5, 3, 1),    # 各パネルの余白（下・左・上・右）
  las   = 1,                   # 軸ラベルを常に水平にする
  oma   = c(0, 0, 1, 0)        # 全体の外側余白（タイトル用に上を確保）
)

plot_panel(
  rec    = rec_male_h,
  sim_df = sim_male_h,
  main   = "男児：身長",
  ylab   = "身長（cm）",
  ylim   = c(100, 185),
  ygrid  = seq(100, 180, 10)
)

plot_panel(
  rec    = rec_female_h,
  sim_df = sim_female_h,
  main   = "女児：身長",
  ylab   = "身長（cm）",
  ylim   = c(100, 185),
  ygrid  = seq(100, 180, 10)
)

plot_panel(
  rec    = rec_male_w,
  sim_df = sim_male_w,
  main   = "男児：体重",
  ylab   = "体重（kg）",
  ylim   = c(10, 90),
  ygrid  = seq(10, 90, 10)
)

plot_panel(
  rec    = rec_female_w,
  sim_df = sim_female_w,
  main   = "女児：体重",
  ylab   = "体重（kg）",
  ylim   = c(10, 90),
  ygrid  = seq(10, 90, 10)
)

par(op)  # グラフィックスパラメータを元に戻す


# ---------------------------------------------------------
# 10. 出力ファイルの保存
# ---------------------------------------------------------
# 出力ファイルの内容説明：
#
# sim_all（simulated_growth_data_from_LMS_only.csv）:
#   スムージング済み LMS モデルから生成された個人レベルの模擬観測データ。
#   各行は 1 サンプルに対応し、sex・item・age・value・L・M・S を含む。
#
# true_centiles_all（true_model_based_centile_curves.csv）:
#   年齢 5〜17 歳を 0.1 年刻みで評価したモデルベースの「真の」パーセンタイル曲線。
#   各行には年齢・性別・項目・スムージング済み L, M, S・
#   p3〜p97 の各パーセンタイル値が含まれる。

setwd(DIR)
write.csv(sim_all,           "simulated_growth_data_from_smoothed_LMS.csv",  row.names = FALSE)
write.csv(true_centiles_all, "true_model_based_centile_curves.csv",       row.names = FALSE)

1.2 シミュレーションの全体設計

　上のRスクリプトは、次の流れで動作します。

① 離散的なLMSテーブルを用意する

　年齢5〜17歳の整数年齢における $\displaystyle{ L }$ 、 $\displaystyle{ M }$ 、 $\displaystyle{ S }$ の値を、性別（男児・女児）×測定項目（身長・体重）の4グループ分あらかじめ用意します。これがシミュレーションの「真の」成長モデルの出発点です。

② スプラインでスムージングする

　用意した離散テーブルを smooth.spline() でスムージングし、任意の年齢（たとえば7.3歳）での $\displaystyle{ L(t) }$ 、 $\displaystyle{ M(t) }$ 、 $\displaystyle{ S(t) }$ を補間できる連続関数として表現します。このとき $\displaystyle{M}$ と $\displaystyle{S}$ は対数スケールでスムージングし、予測値が必ず正になるよう工夫しています（ $\displaystyle{L}$ はそのままのスケール）。スムージングの滑らかさは spar パラメータ（0〜1）で制御でき、値が小さいほど元データに近い曲線に、大きいほど滑らかな曲線になります。

③ 連続関数 $\displaystyle{L(t)}$ 、 $\displaystyle{M(t)}$ 、 $\displaystyle{S(t)}$ を得る

　スムージングの結果として、5〜17歳の任意の年齢点で $\displaystyle{ L }$ 、 $\displaystyle{ M }$ 、 $\displaystyle{ S }$ を予測できる関数オブジェクトが得られます。以降のすべての処理はこの連続関数を通して行われます。

④ 各年齢で確率分布を構成し、乱数を生成する

　 $\displaystyle{n = 10{,}000}$ 人分の年齢を一様分布からランダムに生成し、その年齢ごとに連続LMSモデルから $\displaystyle{ L }$ 、 $\displaystyle{ M }$ 、 $\displaystyle{ S }$ を取得します。次に、標準正規乱数 $\displaystyle{z \sim N(0,1)}$ を発生させ、上述の逆変換式でそれぞれの観測値（身長・体重）を生成します。 $\displaystyle{(1 + L \cdot S \cdot z) \leq 0}$ になるケースはLMSモデルの定義域外であるため、そのような乱数は再サンプリングされますが、通常は極めて稀にしか起きません。

⑤ パーセンタイル曲線を再構成する

　0.1歳刻みの年齢グリッドで $\displaystyle{ L }$ 、 $\displaystyle{ M }$ 、 $\displaystyle{ S }$ を予測し、対象パーセンタイル（3, 10, 25, 50, 75, 90, 97パーセンタイル）に対応する $\displaystyle{z}$ 値を標準正規分布の逆関数で求め、逆変換式で元のスケールに戻します。こうして得られる曲線が「モデルに基づく真のパーセンタイル曲線」です。これは元の離散テーブル値をそのまま繋いだものではなく、スムージング済みLMS関数から導出した理論曲線である点に注意してください。

　全体の流れを図式化すると次のようになります。

離散LMSテーブル → スムージング → 連続LMS関数 → 乱数生成（個体データ）

連続LMS関数 → パーセンタイル逆変換 → 「真の」成長曲線

つまり今回のシミュレーションは、「真の成長曲線」を先に定義し、そこから個体データを逆算する構造になっています。現実のデータ解析では観測データから成長曲線を推定しますが、このシミュレーションではその逆方向の操作を行っています。この仕組みを使うことで、後の推定ステップ（GAMLSSによる解析）で得た曲線が「真の曲線」にどれだけ近いかを定量的に評価できます。

1.3 スクリプトの構成と主な関数

　Rスクリプトは大きく10のセクションに分かれています。各セクションの役割を以下に整理します。

セクション1：グローバル設定

　出力先ディレクトリ、年齢グリッド、パーセンタイルの確率値、サンプルサイズ、スムージングパラメータ、乱数シードなど、スクリプト全体で使う定数をまとめて定義します。パラメータを変更したい場合は、このセクションの値を書き換えるだけで済みます。

セクション2：入力LMSテーブル

　性別×測定項目の4グループそれぞれについて、年齢5〜17歳の整数点での $\displaystyle{L}$ 、 $\displaystyle{M}$ 、 $\displaystyle{S}$ 値を data.frame として定義します。これがシミュレーション全体の入力データです。

セクション3：ユーティリティ関数

　モデル構築・シミュレーション・パーセンタイル計算に使う関数群を定義します。各関数の役割は以下のとおりです。

関数名	役割
`check_lms_table()`	入力LMSテーブルの妥当性を検証する
`qlms()`	LMSモデルの分位点関数（パーセンタイル値を計算する）
`rlms()`	LMSモデルに従う乱数を生成する
`smooth_lms_curve()`	離散LMSテーブルをスプラインでスムージングする
`predict_lms_curve()`	スムージング済みモデルから任意年齢の $\displaystyle{L}$ 、 $\displaystyle{M}$ 、 $\displaystyle{S}$ を予測する
`build_model_from_lms_table()`	入力テーブルから連続LMSモデルを構築する（検証＋スムージング）
`simulate_from_lms_model()`	連続LMSモデルから模擬観測データを生成する
`reconstruct_centiles()`	スムージング済みモデルからパーセンタイル曲線を再構成する
`make_true_centile_df()`	パーセンタイル曲線を `data.frame` 形式に変換する

セクション4：プロット設定

　パーセンタイル曲線の色と線幅を定義します。下位パーセンタイルを青系、上位を橙系、中央値（50パーセンタイル）を黒でより太く描くことで、視覚的に分布の位置が把握しやすい配色を採用しています。

セクション5〜8：モデル構築・データ生成・曲線再構成

　定義した関数を用いて、4グループ分のモデル構築・乱数生成・パーセンタイル曲線の再構成・data.frame への変換を順番に実行します。

セクション9：プロット

　男女×身長・体重の2行×2列パネルで、模擬観測データ（散布図）とパーセンタイル曲線（折れ線）を重ねて描画します。

セクション10：ファイルの保存

　以下の2つのCSVファイルを出力します。

simulated_growth_data_from_smoothed_LMS.csv：個人レベルの模擬観測データ（sex・item・age・value・L・M・S を列として持つ）
true_model_based_centile_curves.csv：モデルに基づく「真の」パーセンタイル曲線（0.1歳刻みで評価した各パーセンタイル値を列として持つ）

後のGAMLSSによる推定では、前者のCSVを観測データとして入力し、後者のCSVを「正解」として推定結果と比較します。

　以上が実際に手元でコードを動かす前に把握しておきたい全体像です。次のセクションからは、生成したデータを使ったGAMLSSによる推定に進みます。

2. BCCG分布による成長曲線推定

　前節で説明したRスクリプトを実行して模擬データと「真の」パーセンタイル曲線のCSVファイルを出力したら、次は今回の本題であるGAMLSSによる成長曲線の推定に進みます。以下のRスクリプトを実行すれば、数値シミュレーションデータ simulated_growth_data_from_smoothed_LMS.csv を読み込み、成長曲線が推定されます。ただし、正しく実行するには、

# 出力先ディレクトリ
DIR <- r"(... ここにファイルパスを指定 ...)"

の部分に、数値シミュレーションデータが存在するフォルダのパスを指定する必要があります（上のスクリプトと共通のパス DIR を指定するということです）。

# =========================================================
# GAMLSS による成長曲線推定（Box-Cox-Cole-Green (BCCG) 分布）
# ---------------------------------------------------------
# BCCG 分布は古典的 LMS 法（Cole & Green, 1992）と対応する。
#   mu    <-> M（中央値）
#   sigma <-> S（変動係数）
#   nu    <-> L（歪度補正べき乗）
#
# 平滑化には P-スプライン pb() を使用し、
# スムージングパラメタは pb() 内で自動推定する。
#
# 【出力】
#   - センタイル曲線プロット（実測値との重ね合わせ）
#   - L / M / S パラメタ推定曲線（真値との比較）
#   - CSV ファイル（推定パラメタ・センタイル値）
# =========================================================

# =========================================================
# 0. パッケージ
# =========================================================
for (p in c("gamlss", "gamlss.dist")) {
  if (!requireNamespace(p, quietly = TRUE)) install.packages(p)
}
library(gamlss)
library(gamlss.dist)

# =========================================================
# 1. 解析設定（ここだけ変更して使う）
# =========================================================

# --- ファイルパス ---
DIR       <- r"(... ここにファイルパスを指定 ...)"
file_sim  <- file.path(DIR, "simulated_growth_data_from_smoothed_LMS.csv")
file_true <- file.path(DIR, "true_model_based_centile_curves.csv")
#   ↑ 真値ファイルが不要なら NULL を指定: file_true <- NULL

# --- 解析対象（"Male"/"Female"、"Height"/"Weight"）---
SEX_TARGET  <- "Male"
ITEM_TARGET <- "Height"

# --- 年齢グリッド（センタイル予測に使用）---
age_grid <- seq(5, 17, by = 0.1)

# --- 描画するセンタイル ---
cent_vals  <- c(0.03, 0.10, 0.25, 0.50, 0.75, 0.90, 0.97) * 100  # パーセント表記
prob_names <- c("p3", "p10", "p25", "p50", "p75", "p90", "p97")

# --- pb() のスムージング設定 ---
# method: "ML", "ML-1", "EM", "GAIC", "GCV" から選択
# k     : GAIC のペナルティ係数（k=2 で AIC 相当）
pb_method <- "GAIC"
pb_k      <- 2

# 各パラメタの自由度上限（大きすぎると過学習）
max_df_mu    <- 10
max_df_sigma <- 8
max_df_nu    <- 6

# --- 真値との比較を描画するか ---
plot_true_curve <- TRUE

# --- CSV 保存 ---
save_csv <- TRUE

# --- 出力ファイル名の接頭辞 ---
tag_out <- paste0(
  tolower(SEX_TARGET),  "_",
  tolower(ITEM_TARGET), "_BCCG_", pb_method
)

# =========================================================
# 2. 補助関数
# =========================================================

# --- データ読み込み ---
# 真値ファイルが存在しない場合は true = NULL を返す
read_analysis_data <- function(file_sim, file_true = NULL) {
  sim_dat  <- read.csv(file_sim, stringsAsFactors = FALSE)
  true_dat <- NULL
  if (!is.null(file_true) && file.exists(file_true)) {
    true_dat <- read.csv(file_true, stringsAsFactors = FALSE)
  }
  list(sim = sim_dat, true = true_dat)
}

# --- sex / item でサブセット（age 昇順）---
subset_group_data <- function(dat, sex, item) {
  out <- dat[dat$sex == sex & dat$item == item, , drop = FALSE]
  if (nrow(out) == 0) stop("指定した sex/item が見つかりません。")
  out <- out[order(out$age), ]
  rownames(out) <- NULL
  out
}

# --- BCCG モデルの当てはめ ---
# mu / sigma / nu をすべて pb(age, ...) でスムージングする
fit_model_bccg <- function(dat_work,
                           pb_method = "GAIC", pb_k = 2,
                           max_df_mu = 10, max_df_sigma = 8, max_df_nu = 6,
                           n_cyc = 200, trace_fit = FALSE) {
  pb_expr <- function(max_df) {
    sprintf("pb(age, method = '%s', k = %g, max.df = %g)",
            pb_method, pb_k, max_df)
  }

  fit <- gamlss(
    formula       = as.formula(paste("value ~", pb_expr(max_df_mu))),
    sigma.formula = as.formula(paste("~",       pb_expr(max_df_sigma))),
    nu.formula    = as.formula(paste("~",       pb_expr(max_df_nu))),
    family  = BCCG(),
    data    = dat_work,
    control = gamlss.control(n.cyc = n_cyc, trace = trace_fit)
  )

  # predict() 呼び出し時の data 再評価エラーを回避
  fit$call$data <- quote(dat_work)
  fit
}

# --- L / M / S パラメタを年齢グリッドで予測 ---
# BCCG: nu = L, mu = M, sigma = S
predict_bccg_params <- function(fit, age_grid) {
  newd <- data.frame(age = age_grid)
  data.frame(
    age = age_grid,
    L   = as.numeric(predict(fit, newdata = newd, what = "nu",    type = "response")),
    M   = as.numeric(predict(fit, newdata = newd, what = "mu",    type = "response")),
    S   = as.numeric(predict(fit, newdata = newd, what = "sigma", type = "response"))
  )
}

# --- センタイル曲線を年齢グリッドで予測 ---
predict_centiles <- function(fit, age_grid, cent_vals, prob_names) {
  mat <- centiles.pred(
    obj     = fit,
    xname   = "age",
    xvalues = age_grid,
    cent    = cent_vals,
    plot    = FALSE,
    data    = fit$data
  )
  # 列名を付けて age 列を上書き（centiles.pred の第1列は xvalues そのもの）
  out        <- as.data.frame(mat)
  names(out) <- c("age", prob_names)
  out$age    <- age_grid
  out
}

# =========================================================
# 3. 描画設定（色・線幅）
# =========================================================

# センタイルごとの色（外側ほど薄色）
centile_cols <- c(
  p3  = "#85B7EB",
  p10 = "#185FA5",
  p25 = "#378ADD",
  p50 = "#001F5B",
  p75 = "#BA7517",
  p90 = "#D85A30",
  p97 = "#EF9F27"
)

# センタイルごとの線幅（中央値を太く）
centile_lwds <- c(1.5, 2.0, 2.0, 3.0, 2.0, 2.0, 1.5)

# =========================================================
# 4. データ読み込みと対象群の抽出
# =========================================================
obj      <- read_analysis_data(file_sim = file_sim, file_true = file_true)
sim_all  <- obj$sim
true_all <- obj$true

dat_work <- subset_group_data(sim_all, sex = SEX_TARGET, item = ITEM_TARGET)

true_target <- NULL
if (!is.null(true_all)) {
  true_target <- subset_group_data(true_all, sex = SEX_TARGET, item = ITEM_TARGET)
}

# =========================================================
# 5. BCCG モデルの当てはめと予測
# =========================================================
cat("BCCG モデルを当てはめています...\n")

fit_bccg <- fit_model_bccg(
  dat_work     = dat_work,
  pb_method    = pb_method,
  pb_k         = pb_k,
  max_df_mu    = max_df_mu,
  max_df_sigma = max_df_sigma,
  max_df_nu    = max_df_nu
)

est_params <- predict_bccg_params(fit_bccg, age_grid)
est_cent   <- predict_centiles(fit_bccg, age_grid, cent_vals, prob_names)

# =========================================================
# 6. CSV 保存
# =========================================================
if (save_csv) {
  write.csv(est_params,
            file.path(DIR, paste0("estimated_params_",    tag_out, ".csv")),
            row.names = FALSE)
  write.csv(est_cent,
            file.path(DIR, paste0("estimated_centiles_", tag_out, ".csv")),
            row.names = FALSE)
  cat("CSV を保存しました:", DIR, "\n")
}

# =========================================================
# 7. 描画（1 画面・2×2 パネル）
#    左上: 成長曲線（実測値 + センタイル）
#    右上: L（歪度）パラメタ
#    左下: M（中央値）パラメタ
#    右下: S（変動係数）パラメタ
# =========================================================
ylab_main <- if (ITEM_TARGET == "Height") "Height (cm)" else "Weight (kg)"

# 真値表示の可否（真値データが存在し、かつ設定が TRUE のとき）
true_for_plot <- if (!is.null(true_target) && plot_true_curve) true_target else NULL

# 真値列を安全に取り出すヘルパー（真値がなければ NULL を返す）
true_col <- function(col_name) {
  if (!is.null(true_for_plot)) true_for_plot[[col_name]] else NULL
}

op <- par(no.readonly = TRUE)
par(mfrow = c(2, 2), mar = c(4, 4.2, 2.8, 1), las = 1)

# ---- パネル 1: 成長曲線 ----
ylim_main <- range(c(dat_work$value,
                     as.matrix(est_cent[, prob_names])), na.rm = TRUE)

plot(dat_work$age, dat_work$value,
     pch = 16, cex = 0.65, col = gray(0.7, alpha = 0.2),
     xlab = "Age (years)", ylab = ylab_main,
     xlim = range(age_grid), ylim = ylim_main,
     main = paste0("BCCG  [", SEX_TARGET, " / ", ITEM_TARGET,
                   "  pb=", pb_method, ", k=", pb_k, "]"))

# 真値センタイル（点線）
if (!is.null(true_for_plot)) {
  for (j in seq_along(prob_names)) {
    lines(true_for_plot$age, true_for_plot[[prob_names[j]]],
          col = 2, lwd = 1.2, lty = 2)
  }
}

# 推定センタイル（実線）
for (j in seq_along(prob_names)) {
  lines(est_cent$age, est_cent[[prob_names[j]]],
        col = centile_cols[j], lwd = centile_lwds[j])
}

legend("topleft",
       legend = c("3%", "10%", "25%", "50%", "75%", "90%", "97%"),
       col = centile_cols, lty = 1, lwd = centile_lwds,
       bg = "white", cex = 0.8)

if (!is.null(true_for_plot)) {
  legend("bottomright",
         legend = c("Estimated centiles", "True centiles"),
         col = c("black", 2),
         lty = c(1, 2), lwd = c(2, 1.5),
         bg = "white", cex = 0.8)
}

# ---- パネル 2–4: L / M / S パラメタ ----
# パネル定義: list(推定値列名, 真値列名, y軸ラベル, タイトル, 凡例位置)
lms_panels <- list(
  list("L", "L", "L  (skewness)",                "L parameter", "bottomright"),
  list("M", "M", paste0("M  (", ylab_main, ")"), "M parameter", "topleft"),
  list("S", "S", "S  (coef. of var.)",            "S parameter", "topright")
)

for (p in lms_panels) {
  y      <- est_params[[p[[1]]]]
  y_true <- true_col(p[[2]])
  ylim   <- range(c(y, y_true), na.rm = TRUE)

  plot(est_params$age, y,
       type = "l", lwd = 2, col = 4,
       xlab = "Age (years)", ylab = p[[3]],
       ylim = ylim, main = p[[4]])

  if (!is.null(y_true)) {
    lines(true_for_plot$age, y_true, lty = 2, lwd = 2, col = 2)
    legend(p[[5]], legend = c("Estimated", "True"), col = c(4, 2),
           lty = c(1, 2), lwd = 2, bg = "white", cex = 0.85)
  }
}

par(op)

2.1 BCCG分布とLMS法の対応

GAMLSSパッケージで使用するBCCG（Box-Cox-Cole-Green）分布は、古典的なLMS法（Cole & Green, 1992）と直接対応しています。具体的には次のとおりです。

LMS法のパラメータ	GAMLSSでの呼称	意味
$L$ （歪度）	`nu`	Box-Cox変換のべき乗パラメータ
$M$ （中央値）	`mu`	分布の中央値
$S$ （変動係数）	`sigma`	分布の広がりを表す変動係数

　LMS法を「GAMLSSの枠組みで実装したもの」がBCCGモデルと理解すると、両者の関係が整理しやすくなります。古典的なLMS法が「3つのパラメータを別々にスムージングして成長曲線を得る」という手順を踏むのに対し、GAMLSSでは尤度に基づく統一的な枠組みのもとで3つのパラメータを同時に推定します。

3.2 推定アルゴリズムの全体像

　このスクリプトは、大きく次の5つのステップで動作します。

① データの読み込みと対象群の抽出

　前節で出力した simulated_growth_data_from_smoothed_LMS.csv（模擬観測データ）と true_model_based_centile_curves.csv（真のパーセンタイル曲線）を読み込みます。解析対象の性別（SEX_TARGET）と測定項目（ITEM_TARGET）を指定することで、4グループ（男女×身長・体重）のうち任意の1グループを抽出して解析できます。

② BCCGモデルの当てはめ

　gamlssパッケージの gamlss() 関数を用いて、次の3つのパラメータを同時に年齢の関数として推定します。

$\displaystyle{ \mu(t) = f_1(t), \quad \sigma(t) = f_2(t), \quad \nu(t) = f_3(t) }$

　ここで $\displaystyle{ f _ 1, f _ 2, f _ 3 }$ はいずれも年齢 $t$ に対するP-スプライン（pb(age)）です。単純な線形回帰や多項式回帰ではなく、データの局所的な構造に柔軟に追従できる非線形スムーザーを使うことで、思春期の急成長など年齢に対して非線形に変化するパターンも適切に捉えられます。

　スムージングの滑らかさは pb() 内で自動推定されます。推定方法は pb_method で指定でき（"GAIC"・"ML" など）、"GAIC" の場合はペナルティ係数 k（デフォルト $k=2$ でAIC相当）のもとで最適な自由度が選ばれます。各パラメータの自由度の上限（max_df_mu、max_df_sigma、max_df_nu）は過学習を防ぐために設定します。

③ パラメータとパーセンタイル曲線の予測

　当てはめたモデルから、0.1歳刻みの年齢グリッドで $L$ 、 $M$ 、 $S$ の推定値と、対象パーセンタイル（3, 10, 25, 50, 75, 90, 97パーセンタイル）の曲線を予測します。パーセンタイル曲線の計算には centiles.pred() 関数を使用します。

④ CSVファイルの保存

　推定されたパラメータ曲線（estimated_params_*.csv）とパーセンタイル曲線（estimated_centiles_*.csv）を指定のディレクトリに保存します。ファイル名には性別・測定項目・スムージング方法が自動的に付加されるため、複数の条件で解析を実行しても出力ファイルが上書きされません。

⑤ 結果の可視化

　2行×2列の4パネルで推定結果を描画します。各パネルの内容は次のとおりです。

左上（成長曲線）：模擬観測データの散布図に、推定されたパーセンタイル曲線（実線）と真のパーセンタイル曲線（点線）を重ね描きします。推定曲線が真の曲線にどれだけ一致しているかを視覚的に確認できます。
右上・左下・右下（L・M・Sパラメータ）：各パラメータの推定曲線（青実線）と真の曲線（赤点線）を年齢に対してプロットします。どのパラメータの推定精度が高く、どのパラメータに乖離が生じやすいかを直接比較できます。

2.3 推定アルゴリズムのポイント

　この推定手法で重要なのは、「分布そのものを年齢の関数として推定している」という点です。

　通常の回帰分析では「平均値がどのように変化するか」を推定しますが、LMS法・BCCGモデルでは中央値（ $M$ ）だけでなく、分布の広がり（ $S$ ）と歪み（ $L$ ）もすべて年齢に依存して変化するものとして同時に推定します。これによって、成長に伴う分布形状の変化——たとえば思春期に体重分布の歪みが大きくなる傾向——も適切にモデルに取り込むことができます。

2.4 スクリプトの主な関数と役割

　スクリプト内で定義されている関数の役割をまとめます。

関数名	役割
`read_analysis_data()`	CSVを読み込み、真値ファイルが存在しない場合は `NULL` を返す
`subset_group_data()`	性別・測定項目でサブセットし、年齢昇順に整列する
`fit_model_bccg()`	GAMLSSでBCCGモデルを当てはめる（P-スプラインで3パラメータを同時推定）
`predict_bccg_params()`	当てはめたモデルから年齢グリッドでL・M・Sを予測する
`predict_centiles()`	`centiles.pred()` を用いてパーセンタイル曲線を年齢グリッドで予測する

2.5 いろいろ試してみよう

　このスクリプトでは、GAMLSSで推定した曲線を「真の曲線」と重ね合わせて比較できるようになっています。せっかくなので、パラメータをいろいろ変えながら、成長曲線推定の感覚を養ってみましょう。

　たとえば、データ生成スクリプトの n_sim を 10,000 から 500 や 100 に減らしてみると、サンプルサイズが小さいときに推定曲線が真値からどれだけ乖離するかを直接体験できます。また、spar_L、spar_M、spar_S の値を変えてスムージングの強弱を調整したり、推定スクリプトの max_df_mu などの自由度上限を変えたりすることで、過学習・過平滑化がパーセンタイル曲線の形状にどう影響するかも確認できます。

　「分布 → 個体データ → 分布の推定」という往復の構造を手を動かして体験することが、LMS法とGAMLSSの理解への一番の近道です。

3. まとめ：分布モデルを仮定するパラメトリック推定は正しい？

　今回は、成長曲線推定においてもっとも広く使われているアプローチを実際に動かしてみました。しかし、問題はここからです。果たしてこの方法はどの程度正しいのか——それが私の率直な疑問です。時間はかかると思いますし、激しい反論もあるかもしれませんが、この方法の限界と欠点を明らかにし、世界に示していきたいと考えています。もちろん実用上は、今回紹介したアプローチで十分に役立つ場面が多いでしょう。しかし、データの構造を理解する科学としては、ここで立ち止まるのではなく、新たな冒険に踏み出す意義があるのです。