子どもの身長や体重の成長曲線を推定する手法としてはLMS法が最も有名であり、近年ではその発展形であるGAMLSS（Generalized Additive Model for Location, Scale and Shape）が主流となっています。しかし、成長曲線の推定アプローチはLMS法だけではありません。
本記事では、データの正規性を仮定するLMS法とは異なるアプローチとして、「分位点回帰 (Quantile regression)」と「P-spline」を組み合わせた手法による推定方法を紹介します。あわせて、実際に動作するRスクリプトを提示し、その具体的な使い方も解説します。
分位点回帰ベースのセンタイル推定アルゴリズムの考え方
LMS法との違い
LMS法では、年齢ごとの分布を
- L：歪度（Box-Cox変換）
- M：中央値
- S：ばらつき
という3つのパラメタで表し、それらを年齢の関数として滑らかに推定します。そして、その分布からセンタイル（パーセンタイル）曲線を求めるのでした（前回、前々回の記事参照）。
一方、ここで扱う方法では、
パーセンタイル曲線そのものを直接推定します
この点が最も重要な違いです。
分位点回帰 (Quantile regression)とは何か
通常の回帰分析では平均値を推定しますが、分位点回帰では中央値や上位・下位の分位点を直接推定します。
たとえば、
- 50% → 中央値
- 90% → 上位の境界
- 10% → 下位の境界
といった曲線を、それぞれ独立に推定できます。
ここで与えるスクリプトでは、デフォルトで、以下の7本のセンタイル曲線を同時に推定します：
- 3%, 10%, 25%, 50%, 75%, 90%, 97%
とはいえ、Rスクリプトの以下の部分を修正すれば、センタイルの位置と数を変更できます。
tau_vec <- c(0.03, 0.10, 0.25, 0.50, 0.75, 0.90, 0.97) tau_names <- c("p3", "p10", "p25", "p50", "p75", "p90", "p97")
実際の実装ではこれに加えて、描画関連の部分も修正してください。
曲線の表現：B-splineとP-spline
成長曲線は単純な関数では表せないため、柔軟な表現が必要です。そのために使われるのが B-spline (B-スプライン) です。B-splineは、複数の基底関数の線形結合として曲線を表現する方法で、滑らかな曲線を柔軟に作ることができます。ただし、そのままだと曲線がくねくねしすぎるかもしれないので（不自然に振動するのは避けたいので）、
係数の2階差分にペナルティを課す（P-spline）
ことで、滑らかさを制御します。PはPenalty（ペナルティ）の"P"です。
非交差制約
複数のパーセンタイル曲線を推定するときに重要なのが、
上位の曲線が下位の曲線を下回らないこと
です。たとえば、90%曲線が50%曲線より下に来てしまうと不自然です。
そこで、このスクリプトでは、
となるように制約を課しています。
単調増加制約（オプション）
身長などでは、一定の年齢範囲では「年齢とともに増加する」という性質があります。
そのため、このスクリプトではオプションとして
年齢に対して単調非減少
という制約を加えることができます（私は、身長ではオン TRUE、体重ではオフ FALSE にしています）。
最適化
以上の要素（分位点回帰、滑らかさ、非交差、単調性）をすべてまとめて、
凸最適化問題として解いています
凸最適化問題とは、目的関数と制約条件がすべて凸（下に凸な形）であるような最適化問題のことです。この性質により、局所解が必ず大域解と一致することが保証されており、効率的かつ確実に最適解を求めることができます。成長曲線の推定では、分位点回帰の損失、滑らかさのペナルティ、非交差制約、単調性制約がすべて凸の条件を満たすため、凸最適化として定式化できます。
実装では CVXR パッケージを使い、ECOSやSCSなどのソルバーで解いています。
CVXR はRで凸最適化問題を記述するためのパッケージで、数式に近い形でモデルを定義できます。実際に問題を解く（数値計算を行う）のは、バックエンドのソルバーです。
- ECOS（Embedded Conic Solver）：等式・不等式制約を含む中規模の問題を高速に解くソルバー。精度と速度のバランスが良く、デフォルトとしてよく使われます。
- SCS（Splitting Conic Solver）：大規模な問題に適したソルバー。ECOSより高速ですが、やや精度が低い場合があります。
CVXR が「問題の定式化」を担い、ECOSやSCSが「実際の数値計算」を担う、という役割分担になっています。
Rでやってみる
以上の処理をRで実装したのが以下のRスクリプトです。このスクリプトを実行するためには、身長、体重データが必要ですので、事前に付録に示したRスクリプトを実行して、数値シミュレーションデータを生成しておいてください。
# ============================================================
# 非交差制約付き P-spline 分位点回帰による成長曲線推定
#
# 概要:
# B-spline 基底と差分ペナルティ（P-spline）を用いた
# 分位点回帰により、成長曲線の複数パーセンタイルを同時推定する。
# 最適化には CVXR パッケージを使用し、以下の制約を課す。
# (A) 異なる分位点曲線が交差しない（非交差制約）
# (B) 各分位点曲線が年齢に対して単調非減少（オプション）
# ============================================================
# ============================================================
# 0. 必要パッケージの自動インストールと読み込み
# ============================================================
load_or_install <- function(pkg) {
if (!require(pkg, character.only = TRUE)) {
install.packages(pkg, dependencies = TRUE)
library(pkg, character.only = TRUE)
}
}
pkgs <- c("splines", "CVXR", "Matrix")
invisible(lapply(pkgs, load_or_install))
# ============================================================
# 1. グローバル設定
# ============================================================
# --- 入出力ディレクトリ ---
DIR <- r"(... データがあるフォルダのパスを指定 ...)"
###
input_file <- "simulated_growth_data_from_smoothed_LMS.csv"
true_file <- "true_model_based_centile_curves.csv"
target_sex <- "Male"
target_item <- "Height"
age_grid <- seq(5, 17, by = 0.1)
tau_vec <- c(0.03, 0.10, 0.25, 0.50, 0.75, 0.90, 0.97)
tau_names <- c("p3", "p10", "p25", "p50", "p75", "p90", "p97")
df_pspline <- 15
degree_pspline <- 3
lambda_pspline <- 5
ngrid_constraint <- 50
intercept_pspline <- TRUE
monotone_increasing_pspline <- TRUE
monotone_grid_n <- 25
# OSQP を最優先：制約数が多い疎問題で ECOS より高速
solver_order <- c("OSQP", "ECOS", "SCS")
col_p97 <- "#EF9F27"; col_p90 <- "#D85A30"; col_p75 <- "#BA7517"
col_p50 <- "#1A1A1A"; col_p25 <- "#378ADD"; col_p10 <- "#185FA5"
col_p3 <- "#85B7EB"
centile_cols <- c(col_p3, col_p10, col_p25, col_p50, col_p75, col_p90, col_p97)
centile_lwds <- c(1.5, 2.5, 2.0, 3.5, 2.0, 2.5, 1.5)
# ============================================================
# 2. 補助関数
# ============================================================
# 2階差分行列（P-spline ペナルティ用）
# diff(diag(K), differences=2) で for ループなしに生成
make_diff2_matrix <- function(K) {
if (K < 3) stop("K must be >= 3.")
diff(diag(K), differences = 2)
}
# 前向き1階差分行列（単調性制約用）
# D1[i,i]=-1, D1[i,i+1]=+1 → D1 %*% v >= 0 ⟺ 単調非減少
make_diff1_matrix <- function(M) {
diff(diag(M)) # (M-1) × M
}
make_bs_design <- function(x, df = 12, degree = 3, intercept = TRUE) {
B <- splines::bs(x, df = df, degree = degree, intercept = intercept)
attr_list <- attributes(B)
list(
B = as.matrix(B),
knots = attr_list$knots,
Boundary.knots = attr_list$Boundary.knots,
degree = attr_list$degree,
intercept = intercept
)
}
predict_bs_design <- function(newx, basis_info) {
as.matrix(splines::bs(
newx,
knots = basis_info$knots,
Boundary.knots = basis_info$Boundary.knots,
degree = basis_info$degree,
intercept = basis_info$intercept
))
}
check_monotone_non_decreasing <- function(x, y, tol = 1e-8) {
all(diff(y[order(x)]) >= -tol)
}
check_all_quantile_curves_monotone <- function(age, qmat, tol = 1e-8) {
apply(qmat, 2, check_monotone_non_decreasing, x = age, tol = tol)
}
# ============================================================
# 3. 非交差制約付き P-spline 分位点回帰（高速化・バグ修正版）
# ============================================================
fit_noncrossing_pspline_qr <- function(
x,
y,
tau_vec = c(0.1, 0.5, 0.9),
df = 12,
degree = 3,
lambda = 5,
ngrid_constraint = 120,
intercept = TRUE,
monotone_increasing = FALSE,
monotone_grid_n = 200,
solver_order = c("OSQP", "ECOS", "SCS")
) {
if (length(x) != length(y)) stop("x と y の長さが一致していません。")
sel <- is.finite(x) & is.finite(y)
x <- as.numeric(x[sel]); y <- as.numeric(y[sel])
if (length(x) < 20) stop("有限な観測点が少なすぎます（20点以上必要）。")
tau_vec <- sort(unique(as.numeric(tau_vec)))
if (any(tau_vec <= 0 | tau_vec >= 1))
stop("tau_vec の各要素は (0, 1) の範囲である必要があります。")
ord <- order(x); x <- x[ord]; y <- y[ord]
ntau <- length(tau_vec)
# --- 計画行列・差分行列（一度だけ計算）---
basis_info <- make_bs_design(x = x, df = df, degree = degree, intercept = intercept)
B <- basis_info$B
K <- ncol(B)
D2 <- make_diff2_matrix(K) # (K-2) × K
# 制約グリッド上の計画行列
xg_cross <- seq(min(x), max(x), length.out = ngrid_constraint)
Bg_cross <- predict_bs_design(xg_cross, basis_info) # ngrid × K
xg_mono <- seq(min(x), max(x), length.out = monotone_grid_n)
Bg_mono <- predict_bs_design(xg_mono, basis_info) # M × K
# --- CVXR 決定変数 ---
beta_list <- lapply(seq_len(ntau), function(j) Variable(K))
# --- 目的関数 ---
obj_terms <- lapply(seq_len(ntau), function(j) {
tau <- tau_vec[j]
resid_j <- y - B %*% beta_list[[j]]
tau * sum_entries(pos(resid_j)) +
(1 - tau) * sum_entries(pos(-resid_j)) +
lambda * sum_squares(D2 %*% beta_list[[j]])
})
objective <- Minimize(Reduce(`+`, obj_terms))
# --- 制約①：非交差（lapply で一括生成）---
constraints_cross <- if (ntau >= 2) {
lapply(seq_len(ntau - 1), function(j)
Bg_cross %*% beta_list[[j]] <= Bg_cross %*% beta_list[[j + 1]]
)
} else {
list()
}
# --- 制約②：単調非減少（差分行列ベクトル化・バグ修正済み）---
# make_diff1_matrix() の diff(diag(M)) で符号を正しく保証。
# D1 %*% Bg_mono をループ外で事前計算し、CVXR にはベクトル不等式
# 1本として渡す（旧: ntau × (M-1) 本のスカラー制約）。
constraints_mono <- if (isTRUE(monotone_increasing)) {
D1 <- make_diff1_matrix(monotone_grid_n) # (M-1) × M
DMB <- D1 %*% Bg_mono # (M-1) × K（定数・事前計算）
lapply(seq_len(ntau), function(j) DMB %*% beta_list[[j]] >= 0)
} else {
list()
}
constraints <- c(constraints_cross, constraints_mono)
# --- 最適化 ---
problem <- Problem(objective, constraints)
result <- NULL
last_error <- NULL
for (sv in solver_order) {
solve_args <- if (sv == "OSQP") {
list(solver = sv, warm_start = TRUE,
eps_abs = 1e-3, eps_rel = 1e-3, max_iter = 5000L)
} else {
list(solver = sv)
}
tmp <- tryCatch(
do.call(solve, c(list(problem), solve_args)),
error = function(e) e
)
if (!inherits(tmp, "error")) { result <- tmp; break }
last_error <- tmp$message
}
if (is.null(result))
stop(paste("すべてのソルバーで最適化に失敗しました。最後のエラー:", last_error))
status_now <- tryCatch(result$status, error = function(e) NA_character_)
if (is.na(status_now) || !(status_now %in% c("optimal", "optimal_inaccurate")))
stop(paste("最適化が正常終了していません。status =", status_now))
# --- 結果収集 ---
coef_list <- lapply(beta_list, function(b) result$getValue(b))
if (any(vapply(coef_list, is.null, logical(1)))) stop("係数の取得に失敗しました。")
coef_mat <- as.matrix(do.call(cbind, coef_list))
colnames(coef_mat) <- paste0("tau_", tau_vec)
# 行列×行列積（BLAS dgemm）で全分位点を一括計算
fitted_mat <- B %*% coef_mat
pred_cross_grid_mat <- Bg_cross %*% coef_mat
pred_mono_grid_mat <- Bg_mono %*% coef_mat
colnames(fitted_mat) <- paste0("tau_", tau_vec)
colnames(pred_cross_grid_mat) <- paste0("tau_", tau_vec)
colnames(pred_mono_grid_mat) <- paste0("tau_", tau_vec)
structure(
list(
x = x,
y = y,
tau_vec = tau_vec,
coef = coef_mat,
fitted = fitted_mat,
xg_cross = xg_cross,
pred_cross_grid = pred_cross_grid_mat,
xg_mono = xg_mono,
pred_mono_grid = pred_mono_grid_mat,
lambda = lambda,
df = df,
degree = degree,
basis_info = basis_info,
monotone_increasing = monotone_increasing,
monotone_grid_n = monotone_grid_n,
result = result
),
class = "noncrossing_pspline_qr"
)
}
# ============================================================
# 4. 予測関数
# ============================================================
predict_noncrossing_pspline_qr <- function(fit, newx) {
Bnew <- predict_bs_design(as.numeric(newx), fit$basis_info)
pred <- Bnew %*% fit$coef
colnames(pred) <- paste0("tau_", fit$tau_vec)
pred
}
# ============================================================
# 5. データの読み込みと対象群の抽出
# ============================================================
setwd(DIR)
dat_all <- read.csv(input_file, stringsAsFactors = FALSE)
req_cols <- c("sex", "item", "age", "value")
if (!all(req_cols %in% names(dat_all)))
stop(sprintf("入力データには %s 列が必要です。", paste(req_cols, collapse = ", ")))
dat_all$sex <- as.character(dat_all$sex)
dat_all$item <- as.character(dat_all$item)
dat_all$age <- as.numeric(dat_all$age)
dat_all$value <- as.numeric(dat_all$value)
dat <- subset(dat_all, sex == target_sex & item == target_item)
if (nrow(dat) == 0)
stop(sprintf("指定した条件のデータが見つかりません: sex=%s, item=%s",
target_sex, target_item))
dat <- dat[is.finite(dat$age) & is.finite(dat$value), ]
dat <- dat[order(dat$age), ]
if (nrow(dat) < 20) stop("対象群のデータ数が少なすぎます（20点以上必要）。")
# 読み込み時点のデータ点数を記録
n_total_loaded <- nrow(dat_all)
n_used <- nrow(dat)
# ============================================================
# 6. 真の曲線の読み込みと補間
# ============================================================
true_all <- read.csv(true_file, stringsAsFactors = FALSE)
req_true_cols <- c("sex", "item", "age", tau_names)
if (!all(req_true_cols %in% names(true_all)))
stop(sprintf("真の曲線データには %s 列が必要です。",
paste(req_true_cols, collapse = ", ")))
true_all$sex <- as.character(true_all$sex)
true_all$item <- as.character(true_all$item)
true_all$age <- as.numeric(true_all$age)
true_curve <- subset(true_all, sex == target_sex & item == target_item)
if (nrow(true_curve) == 0)
stop(sprintf("指定した条件の真の曲線が見つかりません: sex=%s, item=%s",
target_sex, target_item))
true_curve <- true_curve[is.finite(true_curve$age), ]
true_curve <- true_curve[order(true_curve$age), ]
true_interp <- data.frame(sex = target_sex, item = target_item, age = age_grid)
for (nm in tau_names) {
true_interp[[nm]] <- approx(
x = true_curve$age, y = true_curve[[nm]],
xout = age_grid, rule = 2
)$y
}
# ============================================================
# 7. モデルのフィッティング（計算時間を計測）
# ============================================================
fit_time <- system.time({
fit <- fit_noncrossing_pspline_qr(
x = dat$age,
y = dat$value,
tau_vec = tau_vec,
df = df_pspline,
degree = degree_pspline,
lambda = lambda_pspline,
ngrid_constraint = ngrid_constraint,
intercept = intercept_pspline,
monotone_increasing = monotone_increasing_pspline,
monotone_grid_n = monotone_grid_n,
solver_order = solver_order
)
})
pred_age_grid <- predict_noncrossing_pspline_qr(fit, age_grid)
colnames(pred_age_grid) <- tau_names
curve_est <- data.frame(
sex = target_sex, item = target_item, age = age_grid,
pred_age_grid, row.names = NULL
)
# 非交差チェック
is_noncrossing <- all(sapply(seq_len(ncol(pred_age_grid) - 1), function(j)
all(pred_age_grid[, j] <= pred_age_grid[, j + 1] + 1e-8)
))
# 単調性チェック
mono_check <- check_all_quantile_curves_monotone(
age = curve_est$age,
qmat = as.matrix(curve_est[, tau_names]),
tol = 1e-8
)
# ============================================================
# 8. 推定値と真値の比較
# ============================================================
comparison_df <- data.frame(sex = target_sex, item = target_item, age = age_grid)
for (nm in tau_names) {
comparison_df[[paste0(nm, "_est")]] <- curve_est[[nm]]
comparison_df[[paste0(nm, "_true")]] <- true_interp[[nm]]
comparison_df[[paste0(nm, "_diff")]] <- curve_est[[nm]] - true_interp[[nm]]
comparison_df[[paste0(nm, "_absdiff")]] <- abs(curve_est[[nm]] - true_interp[[nm]])
}
rmse_vec <- sapply(tau_names, function(nm)
sqrt(mean((curve_est[[nm]] - true_interp[[nm]])^2, na.rm = TRUE)))
mae_vec <- sapply(tau_names, function(nm)
mean(abs(curve_est[[nm]] - true_interp[[nm]]), na.rm = TRUE))
summary_error_df <- data.frame(
quantile = tau_names,
tau = tau_vec,
RMSE = as.numeric(rmse_vec),
MAE = as.numeric(mae_vec)
)
# ============================================================
# 9. 描画
# ============================================================
op <- par(no.readonly = TRUE)
par(mfrow = c(1, 1), mar = c(4, 4.5, 3, 1), las = 1, oma = c(0, 0, 1, 0))
ylab_text <- if (target_item == "Height") "Height (cm)" else "Weight (kg)"
ylim_all <- range(dat$value, curve_est[, tau_names],
true_interp[, tau_names], na.rm = TRUE)
plot(dat$age, dat$value,
pch = 16, cex = 0.6, col = gray(0.7, 0.30),
xlim = range(age_grid), ylim = ylim_all,
xlab = "Age (years)", ylab = ylab_text,
main = paste0("Estimated vs true growth curves: ",
target_sex, " ", target_item))
qmat_est <- as.matrix(curve_est[, tau_names])
qmat_true <- as.matrix(true_interp[, tau_names])
for (j in seq_len(ncol(qmat_est))) {
lines(curve_est$age, qmat_est[, j], col = centile_cols[j],
lwd = centile_lwds[j], lty = 1)
lines(true_interp$age, qmat_true[, j], col = centile_cols[j],
lwd = 1.5, lty = 2)
}
legend("bottomright", legend = c("Estimated", "True"),
lty = c(1, 2), lwd = c(2, 1.5), col = "black",
bty = "o", bg = "white", cex = 0.85)
legend("topleft",
legend = rev(c("3%", "10%", "25%", "50%", "75%", "90%", "97%")),
col = rev(centile_cols), lty = 1, lwd = rev(centile_lwds),
bty = "o", bg = "white", cex = 0.8)
par(op)
# ============================================================
# 10. 結果の保存
# ============================================================
out_csv_est <- sprintf("estimated_growth_curves_pspline_%s_%s.csv",
target_sex, target_item)
out_csv_cmp <- sprintf("comparison_estimated_vs_true_%s_%s.csv",
target_sex, target_item)
out_csv_err <- sprintf("summary_error_estimated_vs_true_%s_%s.csv",
target_sex, target_item)
write.csv(curve_est, file = out_csv_est, row.names = FALSE)
write.csv(comparison_df, file = out_csv_cmp, row.names = FALSE)
write.csv(summary_error_df, file = out_csv_err, row.names = FALSE)
# ============================================================
# 11. 実行サマリーの出力
# ============================================================
# 精度指標テーブルの行を組み立て
error_rows <- paste0(
sprintf(" %-6s %5.2f %8.4f %8.4f",
summary_error_df$quantile,
summary_error_df$tau,
summary_error_df$RMSE,
summary_error_df$MAE),
collapse = "\n"
)
# 単調性グリッド行（制約が有効な場合のみ）
mono_grid_line <- if (monotone_increasing_pspline) {
sprintf(" 単調性グリッド : %d 点\n", monotone_grid_n)
} else {
""
}
# 文字列を一括組み立て
summary_text <- paste0(
"\n",
strrep("=", 60), "\n",
" 実行サマリー\n",
strrep("=", 60), "\n",
"\n【対象データ】\n",
sprintf(" 入力ファイル : %s\n", input_file),
sprintf(" 対象群 : sex = %s / item = %s\n", target_sex, target_item),
sprintf(" 全データ点数 : %d 点\n", n_total_loaded),
sprintf(" 推定使用点数 : %d 点\n", n_used),
sprintf(" 除外点数 : %d 点（有限値フィルタ）\n", n_total_loaded - n_used),
sprintf(" 年齢範囲（使用） : %.2f 〜 %.2f 歳\n", min(dat$age), max(dat$age)),
"\n【モデル設定】\n",
sprintf(" 推定パーセンタイル : %s\n",
paste(sprintf("%.0f%%", tau_vec * 100), collapse = ", ")),
sprintf(" B-spline 自由度 : %d\n", df_pspline),
sprintf(" B-spline 次数 : %d\n", degree_pspline),
sprintf(" 切片項 : %s\n", intercept_pspline),
sprintf(" ペナルティ強度 λ : %.4g\n", lambda_pspline),
sprintf(" 非交差制約グリッド : %d 点\n", ngrid_constraint),
sprintf(" 単調増加制約 : %s\n",
ifelse(monotone_increasing_pspline, "有効", "無効")),
mono_grid_line,
sprintf(" ソルバー優先順位 : %s\n",
paste(solver_order, collapse = " > ")),
sprintf(" 使用ソルバー : %s\n",
if (!is.null(fit$result$solver)) fit$result$solver else solver_order[1]),
"\n【最適化結果】\n",
sprintf(" ソルバーステータス : %s\n", fit$result$status),
sprintf(" 非交差（age_grid） : %s\n",
ifelse(is_noncrossing, "OK（交差なし）", "NG（交差あり）")),
sprintf(" 単調非減少 : %s\n",
ifelse(all(mono_check), "全曲線 OK",
paste0("NG（", paste(tau_names[!mono_check], collapse = ", "),
" が非単調）"))),
"\n【精度指標（推定 vs 真値）】\n",
sprintf(" %-6s %5s %8s %8s\n", "分位点", "tau", "RMSE", "MAE"),
sprintf(" %s\n", strrep("-", 36)),
error_rows, "\n",
sprintf(" %s\n", strrep("-", 36)),
sprintf(" %-6s %5s %8.4f %8.4f\n", "平均", "",
mean(summary_error_df$RMSE), mean(summary_error_df$MAE)),
"\n【計算時間】\n",
sprintf(" user : %7.2f 秒\n", fit_time["user.self"]),
sprintf(" system : %7.2f 秒\n", fit_time["sys.self"]),
sprintf(" elapsed : %7.2f 秒\n", fit_time["elapsed"]),
"\n【保存ファイル】\n",
paste0(" ", c(out_csv_est, out_csv_cmp, out_csv_err), collapse = "\n"), "\n",
strrep("=", 60), "\n"
)
cat(summary_text)
このスクリプトは、以下の処理を含んでいます：
- B-splineで曲線を表現する
- 分位点回帰で各パーセンタイルを推定する
- 2階差分ペナルティで滑らかさを制御する
- 非交差制約を入れる
- 必要に応じて単調増加制約を入れる
- すべてを同時に最適化する
Rスクリプトの使い方
ここからは、実際に添付のRスクリプトをどのように使うかを説明します。Rスクリプトの実行には、かなりの計算負荷がかかりますので、PCのスペックが貧弱だと、待ちきれないほどの計算時間がかかるかもしれません。ちなみに私のPC（Intel(R) Core(TM) Ultra 9 285, RAM 160GB）では、1万点に対するフィットで、計算時間が 7秒、20万点に対するフィットで 7分20秒でした。100万点でも、85分で処理できました。
0. データの準備とパスの指定
Rスクリプトを実行するには、まず分析対象の身長・体重データ（数値シミュレーションデータ）と、推定結果を真値と比較するためのデータを用意する必要があります。これらのデータは、付録のRスクリプトを実行することで生成されます。
付録のRスクリプトを実行する前に、以下の箇所にデータを保存しているフォルダのパスを指定してください。
# 出力先ディレクトリ DIR <- r"(... ここにファイルパスを指定 ...)"
この DIR は、メインのRスクリプトでも同じパスを指定してください。両スクリプトで DIR を一致させることで、ファイルの読み書きが正しく動作します。
入力データ
付録のRスクリプトを実行すると、以下の2つのCSVファイルが DIR に出力されます。メインのRスクリプトはこれらのファイルを読み込んで動作します。
① 身長・体重データ simulated_growth_data_from_smoothed_LMS.csv
| 列名 | 内容 |
|---|---|
sex |
性別（例："Male", "Female"） |
item |
測定項目（例："Height", "Weight"） |
age |
年齢（数値） |
value |
測定値（身長または体重） |
② 真の成長曲線（比較用） true_model_based_centile_curves.csv
| 列名 | 内容 |
|---|---|
sex, item, age |
上と同様 |
p3, p10, p25, p50, p75, p90, p97 |
各パーセンタイルの真値 |
1. 必要なパッケージ
このスクリプトでは以下の2つのパッケージを使用します。
splines：B-spline基底関数の生成に使用します。CVXR：凸最適化問題の定式化と求解に使用します。
スクリプト内に自動インストールの処理が含まれているため、未インストールの場合でも、そのまま実行すれば自動的にインストールされます。
2. 最初に設定する項目
スクリプト冒頭で以下の項目を設定します。
作業ディレクトリ
# --- 入出力ディレクトリ --- DIR <- r"(... データがあるフォルダのパスを指定 ...)"
「0. データの準備とパスの指定」で指定したパスと同じ値を設定してください。ここに置かれたCSVファイルが読み込まれ、出力結果もここに保存されます。
対象データの指定
target_sex <- "Male" # "Male" または "Female" target_item <- "Height" # "Height" または "Weight"
分析したい性別と測定項目を指定します。
年齢グリッド
age_grid <- seq(5, 17, by = 0.1)
成長曲線を出力する年齢の範囲と刻み幅を指定します。上の例では5歳から17歳まで0.1歳刻みで曲線を推定します。
3. パラメータの設定
推定の挙動を制御する主要なパラメータは以下の3つです。
| パラメータ | 役割 | 目安 |
|---|---|---|
df_pspline |
B-splineの自由度。大きいほど曲線が柔軟になります | まずはデフォルト値で試してください |
lambda_pspline |
滑らかさのペナルティ強度。大きいほど滑らかになります | 大きすぎると過度に平滑化されます |
monotone_increasing_pspline |
単調増加制約の有無（TRUE / FALSE） |
身長など単調に増加する指標には TRUE を推奨します |
初回はデフォルト値のまま実行し、グラフや評価指標を確認しながら調整することをおすすめします。
4. 実行の流れ
スクリプトを実行すると、以下の順に処理が進みます。
- データの読み込み：
DIRからCSVファイルを読み込みます。 - 対象群の抽出：
target_sexとtarget_itemに該当するデータを抽出します。 - 真値の補間：比較用データを
age_gridに合わせて補間します。 - モデルのフィット：凸最適化によって各パーセンタイル曲線を推定します。
- 結果の評価・保存：推定曲線・真値との比較・評価指標をCSVに保存し、グラフを表示します。
5. 出力結果
CSVファイル（DIR に保存）
- 推定されたパーセンタイル曲線
- 真値との比較データ
- 評価指標（RMSE・MAE）
これらにより、推定曲線が真の成長曲線をどの程度再現できているかを定量的に確認できます。
グラフ
グラフには以下の3要素が重ねて表示されます。
- 散布図：元の観測データ
- 実線：推定されたパーセンタイル曲線
- 破線：真のパーセンタイル曲線
グラフを通じて、曲線の滑らかさ・非交差性・フィットの良さを視覚的に確認してください。
おわりに
本記事では、LMS法とは異なる成長曲線の推定法として、分位点回帰 ＋ P-spline ＋ 制約付き最適化 というアプローチを紹介しました。この方法には、以下の特徴があります。
- パーセンタイルの直接推定：LMS法のように正規分布を仮定してパーセンタイルを変換するのではなく、各パーセンタイル曲線を直接データから推定します。
- 非交差の保証：異なるパーセンタイル曲線が交差しないよう、最適化の段階で制約として組み込んでいます。
- 単調性の制御：身長のように年齢とともに単調増加するべき指標に対して、その性質を制約として課すことができます。
LMS法が「分布モデルベース」の推定法であるのに対し、本手法は「曲線直接推定型」の推定法です。一方が他方より常に優れているとは言えず、データの特性や分析の目的に応じて使い分けることが重要です。
今回と前回の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 <- 777
# ---------------------------------------------------------
# 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)
