RでPrincipal axis regressionを実装してみる

r
Principal axis regressionの実装方法について説明します
Published

2026-05-26

Modified

2026-05-26

本記事は、Zhou et al. (2025) の Fig. 4 で用いられている考え方に着想を得て、principal axis regression による trait space の主軸推定と回転角の計算を学ぶためのチュートリアルです。ここで使用するデータはすべて教育目的で作成した架空データであり、Zhou et al. (2025) の元データ、図から読み取った値、または論文の結果を再現するためのデータではありません。解析手順や生態学的解釈の詳細については、原著論文を参照してください。

Principal axis regression (PAR)は、主成分分析と似たような手法であり、データの主要な変動方向を捉えるために使用されます (Bueno et al. 2023)。

Principal axis regressionの概要

Principal axis regressionは、データの共分散行列を計算し、その固有値と固有ベクトルを求めることで、データの主要な変動方向を特定します。これにより、次元削減や特徴抽出が可能になります。

Rでの実装

今回の実装では、次の順番で処理します。

  1. 種間差と種内差をもつダミーデータを作る
  2. 種ごとの代表値を使って、種間レベルの主軸を求める
  3. 各種から1個体ずつ取り出して、種内変異を含む主軸を繰り返し推定する
  4. 種間レベルの主軸と、種内変異を含む主軸の角度差を計算する
  5. 固有値の比を使って、種内変異がtrait spaceの各方向にどれくらい効いているかを可視化する

ここで使う x と y は、実データのtrait名ではなく、説明用の2次元traitです。 x を防御形質、y を栄養獲得形質のようなものとして読み替えると、Fig. 4の構造が理解しやすくなります。

データを作る

rnorm()関数を使用して種ごとの中心位置を作り、そのまわりに種内変異を持つ2次元のデータセットを作成します。 種内変異は、種ごとの中央値が中心位置に戻るように、対称なスコアとして作ります。 このデータセットは、2つのグループ（落葉樹っぽいグループと常緑樹っぽいグループ）を含み、各グループ内で種ごとに個体値がばらつくように設計されています。 ここでは、種平均のばらつきと、同じ種の中での個体差を別々の標準偏差で指定できるようにしています。 さらに、種平均が並ぶ方向と種内変異が伸びる方向を別々に指定します。 このようにすると、Fig. 4cのように「種間軸は斜めだが、種内変異を含む軸は別方向に固まる」データを作れます。

# 乱数を使うため、毎回同じ結果になるようにseedを固定します。
set.seed(123)

# major方向とminor方向の乱数を、指定した角度だけ回転させる補助関数です。
# angle_deg = 0ならx軸方向、angle_deg = 45なら右上がり45度方向になります。
rotate_xy <- function(major, minor, angle_deg) {
  theta <- angle_deg * pi / 180

  cbind(
    x = major * cos(theta) - minor * sin(theta),
    y = major * sin(theta) + minor * cos(theta)
  )
}

# 種内変異の値を、中央値が0になるように対称に作る関数です。
# こうしておくと、種ごとの中央値は種平均をよく反映します。
# その一方で、各種から1個体ずつ選ぶと、種内変異の方向が強く現れます。
balanced_scores <- function(n, target_sd) {
  if (target_sd == 0) {
    return(rep(0, n))
  }

  z <- seq(-1, 1, length.out = n)
  z <- z - mean(z)
  z / sd(z) * target_sd
}

