Rで解析：反復測定や多層データの曲線的な変化をスプラインで当てはめる「MultiSpline」パッケージ

「MultiSpline」パッケージは、反復測定や、学校と生徒のような入れ子の構造を持つデータで、結果が説明変数に対して曲線を描いて変わる様子を、スプラインで当てはめて調べられるパッケージです。渡すのは、1行が1回の測定にあたる縦長のデータフレームで、結果の列、曲線で表したい数値の列、個体や集団を表す列の名前を文字列で指定します。モデルの当てはめには内部でlme4パッケージとmgcvパッケージを使い、結果の図はggplot2パッケージで返ります。ノット数の選択から、予測値、傾きの変化、集団ごとの曲線のばらつきまでを同じモデルから順に取り出せます。成長や経過の曲線を、個体差を考えたうえで読みたいときに役立てられるのではないかと考えます。

パッケージバージョンは0.3.0。Windows 11 x64 (build 26200)のR version 4.6.1で確認しています。

＜おすすめのRに関する書籍です＞

Rユーザのためのtidymodels[実践]入門〜モダンな統計・機械学習モデリングの世界

パッケージごとにバラバラだった機械学習の書き方をtidymodelsで統一する方法を解説。前処理から評価・運用までのパイプラインを組めるようになる。

パッケージのインストール

下記コマンドを実行してください。lmerTestパッケージは係数のp値の計算に、patchworkパッケージは複数の図を並べるときに使います。

# パッケージのインストール
install.packages(c("MultiSpline", "lmerTest", "patchwork"))

# パッケージの読み込み
library("MultiSpline")

コマンド例

詳細はコメント、パッケージのヘルプを確認してください。

スプライン：区間ごとの多項式をなめらかにつないで曲線を表す方法。ノット（つなぎ目）の数は自由度（df）で決まり、多いほど細かく曲がる。クラスター：同じ個体や集団からの測定のまとまり。混合モデルではクラスターごとの違いを変量効果として扱う。

例として、宇治、宇治田原、和束、南山城の茶園24か所で、施肥後0〜12週の新芽の長さを毎週測った架空のデータを作ります。新芽は5週目前後で急に伸び、その後は頭打ちになるとし、伸び始めの長さと伸びの大きさは茶園ごとに変えています。

# 宇治茶の茶園24か所で、施肥後0〜12週の新芽の長さ（cm）を測った架空データを作成
set.seed(1205)
Garden <- paste0(
  rep(c("Uji", "Ujitawara", "Wazuka", "Minamiyamashiro"), each = 6),
  "_", rep(1:6, times = 4)
)
UjiTea <- expand.grid(Week = 0:12, Garden = Garden, stringsAsFactors = FALSE)
UjiTea$Area <- sub("_.*", "", UjiTea$Garden)

# 茶園ごとの伸び始めの長さと伸びの大きさ
Base <- setNames(rnorm(24, mean = 0, sd = 1.2), Garden)
Gain <- setNames(rnorm(24, mean = 1, sd = 0.15), Garden)

# 5週目前後で急に伸び、その後は頭打ちになるS字の伸び方
UjiTea$Shoot <- round(2 + Base[UjiTea$Garden] +
  Gain[UjiTea$Garden] * 9 / (1 + exp(-(UjiTea$Week - 5) / 1.5)) +
  rnorm(nrow(UjiTea), mean = 0, sd = 0.6), 1)

# データの確認
head(UjiTea)
  Week Garden Area Shoot
1    0  Uji_1  Uji   3.0
2    1  Uji_1  Uji   2.4
3    2  Uji_1  Uji   3.1
4    3  Uji_1  Uji   4.2
5    4  Uji_1  Uji   5.7
6    5  Uji_1  Uji   7.0

スプラインの自由度を情報量規準で選ぶ：nl_knotsコマンド

df_rangeオプションに並べた自由度で順にモデルを当てはめ、AICまたはBICが最も小さい自由度を示します。plot = TRUE では、自由度ごとの情報量規準の図も出ます。比較には最尤法を使います。

# 茶園をクラスターとして、自由度2〜6からAICで選ぶ
nl_knots(UjiTea, y = "Shoot", x = "Week", cluster = "Garden", df_range = 2:6)
Best df = 5  (AIC = 950.83)

茶園をクラスターとして、自由度2〜6からAICで選ぶ

スプラインを組み込んだモデルを当てはめる：nl_fitコマンド

xオプションの変数をスプラインに展開してモデルを当てはめます。clusterオプションを指定すると、lme4パッケージの混合モデル（正規分布ならlmer、ほかの分布族ならglmer）になります。randomオプションに lme4 の書き方で変量効果を書くと、cluster、nested、random_slope より優先されます。method = “gam” は単一レベルのモデルだけで使えます。nl_knotsでは自由度5が選ばれましたが、自由度4とのAICの差は1.3ほどなので、ここでは自由度4のモデルを当てはめ、以降の節で使い回します。

