13 Time-series analysis in R-INLA

本章では、INLAで扱える時系列相関を考慮したモデルについて解説する。時系列分析についてより詳しい解説が必要な場合は、他の書籍を参照されたし(e.g., 馬場 2018, 2019; Ravishanker et al. 2022)。

13.1 Simulation study

以下では、シミュレーションデータを用いて話を進める。\(Y_t\)を時系列データとし、以下のモデル式をもとに得られるとする。

\[ \begin{aligned} &Y_t = \rm{Intercept} + \rm{Covariates_i} + \rm{Trend_t} + \epsilon_t \\ &\epsilon_t \sim N(0, \sigma^2_{\epsilon}) \end{aligned} \]

ここでトレンド項(\(\rm{Trend_t}\))が\(Times \times \beta\)のように表せるとすれば、これは通常の線形回帰モデルである。しかしこの場合、時間的な疑似相関が生じてしまう。この問題を解決するため、いわゆるランダムウォークが用いられる。もっとも単純なランダムウォークモデルは以下のように書ける。10

\[ \begin{aligned} &Y_t = \rm{Intercept} + \rm{Covariates_i} + \mu_t + \epsilon_t \\ &\mu_t = \mu_{t-1} + v_t\\ &\epsilon_t \sim N(0, \sigma^2_{\epsilon}) \; and \; v_t \sim N(0, \sigma^2_v) \end{aligned} \]

このモデルでは、互いに独立な残差\(\epsilon_t\)と\(v_t\)がある。以下、これらについて理解するために\(\sigma_{\epsilon}\)と\(\sigma_{v}\)に様々な値を当てはめたときにシミュレートされる\(Y_t\)を図示したのが図13.1である。1行目の3つは\(\sigma_{\epsilon} = 1,\sigma_v = 1\)\(、2行目の3つは\)\(\sigma_{\epsilon} = 1,\sigma_v = 0.1\)\(、3行目の3つは\)\(\sigma_{\epsilon} = 5, \sigma_v = 1\)\(、4行目の3つは\)\(\sigma_{\epsilon} = 1, \sigma_v = 5\)\(、5行目の3つは\)\(\sigma_{\epsilon} = 3, \sigma_v = 3\)である。なお、線は\(\mu_t\)を、点は\(Y_t\)を表す。

sigma_e <- rep(c(1,1,5,1,3), each = 3)
sigma_v <- rep(c(1, 0.1, 1,5,3), each = 3)

set.seed(123)
tmax <- 100
Intercept <-  3.2
y <- matrix(ncol = 15, nrow = tmax)
mu <- matrix(ncol = 15, nrow = tmax)
mu[1, 1:15] <- rnorm(n = 15, mean = 0, sd = 1)
y[1, 1:15] <- Intercept + mu[1,1:15] + rnorm(n = 15, mean = 0, sd = sigma_e)

for(j in seq_along(sigma_e)){
  for(i in 2:100){
     mu[i,j] <- mu[i-1,j] + rnorm(n = 1, 0, sigma_v[j])
     y[i, j] <- Intercept + mu[i,j] + rnorm(n = 1, 0, sigma_e[j]) 
  }
}

y %>% 
  data.frame() %>% 
  mutate(time = 1:n()) %>% 
  pivot_longer(1:15) %>% 
  mutate(name = as.numeric(str_replace_all(name, "X",""))) %>% 
  ggplot(aes(x = time, y = value))+
  geom_point(size = 0.5)+
  geom_line(data = mu %>% 
              data.frame() %>% 
              mutate(time = 1:n()) %>% 
              pivot_longer(1:15) %>% 
              mutate(name = as.numeric(str_replace_all(name, "X",""))),
            aes(y = value + Intercept),
            linewidth = 0.8)+
  facet_rep_wrap(~name, scales = "free",
                 repeat.tick.labels = TRUE,
                 ncol = 3)+
  theme_bw()+
  theme(aspect.ratio = 0.7)+
  labs(y = "Simulated data")