sim_group <- function(
  group,
  n_species = 40,
  n_ind = 7,
  sp_major_sd,
  sp_minor_sd,
  sp_angle,
  ind_major_sd,
  ind_minor_sd,
  ind_angle
) {
  # 種ごとのデータを一時的に入れておくリストです。
  # 最後にrbindして、1つのdata.frameにまとめます。
  out <- list()

  # 1つのグループの中に、n_species個の種を作ります。
  for (i in seq_len(n_species)) {
    sp <- paste0(group, "_sp", i)

    # まず、種ごとの中心位置を決めます。
    # sp_major_sdは主方向の種間差、sp_minor_sdはそれに直交する方向のばらつきです。
    # sp_angleで、種平均がどの方向に並ぶかを決めます。
    sp_mean <- rotate_xy(
      major = rnorm(1, mean = 0, sd = sp_major_sd),
      minor = rnorm(1, mean = 0, sd = sp_minor_sd),
      angle_deg = sp_angle
    )

    # 次に、その種平均のまわりに個体値を発生させます。
    # ind_major_sdは主方向の種内差、ind_minor_sdはそれに直交する方向のばらつきです。
    # ここでは、各種内の値が種平均のまわりで対称になるようにしています。
    # これにより、種中央値は種間軸を保ちつつ、1個体サンプリングでは種内差が現れます。
    ind_major <- sample(balanced_scores(n_ind, ind_major_sd), size = n_ind)
    ind_minor <- sample(balanced_scores(n_ind, ind_minor_sd), size = n_ind)

    # ind_angleで、同じ種の中の個体差がどの方向に伸びるかを決めます。
    ind_dev <- rotate_xy(
      major = ind_major,
      minor = ind_minor,
      angle_deg = ind_angle
    )

    x <- as.numeric(sp_mean[, "x"]) + ind_dev[, "x"]
    y <- as.numeric(sp_mean[, "y"]) + ind_dev[, "y"]

    # group, species, individualを残しておくと、後で
    # 「種ごとの代表値」や「各種から1個体ずつサンプリング」がしやすくなります。
    out[[i]] <- data.frame(
      group = group,
      species = sp,
      individual = seq_len(n_ind),
      x = x,
      y = y
    )
  }

  # 種ごとのdata.frameを縦に結合します。
  do.call(rbind, out)
}

# 落葉樹っぽいグループです。
# 種間差も個体差も主にx方向に出るようにしています。
# そのため、種間軸と種内変異を含む軸はほとんど回転しません。
dat_D <- sim_group(
  group = "D",
  n_species = 40,
  n_ind = 7,
  sp_major_sd = 1.0,
  sp_minor_sd = 0.05,
  sp_angle = 0,
  ind_major_sd = 0.7,
  ind_minor_sd = 0.05,
  ind_angle = 0
)

# 常緑樹っぽいグループです。
# 種平均は斜め方向に並べ、個体差はx方向に強く出るようにしています。
# これにより、灰色の種間軸は青い種内軸からおよそ45度ずれます。
dat_E <- sim_group(
  group = "E",
  n_species = 40,
  n_ind = 7,
  sp_major_sd = 0.75,
  sp_minor_sd = 0.35,
  sp_angle = 55,
  ind_major_sd = 1.8,
  ind_minor_sd = 0.04,
  ind_angle = 0
)

# 2つのグループをまとめて、全データとして扱えるようにします。
dat <- rbind(dat_D, dat_E)

相関を確認する

それぞれのグループ内で、xとyの相関を確認してみましょう。 相関が小さいことを確認できます。 ここで相関を見ているのは、単純な「右上がり・右下がり」の関係が強くない場合でも、分散の向きとして主軸を求められることを示すためです。

tapply(seq_len(nrow(dat)), dat$group, function(i) {
  # tapply()でグループごとの行番号を取り出し、
  # その行だけを使ってxとyの相関係数を計算します。
  cor(dat$x[i], dat$y[i])
})
        D         E 
0.1066764 0.2091585

principal axis を求める関数

この関数では、2次元データから共分散行列を作り、その固有ベクトルをprincipal axisとして取り出します。 第1固有ベクトルはもっとも分散が大きい方向、第2固有ベクトルはそれに直交する方向です。

