目的

間接効果の検証
Multivariate S-mapにおける埋め込み次元が係数推定に与える影響その2(システムのデータの全てではなく一部を用いる場合)
間接効果の推定にtransmitterを含めるかどうかが与える影響

パッケージのロード

library(rEDM); packageVersion("rEDM") #"0.7.1"

## Warning: package 'rEDM' was built under R version 3.4.4

## [1] '0.7.1'

library(dplyr); packageVersion("dplyr") #"0.8.0.1"

## Warning: package 'dplyr' was built under R version 3.4.4

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## [1] '0.8.0.1'

3種food chainモデル

以下のResource, Consumer, Predator3種のfood chainモデルを考える
Rt+1 = Rt + bR×Rt×(1-Rt/K) - aC×Rt×Ct
Ct+1 = Ct + bC×Ct×Rt - aP×Ct×Pt - cC×Ct
Pt+1 = Pt + bP×Ct×Pt - c×Pt
aは捕食効率、bはbirth rate、cはdeath rate
＃Resourceはロジスティック成長を仮定、捕食被食関係は線形を仮定、また、PredatorだけでなくConsumerもひとりでに死ぬ
＃パラメータは系が安定になるように適当に設定

#Rt+1 = Rt + bR*Rt*(1-Rt/K) - aC*Rt*Ct
#Ct+1 = Ct + bC*Ct*Rt - aP*Ct*Pt - cC*Ct
#Pt+1 = Pt + bP*Ct*Pt - c*Pt

bR=3; K=60; aC=0.02
R.f <- function(v){
  v[1, 1] + bR*(1-v[1, 1]/K)*v[1, 1] - aC*v[1, 1]*v[1, 2]
}

bC=1.2; cC=0.7; aP=0.01
C.f <- function(v){
  v[1, 2] + bC*aC*v[1, 2]*v[1, 1] - aP*v[1, 2]*v[1, 3] - cC*v[1, 2]
}

bP=0.25; cP=0.1
P.f <- function(v){
  v[1, 3] + bP*aP*v[1, 3]*v[1, 2] - cP*v[1, 3]
}

関数定義

テスト時系列を作る関数test.ts
Multivariate S-mapをNxtとNytだけで行う関数mvar_smap2
Simplex projectionで埋め込み次元を決定する関数simp_emb
Multivariate S-mapをNxtの時間遅れ埋め込みとNytで行う関数mvar_smapE

##テスト時系列を作る関数
#増加率が正規分布でランダムにばらつく、これの大きさを操作できるようにする
#測定誤差は無いものとする

test.ts <- function(x) {
  g <- matrix(c(30, 30, 30),nrow=1,ncol=3)
  
  eR <- rnorm(300, 0, x)
  eC <- rnorm(300, 0, x)
  eP <- rnorm(300, 0, x)

  for (i in 1:299) {
    gn <- c(R.f(matrix(g[i, ], nrow=1))+eR[i+1], 
            C.f(matrix(g[i, ], nrow=1))+eC[i+1], 
            P.f(matrix(g[i, ], nrow=1))+eP[i+1])
    g <- rbind(g, gn)
  } # g:1列目にresource, 2列目にconsumer, 3列目にpredator, 300 時間ステップ分
  
  g <- cbind(g,rnorm(300,25,x)) #4列目に乱数の列であるzの動態を挿入
  dimnames(g) <- NULL
  g
}


##simplex projectionで埋め込み次元を決定する関数
simp_emb <- function(x){
  simp_x_res0 <- 
    simplex(x, E=1:10, lib=c(1, 50), pred=c(1, 50), silent=TRUE)
  simp_x_res0[which.min(simp_x_res0$rmse), "E"]
}

##Multivariate S-mapを任意の種xと種y2種で行う関数
mvar_smap2 <- function(v, x, y) {
  mlt_smap_res0 <- block_lnlp(cbind(v[,x], v[,y]), tp = 1, target_column = 1,
                              method = "s-map", theta = seq(1,10,by=0.1), silent = T)
  Tm <- mlt_smap_res0[which.min(mlt_smap_res0$rmse), "theta"]
  block_lnlp(cbind(v[,x], v[,y]), tp = 1, target_column = 1,
             method = "s-map", theta = Tm, save_smap_coefficients = T, silent = T)
}

##Multivariate S-mapを任意の種x,y,z3種で行う関数
mvar_smap3 <- function(v, x, y, z) {
  mlt_smap_res0 <- block_lnlp(cbind(v[,x], v[,y], v[,z]), tp = 1, target_column = 1,
                              method = "s-map", theta = seq(1,10,by=0.1), silent = T)
  Tm <- mlt_smap_res0[which.min(mlt_smap_res0$rmse), "theta"]
  block_lnlp(cbind(v[,x], v[,y], v[,z]), tp = 1, target_column = 1,
             method = "s-map", theta = Tm, save_smap_coefficients = T, silent = T)
}