Examples of simulated data sets with different random walk trends. Each of the three panels on the same row represents a different simulation from the same model. The number above a panel is the simulation number. In simulations 1–3 (top row) we used σ ε = 1 and σ v = 1; in simulations 4–6 (second row) we used σ ε = 1 and σ v = 0.1; in simulations 7–9 we used σ ε = 5 and σ v = 1; in simulations 10–12 we used σ ε = 1 and σ v = 5; and in simulations 13–15 we used σ ε = 3 and σ v = 3.

図13.1: Examples of simulated data sets with different random walk trends. Each of the three panels on the same row represents a different simulation from the same model. The number above a panel is the simulation number. In simulations 1–3 (top row) we used σ ε = 1 and σ v = 1; in simulations 4–6 (second row) we used σ ε = 1 and σ v = 0.1; in simulations 7–9 we used σ ε = 5 and σ v = 1; in simulations 10–12 we used σ ε = 1 and σ v = 5; and in simulations 13–15 we used σ ε = 3 and σ v = 3.


この図からは、\(\sigma_v\)が小さいとトレンド項(\(\mu_t\))の変動も小さく、なめらかであることが分かる(e.g., 2行目)。反対に、\(\sigma_v\)が大きいと\(\mu_t\)の変動が大きくなる。一方で、\(\sigma_{\epsilon}\)は大きいほど\(Y_t\)がトレンド項\(\mu_t\)の周りで大きくばらつく。ランダムウォークについてさらに詳しい解説は Durbin and Koopman (2012) や 馬場 (2018) 、馬場 (2019) などを参照。

13.5 Multivariate time series for Hawaiian birds

13.5.1 Importing and preparing the data

本節では、 Zuur et al. (2007) と Reed et al. (2011) で分析されている、ハワイに生息する3種の固有種を1年に2回測定したデータを用いる。

データは以下の通り。Birdsは確認された羽数を、Yearは観測年を、Speciesは種のID(1: stilts, 2: coots, 3: moorhen)を、Islandはデータを取得した島(1: オアフ, 2: マウイ, 3: カウアイとニイハウ)を示す。SpeciesとIslandは実際の名前を記した変数を作成する。

hb <- read_delim("data/HawaiiBirdsV2.txt")

speciesname <- c("Stilts", "Coots", "Moorhens")
islandname <- c("Oahu", "Maui", "Kauai and Niihau")

hb2 <- hb %>% 
  mutate(fSpecies = speciesname[hb$Species],
         fIsland = islandname[hb$Island]) %>% 
  mutate(species_island = str_c(fSpecies, ".",fIsland)) %>% 
  mutate(fSpecies = fct_relevel(fSpecies, "Stilts", "Moorhens"),
         fIsland = fct_relevel(fIsland, "Maui", "Oahu")) 

datatable(hb2,
          options = list(scrollX = 30),
          filter = "top")

13.5.2 Data exploration

各島における各種の羽数の推移を示したのが図13.5である。島ごとにトレンドは異なっていそうである。一方で、オアフ島のMoorhensとCootsは似たトレンドを持つ。

hb2 %>% 
  ggplot(aes(x = Year, y = Birds))+
  geom_point()+
  geom_line()+
  facet_rep_grid(fSpecies ~ fIsland,
                 repeat.tick.labels = TRUE,
                 scales = "free")+
  theme_bw()+
  labs(y = "Bird abundance")
Time-series plot of the bird abundance at three islands and three species.

図13.5: Time-series plot of the bird abundance at three islands and three species.


13.5.3 Model formulation

ここでは、負の二項分布モデルを用いる。考えなければいけないのは、トレンド項をどのように扱うかである。以下の4つのパターンが考えられる。

  1. 全ての島のすべての種が同じトレンドを持つ
  2. 島ごとに異なるトレンドを持つが、種によるトレンドの違いはない
  3. 種ごとに異なるトレンドを持つが、島によるトレンドの違いはない
  4. 各島の各種ごとに異なるトレンドを持つ