principal_axis <- function(x, y) {
  # xとyを2列の行列にまとめます。
  # この行列の各行が1個体または1種、各列が1つのtraitです。
  X <- cbind(x, y)

  # 共分散行列です。
  # 対角成分はxとyそれぞれの分散、非対角成分はx-yの共分散です。
  S <- cov(X)

  # 固有値分解を行います。
  # eigen()は、分散が大きい順に固有値と固有ベクトルを返します。
  eg <- eigen(S)

  vectors <- eg$vectors

  # 表示を安定させるため、x方向が正になるようにそろえる
  # 固有ベクトルは向きを反転しても同じ軸を表すため、
  # 図示のたびに線の向きが反転しないように符号をそろえます。
  for (j in seq_len(ncol(vectors))) {
    if (vectors[1, j] < 0) {
      vectors[, j] <- -vectors[, j]
    }
  }

  v1 <- vectors[, 1]

  # 後の計算と作図で使う値をまとめて返します。
  list(
    covariance = S,
    eigenvalues = eg$values,
    vectors = vectors,
    first_vector = v1,
    first_eigenvalue = eg$values[1],
    center = colMeans(X),
    # x軸から見た第1主軸の角度です。
    # atan2()を使うと、ベクトルのx成分・y成分から角度を直接計算できます。
    angle_from_x = atan2(v1[2], v1[1]) * 180 / pi
  )
}

angle_between_axes <- function(v1, v2) {
  # 2つのベクトルのなす角を内積から計算します。
  # abs()を使うことで、向きが反対の同じ軸を0度として扱います。
  cos_angle <- abs(sum(v1 * v2) / sqrt(sum(v1^2) * sum(v2^2)))
  acos(cos_angle) * 180 / pi
}

principal_axis() は、要するに以下の処理を行っています。

eigen(cov(cbind(x, y)))

種平均の trait space を作る

論文では、mean–variance effect を避けるために log変換後の中央値を使っているようですが、今回は単純に中央値をとってみます。 ここで作る species_med は、種ごとに1点だけを持つデータです。 したがって、このデータから求める主軸は、種間差を反映した interspecific axis とみなせます。

species_med <- aggregate(
  # xとyを同時に集計します。
  cbind(x, y) ~ group + species,
  data = dat,
  # 平均ではなく中央値を使います。
  # 外れ値に少し強い代表値として扱えます。
  FUN = median
)
head(species_med)
  group species           x           y
1     D   D_sp1 -0.56047565 -0.01150887
2     D  D_sp10  0.54839696  0.01193659
3     D  D_sp11  0.03768285  0.05881101
4     D  D_sp12 -0.05556197  0.02597036
5     D  D_sp13  1.87141801  0.03014352
6     D  D_sp14 -0.42249683 -0.10266236

種平均ベースの主軸 interspecific axis

まずは落葉樹っぽいグループと常緑樹っぽいグループを分けて、種平均ベースの主軸を求めます。 この時点では、各種の中に複数個体があることは無視され、各種が1点で表されています。

pa_inter_D <- with(
  # 落葉樹っぽいグループの種中央値だけを使います。
  subset(species_med, group == "D"),
  principal_axis(x, y)
)

pa_inter_E <- with(
  # 常緑樹っぽいグループの種中央値だけを使います。
  subset(species_med, group == "E"),
  principal_axis(x, y)
)

# x軸から見た第1主軸の角度を確認します。
pa_inter_D$angle_from_x
[1] 0.7338449
pa_inter_E$angle_from_x
[1] 47.56799

これが、論文でいう interspecific eigenvector に相当します。

各種から1個体ずつランダムに選ぶ

論文では、種ごとのサンプル数の不均衡を避けるため、各反復で「各種から1個体」を選び、これを50回繰り返しています。 この処理により、種数はそのままにしつつ、種内変異を反映した個体レベルのtrait spaceを繰り返し作れます。

sample_one_per_species <- function(d) {
  # speciesごとにdata.frameを分割します。
  sp_list <- split(d, d$species)

  sampled <- lapply(sp_list, function(z) {
    # 各種の中から、ランダムに1行だけ選びます。
    z[sample(seq_len(nrow(z)), 1), ]
  })

  # 種ごとに選ばれた1個体を、1つのdata.frameに戻します。
  do.call(rbind, sampled)
}

