4.3 Pre-network permutation
4.3.1 Pre-network permutationの手法
Node permutationとともによく用いられるパーミュテーションの方法が、pre-network (またはdata stream) permutationである。この方法では、隣接行列(ネットワーク)を作成する前の生データに対してパーミュテーションを行う(Farine, 2017)。
Pre-network permutationは、もともと“group by individual”(第2.2節参照)形式のデータに適用するために考案された手法である(Bejder et al., 1998; Whitehead et al., 1999)。このような形式のデータとしては、下のようなものがあげられる。これは、金華山島\(B_1\)群で2021年交尾期の各観察日(date)にそれぞれのメスが群れ内で観察されたか否かを示したものである。1が観察されたこと、0が観察されなかったことを表す。
Pre-network permutationのプロセスを模式的に表したのが図4.11である(Bejder et al., 1998; Whitehead et al., 1999)。以下のプロセスで行う。
“group by individual”形式のデータからランダムに2つのグループ(行)と2個体を選ぶ。ただし、各個体はいずれか一方のグループ(行)のみで確認され、かつ各グループ(行)ではいずれか一方の個体のみしか確認されていないものの中から選ぶ。例えば、図4.11の左の行列の下線の組み合わせがこのような条件を満たす。
2つのグループ(行)における各個体の数値を入れ替える(図4.11の右の行列)。
入れ替えた後のデータをもとにassociation index(e.g., SRIなど)を算出し、隣接行列(ネットワーク)を作成する。
2~3の過程を少なくとも1000回繰り返す。ただし、各ステップで入れ替えを行うデータは、前のステップで入れ替えを行った後のデータである。例えば図4.11では、次のステップにおいては右の行列(= 前のステップで入れ替えた後のデータ)に対してランダムな入れ替えを同様に行う。
図4.11: Pre-network permutationのプロセス。Whitehead et al. (1999)から引用。
この方法を用いることで、各グループ(行)で確認された個体数と、各個体が確認された回数を保ったままランダムなネットワークを無数に生成することができる(Farine, 2017)。また、ランダムに選ぶ2個体を同じ属性(e.g., 性別)の個体に限定したり、観察が複数の調査期間/場所で行われた際にそれぞれの調査期間/場所内で入れ替えを行ったりすることで、node permutationで問題となっていた観察バイアスや観察場所・期間の影響を排除することができる(Whitehead et al., 2005)。
検定を行う際には、node permutationのときと同様に、実データをもとに算出した統計検定量と、各ステップでランダムな入れ替えによって生成されたデータから算出された統計検定量の分布を比較することでP値を算出する。
なお、pre-network permutationは近年適用範囲が広がり、“group by individual”形式のデータだけでなく、GPSデータや個体追跡データに対しても適用できるようになってきている。詳しい方法については、Farine (2017) を参照。
4.3.2 分析例
Pre-network permutationを用いた分析は、node permutationで行ったすべての分析(線形モデル、相関分析、MRQAP検定など)に対して適用することができる。
ここでは、簡単な例としてGLMによる分析を行う。
下のデータは、金華山島\(B_1\)群で2019～2021年において、各観察日(date)に各メスが観察されたか否かを”group by individual”形式で示したものである(presence_multi)。
このデータを用いてSRI(第2.2節参照)を算出し、SRIネットワークにおける各個体の固有ベクトル中心性(eigen)が順位(rank)と年齢(age)と関連するかを調べたいとする。 なお、pre-network permutationは、各観察年(period)内で行うものとする。
ここでは、以下のような一般化線形モデルを考える。
- 分布: ガンマ分布
- リンク関数: log関数
- 応答変数: 固有ベクトル中心性(
eigen)
- 説明変数: 順位(
rank)、年齢(age)
4.3.2.1 snaパッケージを用いる方法
まずは実データからSRIを算出して隣接行列を作成する(第2.2参照)。
##
presence_mat_multi <- get_network(presence_multi
## 日付と調査期間は除く
%>% select(-date, -period),
data_format = "GBI",
association_index = "SRI")## Generating 15 x 15 matrixそこからANTsパッケージで固有ベクトル中心性を算出し(第3.1節参照)、メスの属性データに追加する。
## 属性データ
att_females <- data.frame(femaleID = colnames(prox_mat),
rank = c(8,9,11,12,1,2,3,5,4,13,15,14,6,7,10),
age = c(12,9,10,10,14,7,12,8,6,8,10,15, 10, 10, 10))
att_eigen <- met.eigen(presence_mat_multi,
df = att_females,
dfid = 1)算出したものがこちら。
その後、node permutationの時と同様に実データに対してモデリングを行い、統計検定量(= 係数の推定値)を求める。
推定結果は以下のようになる。
##
## Call: glm(formula = eigen ~ rank + age, family = Gamma("log"), data = att_eigen)
##
## Coefficients:
## (Intercept) rank age
## 0.237425 -0.006943 -0.027929
##
## Degrees of Freedom: 14 Total (i.e. Null); 12 Residual
## Null Deviance: 0.3288
## Residual Deviance: 0.247 AIC: -14.21係数の推定値を抽出する。
それでは、pre-network permutationを行い、それによって生成されたランダムなネットワークに対しても同様の分析を行って係数の推定値を抽出していく。パーミュテーションは10000回行うものとする。
snaパッケージでは、random_permutation関数を用いてpre-network permutationを行うことができる。関数では、locations =、days =にデータの各行の観察場所/時間を、classes =に各個体の属性を指定できる。そのうえで、within_location = TRUE、within_day = TRUE、within_class = TRUEとすることで、各観察場所/期間内、あるいは同じ属性の個体同士でしかデータの入れ替えが行われないようにすることができる。今回は各観察年内で入れ替えを行うので、days = periodとし、within_period = TRUEとする。
period <- presence_multi$period
set.seed(123)
presence_mat_rand <- network_permutation(presence_multi
%>% select(-date,-period),
data_format = "GBI",
association_index = "SRI",
#パーミュテーション回数
permutations = 10000,
days = period,
within_day = TRUE)## No association matrix provided, generating 15 x 15 matrix
## Generating 15 x 15 matrix
## Starting permutations, generating 10000 x 15 x 15 matrixこれで、ランダムなネットワークが10000個生成された。例えば、120個目のネットワークの隣接行列は以下のようになる(表4.7)。
presence_mat_rand[120, , ] %>%
kable(digits = 2, align = "c", caption = "120個目の隣接行列") %>%
kable_styling(font_size = 6, full_width = FALSE) | 0.00 | 0.63 | 0.46 | 0.44 | 0.47 | 0.47 | 0.53 | 0.45 | 0.45 | 0.47 | 0.45 | 0.54 | 0.45 | 0.46 | 0.46 |
| 0.63 | 0.00 | 0.65 | 0.61 | 0.65 | 0.67 | 0.74 | 0.66 | 0.65 | 0.65 | 0.65 | 0.52 | 0.63 | 0.66 | 0.68 |
| 0.46 | 0.65 | 0.00 | 0.89 | 0.90 | 0.92 | 0.82 | 0.93 | 0.95 | 0.84 | 0.95 | 0.57 | 0.86 | 0.91 | 0.85 |
| 0.44 | 0.61 | 0.89 | 0.00 | 0.83 | 0.85 | 0.75 | 0.87 | 0.88 | 0.80 | 0.88 | 0.54 | 0.79 | 0.84 | 0.81 |
| 0.47 | 0.65 | 0.90 | 0.83 | 0.00 | 0.94 | 0.86 | 0.90 | 0.90 | 0.88 | 0.92 | 0.63 | 0.81 | 0.90 | 0.81 |
| 0.47 | 0.67 | 0.92 | 0.85 | 0.94 | 0.00 | 0.86 | 0.92 | 0.94 | 0.86 | 0.94 | 0.60 | 0.84 | 0.90 | 0.84 |
| 0.53 | 0.74 | 0.82 | 0.75 | 0.86 | 0.86 | 0.00 | 0.85 | 0.83 | 0.85 | 0.83 | 0.57 | 0.75 | 0.84 | 0.80 |
| 0.45 | 0.66 | 0.93 | 0.87 | 0.90 | 0.92 | 0.85 | 0.00 | 0.97 | 0.82 | 0.96 | 0.58 | 0.85 | 0.93 | 0.87 |
| 0.45 | 0.65 | 0.95 | 0.88 | 0.90 | 0.94 | 0.83 | 0.97 | 0.00 | 0.83 | 0.96 | 0.57 | 0.87 | 0.94 | 0.87 |
| 0.47 | 0.65 | 0.84 | 0.80 | 0.88 | 0.86 | 0.85 | 0.82 | 0.83 | 0.00 | 0.83 | 0.54 | 0.75 | 0.81 | 0.80 |
| 0.45 | 0.65 | 0.95 | 0.88 | 0.92 | 0.94 | 0.83 | 0.96 | 0.96 | 0.83 | 0.00 | 0.58 | 0.88 | 0.94 | 0.87 |
| 0.54 | 0.52 | 0.57 | 0.54 | 0.63 | 0.60 | 0.57 | 0.58 | 0.57 | 0.54 | 0.58 | 0.00 | 0.57 | 0.60 | 0.56 |
| 0.45 | 0.63 | 0.86 | 0.79 | 0.81 | 0.84 | 0.75 | 0.85 | 0.87 | 0.75 | 0.88 | 0.57 | 0.00 | 0.84 | 0.86 |
| 0.46 | 0.66 | 0.91 | 0.84 | 0.90 | 0.90 | 0.84 | 0.93 | 0.94 | 0.81 | 0.94 | 0.60 | 0.84 | 0.00 | 0.84 |
| 0.46 | 0.68 | 0.85 | 0.81 | 0.81 | 0.84 | 0.80 | 0.87 | 0.87 | 0.80 | 0.87 | 0.56 | 0.86 | 0.84 | 0.00 |
それでは、各ランダムネットワークに対してモデリングを行い、係数の推定値を抽出していく。b_rank_randとb_rank_ageという空のベクトルを作成し、そこにランダムネットワークから算出された係数の推定値を入れている。
## 空の行列を作成
b_rank_rand <- rep(0,10000)
b_age_rand <- rep(0,10000)
for(i in 1:10000){
att_eigen_rand <- met.eigen(presence_mat_rand[i,,], df = att_females)
r_glm_rand <- glm(data = att_eigen_rand,
formula = eigen ~ rank + age,
family = Gamma("log"))
b_rank_rand[i] <- coef(r_glm_rand)[[2]]
b_age_rand[i] <- coef(r_glm_rand)[[3]]
}それでは、最後に実データの回帰係数をランダムネットワークの回帰係数の分布と比較する(図4.12)。灰色のヒストグラムがランダムネットワークの回帰係数の分布を、赤い線が実データの回帰係数を表している。
図4.12: 実データの回帰係数とランダムネットワークの回帰係数の分布の比較
そこで、実データの回帰係数がランダムネットワークの回帰係数の何%より大きいか/小さいかを計算することで、\(P\)値を算出する。
算出は以下のように行える。
順位の係数については、P = 0.067なので、有意な関連はないよう(有意傾向)。
## [1] 0.0674一方で、年齢の係数(\(\beta_2\))については、P = 0.000なので、有意である。
## [1] 04.3.2.2 ANTsパッケージを用いる方法
ANTsパッケージでも同じ同様の分析を行うことができる。ANTsパッケージではpre-network permutationを行う前に”group by individual”形式のデータをgbi.to.df関数で以下のように縦長のデータフレーム(presence_df)に直す必要がある。
## まずはマトリックスに直す。
presence_multi %>%
column_to_rownames(var = "date") %>%
select(-period) %>%
as.matrix()-> presence_multi_b
## gbi.to.df関数でデータフレーム形式に
presence_df <- gbi.to.df(presence_multi_b)
## periodを追加
presence_df %>%
mutate(period = ifelse(year(scan) == "2019",1,
ifelse(year(scan) == "2020",2,3))) -> presence_dfこのように、データフレームには観察した日付(scan)と観察された個体の名前(ID)が記されている。
このデータフレームをもとに、perm.ds.grp関数でパーミュテーションを行うことができる。この際、ctrlf = periodと指定することで、各観察年内で入れ替えを行うことができる。
set.seed(123)
perm_ds <- perm.ds.grp(presence_df_b,
## 観察日
scan = "scan",
## パーミュテーション回数
nperm = 10000,
ctrlf = "period",
index = "sri",
progress = FALSE)これで、実データにもとづくネットワーク(隣接行列)1つと、ランダムネットワークが10000個生成された。各ネットワークにおける固有ベクトル中心性を算出し、それぞれをメスの属性データに追加していく。
## メスの属性データ
att_females <- data.frame(femaleID = colnames(prox_mat),
rank = c(8,9,11,12,1,2,3,5,4,13,15,14,6,7,10),
age = c(12,9,10,10,14,7,12,8,6,8,10,15, 10, 10, 10))
att_eigen_list <- met.eigen(perm_ds,
df = att_females,
dfid = 1)続いてstat.glm関数を用い、実データとランダムネットワークにおけるデータに対してモデリングを行い、パラメータの推定値を算出する。
glm_ants <- stat.glm(ant = att_list,
formula = eigen ~ age + rank,
family = Gamma(link = "log"),
## 実データのデータフレーム
oda = presence_df,
progress = FALSE)## Original model :
##
##
## Call:
## c("eigen ~ age + rank , family = Gamma", "eigen ~ age + rank , family = log",
## "eigen ~ age + rank , family = function (mu) \nlog(mu)", "eigen ~ age + rank , family = function (eta) \npmax(exp(eta), .Machine$double.eps)",
## "eigen ~ age + rank , family = function (mu) \nmu^2", "eigen ~ age + rank , family = function (y, mu, wt) \n-2 * wt * (log(ifelse(y == 0, 1, y/mu)) - (y - mu)/mu)",
## "eigen ~ age + rank , family = function (y, n, mu, wt, dev) \n{\n n <- sum(wt)\n disp <- dev/n\n -2 * sum(dgamma(y, 1/disp, scale = mu * disp, log = TRUE) * wt) + 2\n}",
## "eigen ~ age + rank , family = function (eta) \npmax(exp(eta), .Machine$double.eps)",
## "eigen ~ age + rank , family = expression({\n if (any(y <= 0)) \n stop(\"non-positive values not allowed for the 'Gamma' family\")\n n <- rep.int(1, nobs)\n mustart <- y\n})",
## "eigen ~ age + rank , family = function (mu) \nall(is.finite(mu)) && all(mu > 0)",
## "eigen ~ age + rank , family = function (eta) \nTRUE", "eigen ~ age + rank , family = function (object, nsim) \n{\n wts <- object$prior.weights\n if (any(wts != 1)) \n message(\"using weights as shape parameters\")\n ftd <- fitted(object)\n shape <- MASS::gamma.shape(object)$alpha * wts\n rgamma(nsim * length(ftd), shape = shape, rate = shape/ftd)\n}",
## "eigen ~ age + rank , family = NA")
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.237425 0.159228 1.491
## age -0.027929 0.014937 -1.870
## rank -0.006943 0.008129 -0.854
##
## (Dispersion parameter for Gamma family taken to be 0.01815169)
##
## Null deviance: 0.32879 on 14 degrees of freedom
## Residual deviance: 0.24701 on 12 degrees of freedom
## AIC: -14.212
##
## Number of Fisher Scoring iterations: 5
##
## in分析結果は、ant関数を用いることで確認できる。
モデルの診断プロット（残差 vs 予測値、QQプロット）および実データの回帰係数とランダムネットワークの回帰係数の分布の比較は以下のように可視化できる。
図4.13: 実データの回帰係数とランダムネットワークの回帰係数の分布の比較
図4.14: 実データの回帰係数とランダムネットワークの回帰係数の分布の比較
検定結果は以下の通り。p.leftとp_rightはランダムネットワークの係数が実際の係数よりも小さい/大きい確率を表す。snaパッケージを用いた場合と少し異なる結果が得られており、順位の関連も有意である。この違いはパッケージごとのアルゴリズムの違いに起因？
## Estimate Std. Error t value p.left p.rigth p.two.sides
## (Intercept) 0.237424995 0.159227844 1.4911022 0.9999 0.0001 2e-04
## age -0.027928830 0.014936798 -1.8698004 0.0000 1.0000 0e+00
## rank -0.006942727 0.008129111 -0.8540574 0.0150 0.9850 3e-02
## lower.ci uper.ci mean
## (Intercept) 0.198183854 0.221152612 0.209428988
## age -0.026096189 -0.023854130 -0.024945476
## rank -0.006744554 -0.005549161 -0.0061203014.3.3 Pre-network permutationの欠点
以上のように、pre-network permutationは観察バイアスや観察場所/期間の影響を排除することのできる強力な手法であり、多くの研究で用いられている(Farine, 2017; Farine & Carter, 2022; Farine & Whitehead, 2015)。一方で、ここ数年の研究によってpre-network permutationは相関分析や回帰分析の検定を行う手法としては適切でないことが指摘されるようになっている(Farine & Carter, 2022; Puga‐Gonzalez et al., 2021; Weiss et al., 2021)。
回帰分析を例にとって考えてみよう。回帰分析の目的は説明変数(X)と応答変数(Y)に関連があるかを調べることなので、その帰無仮説は「両者の間に関連がない(\(\beta = 0\))」である。Node permutationでは隣接行列の行と列に対してパーミュテーションを行うことで、ネットワークの構造自体は保ちつつ、XとYの関連がランダムなネットワークを作り出すことができていた(第4.1節と第4.2節を参照)。
一方で、pre-network permutationは生のassociation/交渉データに対してパーミュテーションを行うので、ネットワークの構造自体に大きな変化が生じてしまい、辺の重みや中心性指標の分布が大きく変わってしまう(Weiss et al., 2021)。その結果、現実世界における辺の重みや中心性指標の分布からは乖離してしまい、適切な帰無分布にならなくなってしまう可能性がある(Farine & Carter, 2022)。一般的に、pre-network permutationによって辺の重みや中心性指標の分布のばらつきは小さくなってしまうので(つまり、実際のデータの値のほうが極端になる傾向があるので)、第一種の過誤を犯す確率が高くなってしまう(図4.15)。
図4.15: Pre-network permutationを繰り返すことで値の分布のばらつきが小さくなる。Weiss et al. (2021)から引用。
この問題を解消するため、Farine & Carter (2022) はnode permutationとpre-network permutationの手法を組み合わせたdouble permutation法を提唱し、この手法を用いることでそれぞれの手法を単独で用いるよりも誤った結論を導いてしまう可能性が低くなることを指摘している。一方で Weiss et al. (2021) は、①Generalized affiliation indicesなどの交絡要因をはじめから除外できる指標を用いること、②観察バイアスや観察場所や期間の情報を明示的にモデルに取り込むこと(e.g., 観察回数を説明変数に入れたり、ランダム効果に場所や期間を入れるなど)(c.f. Franks et al., 2021)でnode permutationだけでも適切な検定ができることを指摘している。
いずれにしても、ここ数年だけでもパーミュテーション検定の妥当性を検証する論文が多数出ていることからもわかるように(Evans et al., 2020; Farine & Carter, 2022; Franks et al., 2021; Hart et al., 2022; Puga‐Gonzalez et al., 2021; Weiss et al., 2021)、どのような分析方法が最も適切かについてはまだ結論が出ていないのが現状だ。今後新たにどのような検証結果が出るかに注視していく必要があるだろう。