いずれのモデルでもモデルの以下の部分は変わらない。なお、\(t\)は年を、\(i\)は種を、\(j\)は島を表す添え字である。

\[ \begin{aligned} &Birds_{tij} \sim NB(\mu_{tij}, \theta)\\ &E(Birds_{tij}) = \mu_{tij} \\ &var(Birds_{tij}) = \mu_{tij} + \mu_{tij}^2/\theta \end{aligned} \]

一方で、\(\mu_{tij}\)をどうモデリングするかが4つの場合で異なる。\(a_{tij}\)は8つの時系列ごとに異なる切片である。

\[ \begin{aligned} log(\mu_{tij})= \begin{cases} a_{tij} + Trend_t \;\; &(1)\\ a_{tij} + Trend_{tj} \;\; &(2)\\ a_{tij} + Trend_{ti} \;\; &(3)\\ a_{tij} + Trend_{tij} \;\; &(4) \end{cases} \end{aligned} \]

13.5.4 Executing the models

13.5.4.1 One trend for all series

それでは、まず1の場合のモデルを実行する。データから応答変数が欠損値である行を削除する。

hb3 <- hb2 %>% 
  drop_na(Birds)

ハイパーパラメータの事前分布には、\(U=1\)のPC事前分布を用いる。

hyper.prec <- list(theta = list(prior = "pc.prec",
                                param = c(1, 0.01)))

f14_9a <- Birds ~ species_island + f(Year, model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec)

m14_9a <- inla(f14_9a,
               control.compute = list(dic = TRUE),
               family = "nbinomial",
               data = hb3)

モデルの事後平均から算出したピアソン残差は以下のように求められる。

mu <- m14_9a$summary.fitted.values$mean
theta <- m14_9a$summary.hyperpar[1,1]
E1 <- (hb3$Birds - mu)/sqrt(mu + mu^2/theta)

モデルから求められた\(\mu_{tij}\)の事後平均とその95%確信区間を図示すると図13.6のようになる。これらの曲線は、切片以外はすべて同じである。

m14_9a$summary.fitted.values %>% 
  bind_cols(hb3) %>% 
  ggplot(aes(x = Year, y = mean))+
  geom_line()+
  geom_ribbon(aes(ymax = `0.025quant`, ymin = `0.975quant`),
              alpha = 0.2)+
  geom_point(aes(y = Birds))+
  theme_bw()+
  theme(aspect.ratio = 0.8)+
  facet_rep_grid(fSpecies~fIsland, repeat.tick.labels = TRUE,
                 scales = "free")+
  labs(y = "Bird abundance")
Fitted values of the model m14_9a

図13.6: Fitted values of the model m14_9a


モデルの事後平均から算出したピアソン残差を用いて分散パラメータを算出すると、やや過少分散？

mu_a <- m14_9a$summary.fitted.values$mean
theta_a <- m14_9a$summary.hyperpar[1,1]
E1_a <- (hb3$Birds - mu_a)/sqrt(mu_a + mu^2/theta_a)
p_a <- nrow(m14_9a$summary.fixed)
N <- nrow(hb3)

phi_a <- sum(E1^2)/(N - p_a)

phi_a
## [1] 0.8425054

最後に、残差の自己相関を調べる。

hb3 %>% 
  mutate(resid = E1_a) %>% 
  group_by(species_island) %>% 
  arrange(Year, .by_group = TRUE) -> hb3b

species_island <- unique(hb3b$species_island)
all.out_a <- NULL