##2変数multivariate S-mapをフル埋め込み(末尾入れ替え)で行う関数
mvar_smapE <- function(x, E){
  mvarts <- cbind(x[,1], x[,2])
  for (i in 1:(E-2)){
    mvarts <- cbind(mvarts, lag(x[,1], n=as.integer((1:(E-2))[i])))
  }
  mlt_smap_res0 <- block_lnlp(mvarts, tp = 1, target_column = 1,
                              method = "s-map", theta = seq(1,10,by=0.1), silent = T)
  Tm <- mlt_smap_res0[which.min(mlt_smap_res0$rmse), "theta"]
  block_lnlp(mvarts, tp = 1, target_column = 1,
             method = "s-map", theta = Tm, save_smap_coefficients = T, silent = T)
}

時系列の例

図示した時系列について、PredatorがResourceに与える間接効果を検証した。1.ResourceとPredatorの時系列のみを用いる方法、2.Resource埋め込み後にPredatorの時系列に置換する方法、3.全ての時系列を使う方法(fully multivariate)を比較した。
Consumer(間接効果のtransmitter)が入るか入らないかで結果が異なった。
Consumerがいない場合では捕食者が資源に与える間接効果は正として推定された。
一方、Consumerを含めた方法(ここではfully multivariate)では、ゼロに近い負の値が推定された。

set.seed(1)
g <- test.ts(1)
subset.g <- g[201:300, ] #最初の200ステップは捨てる

par(ps=14, lwd=1.5, family="sans")
plot(subset.g[, 1], type="l", ylim=c(10, 50), ylab="Abundance", xlab="Time", xaxt="n")
lines(subset.g[, 2], type="l", ylim=c(10, 50), col="blue")
lines(subset.g[, 3], type="l", ylim=c(10, 50), col="red")
lines(subset.g[, 4], type="l", ylim=c(10, 50), lty=3)
axis(1, labels=c(200, 220, 240, 260, 280, 300), at=c(0, 20, 40, 60, 80, 100))
legend("bottomleft", 
       legend=c("Resource", "Consumer", "Predator", "Gaussian noise"), 
       col=c("black", "blue", "red", "black"), lty=c(1, 1, 1, 3))

#この時系列について、PredatorがResourceに与える間接効果を検証したい
#そこで、1.ResourceとPredatorの時系列のみを用いる方法、
#2.Resource埋め込み後にPredatorの時系列に置換する方法、
#3.全ての時系列を使う方法(fully multivariate)を比較

mlt_smap_res1 <- mvar_smap2(subset.g, 1, 3) #方法1

E <- simp_emb(subset.g[,1]) #方法2
E

## [1] 7

mlt_smap_res2 <- mvar_smapE(subset.g[, c(1, 3)], E)

mlt_smap_res3 <- mvar_smap3(subset.g, 1, 2, 3) #方法3


par(mfrow=c(4, 1), mar=c(3, 4, 1, 1), ps=18, lwd=1.5, family="sans")

plot(subset.g[, 1], type="l", ylim=c(25, 50), ylab="Abundance", xlab="Time", xaxt="n", lty=3)
lines(subset.g[, 3], ylim=c(10, 60), col="red", lty=3)
lines(pred ~ time, mlt_smap_res1$model_output[[1]], lwd=1, col=4)
lines(pred ~ time, mlt_smap_res2$model_output[[1]], lwd=1, col=5)
lines(pred ~ time, mlt_smap_res3$model_output[[1]], lwd=1, col=6)
axis(1, labels=c(200, 220, 240, 260, 280, 300), at=c(0, 20, 40, 60, 80, 100))
legend("bottomleft", 
       legend=c("Resource", "Predator", "Predict 1", "Predict 2", "Predict 3"), 
       col=c(1, 2, 4, 5, 6), lty=c(3, 3, 1, 1, 1))

plot(mlt_smap_res1$smap_coefficients[[1]]$c_2, type="l", ylim=c(-1, 1), 
     ylab="Predator -> Resource", xlab="Time", xaxt="n")
axis(1, labels=c(200, 220, 240, 260, 280, 300), at=c(0, 20, 40, 60, 80, 100))
abline(h=0, lty=2)
legend("topleft", legend=c("Predict 1"), bty="n")

plot(mlt_smap_res2$smap_coefficients[[1]]$c_2, type="l", ylim=c(-1, 1), 
     ylab="Predator -> Resource", xlab="Time", xaxt="n")
