読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺　頭角四蹄都過了　因甚麼尾巴過不得

時系列データで相関を出してはいけないのか（失業と自殺は関係あるのか2）

R 時系列 仮説検定 乱数

経緯

男性に関して言えば、失業率と自殺率は強い相関を持つことが舞田敏彦らによって指摘されている。相関係数は0.7224。

データえっせい: 失業率が１％上がると，自殺者は何人増えるか？


一方で時系列データでは相関を出してはいけない、というような懐疑的な声も聞く。

ここでは自殺率と失業率の相関が「有意」といえるのかシミュレーションを用いて調べる。

結論としては、有意といえるとぼくは思う。

見せかけの回帰

時系列データで相関を出してはいけない、と言われている理由はおそらく「ランダムウォークする系列どうしで相関を出すと高い相関が出やすい」からだろう（違ったらごめん）。

ランダムウォークとは一歩前の値に足し算（引き算でもいい）していくことで作られた系列のこと。

これはまったく別々に作られた系列どうしで相関を算出しても高い値がでやすい。

1万回のシミュレーションの結果、2つのランダムウォーク系列の相関係数のヒストグラムはこんな感じだ。

f:id:abrahamcow:20160719085920p:plain

シミュレーションに用いた R のコードはこちら。

sim_rho0<-numeric(10000)
for(i in 1:10000){
  sim1 <- cumsum(rnorm(100))
  sim2 <- cumsum(rnorm(100))
  sim_rho0[i] <-cor(sim1,sim2)
}
hist(sim_rho0)

ランダムウォークする試行回数の二項分布

ただ、舞田の指摘した相関は失業「率」と、自殺「率」の相関だった。

日本の人口がランダムウォークを整数に直したもの、自殺者数と失業者数がそれを試行回数パラメータにした二項分布で決まるとすると、失業率、自殺率は二項分布の成功率パラメータの推定値である。

これはもちろんランダムウォークではない。

1万回のシミュレーションの結果、相関係数のヒストグラムはこんな感じで、正規分布に近いような形をしている。

f:id:abrahamcow:20160719092017p:plain

相関係数を求めるのは、十分意味があると思う。

シミュレーションに用いた R のコードはこちら。

sim_rho1<-numeric(10000)
for(i in 1:10000){
  sim_den <- round(cumsum(c(10000,rnorm(99,0,10))))
  sim1 <- rbinom(100,sim_den,0.01)
  sim2 <-  rbinom(100,sim_den,0.04)
  sim_rho1[i] <-cor(sim1/sim_den,sim2/sim_den)
}
hist(sim_rho1)

ロジットがランダムウォーク

難しいのは、失業率、自殺率をロジット変換したものがランダムウォークしている場合である。

1万回のシミュレーションの結果、相関係数のヒストグラムはこんな感じになる。

f:id:abrahamcow:20160719093243p:plain

単純に相関（だけ）を出すのは危険な感じがする。

しかし、一応62期分（1953～2014年）の系列を1万回シミュレーションした結果、0.7224以上の高い相関が出る割合は約0.04であり、5%水準で有意といえると思う。（追記：ちょっとまてよ。この場合相関係数は正規分布の平均と分散に依存するかもしれないな。後で調べる）

シミュレーションに用いた R のコードはこちら。

library(boot)
sim_rho2<-numeric(10000)
for(i in 1:10000){
  sim1 <- inv.logit(cumsum(rnorm(62)))
  sim2 <- inv.logit(cumsum(rnorm(62)))
  sim_rho2[i] <-cor(sim1,sim2)
}
hist(sim_rho2)
round(sum(sim_rho2 >= 0.7224)/10000,2)

無相関検定について

相関係数の検定（無相関検定）では、r を相関係数として、n をサンプルサイズとして、

t_0=\frac{r\sqrt{n-2}}{\sqrt{1-r^2}}

が、帰無仮説（相関係数=0）の下で自由度 n-2 の t 分布に従うことを利用しておこなう。

パターン1「ランダムウォークする試行回数の二項分布」の場合、t_0 は t 分布に従うようなので、ふつうに相関係数の検定をしてもよいだろう。

パターン1の場合について、t_0 のヒストグラムと t 分布を重ねて描いた図は以下。

f:id:abrahamcow:20160719112110p:plain

R のコードはこちら。

t1 <-sim_rho1*sqrt((100-2))/sqrt(1-sim_rho1^2)
hist(t1,freq = FALSE)
curve(dt(x,100-2),add = TRUE)

パターン2「ロジットがランダムウォーク」の場合、t_0 の分布は t 分布とはだいぶ異なる。

パターン2の場合について、t_0 のヒストグラムと t 分布を重ねて描いた図は以下。

f:id:abrahamcow:20160719112604p:plain

R のコードはこちら。

t2 <-sim_rho2*sqrt((62-2))/sqrt(1-sim_rho2^2)
hist(t2,freq = FALSE)
curve(dt(x,62-2),add = TRUE)

パターン2の場合はふつうの相関係数の検定はやってはいけない。

上でやったように「ロジットがランダムウォーク」するような系列をシミュレーションして相関係数の分布を自分で出してやる必要があるだろう。

まとめ

のふたつの場合について、相関係数の分布を検討したが、見落としていることがあったら教えてくれるとうれしいです。

関連エントリ

abrahamcow.hatenablog.com