for(i in seq_along(species_island)){
  data <- hb3b %>% filter(species_island == species_island[i])
  
  out.acf <- acf(data$resid,
                 lag.max = 15,
                 plot = FALSE)
  
  out.df <- data.frame(Timelag = out.acf$lag,
                       Acf = out.acf$acf,
                       SE = qnorm(0.975)/sqrt(out.acf$n.used),
                       ID = species_island[i])
  
  all.out_a <- bind_rows(all.out_a, out.df)
}

自己相関はまだかなり残っている。

all.out_a %>% 
  ggplot(aes(x = Timelag, y = 0))+
  geom_segment(aes(xend = Timelag, yend = Acf))+
  geom_ribbon(aes(ymax = SE, ymin = -SE),
              alpha = 0.3)+
  theme_bw()+
  theme(aspect.ratio = 0.8)+
  facet_rep_wrap(~ID, repeat.tick.labels = TRUE)+
  labs(y = "Auto-correlation")

13.5.4.2 One trend for each island

続いて、2のモデルを実行する。このモデルでは、各島ごとに異なるトレンドがあると仮定する。

このようなモデルを実行するとき、まず複数のランダムウォークの事前分布を標準化する必要がある。

inla.setOption(scale.model.default = TRUE)

また、新たな変数として、各島を表すダミー変数を作成する。また、Yearを3つのモデルに使用するので観察年を含むYear2、Year3という変数も作成する。

hb3 %>% 
  mutate(Oahu = as.numeric(fIsland == "Oahu"),
         Maui = as.numeric(fIsland == "Maui"),
         KN = as.numeric(fIsland == "Kauai and Niihau")) %>% 
  mutate(Year2 = Year,
         Year3 = Year)-> hb4

モデルは以下のように実行できる。

f14_9b <- Birds ~ species_island + f(Year, Oahu,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year2, Maui,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year3, KN,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec)
  

m14_9b<- inla(f14_9b,
               control.compute = list(dic = TRUE),
               family = "nbinomial",
               data = hb4)

3つのランダムウォークごとに\(\sigma_v\)が推定されている。

m14_9b$summary.hyperpar

モデルから求められた\(\mu_{tij}\)の事後平均とその95%確信区間を図示すると図13.7のようになる。

m14_9b$summary.fitted.values %>% 
  bind_cols(hb4) %>% 
  ggplot(aes(x = Year, y = mean))+
  geom_line()+
  geom_ribbon(aes(ymax = `0.025quant`, ymin = `0.975quant`),
              alpha = 0.2)+
  geom_point(aes(y = Birds))+
  theme_bw()+
  theme(aspect.ratio = 0.8)+
  facet_rep_grid(fSpecies~fIsland, repeat.tick.labels = TRUE,
                 scales = "free")+
  labs(y = "Bird abundance")
Fitted values of the model m14_9b

図13.7: Fitted values of the model m14_9b


先ほどのモデルと比較すると、トレンドを島ごとにした方がはるかに良い。

c(m14_9a$dic$dic, m14_9b$dic$dic)
## [1] 3822.675 3783.597

モデルの事後平均から算出したピアソン残差を用いて分散パラメータを算出すると、やや過少分散？

mu_b <- m14_9b$summary.fitted.values$mean
theta_b <- m14_9b$summary.hyperpar[1,1]
E1_b <- (hb3$Birds - mu_b)/sqrt(mu_b + mu_b^2/theta_b)
p_b <- nrow(m14_9b$summary.fixed)
N <- nrow(hb4)

phi_b <- sum(E1_b^2)/(N - p_b)

phi_b
## [1] 0.7737521

最後に、残差の自己相関を調べる。

hb4 %>% 
  mutate(resid = E1_b) %>% 
  group_by(species_island) %>% 
  arrange(Year, .by_group = TRUE) -> hb4b

species_island <- unique(hb4b$species_island)
all.out_b <- NULL