50回リサンプリングして rotation angle を計算する

run_resampling() は、Fig. 4の上段と下段に使う値をまとめて計算する関数です。 入力は、個体レベルのデータ dat_group と、種代表値のデータ species_med_group です。 出力には、種間軸、リサンプリングごとの種内軸、軸の角度差、固有値比が入ります。

run_resampling <- function(dat_group, species_med_group, n_iter = 50) {
  # 種代表値から求める主軸です。
  # これを基準となるinterspecific axisとして使います。
  pa_inter <- principal_axis(
    species_med_group$x,
    species_med_group$y
  )

  # 反復ごとの結果を入れるベクトルとリストを用意します。
  angles <- numeric(n_iter)
  eig_ratio_dim1 <- numeric(n_iter)
  eig_ratio_dim2 <- numeric(n_iter)
  axes <- vector("list", n_iter)
  samples <- vector("list", n_iter)

  for (i in seq_len(n_iter)) {
    # 各種から1個体ずつ選びます。
    # この1セットが、種内変異を反映したtrait spaceになります。
    sampled <- sample_one_per_species(dat_group)

    # サンプリングされた個体データから主軸を求めます。
    # これをintraspecific variationを含んだ軸として扱います。
    pa_intra <- principal_axis(
      sampled$x,
      sampled$y
    )

    # 種間軸と、種内変異を含む軸の角度差です。
    # Fig. 4上段の角度に対応する量です。
    angles[i] <- angle_between_axes(
      pa_inter$first_vector,
      pa_intra$first_vector
    )

    # 第1軸方向の分散が、種間レベルに対してどれくらい大きいかを見ます。
    # 1より大きければ、種内変異を含めたときにその方向の分散が増えています。
    eig_ratio_dim1[i] <- pa_intra$eigenvalues[1] /
      pa_inter$eigenvalues[1]

    # 第2軸方向についても同じ比を計算します。
    eig_ratio_dim2[i] <- pa_intra$eigenvalues[2] /
      pa_inter$eigenvalues[2]

    # 後で作図できるように、軸そのものとサンプルも保存します。
    axes[[i]] <- pa_intra
    samples[[i]] <- sampled
  }

  list(
    pa_inter = pa_inter,
    angles = angles,
    eig_ratio = eig_ratio_dim1,
    eig_ratio_dim1 = eig_ratio_dim1,
    eig_ratio_dim2 = eig_ratio_dim2,
    axes = axes,
    samples = samples
  )
}

res_all <- run_resampling(
  # 全データをまとめた場合です。
  dat_group = dat,
  species_med_group = species_med,
  n_iter = 50
)

res_D <- run_resampling(
  # 落葉樹っぽいグループだけで同じ処理を行います。
  dat_group = subset(dat, group == "D"),
  species_med_group = subset(species_med, group == "D"),
  n_iter = 50
)

res_E <- run_resampling(
  # 常緑樹っぽいグループだけで同じ処理を行います。
  dat_group = subset(dat, group == "E"),
  species_med_group = subset(species_med, group == "E"),
  n_iter = 50
)

結果の確認

まず、角度差と固有値比の分布を数値で確認します。 angles は種間軸からどれくらい回転したかを表し、eig_ratio_dim1 と eig_ratio_dim2 は第1軸・第2軸それぞれの分散の増え方を表します。

# 主軸の回転角です。
summary(res_all$angles)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  11.13   13.44   14.81   15.23   16.89   21.29 
summary(res_D$angles)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.03909 0.30727 0.50888 0.53914 0.69654 1.46247 
summary(res_E$angles)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  37.12   42.14   43.19   43.40   44.84   50.02 
# 第1軸方向の固有値比です。
summary(res_all$eig_ratio_dim1)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  3.004   3.585   3.886   3.929   4.267   4.756 
summary(res_D$eig_ratio_dim1)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.087   1.465   1.669   1.659   1.792   2.217 
summary(res_E$eig_ratio_dim1)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  3.157   4.630   5.193   5.183   5.775   6.892 
# 第2軸方向の固有値比です。
summary(res_all$eig_ratio_dim2)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.043   1.173   1.217   1.217   1.276   1.339 
summary(res_D$eig_ratio_dim2)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.107   1.366   1.597   1.589   1.771   2.097 
summary(res_E$eig_ratio_dim2)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.245   2.807   2.867   2.840   2.952   3.083