axis(1, labels=c(200, 220, 240, 260, 280, 300), at=c(0, 20, 40, 60, 80, 100))
abline(h=0, lty=2)
legend("topleft", legend=c("Predict 2"), bty="n")

plot(mlt_smap_res3$smap_coefficients[[1]]$c_3, type="l", ylim=c(-1, 1), 
     ylab="Predator -> Resource", xlab="Time", xaxt="n")
axis(1, labels=c(200, 220, 240, 260, 280, 300), at=c(0, 20, 40, 60, 80, 100))
abline(h=0, lty=2)
legend("topleft", legend=c("Predict 3"), bty="n")

dev.off()

## null device 
##           1

時系列100個によるシミュレーション

Fully multivariateの予測精度が最も高かった。
Consumerが入っているfully multivariateでは、間接効果はほぼゼロ(やや負？？)として推定された。
(なお、Consumer⇒ResourceのS-map係数は負になった)
こころもち、時間遅れ時系列も含めた方がResourceとPredatorのみより予測精度が高い（しかしあまり大きな差ではない）
モデルの係数aC×aP=0.0002は、Rt+1 =の式中のPt-1にPt = Pt-1 + …の式変形を強引に代入して導いたものなので意味ないかもしれない

set.seed(1)
check.par <- NULL
for (i in 1:100){
  subset.g <- test.ts(1)[201:300, ]
  mlt_smap_res1 <- mvar_smap2(subset.g, 1, 3)
  rho_obs <- mlt_smap_res1$rho
  rmse_obs <- mlt_smap_res1$rmse
  IS_obs <- mean(mlt_smap_res1$smap_coefficients[[1]]$c_2, na.rm=TRUE)
  par_obs <- mean(mlt_smap_res1$smap_coefficients[[1]]$c_2*(1+bP*aP*lag(subset.g[,2])-cP)/
                    (subset.g[,1]*lag(subset.g[,2])), na.rm=TRUE)
  check.par <- rbind(check.par, data.frame(rho_obs, rmse_obs, IS_obs, par_obs))
}
check.par_2 <- check.par 

set.seed(1)
check.par <- NULL
for (i in 1:100){
  subset.g <- test.ts(1)[201:300, ]
  Ex <- simp_emb(subset.g)
  if(Ex<3){
  } else {
    mlt_smap_res1 <- mvar_smapE(subset.g[, c(1, 3)], Ex)
    rho_obs <- mlt_smap_res1$rho
    rmse_obs <- mlt_smap_res1$rmse
    IS_obs <- mean(mlt_smap_res1$smap_coefficients[[1]]$c_2, na.rm=TRUE)
    par_obs <- mean(mlt_smap_res1$smap_coefficients[[1]]$c_2*(1+bP*aP*lag(subset.g[,2])-cP)/
                      (subset.g[,1]*lag(subset.g[,2])), na.rm=TRUE)
    check.par <- rbind(check.par, data.frame(rho_obs, rmse_obs, IS_obs, par_obs))
  }
}
check.par_E <- check.par 

set.seed(1)
check.par <- NULL
for (i in 1:100){
  subset.g <- test.ts(1)[201:300, ]
  mlt_smap_res1 <- mvar_smap3(subset.g, 1, 2, 3)
  rho_obs <- mlt_smap_res1$rho
  rmse_obs <- mlt_smap_res1$rmse
  IS_obs <- mean(mlt_smap_res1$smap_coefficients[[1]]$c_3, na.rm=TRUE)
  par_obs <- mean(mlt_smap_res1$smap_coefficients[[1]]$c_3*(1+bP*aP*lag(subset.g[,2])-cP)/
                    (subset.g[,1]*lag(subset.g[,2])), na.rm=TRUE)
  check.par <- rbind(check.par, data.frame(rho_obs, rmse_obs, IS_obs, par_obs))
}
check.par_full <- check.par 


check.par <- rbind(check.par_2, check.par_E, check.par_full)
check.par$method <- as.factor(c(rep("2", times=nrow(check.par_2)), 
                                rep("E", times=nrow(check.par_E)), 
                                rep("Fully", times=nrow(check.par_full))))

par(mfrow=c(1, 4), ps=15, family="sans", lwd=1.5)
plot(rho_obs~method, check.par, ylab="rho")
plot(rmse_obs~method, check.par, ylab="RMSE")
plot(IS_obs~method, check.par, ylab="S-map coef (Predator -> Resource)")
abline(h=0, lty=3)
plot(par_obs~method, check.par, ylab="aC x aP (?)")
abline(h=0.0002, lty=3)

MultivariateSmap.test2

Koya Hashimoto

2020年7月15日

目的

パッケージのロード

3種food chainモデル

関数定義

時系列の例

時系列100個によるシミュレーション