for(i in seq_along(species_island)){
  data <- hb4b %>% filter(species_island == species_island[i])
  
  out.acf <- acf(data$resid,
                 lag.max = 15,
                 plot = FALSE)
  
  out.df <- data.frame(Timelag = out.acf$lag,
                       Acf = out.acf$acf,
                       SE = qnorm(0.975)/sqrt(out.acf$n.used),
                       ID = species_island[i])
  
  all.out_b <- bind_rows(all.out_b, out.df)
}

自己相関はまだかなり残っているように見える。

all.out_b %>% 
  ggplot(aes(x = Timelag, y = 0))+
  geom_segment(aes(xend = Timelag, yend = Acf))+
  geom_ribbon(aes(ymax = SE, ymin = -SE),
              alpha = 0.3)+
  theme_bw()+
  theme(aspect.ratio = 0.8)+
  facet_rep_wrap(~ID, repeat.tick.labels = TRUE)+
  labs(y = "Auto-correlation")

13.5.4.3 One trend for each species

続いて、3のモデルを実行する。このモデルでは、各種ごとに異なるトレンドがあると仮定する。

新たな変数として、各島を表すダミー変数を作成する。また、Yearを3つのモデルに使用するので観察年を含むYear2、Year3という変数も作成する。

hb3 %>% 
  mutate(Stilts = as.numeric(fSpecies == "Stilts"),
         Coots = as.numeric(fSpecies == "Coots"),
         Moorhens = as.numeric(fSpecies == "Moorhens")) %>% 
  mutate(Year2 = Year,
         Year3 = Year)-> hb5

モデルは以下のように実行できる。

f14_9c <- Birds ~ species_island + f(Year, Stilts,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year2, Coots,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year3, Moorhens,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec)
  

m14_9c <- inla(f14_9c,
               control.compute = list(dic = TRUE),
               family = "nbinomial",
               data = hb5)

モデルから求められた\(\mu_{tij}\)の事後平均とその95%確信区間を図示すると図13.8のようになる。

m14_9c$summary.fitted.values %>% 
  bind_cols(hb4) %>% 
  ggplot(aes(x = Year, y = mean))+
  geom_line()+
  geom_ribbon(aes(ymax = `0.025quant`, ymin = `0.975quant`),
              alpha = 0.2)+
  geom_point(aes(y = Birds))+
  theme_bw()+
  theme(aspect.ratio = 0.8)+
  facet_rep_grid(fSpecies~fIsland, repeat.tick.labels = TRUE,
                 scales = "free")+
  labs(y = "Bird abundance")
Fitted values of the model m14_9c

図13.8: Fitted values of the model m14_9c


モデルの事後平均から算出したピアソン残差を用いて分散パラメータを算出すると、やや過少分散？

mu_c <- m14_9c$summary.fitted.values$mean
theta_c <- m14_9c$summary.hyperpar[1,1]
E1_c <- (hb3$Birds - mu_c)/sqrt(mu_c + mu_c^2/theta_c)
p_c <- nrow(m14_9c$summary.fixed)
N <- nrow(hb5)

phi_c <- sum(E1_c^2)/(N - p_c)

phi_c
## [1] 0.7466784

最後に、残差の自己相関を調べる。

hb5 %>% 
  mutate(resid = E1_c) %>% 
  group_by(species_island) %>% 
  arrange(Year, .by_group = TRUE) -> hb5b

species_island <- unique(hb5b$species_island)
all.out_c <- NULL

for(i in seq_along(species_island)){
  data <- hb5b %>% filter(species_island == species_island[i])
  
  out.acf <- acf(data$resid,
                 lag.max = 15,
                 plot = FALSE)
  
  out.df <- data.frame(Timelag = out.acf$lag,
                       Acf = out.acf$acf,
                       SE = qnorm(0.975)/sqrt(out.acf$n.used),
                       ID = species_island[i])
  
  all.out_c <- bind_rows(all.out_b, out.df)
}

自己相関はまだかなり残っているように見える。

