はじめに

こちらの書籍を読んでいます。

Rによるベイジアン動的線形モデル (統計ライブラリー)

作者:G.ペトリス,S.ペトローネ,P.カンパニョーリ
発売日: 2013/05/08
メディア: 単行本

その中で、動的因子分析について一部触れていたので、packageとRstanでの実行を行いたいと思います。

動的因子分析

m個の時系列がp個の共通因子によって変動すると考えることができる。

例えば、観測系列 $Y_t$ がp個の相関を持つランダムウォーク $\mu_t$ に線形従属していると、次のようなモデルとして書くことができる。

$Y_t = A \mu_t + v_t, v_t \sim N(0,V)$

$\mu_t = \mu_{t-1}+ w_t, w_t \sim N(0,W)$

ここで、 $A$ は因子負荷行列、 $\mu_t$ は時間と共に遷移する因子時系列要素である。

$V$ はm×1の観測系列の分散、 $W$ は因子時系列の分散を表現するp×pの対角行列である。

また、未知パラメータの推定を行うためには、 $A$ の $(i,j)$ 要素 $A_{(i,j)}$ が $j>i$ において0となるような制約をかける必要がある。

これは次のように書くこともできる。

$A = \begin{eqnarray} \left[ \begin{array}{1} T \\ B \end{array} \right] \end{eqnarray}$

ここで、 $T$ はp×pの下三角行列、 $B$ は(m-p)×pの矩形行列である。

MARSS packageを用いた場合

動的因子モデルは、MARSS packageで実行ができます。基本的にはこちらのコードを用いました。

用いたデータはMARSS packageに含まれている、ワシントン湖のプランクトン量のデータであるlakeWAplanktonデータセットをもちいます。

このデータセットから月ごとの5種類の植物プランクトンの計測値を用いました。

この時の値は、1mLあたりの細胞数を標準化したものです。

具体的なデータを見てみます。

data(lakeWAplankton, package = "MARSS")

all_dat <- lakeWAplanktonTrans 

dat_1980 <- 
  all_dat %>% 
  as.data.frame() %>% 
  filter(Year >= 1980, Year <= 1989) %>% 
  dplyr::select(Year, Month, Cryptomonas, Diatoms, Greens, Unicells, Other.algae) 

dat_1980 %>%
  rowid_to_column("index") %>% 
  pivot_longer(cols = -c("index", "Year", "Month"), names_to = "type", values_to = "value") %>% 
  ggplot() + 
  aes(x = index, y = value) +
  geom_point(aes(col = type), size = 0.5) + 
  geom_line(aes(col = type), size = 0.4) + 
  facet_grid(type ~ .) +
  theme_light() +
  ylab("Abundance index")+xlab("")+ 
  theme(legend.position = 'none')

f:id:saltcooky:20201128013807p:plain

モデルの推定のために色々変数等を準備します。

それらを定義は次のようにしていきます。

N_ts <- ncol(dat_1980)-2
TT <- nrow(dat_1980)

dat_1980 <- t(dat_1980)
y_bar <- apply(dat_1980, 1, mean, na.rm = TRUE)
dat <- dat_1980 - y_bar
rownames(dat) <- rownames(dat_1980)
dat <- dat[-c(1,2),]

# 因子負荷行列の定義
Z_vals <- list("z11", 0, 0, "z21", "z22", 0, "z31", "z32", "z33", 
               "z41", "z42", "z43", "z51", "z52", "z53")
ZZ <- matrix(Z_vals, nrow = N_ts, ncol = 3, byrow = TRUE)
ZZ

aa <- "zero"
DD <- "zero"  
dd <- "zero"  
RR <- "diagonal and unequal"

mm <- 3
BB <- "identity"  
uu <- "zero"  
CC <- "zero"  
cc <- "zero"  
QQ <- "identity" 

# モデルに与えるデータをまとめる
mod_list <- list(Z = ZZ, A = aa, D = DD, d = dd, R = RR, B = BB, 
                 U = uu, C = CC, c = cc, Q = QQ)

# 初期値の設定
init_list <- list(x0 = matrix(rep(0, mm), mm, 1))

# コントロールパラメータの設定
con_list <- list(maxit = 3000, allow.degen = TRUE)

このモデルを推定します。

# モデル推定
dfa_1 <- MARSS(y = dat,
               model = mod_list,
               inits = init_list, 
               control = con_list)

推定された因子時系列を見てみましょう。

# 因子負荷行列の抽出
Z_est <- coef(dfa_1, type = "matrix")$Z