# 茶園ごとの変量切片を入れた、自由度4の自然3次スプラインのモデル
fit <- nl_fit(UjiTea, y = "Shoot", x = "Week", cluster = "Garden", df = 4)

係数の表を出す：nl_summaryコマンド

固定効果の係数、標準誤差、自由度、検定統計量、p値を表で返します。スプラインの係数は1つずつでは解釈しにくいので、曲線の形は後の節の予測値と傾きで読みます。

# 係数の表を表示
nl_summary(fit)
# A tibble: 5 × 6
  Term                       Estimate Std.Error    df statistic p.value
  <chr>                         <dbl>     <dbl> <dbl>     <dbl>   <dbl>
1 (Intercept)                    1.86     0.41   31.7      4.54 0.00008
2 splines::ns(Week, df = 4)1     5.96     0.257 284       23.2  0
3 splines::ns(Week, df = 4)2     7.83     0.241 284       32.6  0
4 splines::ns(Week, df = 4)3     8.5      0.443 284       19.2  0
5 splines::ns(Week, df = 4)4     8.09     0.211 284       38.3  0

スプラインの項をまとめて検定する：nl_waldコマンド

スプラインに展開した項のすべての係数が0かどうかを、Wald検定でまとめて調べます。説明変数が結果に関係しているかを1つのp値で確かめられます。

# 施肥後の週の項をまとめて検定
nl_wald(fit)
  term    chisq df p.value
1 Week 3547.585  4       0

＜おすすめのRに関する書籍です＞

Rクックブック 第2版 | J.D.Long, Paul Teetor, 大橋真也

「これをやりたい」から引ける約300本のレシピ集。手元に置いて困ったときに開く、リファレンス的な使い方が合う一冊。

級内相関係数を求める：nl_iccコマンド

変量効果の分散が、全体の分散（固定効果の分を除く）のうちどれだけを占めるかを、クラスターごとに返します。この例では、新芽の長さのばらつきのおよそ8割が茶園の違いによるものです。

# 茶園と残差の割合
nl_icc(fit)
   Garden  Residual
0.7901774 0.2098226

決定係数と分散の内訳を出す：nl_r2コマンド

固定効果だけで説明できる割合（周辺R2）と、変量効果まで含めた割合（条件付きR2）を出し、全体の分散を固定効果、クラスター、残差に分けて示します。

# 決定係数と分散の内訳
nl_r2(fit)

=== MultiSpline R-squared Decomposition ===
Model type: LMM

  Marginal  R2m = 0.7053   (fixed effects only)
  Conditional R2c = 0.9382  (fixed + all random effects)

Variance Partition (r2_mlm-style):
     Component Variance Proportion
 Fixed effects  10.1605     0.7053
        Garden   3.3544     0.2329
      Residual   0.8907     0.0618

クラスターごとの曲線のばらつきを調べる：nl_hetコマンド

変量切片だけのモデルを、スプラインの変量傾きを入れたモデルに当てはめ直し、尤度比検定で比べます。曲線の形そのものがクラスターごとに違うかを調べられます。plot = TRUE では、クラスターごとの曲線と全体の平均の曲線を重ねた図が出ます。分散が0の境界で検定するため、p値は保守的です。

# 茶園ごとの伸び方の違いを検定し、12園の曲線をプロット
het <- nl_het(fit, n_clusters_plot = 12)
Refitting with random spline slopes to assess heterogeneity ...
boundary (singular) fit: see help('isSingular')

--- LRT (ML): random spline slopes vs random intercepts ---
            model npar      AIC    logLik    Chisq Df     p.value
 random intercept    7 952.1294 -469.0647       NA NA          NA
    random slopes   21 768.9839 -363.4920 211.1455 14 2.88049e-37
p-value is conservative: variances are tested on the boundary.

茶園ごとの伸び方の違いを検定し、12園の曲線をプロット

直線や多項式のモデルと比べる：nl_compareコマンド

同じデータで直線と多項式のモデルを当てはめ、スプラインのモデルとAIC、BIC、対数尤度で比べます。直線のモデルとの尤度比検定も付きます。曲線で表す必要があるかを確かめるときに使います。

# 直線、2次、3次の多項式と比較
nl_compare(fit, polynomial_degrees = 2:3)