all.out_c %>% 
  ggplot(aes(x = Timelag, y = 0))+
  geom_segment(aes(xend = Timelag, yend = Acf))+
  geom_ribbon(aes(ymax = SE, ymin = -SE),
              alpha = 0.3)+
  theme_bw()+
  theme(aspect.ratio = 0.8)+
  facet_rep_wrap(~ID, repeat.tick.labels = TRUE)+
  labs(y = "Auto-correlation")

13.5.4.4 One trend for each species and island combination

最後に、4のモデルを実行する。このモデルでは、各島の各種ごとに異なるトレンドがあると仮定する。

新たな変数として、島と種の組み合わせを表すダミー変数を作成する。また、Yearを8つのモデルに使用するので観察年を含むYear2からYear8という変数も作成する。

hb3 %>% 
  mutate(Stilts_Oahu = as.numeric(species_island == "Stilts.Oahu"),
         Stilts_Maui = as.numeric(species_island == "Stilts.Maui"),
         Stilts_KN = as.numeric(species_island == "Stilts.Kauai and Niihau"),
         Coots_Oahu = as.numeric(species_island == "Coots.Oahu"),
         Coots_Maui = as.numeric(species_island == "Coots.Maui"),
         Coots_KN = as.numeric(species_island == "Coots.Kauai and Niihau"),
         Moorhens_Oahu = as.numeric(species_island == "Moorhens.Oahu"),
         Moorhens_Maui = as.numeric(species_island == "Moorhens.Maui")) %>% 
  mutate(Year2 = Year,
         Year3 = Year,
         Year4 = Year,
         Year5 = Year,
         Year6 = Year,
         Year7 = Year,
         Year8 = Year) -> hb6

モデルは以下のように実行できる。

f14_9d <- Birds ~ species_island + f(Year, Stilts_Oahu,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year2, Stilts_Maui,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year3, Stilts_KN,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year4, Coots_Oahu,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year5, Coots_Maui,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year6, Coots_KN,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year7, Moorhens_Oahu,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year8, Moorhens_Maui,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) 
  

m14_9d <- inla(f14_9d,
               control.compute = list(dic = TRUE),
               family = "nbinomial",
               data = hb6)

モデルから求められた\(\mu_{tij}\)の事後平均とその95%確信区間を図示すると図13.9のようになる。より実データに近い予測値になっていることが分かる。

m14_9d$summary.fitted.values %>% 
  bind_cols(hb4) %>% 
  ggplot(aes(x = Year, y = mean))+
  geom_line()+
  geom_ribbon(aes(ymax = `0.025quant`, ymin = `0.975quant`),
              alpha = 0.2)+
  geom_point(aes(y = Birds))+
  theme_bw()+
  theme(aspect.ratio = 0.8)+
  facet_rep_grid(fSpecies~fIsland, repeat.tick.labels = TRUE,
                 scales = "free")+
  labs(y = "Bird abundance")
Fitted values of the model m14_9d

図13.9: Fitted values of the model m14_9d


モデルの事後平均から算出したピアソン残差を用いて分散パラメータを算出すると、かなり過少分散？

mu_d <- m14_9d$summary.fitted.values$mean
theta_d <- m14_9d$summary.hyperpar[1,1]
E1_d <- (hb6$Birds - mu_d)/sqrt(mu_d + mu_d^2/theta_d)
p_d <- nrow(m14_9d$summary.fixed)
N <- nrow(hb6)

phi_d <- sum(E1_d^2)/(N - p_d)

phi_d
## [1] 0.6304094

最後に、残差の自己相関を調べる。

hb6 %>% 
  mutate(resid = E1_d) %>% 
  group_by(species_island) %>% 
  arrange(Year, .by_group = TRUE) -> hb6b

species_island <- unique(hb6b$species_island)
all.out_d <- NULL