Fig. 4 のようにまとめて描く

Zhou et al. (2025) のFig. 4では、上段でtrait spaceの回転、下段でinterspecificな固有値に対するintraspecificな固有値の比を示しています。 ここでも同じ構成で、全データ、落葉樹っぽいグループ、常緑樹っぽいグループを並べてみます。

# 図全体で使う色です。
# all, D, Eは点と輪郭の色、intraとinterは主軸の色です。
cols <- c(
  all = "#e41a1c",
  D = "#e6ab02",
  E = "#1b9e77",
  intra = "#6bb7de",
  inter = "#6b6b6b"
)

# principal_axis()の結果から、指定した固有ベクトルの方向に線を描く関数です。
# dim = 1なら第1軸、dim = 2なら第2軸を描きます。
draw_axis <- function(pa, dim = 1, len = 2, col = "black", lwd = 3, lty = 1) {
  # 線の中心は、データ点の中心に置きます。
  c0 <- pa$center

  # 固有ベクトルは長さ1の方向ベクトルです。
  # lenをかけることで、図上で見やすい長さに伸ばします。
  v <- pa$vectors[, dim]

  segments(
    x0 = c0[1] - len * v[1],
    y0 = c0[2] - len * v[2],
    x1 = c0[1] + len * v[1],
    y1 = c0[2] + len * v[2],
    col = col,
    lwd = lwd,
    lty = lty
  )
}

# 散布図の点がどこに集中しているかを、輪郭線で示す関数です。
# MASS::kde2d()で2次元カーネル密度を推定し、その等高線を重ねます。
draw_density <- function(x, y, col, xlim, ylim) {
  # MASSが使えない環境でも本文の主要な計算は止めないようにします。
  if (!requireNamespace("MASS", quietly = TRUE)) {
    return(invisible(NULL))
  }

  # limsを指定して、3つのtrait panelで同じ座標範囲の密度を計算します。
  dens <- MASS::kde2d(x, y, n = 80, lims = c(xlim, ylim))
  z <- as.vector(dens$z)
  z <- z[z > 0]

  # 密度の高い部分だけを輪郭線として描きます。
  # Fig. 4の分布のかたまりを示す線に近い役割です。
  contour(
    dens,
    add = TRUE,
    levels = quantile(z, probs = c(0.75, 0.9), names = FALSE),
    drawlabels = FALSE,
    col = col,
    lwd = c(1, 1.4)
  )
}

# 各trait panelの左下に、Dim.1とDim.2の方向を示す小さな矢印を描きます。
# これは実際の推定軸ではなく、図の読み方を助けるための模式的な矢印です。
draw_dim_arrows <- function(xlim, ylim) {
  # パネルの座標範囲から、矢印の開始位置と長さを相対的に決めます。
  x0 <- xlim[1] + 0.08 * diff(xlim)
  y0 <- ylim[1] + 0.13 * diff(ylim)
  x1 <- x0 + 0.20 * diff(xlim)
  y1 <- y0 + 0.20 * diff(ylim)

  arrows(x0, y0, x1, y0, length = 0.08, lwd = 1.4)
  arrows(x0, y0, x0, y1, length = 0.08, lwd = 1.4)
  text(x1, y0 + 0.06 * diff(ylim), "Dim.1", cex = 0.85, adj = c(1, 0))
  text(x0 + 0.06 * diff(xlim), y1, "Dim.2", cex = 0.85, srt = 90, adj = c(1, 1))
}

