今回は、下図のような身長や体重の成長曲線の推定についてのお話です。
余談ですが、最近、医学・生理学の分野で気になることがあります。数理的な手法が、いつの間にか「魔術」扱いされているという現象です。私自身にとって、数式は言葉の一つにすぎません。丁寧に読めば、むしろ理解を助けてくれる便利な道具です。ところが、数式が苦手な方の目には「怪しげな呪文」に映るようで、どれだけ丁寧な分析結果を示しても「なんか騙されてる気がする…」と疑いから入られてしまうことがあります。あるいは逆に、中身を理解しないまま「呪文（指標）」を盲目的に信じてしまう。この残念な状況を、私は皮肉を込めて「数理の魔術化」と呼んでいます。 もちろん、「数学を完璧にマスターしてください」などと言うつもりはありません。それでも——お互いが少しずつ歩み寄ることはできるはずです。数理の側が「わかりやすく伝える努力」をし、医学・生理学の側が「とりあえず話を聞いてみる姿勢」を持つ。その小さな積み重ねが、思いのほか大きな進歩につながると感じています。「なんとなくわかった、面白そう」 ——そのくらいの距離感から始まる共同作業が、案外、分野全体を動かす力になるものです。ということで、今回からしばらく、子どもの身体計測データの数理を「魔術」ではなく「科学」として見直す試みをしていきたいと思います。
子どもの身長、体重、BMI（Body Mass Index：体重[kg]÷身長[m]²で計算される肥満度の指標）は、年齢とともに変化するだけでなく、そのデータの分布構造（散らばり方）も変わっていきます。特に体重やBMIのデータ、あるいは、成長が急加速する「成長スパート」の始まりや終わりの時期の身長データでは、分布が左右非対称にゆがむ（skew）ことが知られています。たとえば子どもの体重やBMIの分布では、低い値の側には明確な下限（0以下にはなれない、生命維持に最低限必要な下限）がある一方、肥満側の高い値には上限がなく、極端に高い値を持つ例が存在する可能性があります。そのため、分布のすそ野は右側に長く引っ張られた形、つまり右にゆがんだ分布になります。
統計の教科書では正規分布（左右対称な釣り洋鐘型の分布）がよく登場し、その性質は多くの研究者になじみのあるものです。そこで、ゆがんだ分布を正規分布に近づけられれば、正規分布の解析法を使えるので便利です。そのようなデータ分布の形状を変換するための有名な方法としてBox-Cox変換があります。これは、データに対数（log）やべき乗（累乗）などの数学的操作を施すことで、ゆがみを補正し、正規分布に近い形に整える手法です。
LMS法（Lambda-Mu-Sigma法）は、このBox-Cox変換と類似した方法です。子どもの成長データのように「年齢によって分布のゆがみ方が変化する」場合にも対応できるよう工夫されており、成長曲線のパーセンタイル（「同じ年齢の子どもの中で上から・下から何％に位置するか」を示す値）を推定する際に広く使われています。
本記事では、このLMS法の基本的な考え方を解説するとともに、その限界や注意点についても批判的に考察します。
1. LMS法の基本アイデア
LMS法は、年齢ごとに分布の形を以下の3つのパラメタで表現します：
- L（Lambda）：歪度（skewness）を補正するためのBox-Cox変換パラメタ
- M（Mu）：中央値（median）
- S（Sigma）：変動係数（coefficient of variation）
つまり、ある年齢 ( t ) におけるデータ分布を
という3つの関数で表現します。
そして、LMS法は、次の2つの要素から構成されます。
(1) 固定された年齢における分布の正規化（Box-Cox的側面） (2) 年齢依存のパラメタを滑らかな関数として推定（スプライン的側面）
まずは、これらの側面を個別に説明します。
(1) 固定された年齢における分布の正規化（Box-Cox的側面）
まず、ある特定の年齢 を固定して考えます。このとき、データ
に対して
という変換を行います。 この変換は、Box-Cox変換と同じべき乗変換（対数変換を含む）に基づいており、 分布の歪みを補正し、標準正規分布に近づけることを目的としています。 したがって、この段階では
各年齢ごとに分布を正規分布に近づける変換を行っている
と理解できます。
分布の正規化をRで体験
分布の正規化を理解するために、以下のRスクリプトを実行してみてください。このスクリプトでは、LMSモデルに適合する歪んだデータを生成し、そのデータからLMSパラメタ を最尤法で推定します。その後、推定されたパラメタを用いてデータを正規化し、ヒストグラムとWorm plotで結果を確認します。一般的なLMS法では、LMSパラメタの年齢依存性も同時に推定しますが、ここでは、年齢を固定（L, M, Sを定数）としてあります。その意味で、このデモは、Box-Cox型の変換に標準化を組み合わせた変換を行っています。
############################################################
# LMS Method Demo (Fixed Age, Parameter Estimation)
#
# Purpose:
# - Generate skewed data using LMS model
# - Estimate LMS parameters via maximum likelihood
# - Normalize data using estimated LMS parameters
# - Diagnose normality using Worm plot
#
# Output:
# - 2×2 plots (Histogram + Worm plot)
# - Summary of parameter estimation and statistics
############################################################
rm(list = ls())
# ==========================================================
# 1. True LMS parameters (for simulation)
# ==========================================================
L_true <- -0.8 # skewness parameter
M_true <- 60 # median
S_true <- 0.2 # coefficient of variation
n <- 1000
set.seed(1)
# ==========================================================
# 2. Generate skewed data
# ==========================================================
z_true <- rnorm(n)
y <- if (abs(L_true) > 1e-8) {
M_true * (1 + L_true * S_true * z_true)^(1 / L_true)
} else {
M_true * exp(S_true * z_true)
}
stopifnot(all(y > 0)) # LMS requires positive data
# ==========================================================
# 3. Utility functions
# ==========================================================
skewness <- function(x) {
mean((x - mean(x))^3) / sd(x)^3
}
lms_transform <- function(y, L, M, S) {
if (abs(L) > 1e-8) {
((y / M)^L - 1) / (L * S)
} else {
log(y / M) / S
}
}
# ==========================================================
# 4. Maximum likelihood estimation of LMS parameters
# ==========================================================
negloglik_lms <- function(par, y) {
L <- par[1]
M <- exp(par[2])
S <- exp(par[3])
if (!is.finite(L) || !is.finite(M) || !is.finite(S)) return(1e12)
z <- if (abs(L) > 1e-8) {
((y / M)^L - 1) / (L * S)
} else {
log(y / M) / S
}
if (any(!is.finite(z))) return(1e12)
log_jacobian <- -log(S) + (L - 1) * log(y) - L * log(M)
log_lik <- sum(dnorm(z, log = TRUE) + log_jacobian)
return(-log_lik)
}
fit <- optim(
par = c(-0.5, log(median(y)), log(sd(y) / mean(y))),
fn = negloglik_lms,
y = y,
method = "BFGS",
control = list(maxit = 2000)
)
L_hat <- fit$par[1]
M_hat <- exp(fit$par[2])
S_hat <- exp(fit$par[3])
# ==========================================================
# 5. LMS transformation
# ==========================================================
z_lms <- lms_transform(y, L_hat, M_hat, S_hat)
# ==========================================================
# 6. Worm plot function
# ==========================================================
worm_components <- function(x, standardized = FALSE) {
x_std <- if (standardized) x else as.numeric(scale(x))
x_std <- sort(x_std)
n <- length(x)
qt <- qnorm((seq_len(n) - 0.5) / n)
sd_x <- sd(x_std)
dev <- x_std - sd_x * qt
z_plot <- seq(-4, 4, length.out = 500)
SE <- sd_x * sqrt(pnorm(z_plot) * (1 - pnorm(z_plot)) / n) / dnorm(z_plot)
list(qt = qt, dev = dev, z_plot = z_plot, SE = SE)
}
wc_before <- worm_components(y, FALSE)
wc_after <- worm_components(z_lms, TRUE)
# ==========================================================
# 7. Visualization
# ==========================================================
par(mfrow = c(2, 2), mar = c(5, 5, 3, 2))
# Histogram (before)
hist(y, breaks = 30, col = 2, border = "white",
main = "Before LMS transform", xlab = "y")
# Histogram (after)
hist(z_lms, breaks = 30, col = 4, border = "white",
main = "After LMS transform", xlab = "z")
# Worm plot (before)
plot(wc_before$qt, wc_before$dev,
type = "n", xlim = c(-4, 4), ylim = c(-2, 2),
main = "Worm Plot (Before)",
xlab = "Normal Quantile", ylab = "Deviation")
grid()
lines(wc_before$z_plot, 1.96 * wc_before$SE, lty = 2)
lines(wc_before$z_plot, -1.96 * wc_before$SE, lty = 2)
abline(h = 0); abline(v = 0, lty = 3)
points(wc_before$qt, wc_before$dev, col = 2, pch = 16)
# Worm plot (after)
plot(wc_after$qt, wc_after$dev,
type = "n", xlim = c(-4, 4), ylim = c(-2, 2),
main = "Worm Plot (After)",
xlab = "Normal Quantile", ylab = "Deviation")
grid()
lines(wc_after$z_plot, 1.96 * wc_after$SE, lty = 2)
lines(wc_after$z_plot, -1.96 * wc_after$SE, lty = 2)
abline(h = 0); abline(v = 0, lty = 3)
points(wc_after$qt, wc_after$dev, col = 4, pch = 16)
# ==========================================================
# 8. Summary (single output)
# ==========================================================
param_table <- data.frame(
Parameter = c("L", "M", "S"),
True = round(c(L_true, M_true, S_true), 4),
Estimated = round(c(L_hat, M_hat, S_hat), 4)
)
result_table <- data.frame(
Condition = c("Before", "After"),
Mean = round(c(mean(y), mean(z_lms)), 4),
SD = round(c(sd(y), sd(z_lms)), 4),
Skewness = round(c(skewness(y), skewness(z_lms)), 4)
)
text1 <- paste(capture.output(print(param_table, row.names = FALSE)), collapse = "\n")
text2 <- paste(capture.output(print(result_table, row.names = FALSE)), collapse = "\n")
full_output <- paste0(
"\n========================================\n",
"LMS Method Summary\n",
"========================================\n\n",
"[Parameter estimation]\n", text1, "\n\n",
"[Distribution summary]\n", text2, "\n\n",
"Optimization status: ", ifelse(fit$convergence == 0, "OK", "Check"), "\n",
"========================================\n"
)
cat(full_output)このスクリプトを実行すれば、下のような図が出力されます。
この図のように、変換前の分布（図左上）は非対称ですが、LMS法を用いた変換後は左右対称に近い形（図右上）になります。この変化は、尺度の変換だけでなく、分布の歪みが補正されていることを意味します。図の下段のWorm plotは、正規分布からのずれを確認するための図です。変換前のプロットでは、点が系統的な曲線パターンを示し、分布に歪みがあることが示されます。変換後では、点が0付近に沿って並び、正規分布に近づいていることが確認できます。
さらに、Rスクリプトの実行結果として出力される Parameter estimation の表（以下）では、データ生成に用いた真のパラメタと推定値が比較されます。
======================================== LMS Method Summary ======================================== [Parameter estimation] Parameter True Estimated L -0.8 -0.7724 M 60.0 59.8964 S 0.2 0.2067 [Distribution summary] Condition Mean SD Skewness Before 62.4104 14.5338 1.7138 After 0.0000 1.0005 -0.0014 Optimization status: OK ========================================
このように、推定値は真値と近い値になります。このことから、観測データがべき乗変換によって正規分布に従うデータに変換できるのであれば、 を正しく推定できることが確認できます。ここで、
は分布の中心、
はばらつき、
は歪みの補正に対応しています。
次に、Distribution summary の表を見てみると、変換前のデータ では、歪度が0から離れており、分布が対称ではありません。一方、LMS変換後の
では、平均が0付近、標準偏差が1付近となり、歪度も0に近づきます。これは、LMS変換によって分布が標準正規分布に近づいたことを示しています。
この例では年齢を固定していますが、LMS法では通常、パラメタは年齢の関数
として推定されます。このスクリプトは、その基本的な変換と推定の流れを単一の分布で示したものです。
以上より、LMS法では、歪んだ分布を3つのパラメタで表現し、それを用いて正規分布に近い尺度へ変換できることが確認できます。
(2) 年齢依存のパラメタを滑らかな関数として推定（スプライン的側面）
次に重要なのは、LMS法では
を単に年齢ごとに独立に推定するのではなく、
年齢
に依存した 滑らかな関数 として推定する
という点です。
具体的には、例えば は
のように、基底関数 の線形結合として表されます。同様に
,
も滑らかな関数として表現されます。
このようにすることで、年齢ごとのばらつきを直接追うのではなく、分布パラメタの滑らかな変化を仮定することになります。この「滑らかさ」の導入方法には、実装上の違いがあります。
■ Cole–Green法（元のLMS法）
ColeとGreenが1992年に提案した方法では、 を年齢の関数として、主に3次スプラインで表現し、ペナルティ付き尤度により滑らかさを制御します。
■ GAMLSS（Generalized Additive Models for Location, Scale and Shape）
最近では、GAMLSSという方法が使われるようになっているようです。この方法では、（あるいは
）を説明変数の関数として扱い、Pスプラインなどを用いた加法モデルとして記述します。
滑らかな関数表現をRで体験
上の理論的説明を具体的に理解するために、Rを用いて「滑らかな関数として推定する」という考え方を実際に体験します。ここでは、LMSの推定そのものではなく、単一の応答変数 (y) を説明変数 (x) の滑らかな関数として推定する問題に単純化しています。
まず、不等間隔の数値 を生成し、それに対して成長曲線のような単調増加関数を真の関数として設定し、観測誤差を加えたデータを作成します。これは、実際の身長や体重のように「年齢に対して滑らかに変化する量」を模したものです。
############################################################
# smoothing_spline_comparison.R
#
# 不等間隔の小数点データに対する平滑化スプラインの比較
#
# 比較手法：
# 1) 3次平滑化スプライン（Cubic Smoothing Spline）
# smooth.spline() を使用。Cole & Green (1992) の
# LMS法で用いられる3次スプライン型の平滑化と
# 本質的に同じ考え方を持つ古典的手法
#
# 2) Pスプライン（P-spline）
# Eilers & Marx (1996) の方法に基づく自作実装。
# Bスプライン基底 + 差分ペナルティによる正則化最小二乗法。
#
# 参考文献：
# Cole, T.J. and Green, P.J. (1992). Smoothing reference
# centile curves. Statistics in Medicine, 11, 1305-1319.
# Eilers, P.H.C. and Marx, B.D. (1996). Flexible smoothing
# with B-splines and penalties. Statistical Science, 11(2),
# 89-102.
#
# 動作確認済み環境：R >= 4.0、パッケージ splines（標準添付）
############################################################
# ==========================================================
# 0. データ生成
# ==========================================================
set.seed(123)
n <- 30
# 不等間隔の x を一様分布から生成してソート
x <- sort(runif(n, min = 0.2, max = 18.0))
# 真の滑らかな関数（成長曲線らしい単調増加関数）
f_true <- function(x) {
50 + 120 / (1 + exp(-(x - 9) / 2.2))
}
y_true <- f_true(x)
y <- y_true + rnorm(n, sd = 3.0) # 誤差は少し大きめでも自然
# 描画・予測用の等間隔グリッド
xg <- seq(min(x), max(x), length.out = 500)
yg_true <- f_true(xg)
# ==========================================================
# 1. 3次平滑化スプライン（Cubic Smoothing Spline）
# ==========================================================
#
# smooth.spline() は3次スプラインにクロスバリデーション
# またはパラメータ spar で滑らかさを制御する。
# spar = 0 に近いほど補間的、1 に近いほど過度に平滑化される。
# ----------------------------------------------------------
fit_cubic <- smooth.spline(x = x, y = y, spar = 0.7)
pred_cubic <- predict(fit_cubic, x = xg)$y
# ==========================================================
# 2. Pスプライン（P-spline）の自作実装
# ==========================================================
#
# モデル： y ≈ B %*% a
# B : Bスプライン基底行列（splines::bs() で構築）
# a : 係数ベクトル
#
# 目的関数：
# min_a || y - B a ||^2 + lambda * || D a ||^2
# ^^^^^^^^^^^^ ^^^^^^^^^^^^^^^^^
# 残差二乗和 2階差分ペナルティ
#
# ペナルティ項 lambda * ||Da||^2 が隣接係数の滑らかさを強制する。
# lambda が大きいほど滑らかな推定曲線が得られる。
#
# 解析解（正規方程式）：
# a_hat = (B'B + lambda * D'D)^{-1} B'y
# ----------------------------------------------------------
library(splines) # Bスプライン基底生成に使用（標準パッケージ）
#' Pスプライン平滑化
#'
#' @param x 説明変数（数値ベクトル）
#' @param y 目的変数（数値ベクトル）
#' @param x_new 予測を行いたい新しい x 値（NULL の場合は予測なし）
#' @param nseg 内部ノット間隔数の目安（基底関数の個数 ≈ nseg + degree）
#' @param degree Bスプラインの次数（3 = 3次スプライン）
#' @param lambda 差分ペナルティの強さ（大きいほど滑らか）
#'
#' @return 係数・あてはめ値・予測値などを格納したリスト
fit_pspline <- function(x, y,
x_new = NULL,
nseg = 16,
degree = 3,
lambda = 1) {
xr <- range(x)
# 内部ノットを等間隔に配置（Pスプラインでは密なノットが推奨）
knots_inner <- seq(xr[1], xr[2], length.out = nseg + 1)
knots_inner <- knots_inner[-c(1, length(knots_inner))] # 端点を除く
# --- Bスプライン基底行列の構築 ---
B <- bs(
x,
knots = knots_inner,
degree = degree,
intercept = TRUE, # 切片項を含む
Boundary.knots = xr
)
K <- ncol(B) # 基底関数の総数
# --- 2階差分ペナルティ行列 D（サイズ: (K-2) × K） ---
# D %*% a は係数の2階差分を返す。
# ペナルティ項 D'D は係数の2次変化（曲率）を罰する。
D <- diff(diag(K), differences = 2)
# --- 正則化付き正規方程式を解く ---
BtB <- crossprod(B) # B'B
Bty <- crossprod(B, y) # B'y
P <- crossprod(D) # D'D（ペナルティ行列）
a_hat <- solve(BtB + lambda * P, Bty) # 最小化の解
fitted <- as.vector(B %*% a_hat)
# --- 出力リストの作成 ---
out <- list(
coef = a_hat,
fitted = fitted,
x = x,
y = y,
knots_inner = knots_inner,
degree = degree,
Boundary.knots = xr,
lambda = lambda,
B = B
)
# x_new が与えられた場合は新しい x での予測値も計算
if (!is.null(x_new)) {
B_new <- bs(
x_new,
knots = knots_inner,
degree = degree,
intercept = TRUE,
Boundary.knots = xr
)
out$y_new <- as.vector(B_new %*% a_hat)
out$x_new <- x_new
out$B_new <- B_new
}
return(out)
}
# --- Pスプライン適合 ---
fit_p <- fit_pspline(x = x, y = y, x_new = xg,
nseg = 12, degree = 3, lambda = 0.8)
pred_p <- fit_p$y_new
# ==========================================================
# 3. Bスプライン基底の可視化ユーティリティ
# ==========================================================
#
# Pスプラインの動作原理を理解するための補助的プロット。
# 各基底関数は局所的なサポートを持ち、その線形結合として
# 推定曲線が表現されることを確認できる。
# ----------------------------------------------------------
#' Bスプライン基底関数を描画する
#'
#' @param x データ点（range の決定に使用）
#' @param nseg 区間数
#' @param degree Bスプラインの次数
#' @param ngrid 描画用グリッド点の数
plot_bspline_basis <- function(x,
nseg = 16,
degree = 3,
ngrid = 400) {
xr <- range(x)
knots_inner <- seq(xr[1], xr[2], length.out = nseg + 1)
knots_inner <- knots_inner[-c(1, length(knots_inner))]
xg <- seq(xr[1], xr[2], length.out = ngrid)
Bg <- bs(
xg,
knots = knots_inner,
degree = degree,
intercept = TRUE,
Boundary.knots = xr
)
matplot(
xg, Bg,
type = "l", lty = 1, lwd = 1,
xlab = "x", ylab = "基底関数の値",
main = sprintf("Bスプライン基底関数 (次数 %d, %d 基底)",
degree, ncol(Bg))
)
abline(v = knots_inner, lty = 3, col = "gray50")
mtext("縦破線：内部ノット位置", side = 1, line = 3, cex = 0.8)
}
# ==========================================================
# 4. 結果の描画（2×2 レイアウト）
# ==========================================================
op <- par(no.readonly = TRUE) # グラフィクス設定を退避
par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
# ---- (1) サンプルデータと真の関数 ----------------------
plot(
x, y,
pch = 16, cex = 1, col = gray(0.6),
xlab = "x", ylab = "y",
main = "サンプルデータ"
)
lines(xg, yg_true, lwd = 3, col = 4)
legend("topleft",
legend = c("真の関数", "観測値"),
lty = c(1, NA), pch = c(NA, 16),
lwd = c(3, NA), col = c(4, gray(0.6)),
bty = "o")
# ---- (2) 3次平滑化スプライン ---------------------------
plot(
x, y,
pch = 16, cex = 1, col = gray(0.6),
xlab = "x", ylab = "y",
main = "3次平滑化スプライン (spar = 0.7)"
)
lines(xg, yg_true, lwd = 3, lty = 2, col = 4)
lines(xg, pred_cubic, lwd = 3, col = 2)
legend("topleft",
legend = c("真の関数", "平滑化曲線"),
lty = c(2, 1), lwd = c(3, 3), col = c(4, 2),
bty = "o")
# ---- (3) Pスプライン -----------------------------------
plot(
x, y,
pch = 16, cex = 1, col = gray(0.6),
xlab = "x", ylab = "y",
main = "Pスプライン (nseg = 12, lambda = 0.8)"
)
lines(xg, yg_true, lwd = 3, col = 4, lty = 2)
lines(xg, pred_p, lwd = 3, col = 2)
abline(v = fit_p$knots_inner, lty = 3, col = "gray60")
legend("topleft",
legend = c("真の関数", "Pスプライン", "内部ノット"),
lty = c(2, 1, 3),
lwd = c(3, 3, 1),
col = c(4, 2, "gray60"),
bty = "o",
bg = "white")
# ---- (4) Bスプライン基底関数 ---------------------------
plot_bspline_basis(x, nseg = 12, degree = 3)
par(op) # グラフィクス設定を復元
上の図に示したように、本コードでは、このデータに対して次の2種類の平滑化を行っています。
(1) 3次平滑化スプライン（Cubic smoothing spline）
smooth.spline() により、3次スプライン関数を用いてデータを近似します。この方法では、
のような目的関数を最小化しており、曲線の曲率（2階微分）に対するペナルティによって滑らかさが制御されます。
これは、Cole–Green法において「3次スプラインで (L(t), M(t), S(t)) を表現し、ペナルティ付き尤度で滑らかさを制御する」という枠組みと本質的に類似した構造を持っています。
(2) Pスプライン（P-spline）
一方、Pスプラインでは、関数を
と Bスプライン基底の線形結合で表現し、係数ベクトル に対して
という形のペナルティを課します。ここで は2階差分行列であり、隣接する係数の変化（離散的な曲率）を抑える役割を持ちます。
このように、Pスプラインは「多くの局所基底を用いて関数を表現し、その係数に対して滑らかさ制約を課す」という構造になっており、GAMLSSで広く用いられている方法です。一方で、GAMLSSパッケージではこのようなPスプラインだけでなく、Cole–Green法に対応する3次スプライン型の平滑化（例：cs() など）も選択することができ、目的やデータ構造に応じて平滑化手法を柔軟に使い分けることが可能です。
以下のように続けると、理論 → 実装 → アルゴリズム、という流れが自然につながります。
2. 実際のアルゴリズムの考え方
実際のLMS法のアルゴリズムでは、正規化とパラメタの平滑化を個別に処理するのではなく、すべてを同時に最適化します。具体的には、LMS法やGAMLSSでは、パラメタ をそれぞれスプラインで表現し、それらを用いて定義される全体の対数尤度（あるいはペナルティ付き尤度）を反復的に最大化します。すなわち、
を同時に考え、その合計が最もよくなるように、すべてのスプライン係数を含むパラメタを反復計算により同時に最適化します。
アルゴリズムのポイント
実際の計算の流れは、次のように整理できます。
① モデルの設定
各パラメタ（例：）をスプライン基底で表現する。
② 目的関数の構築
データに対する対数尤度に、滑らかさのペナルティを加える。
（例：3次スプラインでは 、Pスプラインでは差分）
③ 同時最適化
すべての係数（ など）をまとめて最適化するために、逐次更新（例：IRLSやbackfitting）による反復計算によって解を求める。
④ 収束判定
パラメタの変化が十分小さくなった時点で終了。
以上のように、LMS法やGAMLSSのアルゴリズムでは、単に「各年齢で別々に推定する」のではなく、
分布の形状とその滑らかさを同時に最適化する
という点が特徴になります。その結果として、全体として一貫した滑らかな成長曲線（あるいは分布パラメタの変化）が得られます。
3. まとめ
ここまで、LMS法およびGAMLSSのアルゴリズムの考え方を整理してきました。これらの手法は、成長データのような年齢依存的な分布の変化を柔軟にモデル化できる点で優れており、数千から数万例程度のデータを用いて小児の成長曲線を作成する場面など、多くの場合に有効に機能します。
しかし、LMS法やGAMLSSが「最善の方法である」と主張したいわけではありません。どのような統計手法にも前提・仮定があり、その限界を正しく理解したうえで使用することが重要です。
LMS法の中心的な仮定は、べき乗変換（Box-Cox変換）によってデータを正規分布に近似できることです。データ数が限られている場合にはこの仮定が実用上許容されることも多いですが、数百万人規模の大規模データを用いると、わずかな分布の歪みや裾の厚さも統計的に明確に検出されるようになります。たとえば、身長分布は成長スパートの開始期・終了期において非対称性が生じやすく、また一般集団では裾がやや厚い分布となる傾向があります。体重やBMIにいたっては、変換によって正規分布に十分近似することはさらに困難です。このような状況では、正規分布への変換を前提とするLMS法の信頼性には、注意が必要な制約があります。。
GAMLSSはより柔軟な分布族を扱える点でLMS法を拡張していますが、適切な分布族の選択が結果を大きく左右するため、大規模データにおける複雑な分布形状すべてに対応できるとは限りません。
統計手法を適切に使うためには、「何を仮定しているか」「どのような場合に限界があるか」を常に意識することが求められます。データの規模や性質、そして研究目的に照らして、手法の有効性と限界を見極めたうえで適切に選択・適用することが大切です。
付録：成長曲線を描画するRスクリプト
以下のRスクリプトを実行すれば、この記事の最初に掲載した成長曲線が描けます。
# =========================================================
# 日本人男女の身長・体重パーセンタイル成長曲線
#
# 【概要】
# 文部科学省・厚生労働省の学校保健統計調査などに基づく
# 日本人小児（5〜17歳）の身長・体重パーセンタイル値を
# 三次スプライン補間で滑らかな成長曲線として描画する。
#
# 【出力レイアウト】
# 2行2列の4パネル構成
# 左上 : 男性 身長 右上 : 女性 身長
# 左下 : 男性 体重 右下 : 女性 体重
#
# 【パーセンタイルの読み方】
# 例）50パーセンタイル（p50）＝同年齢の子どもを身長順に
# 並べたとき、ちょうど真ん中（中央値）に位置する値。
# 3〜97パーセンタイルの範囲が正常域の目安として使われる。
#
# 【色の凡例】
# オレンジ系（暖色）: 高身長・高体重側（75〜97%）
# 黒 : 中央値（50%）
# ブルー系（寒色） : 低身長・低体重側（3〜25%）
# =========================================================
# =========================================================
# 1. 年齢ベクトル
# =========================================================
# 実測データの年齢（整数値：5〜17歳）
age <- 5:17
# スプライン補間後の滑らかな描画用年齢軸
# （5〜17歳を500等分した細かい刻み → 曲線が滑らかに見える）
fine_age <- seq(5, 17, length.out = 500)
# =========================================================
# 2. 身長データ（男性） 単位: cm
#
# 各ベクトルは age（5〜17歳）に対応する13個の値。
# p3〜p97 はそれぞれ3・10・25・50・75・90・97パーセンタイル。
# =========================================================
male_h_p3 <- c(101.6,107.4,113.1,118.1,123.3,128.2,133.1,138.7,146.1,153.3,157.5,158.8,159.4)
male_h_p10 <- c(104.3,110.2,116.0,121.3,126.5,131.5,136.9,143.4,151.2,157.5,161.0,162.2,162.9)
male_h_p25 <- c(107.1,113.1,119.0,124.5,129.9,135.1,140.9,148.3,156.1,161.7,164.5,165.8,166.5)
male_h_p50 <- c(110.3,116.4,122.4,128.2,133.7,139.2,145.6,153.7,161.2,166.1,168.5,169.7,170.5)
male_h_p75 <- c(113.5,119.8,126.0,131.9,137.6,143.5,150.7,159.1,166.1,170.4,172.5,173.7,174.5)
male_h_p90 <- c(116.5,122.9,129.2,135.4,141.3,147.6,155.5,164.1,170.3,174.1,176.1,177.3,178.1)
male_h_p97 <- c(119.4,126.1,132.6,138.9,145.0,151.8,160.6,169.0,174.3,177.7,179.7,180.9,181.6)
# =========================================================
# 3. 身長データ（女性） 単位: cm
#
# 女性は思春期の身長スパートが男性より早く（およそ10〜12歳頃）、
# 14歳以降は成長がほぼ止まる傾向がデータに現れている。
# =========================================================
female_h_p3 <- c(100.7,106.5,112.0,117.2,122.5,128.0,134.3,140.9,144.6,146.2,146.8,147.5,147.6)
female_h_p10 <- c(103.4,109.3,114.9,120.3,125.9,132.0,138.7,144.6,147.9,149.4,149.9,150.5,150.8)
female_h_p25 <- c(106.2,112.2,118.0,123.6,129.4,136.0,142.9,148.4,151.2,152.6,153.2,153.7,154.0)
female_h_p50 <- c(109.4,115.5,121.4,127.3,133.6,140.6,147.5,152.4,154.8,156.2,156.8,157.2,157.7)
female_h_p75 <- c(112.6,118.9,125.0,131.3,137.9,145.3,151.9,156.2,158.5,159.9,160.5,160.9,161.3)
female_h_p90 <- c(115.6,122.0,128.3,134.9,142.0,149.6,155.7,159.6,161.7,163.1,163.8,164.3,164.7)
female_h_p97 <- c(118.5,125.2,131.6,138.7,146.3,153.9,159.3,162.9,165.0,166.4,167.0,167.8,168.0)
# =========================================================
# 4. 体重データ（男性） 単位: kg
#
# 高パーセンタイル（p75〜p97）は思春期以降に急増する傾向があり、
# 肥満評価の際は同年齢・同性の分布との比較が重要。
# =========================================================
male_w_p3 <- c(15.0,16.5,18.4,20.2,22.2,24.4,26.9,30.5,35.0,39.7,43.4,45.5,46.8)
male_w_p10 <- c(16.0,17.7,19.7,21.8,24.2,26.8,29.8,33.9,38.7,43.4,47.1,49.1,50.5)
male_w_p25 <- c(17.2,19.0,21.3,23.8,26.6,29.6,33.2,38.0,43.1,47.8,51.5,53.3,54.9)
male_w_p50 <- c(18.6,20.9,23.5,26.5,29.9,33.6,38.1,43.6,48.9,53.5,57.4,58.9,60.7)
male_w_p75 <- c(20.4,23.1,26.3,30.0,34.3,38.9,44.3,50.6,56.0,60.6,64.8,65.9,67.8)
male_w_p90 <- c(22.3,25.7,29.6,34.3,39.8,45.4,51.6,58.6,63.9,68.4,73.2,73.7,75.8)
male_w_p97 <- c(24.6,28.9,34.0,40.1,47.4,54.4,61.3,68.6,73.4,77.9,83.9,83.4,85.8)
# =========================================================
# 5. 体重データ（女性） 単位: kg
#
# 女性の体重は身長と同様、12〜13歳頃にスパートのピークを迎え、
# その後の増加は男性より緩やかになる。
# =========================================================
female_w_p3 <- c(14.7,16.1,17.9,19.7,21.8,24.2,27.5,31.7,35.6,38.1,39.1,40.2,40.5)
female_w_p10 <- c(15.7,17.3,19.2,21.3,23.7,26.6,30.5,34.9,38.6,41.0,42.1,43.1,43.5)
female_w_p25 <- c(16.9,18.7,20.8,23.3,26.1,29.6,34.1,38.7,42.1,44.4,45.6,46.4,47.0)
female_w_p50 <- c(18.3,20.5,23.0,25.9,29.4,33.7,38.8,43.5,46.6,48.8,50.1,50.9,51.6)
female_w_p75 <- c(20.1,22.7,25.7,29.4,33.6,38.8,44.4,49.1,52.0,54.1,55.7,56.3,57.1)
female_w_p90 <- c(22.0,25.3,28.8,33.4,38.6,44.6,50.5,55.1,57.7,59.8,61.8,62.2,63.3)
female_w_p97 <- c(24.2,28.5,32.7,38.7,45.2,52.0,57.7,62.1,64.4,66.5,69.2,69.5,70.8)
# =========================================================
# 6. 三次スプライン補間
#
# 【なぜスプラインか】
# 元データは1歳刻みの離散値のため、直線で結ぶと折れ線に
# なってしまう。三次スプライン補間を用いることで、各区間を
# 滑らかな三次多項式でつなぎ、生物学的に自然な成長曲線を
# 再現できる。
#
# 【spline() vs smooth.spline()】
# spline(method = "fmm") : 補間（全データ点を厳密に通過）
# smooth.spline() : 平滑化（データ点を近似するが通過しない）
# → 実測パーセンタイル値を正確に保持するため spline() を使用。
#
# 【method = "fmm" とは】
# Forsythe-Malcolm-Moler 法。端点条件（not-a-knot）を使った
# 標準的な三次スプライン。端点でも自然な曲がりを保つ。
# =========================================================
# 補間を行うヘルパー関数
# 引数 y : 13個の実測パーセンタイル値（age に対応）
# 返値 : fine_age（500点）に対応する補間後の値ベクトル
spl <- function(y) spline(age, y, xout = fine_age, method = "fmm")$y
# ---- 男性身長の補間 ----
male_h_s3 <- spl(male_h_p3)
male_h_s10 <- spl(male_h_p10)
male_h_s25 <- spl(male_h_p25)
male_h_s50 <- spl(male_h_p50)
male_h_s75 <- spl(male_h_p75)
male_h_s90 <- spl(male_h_p90)
male_h_s97 <- spl(male_h_p97)
# ---- 女性身長の補間 ----
female_h_s3 <- spl(female_h_p3)
female_h_s10 <- spl(female_h_p10)
female_h_s25 <- spl(female_h_p25)
female_h_s50 <- spl(female_h_p50)
female_h_s75 <- spl(female_h_p75)
female_h_s90 <- spl(female_h_p90)
female_h_s97 <- spl(female_h_p97)
# ---- 男性体重の補間 ----
male_w_s3 <- spl(male_w_p3)
male_w_s10 <- spl(male_w_p10)
male_w_s25 <- spl(male_w_p25)
male_w_s50 <- spl(male_w_p50)
male_w_s75 <- spl(male_w_p75)
male_w_s90 <- spl(male_w_p90)
male_w_s97 <- spl(male_w_p97)
# ---- 女性体重の補間 ----
female_w_s3 <- spl(female_w_p3)
female_w_s10 <- spl(female_w_p10)
female_w_s25 <- spl(female_w_p25)
female_w_s50 <- spl(female_w_p50)
female_w_s75 <- spl(female_w_p75)
female_w_s90 <- spl(female_w_p90)
female_w_s97 <- spl(female_w_p97)
# =========================================================
# 7. 色設定
#
# 【設計方針】
# 成長曲線の「高い側（肥満傾向）＝暖色」「低い側（やせ傾向）＝寒色」
# という直感的な色対応を採用。中央値は最も目立つ黒で強調。
#
# オレンジ系グラデーション（高パーセンタイル側）:
# p97（明るいオレンジ） → p90（濃いオレンジ赤） → p75（琥珀）
# ブルー系グラデーション（低パーセンタイル側）:
# p25（明るい青） → p10（濃い青） → p3（薄い水色）
# =========================================================
col_p97 <- "#EF9F27" # 明るいオレンジ （最高パーセンタイル）
col_p90 <- "#D85A30" # 濃いオレンジ赤
col_p75 <- "#BA7517" # 琥珀・ダークオレンジ
col_p50 <- "#1A1A1A" # ほぼ黒 （中央値：最重要線）
col_p25 <- "#378ADD" # 明るい青
col_p10 <- "#185FA5" # 濃い青
col_p3 <- "#85B7EB" # 薄い水色 （最低パーセンタイル）
# =========================================================
# 8. パネル描画関数
#
# 【引数の説明】
# main_title : パネルタイトル（例："男性 身長"）
# ylab : y軸ラベル（例："身長（cm）"）
# s3〜s97 : 各パーセンタイルの補間済みベクトル（長さ500）
# ylim : y軸の表示範囲 c(下限, 上限)
# y_grid : 水平グリッド線を引くy値のベクトル（NULLなら省略）
# legend_pos : 凡例の位置（"bottomright" など）
#
# 【描画順序】
# ① 背景グリッド（最背面）
# ② 軸・枠
# ③ 各パーセンタイル曲線（上位から順に描画）
# → 後から描いた線が上に重なるため、中央値（p50）を
# 最後に描いて他の線の上に表示
# ④ 凡例
# =========================================================
plot_growth_panel <- function(main_title,
ylab,
s3, s10, s25, s50, s75, s90, s97,
ylim,
y_grid = NULL,
legend_pos = "bottomright") {
# ---- ベースプロット（type = "n" で軸のみ、データは描かない） ----
plot(fine_age, s50,
type = "n", # データ点・線を描かずに枠だけ作る
ylim = ylim,
xlim = c(5, 17),
xlab = "年齢（歳）",
ylab = ylab,
main = main_title,
axes = FALSE) # 軸は後で手動描画（目盛り位置を制御するため）
# ---- 補助グリッド線（データより先に描いて背面に置く） ----
if (!is.null(y_grid)) {
abline(h = y_grid, col = "gray92", lty = 1) # 水平グリッド
}
abline(v = 5:17, col = "gray96", lty = 1) # 垂直グリッド（1歳ごと）
# ---- 軸・枠 ----
axis(1, at = 5:17) # x軸：5〜17歳を1刻みで表示
axis(2, las = 1) # y軸：ラベルを水平表示（las=1）
box() # プロット枠を描画
# ---- パーセンタイル曲線の描画 ----
# 線幅（lwd）の設計:
# 外側の補助線（3%, 97%）: 細め（lwd = 1.5）
# 中間の参照線（10%, 90%, 25%, 75%）: やや太め（lwd = 2.0）
# 中央値（50%）: 最も太く強調（lwd = 3.0）
lines(fine_age, s97, col = col_p97, lwd = 1.5) # 97パーセンタイル
lines(fine_age, s90, col = col_p90, lwd = 2.0) # 90パーセンタイル
lines(fine_age, s75, col = col_p75, lwd = 2.0) # 75パーセンタイル
lines(fine_age, s50, col = col_p50, lwd = 3.0) # 50パーセンタイル（中央値）
lines(fine_age, s25, col = col_p25, lwd = 2.0) # 25パーセンタイル
lines(fine_age, s10, col = col_p10, lwd = 2.0) # 10パーセンタイル
lines(fine_age, s3, col = col_p3, lwd = 1.5) # 3パーセンタイル
# ---- 凡例 ----
# 凡例はプロット上の曲線の並び順（上＝97% → 下＝3%）に合わせる
legend(legend_pos,
legend = c("97%", "90%", "75%", "50%（中央値）", "25%", "10%", "3%"),
col = c(col_p97, col_p90, col_p75, col_p50, col_p25, col_p10, col_p3),
lty = 1, # すべて実線（破線なし）
lwd = c(1.5, 2.0, 2.0, 3.0, 2.0, 2.0, 1.5),
bty = "o", # 凡例ボックスあり
bg = "white", # 凡例の背景色
cex = 0.82) # 凡例の文字サイズ
}
# =========================================================
# 9. 4パネル描画
#
# 【par() パラメータの説明】
# mfrow = c(2, 2) : 2行2列のパネル配置（左上から右下へ順番に埋める）
# las = 1 : 軸ラベルを常に水平表示
# mar = c(4,4.5,3,1): 各パネルの余白（下・左・上・右）単位: 行数
# oma = c(0,0,2,0): 全体の外側余白（上側2行をタイトル用に確保）
#
# 【パネル配置】
# [男性身長] [女性身長]
# [男性体重] [女性体重]
# =========================================================
# 現在の par 設定を保存（関数終了後に元に戻すため）
op <- par(no.readonly = TRUE)
par(mfrow = c(2, 2),
las = 1,
mar = c(4, 4.5, 3, 1),
oma = c(0, 0, 2, 0))
# ---- パネル1: 男性 身長 ----
plot_growth_panel(
main_title = "男性 身長",
ylab = "身長（cm）",
s3 = male_h_s3, s10 = male_h_s10, s25 = male_h_s25,
s50 = male_h_s50, s75 = male_h_s75, s90 = male_h_s90, s97 = male_h_s97,
ylim = c(100, 185),
y_grid = seq(100, 180, 10),
legend_pos = "bottomright"
)
# ---- パネル2: 女性 身長 ----
plot_growth_panel(
main_title = "女性 身長",
ylab = "身長（cm）",
s3 = female_h_s3, s10 = female_h_s10, s25 = female_h_s25,
s50 = female_h_s50, s75 = female_h_s75, s90 = female_h_s90, s97 = female_h_s97,
ylim = c(100, 185),
y_grid = seq(100, 180, 10),
legend_pos = "bottomright"
)
# ---- パネル3: 男性 体重 ----
plot_growth_panel(
main_title = "男性 体重",
ylab = "体重（kg）",
s3 = male_w_s3, s10 = male_w_s10, s25 = male_w_s25,
s50 = male_w_s50, s75 = male_w_s75, s90 = male_w_s90, s97 = male_w_s97,
ylim = c(10, 90),
y_grid = seq(10, 90, 10),
legend_pos = "topleft"
)
# ---- パネル4: 女性 体重 ----
plot_growth_panel(
main_title = "女性 体重",
ylab = "体重（kg）",
s3 = female_w_s3, s10 = female_w_s10, s25 = female_w_s25,
s50 = female_w_s50, s75 = female_w_s75, s90 = female_w_s90, s97 = female_w_s97,
ylim = c(10, 90),
y_grid = seq(10, 90, 10),
legend_pos = "topleft"
)
# ---- 全体タイトル（outer = TRUE で4パネルの外側上部に表示） ----
mtext("日本人男女の身長・体重パーセンタイル成長曲線（5〜17歳）",
outer = TRUE, cex = 1.3)
# par 設定を元に戻す
par(op)