for(i in seq_along(species_island)){
  data <- hb6b %>% filter(species_island == species_island[i])
  
  out.acf <- acf(data$resid,
                 lag.max = 15,
                 plot = FALSE)
  
  out.df <- data.frame(Timelag = out.acf$lag,
                       Acf = out.acf$acf,
                       SE = qnorm(0.975)/sqrt(out.acf$n.used),
                       ID = species_island[i])
  
  all.out_d <- bind_rows(all.out_d, out.df)
}

自己相関はほとんど消えている。

all.out_d %>% 
  ggplot(aes(x = Timelag, y = 0))+
  geom_segment(aes(xend = Timelag, yend = Acf))+
  geom_ribbon(aes(ymax = SE, ymin = -SE),
              alpha = 0.3)+
  theme_bw()+
  theme(aspect.ratio = 0.8)+
  facet_rep_wrap(~ID, repeat.tick.labels = TRUE)+
  labs(y = "Auto-correlation")

13.5.5 Mixing Poisson and negative binomial distributions

データを見てみると、時系列によってばらつきの大きさに差があることが分かる。このことから、時系列ごとにポワソン分布を適用したり負の二項分布を適用したらよりデータに合った分析ができ、先ほどのモデルで見られた過少分散も解決できる可能性があることが分かる。INLAではこれを行うことができる。

この分析を行うためには、以下のように各列に各時系列のデータが入った以下のような行列を作成する必要がある。

Y <- matrix(NA, nrow = nrow(hb6), ncol = 8)

hb6 %>% 
  mutate(species_island = fct_relevel(species_island, unique(hb6$species_island))) -> hb6
  
hb6 %>% 
  select(species_island, Birds, Year) %>% 
  group_by(species_island) %>% 
  arrange(Year, .by_group = TRUE) %>% 
  ungroup() %>% 
  pivot_wider(names_from = species_island,
              values_from = Birds) %>% 
  select(-Year) -> hb7

hb6 %>% 
  group_by(species_island) %>% 
  summarise(N = n()) %>% 
  ungroup() %>% 
  mutate(cum = cumsum(N))-> n_cat

for(i in seq_along(species_island)){
  if(i == 1){
    Y[1:n_cat[[3]][i], i] <- na.omit(hb7[[i]])
  }else{
    Y[(n_cat[[3]][i-1]+1):n_cat[[3]][i], i] <- na.omit(hb7[[i]])
  }
}

Y %>% 
  data.frame()


そして、この行列に加えて観察年などを含むリストを作成する。

data_list <- list(Birds = Y,
                  Stilts_Oahu = hb6$Stilts_Oahu,
                  Stilts_Maui = hb6$Stilts_Maui,
                  Stilts_KN = hb6$Stilts_KN,
                  Coots_Oahu = hb6$Coots_Oahu,
                  Coots_Maui = hb6$Coots_Maui,
                  Coots_KN = hb6$Coots_KN,
                  Moorhens_Oahu = hb6$Moorhens_Oahu,
                  Moorhens_Maui = hb6$Moorhens_Maui,
                  species_island = hb6$species_island,
                  Year = hb6$Year,
                  Year2 = hb6$Year,
                  Year3 = hb6$Year,
                  Year4 = hb6$Year,
                  Year5 = hb6$Year,
                  Year6 = hb6$Year,
                  Year7 = hb6$Year,
                  Year8 = hb6$Year)

モデルは、以下のように実行できる。

f14_9e <- Birds ~ species_island + f(Year, Stilts_Oahu,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year2, Stilts_Maui,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year3, Stilts_KN,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year4, Coots_Oahu,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year5, Coots_Maui,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year6, Coots_KN,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year7, Moorhens_Oahu,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) +
                                   f(Year8, Moorhens_Maui,
                                     model = "rw1",
                                     scale.model = TRUE,
                                     hyper = hyper.prec) 
  

m14_9e <- inla(f14_9e,
               control.compute = list(dic = TRUE),
               family = c("nbinomial", "nbinomial", "nbinomial", "poisson", "poisson", "nbinomial", "nbinomial","nbinomial"),
               data = data_list)