# リサンプリングで得た角度分布を、Fig. 4風のラベルに整形します。
# ここでは簡易的に、2.5%点と97.5%点を信頼区間のように表示しています。
angle_label <- function(angles) {
  ci <- quantile(angles, probs = c(0.025, 0.975), names = FALSE)

  # 下限が0より大きければ、0度からずれているとみなしてP < 0.05と表示します。
  # 厳密な検定ではなく、このシミュレーション図用の簡易表示です。
  p_text <- if (ci[1] > 0) "P < 0.05" else "n.s."

  sprintf(
    "%.2f°  %s\n(%.2f - %.2f)",
    mean(angles),
    p_text,
    ci[1],
    ci[2]
  )
}

# 上段のtrait space panelを1枚描く関数です。
# 点、密度輪郭、種内軸、種間軸、角度ラベルをまとめて描きます。
plot_trait_panel <- function(
  dat_group,
  res,
  main,
  panel,
  col,
  xlim,
  ylim,
  show_y_label = FALSE,
  show_legend = FALSE
) {
  # まず個体値の散布図を描きます。
  # asp = 1にして、x方向とy方向の1単位が同じ長さになるようにしています。
  plot(
    dat_group$x,
    dat_group$y,
    pch = 16,
    cex = 0.65,
    col = adjustcolor(col, alpha.f = 0.8),
    xlim = xlim,
    ylim = ylim,
    asp = 1,
    xlab = "Defense-like trait",
    ylab = if (show_y_label) "Nutrient-like trait" else "",
    main = paste0("(", panel, ")  ", main),
    bty = "l",
    las = 1
  )

  # 点の分布のかたまりを輪郭線で重ねます。
  draw_density(
    dat_group$x,
    dat_group$y,
    col = col,
    xlim = xlim,
    ylim = ylim
  )

  # 主軸の線の長さです。
  # 座標範囲に対する相対値として決めると、パネル間で見た目がそろいます。
  axis_len <- 0.42 * min(diff(xlim), diff(ylim))

  # 各リサンプリングで推定された種内の第1軸を薄い青で重ねます。
  # 50本重ねることで、種内変異を含む主軸がどの方向に固まるかが見えます。
  for (pa in res$axes) {
    draw_axis(
      pa,
      dim = 1,
      len = axis_len,
      col = adjustcolor(cols["intra"], alpha.f = 0.16),
      lwd = 1
    )
  }

  # 種代表値から求めた種間軸を灰色の太線で描きます。
  # dim = 1がrotation angleの比較対象で、dim = 2はそれに直交する軸です。
  # 薄青の第1軸と灰色の第1軸のずれが、rotation angleの直感的な意味です。
  draw_axis(res$pa_inter, dim = 1, len = axis_len, col = cols["inter"], lwd = 3)
  draw_axis(res$pa_inter, dim = 2, len = axis_len, col = cols["inter"], lwd = 3)

  # 左下のDim.1 / Dim.2の補助矢印を描きます。
  draw_dim_arrows(xlim, ylim)

  # 左上に角度と範囲を表示します。
  text(
    xlim[1] + 0.03 * diff(xlim),
    ylim[2] - 0.06 * diff(ylim),
    angle_label(res$angles),
    adj = c(0, 1),
    cex = 0.85,
    font = 2
  )

  # 最初のpanelだけ凡例を表示します。
  if (show_legend) {
    legend(
      "bottomright",
      legend = c("Intraspecies", "Interspecies"),
      col = c(cols["intra"], cols["inter"]),
      lwd = c(3, 3),
      bty = "n",
      cex = 0.8
    )
  }
}