# バリマックス回転の逆行列をえる
H_inv <- varimax(Z_est)$rotmat

# 回転後の因子行列を取得
Z_rot <-
  Z_est %*% H_inv %>% 
  as.data.frame() %>% 
  mutate(type = rownames(dat_1980) %>% tail(5)) %>% 
  rename(facter1=V1, facter2=V2, facter3=V3)

# 各因子時系列の取得
proc_rot <- 
  solve(H_inv) %*% dfa_1$states %>% 
  t() %>% 
  as.data.frame() %>% 
  rename(facter1=V1, facter2=V2, facter3=V3)

# 各因子時系列の描写
proc_rot %>% 
  pivot_longer(cols = everything(), names_to = "facter", values_to = "moment") %>% 
  rowid_to_column("time") %>% 
  ggplot+
  aes(x = time, y = moment) +
  geom_point(aes(col = facter), size = 0.5) + 
  geom_line(aes(col = facter), size = 0.4) + 
  facet_grid(facter ~ .)+
  theme_light()+
  ylab("")+xlab("")+ theme(legend.position = 'none')

f:id:saltcooky:20201128014712p:plain

次に因子負荷量を見てみます。

# 因子負荷量の描写
Z_rot %>% 
  pivot_longer(cols = -type, names_to = "facter", values_to = "moment") %>% 
  ggplot+
  aes(x = facter, y = moment, col = facter, fill=facter) +
  geom_bar(stat = "identity") + 
  facet_grid(type ~ .)+
  theme_light()+
  ylab("")+xlab("")+ theme(legend.position = 'none')

f:id:saltcooky:20201128014742p:plain

Rstanを用いた場合

次にRstanで実装してみます。

stanコードは次のような感じです。合っているのか不安です。間違っている気がします。

data {
  int T; // 時間ステップ数
  int M; // 観測変数の数
  int P; // 因子数
  matrix[M,T] Y; //観測データ
}

parameters {
  cholesky_factor_cov[P] D; //下三角行列
  matrix[M-P,P] B; // 短形行列
  
  matrix[P,T] mu; // 推定される因子時系列状態
  vector<lower = 0>[M] V; // 観測状態の分散 
}

transformed parameters {
  matrix[M,T] Yhat; // 推定される状態
  cov_matrix[P] W; // 因子状態の分散
  matrix[M,P] A; // 因子負荷行列
  
  for(m in 1:M){
      if(m<=P){
        A[m,] = D[m,];
      }
      else{
        A[m,] = B[m-P,];
      }
  }
  
  for(p1 in 1:P){
    for(p2 in 1:P){
      if(p1==p2){
        W[p1,p2] = 1;
      }else{
        W[p1,p2] = 0;
      }
    }
  }
  
  for(t in 1:T){
    for(m in 1:M){
      Yhat[m,t] = A[m,] * mu[,t];
    }
  }
}

model {

  for(t in 1:T) {
    for(p in 1:P) { 
      if(t > 2){
        mu[,t] ~ normal(2*mu[,t-1] -mu[,t-1], W[p,p]);
      }
      else{
        mu[,t] ~ normal(0, W[p,p]);
      }
    }
        
    for(m in 1:M){
      if(Y[m,t] != -999) {
        Y[m,t] ~ normal(Yhat[m,t], V[m]);
      }
    }
  }
  
  for(p in 1:P){
    W[p,p] ~ gamma(2, 2);
    for(m in 1:M){
      A[m,p] ~ normal(0,5);
    }
  }
  
  V ~ gamma(2, 2);
}

データの整形は次のように行いました。

cols_list <- c("Cryptomonas", "Diatoms", "Greens", "Unicells", "Other.algae")

dat_1980 <- all_dat %>% 
  as.data.frame() %>% 
  filter(Year >= 1980, Year <= 1989) %>% 
  dplyr::select(Year, Month, one_of(cols_list)) %>% 
  mutate(across(all_of(cols_list), ~if_else(is.na(.), -999, .))) # 未観測のポイントは-999に置換

run_df <-
  dat_1980 %>% 
  dplyr::select(-Year, -Month) %>% 
  t()

そして、モデルをキックしMCMCを走らせます。

dfm_model <- stan(file = "./model/DFA_Stan.stan", 
                  data = list(T = ncol(run_df), 
                              M = nrow(run_df),
                              P = 3,
                              Y = run_df), 
                  chains = 4,
                  warmup = 2000,
                  iter = 4000,
                  control=list(adapt_delta=0.99)
                  )