モデルの事後平均から算出したピアソン残差を用いて分散パラメータを算出すると、概ね過少分散が解消されている。

mu_e <- m14_9e$summary.fitted.values$mean
theta_e <- m14_9e$summary.hyperpar[1,1]
E1_e <- (hb6$Birds - mu_e)/sqrt(mu_e + mu_e^2/theta_e)
p_e <- nrow(m14_9e$summary.fixed)
N <- nrow(hb6)

phi_e <- sum(E1_e^2)/(N - p_e)

phi_e
## [1] 1.302552

モデルから求められた\(\mu_{tij}\)の事後平均とその95%確信区間を図示すると図13.10のようになる。より実データに近い予測値になっていることが分かる。

m14_9e$summary.fitted.values %>% 
  bind_cols(hb4) %>% 
  ggplot(aes(x = Year, y = mean))+
  geom_line()+
  geom_ribbon(aes(ymax = `0.025quant`, ymin = `0.975quant`),
              alpha = 0.2)+
  geom_point(aes(y = Birds))+
  theme_bw()+
  theme(aspect.ratio = 0.8)+
  facet_rep_grid(fSpecies~fIsland, repeat.tick.labels = TRUE,
                 scales = "free")+
  labs(y = "Bird abundance")
Fitted values of the model m14_9d

図13.10: Fitted values of the model m14_9d


DICも先ほどより改善している。

c(m14_9d$dic$dic, m14_9e$dic$dic)
## [1] 3613.067 3298.294

References

Blangiardo M, Cameletti M (2015) Spatial and spatio-temporal bayesian models with R - INLA. John Wiley & Sons
Carroll R, Lawson AB, Faes C, et al (2015) Comparing INLA and OpenBUGS for hierarchical poisson modeling in disease mapping. Spat Spatiotemporal Epidemiol 14-15:45–54
Crozier LG, Scheuerell MD, Zabel RW (2011) Using time series analysis to characterize evolutionary and plastic responses to environmental change: A case study of a shift toward earlier migration date in sockeye salmon. Am Nat 178:755–773
Durbin J, Koopman SJ (2012) Time series analysis by state space methods. OUP Oxford
Etheridge DM, Steele LP, Francey RJ, et al (1998) Atmospheric methane between 1000 AD and present: Evidence of anthropogenic emissions and climatic variability. Journal of
Mauritzen M, Belikov SE, Boltunov AN, et al (2003) Functional responses in polar bear habitat selection. Oikos 100:112–124
Pierce GJ, Santos MB, Smeenk C, et al (2007) Historical trends in the incidence of strandings of sperm whales (physeter macrocephalus) on north sea coasts: An association with positive temperature anomalies. Fish Res 87:219–228
Ravishanker N, Raman B, Soyer R (2022) Dynamic time series models using R-INLA: An applied perspective. https://ramanbala.github.io/dynamic-time-series-models-R-INLA/; CRC Press
Reed JM, Elphick CS, Ieno EN, Zuur AF (2011) Long-term population trends of endangered hawaiian waterbirds. Popul Ecol 53:473–481
Simpson DP, Rue H, Martins TG, et al (2014) Penalising model component complexity: A principled, practical approach to constructing priors. Stat Sci 32:1–28
Smeenk C (1997) Strandings of sperm whales physeter macrocephalus in the north sea: History and patterns
Zuur AF, Ieno EN, Smith GM (2007) Analyzing ecological data. Springer New York
馬場真哉 (2019) R と stan ではじめるベイズ統計モデリングによるデータ分析入門. 講談社
馬場真哉 (2018) 時系列分析と状態空間モデルの基礎 : RとStanで学ぶ理論と実装. プレアデス出版

  1. 基本的には、 馬場 (2019) でローカルレベルモデルといわれているモデルである。↩︎

  2. 馬場 (2019) で平滑化トレンドレベルモデルといわれているモデルである。↩︎