# 下段の固有値比のboxplotを描く関数です。
# valuesは、All data, Deciduous, Evergreenの3つの数値ベクトルを入れたlistです。
plot_ratio_panel <- function(values, panel, main, colors, ylim = NULL) {
  # y軸範囲が指定されていない場合は、データと基準線1が入る範囲を自動で決めます。
  if (is.null(ylim)) {
    ylim <- range(c(unlist(values), 1))
    ylim <- ylim + c(-0.12, 0.12) * diff(ylim)
  }

  # 箱ひげ図で、固有値比の分布をグループ間で比較します。
  boxplot(
    values,
    names = c("All data", "Deciduous", "Evergreen"),
    col = adjustcolor(colors, alpha.f = 0.85),
    border = colors,
    outline = FALSE,
    ylim = ylim,
    ylab = "Eigenvalue intraspecies / interspecies",
    main = paste0("(", panel, ")  ", main),
    las = 1
  )

  # 箱ひげ図の上に、個々のリサンプリング結果を点で重ねます。
  # ばらつきの大きさを見やすくするため、横方向に少しjitterします。
  for (i in seq_along(values)) {
    stripchart(
      values[[i]],
      at = i,
      vertical = TRUE,
      method = "jitter",
      jitter = 0.12,
      pch = 16,
      cex = 0.6,
      col = adjustcolor(colors[i], alpha.f = 0.35),
      add = TRUE
    )
  }

  # ratio = 1は、種内変異を入れても分散が増えていない基準線です。
  # 1より上に分布すれば、その軸方向で種内変異の寄与が大きいと読めます。
  abline(h = 1, lty = 2, lwd = 1.4)

  # Fig. 4の見た目に合わせて、基準線の近くにアスタリスクを置いています。
  text(seq_along(values), rep(1, length(values)), "*", pos = 3, cex = 1.2)
}

# 下段の3グループに対応する色です。
plot_cols <- c(cols["all"], cols["D"], cols["E"])

# 上段3panelで同じxlim/ylimを使うため、全データの範囲から共通範囲を作ります。
lim <- range(c(dat$x, dat$y))
lim <- lim + c(-0.08, 0.08) * diff(lim)

# par()を変更するので、最後に元へ戻せるように保存しておきます。
old_par <- par(no.readonly = TRUE)

# layout()で5枚のpanel配置を作ります。
# 上段は3枚、下段は中央寄せで2枚です。
layout(
  matrix(
    c(1, 1, 2, 2, 3, 3, 0, 4, 4, 5, 5, 0),
    nrow = 2,
    byrow = TRUE
  ),
  heights = c(1, 1.1)
)

# 余白をやや詰めて、複数panelを1枚の図として見やすくします。
par(mar = c(4, 4, 2, 1), oma = c(0, 0, 0, 0))

# 上段: trait spaceの回転を示す3panelです。
plot_trait_panel(
  dat,
  res_all,
  "All data",
  "a",
  cols["all"],
  lim,
  lim,
  TRUE,
  TRUE
)
plot_trait_panel(
  subset(dat, group == "D"),
  res_D,
  "Deciduous",
  "b",
  cols["D"],
  lim,
  lim
)
plot_trait_panel(
  subset(dat, group == "E"),
  res_E,
  "Evergreen",
  "c",
  cols["E"],
  lim,
  lim
)

# 下段左: 第1軸方向の固有値比です。
plot_ratio_panel(
  list(res_all$eig_ratio_dim1, res_D$eig_ratio_dim1, res_E$eig_ratio_dim1),
  "d",
  "Dim1",
  plot_cols
)

# 下段右: 第2軸方向の固有値比です。
plot_ratio_panel(
  list(res_all$eig_ratio_dim2, res_D$eig_ratio_dim2, res_E$eig_ratio_dim2),
  "e",
  "Dim2",
  plot_cols
)

# 変更した作図設定を元に戻します。
par(old_par)

References

Bueno, C. Guillermo, Aurele Toussaint, Sabrina Träger, et al. 2023. “Reply to: The Importance of Trait Selection in Ecology.” Nature 618 (7967): E31–34. https://doi.org/10.1038/s41586-023-06149-7.
Zhou, Guangkai, Daniel F. Petticord, Xingchang Wang, et al. 2025. “Intraspecific Variation in the Growth–Defense Trade-Off Among Deciduous and Evergreen Broadleaf Woody Plants.” New Phytologist, ahead of print, November 23. https://doi.org/10.1111/nph.70781.