=== MultiSpline Model Comparison ===

   Model      AIC      BIC   LogLik npar Deviance LRT_vs_linear  LRT_p
  Linear 1085.408 1100.380 -538.704    4 1077.408            NA   <NA>
 Poly(2) 1055.292 1074.007 -522.646    5 1045.292        32.117 <0.001
 Poly(3)  974.230  996.688 -481.115    6  962.230       115.178 <0.001
  Spline  952.129  978.330 -469.065    7  938.129       139.279 <0.001

  Best model by AIC: Spline

＜おすすめのRに関する書籍です＞

StanとRでベイズ統計モデリング Wonderful R 2 | 松浦健太郎, 石田基広

Stanによるベイズ統計モデリングの実践書。現実のデータ解析を念頭に、モデリングの考え方から丁寧に解説する。

当てはめた曲線の予測値を求める：nl_predictコマンド

xの値ごとに、当てはめた曲線の予測値と標準誤差、信頼区間を返します。type = “conditional” は変量効果を0とした典型的なクラスターの曲線、”marginal” は母集団で平均した曲線です。x_seqを指定しないときは、xの範囲を細かく区切った値で求めます。

# 2週ごとの予測値と95%信頼区間
pred <- nl_predict(fit, x_seq = seq(0, 12, by = 2))
pred
# A tibble: 7 × 5
   Week   fit se.fit   lwr   upr
  <dbl> <dbl>  <dbl> <dbl> <dbl>
1     0  1.86  0.410  1.06  2.66
2     2  2.37  0.389  1.61  3.13
3     4  4.33  0.387  3.58  5.09
4     6  7.22  0.393  6.45  7.99
5     8  9.04  0.387  8.28  9.80
6    10  9.84  0.389  9.08 10.6
7    12 10.2   0.410  9.36 11.0

曲線の傾きと曲がり方を求める：nl_derivativesコマンド

当てはめた曲線の一次微分（傾き、d1）と二次微分（曲がり方、d2）を、信頼区間付きで返します。この例では、新芽が最も速く伸びるのは5週目で、1週あたり約1.5cmです。0週目と11週目以降は、傾きの信頼区間が0をまたいでいます。

# 1週ごとの傾きと信頼区間
deriv <- nl_derivatives(fit, x_seq = 0:12)
deriv[, c("Week", "d1", "d1_lwr", "d1_upr")]
# A tibble: 13 × 4
    Week    d1  d1_lwr d1_upr
   <dbl> <dbl>   <dbl>  <dbl>
 1     0 0.126 -0.0890  0.341
 2     1 0.223  0.0305  0.415
 3     2 0.512  0.384   0.639
 4     3 0.994  0.910   1.08
 5     4 1.41   1.26    1.55
 6     5 1.50   1.36    1.64
 7     6 1.26   1.18    1.35
 8     7 0.896  0.757   1.04
 9     8 0.606  0.461   0.751
10     9 0.390  0.306   0.474
11    10 0.240  0.112   0.367
12    11 0.150 -0.0422  0.342
13    12 0.120 -0.0952  0.335

極大・極小と変曲点の位置を探す：nl_turning_pointsコマンド

nl_derivativesの結果から、傾きが0になる点（極大・極小）、曲がる向きが変わる区間（変曲点）、増加・減少が続く区間を探します。この例は単調に伸び続けるS字なので極大・極小は無く、4〜5週目に下に凸から上に凸へ変わる区間が見つかります。

# 極大・極小、変曲点、増減の区間
nl_turning_points(deriv)
$turning_points
[1] Week fit  type
<0 rows> (or 0-length row.names)

$inflection_regions
  x_start x_end         direction
1       4     5 convex_to_concave

$slope_regions
  x_start x_end  direction
1       0    12 increasing

予測値の曲線と傾きをプロット：nl_plotコマンド

nl_predictとnl_derivativesの結果を、ggplot2パッケージのプロットにします。type = “trajectory” は予測値の曲線、”slope” は傾き、”curvature” は曲がり方、”combo” は予測値と傾きを横に並べた図です。

# 細かく区切った週で予測値と傾きを求める
pred_all <- nl_predict(fit)
deriv_all <- nl_derivatives(fit)

# 予測値の曲線と傾きを並べてプロット
nl_plot(pred_all, deriv_all,
  type = "combo",
  x_lab = "Weeks after fertilization", y_lab = "Shoot length (cm)"
)

# 曲がり方（二次微分）をプロット
nl_plot(deriv_df = deriv_all, type = "curvature")

予測値の曲線と傾きを並べてプロット

曲がり方（二次微分）をプロット

＜おすすめのRに関する書籍です＞

データに触れながら学ぶ統計学 R言語のUIツール「Exploratory」で身につくデータサイエンスの手法 | 西田勘一郎

RベースのUIツール「Exploratory」を使い、コードを書かずに統計とデータサイエンスの手法を体験できる。

この記事が誰かの役に立ちますように。