裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

2 個の時系列データの相関を考えるときは...(その2)

2014年09月12日 | 統計学

元データがない(?)ようなので,グラフから読み取って分析してみる。

> d = read.table("takahashi.dat", header=TRUE) # データは末尾に掲載
> plot(x ~ t, type="o", col="blue", pch=16, ylim=c(-1.5, 2), ylab="x", xlab="t")
> points(y ~ t, type="o", col="red", pch=16, yaxt="n", ylab="y")
> axis(4, at=seq(-1.5, 2, length=11), labels=seq(90, 140, by=5))
> mtext("y", 4, 1.8)


> (r_xy = cor(x, y)) # 単相関係数
[1] 0.8990529
グラフから読み取ったデータなので,ちょっと誤差がある

> (r_tx = cor(t, x)) # 時間との単相関係数
[1] 0.9478756

> (r_ty = cor(t, y)) # 時間との単相関係数
[1] 0.9603337

> (r_xy.t = (r_xy-r_tx*r_ty)/sqrt(1-r_tx^2)/sqrt(1-r_ty^2)) # 偏相関係数
[1] -0.1263188     ★1★ この結果を後で参照する
このデータでは,偏相関係数は弱いながらも負の相関ということになった

> n = length(t)
> (t.stat = abs(r_xy.t)*sqrt(n-3)/sqrt(1-r_xy.t^2)) # 母偏相関係数=0 の検定統計量
[1] 0.7745721
> (p.value = pt(t.stat, df=n-3, lower.tail=FALSE)*2) # P 値を求める(両側検定)
[1] 0.4435144
有意な偏相関係数ではないということである

なお,偏相関係数は,以下のような原理で求められる。
> resid.x = lm(x ~ t)$residuals # resid.x は,時間 t で予測できない部分(時間 t に影響されない部分)
> resid.y = lm(y ~ t)$residuals # resid.y は,時間 t で予測できない部分(時間 t に影響されない部分)
> cor(resid.x, resid.y) # resid.x と resid.y の単相関係数が,すなわち,時間 t の影響を除いた,正味の x, y の相関係数
[1] -0.1263188   ★2★ この値が,上の ★1★ で示したのと同じになる

ちなみに,resid.x で resid.y を予測するということは,時間の影響を取り除いて x で y を予測することになる
> summary(lm(resid.y ~ resid.x))

Call:
lm(formula = resid.y ~ resid.x)

Residuals:
     Min       1Q   Median       3Q      Max
-0.29687 -0.16963 -0.03915  0.13518  0.52168

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   0.0000     0.0314   0.000    1.000
resid.x      -0.1010     0.1286  -0.785    0.437

Residual standard error: 0.1986 on 38 degrees of freedom
Multiple R-squared:  0.01596,    Adjusted R-squared:  -0.009939
F-statistic: 0.6162 on 1 and 38 DF,  p-value: 0.4373

resid.x と resid.y を,平均値=0,分散=1に正規化して,予測式を作ると以下のようになる
> summary(lm(scale(resid.y) ~ scale(resid.x)))

Call:
lm(formula = scale(resid.y) ~ scale(resid.x))

Residuals:
    Min      1Q  Median      3Q     Max
-1.5025 -0.8585 -0.1982  0.6842  2.6403

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)
(Intercept)     7.687e-17  1.589e-01   0.000    1.000
scale(resid.x) -1.263e-01  1.609e-01  -0.785    0.437

Residual standard error: 1.005 on 38 degrees of freedom
Multiple R-squared:  0.01596,    Adjusted R-squared:  -0.009939
F-statistic: 0.6162 on 1 and 38 DF,  p-value: 0.4373

予測式の係数 scale(resid.x) が -1.263e-01 になっている。これは,★2★ の相関係数と同じ(すなわち★1★の偏相関係数と同じ)。これは,直線回帰の回帰係数と相関係数の関係からそうなる。

以上で~~すっ。

> d
    t      x      y
1   0 -1.282 -0.828
2   1 -1.273 -1.033
3   2 -1.333 -0.894
4   3 -1.308 -0.760
5   4 -1.339 -1.237
6   5 -1.020 -1.012
7   6 -0.949 -0.952
8   7 -1.096 -0.961
9   8 -1.055 -0.837
10  9 -0.920 -0.786
11 10 -0.914 -0.630
12 11 -1.001 -0.447
13 12 -0.897 -0.309
14 13 -0.915 -0.275
15 14 -0.957 -0.277
16 15 -0.887 -0.180
17 16 -0.859 -0.056
18 17 -0.531 -0.036
19 18 -0.562 -0.244
20 19 -0.718 -0.339
21 20 -0.027 -0.257
22 21  0.678 -0.245
23 22  0.235 -0.170
24 23  0.149 -0.151
25 24  0.159  0.010
26 25  0.181  0.084
27 26  0.181  0.320
28 27  0.291  0.507
29 28  0.487  0.684
30 29  0.464  0.610
31 30  0.519  0.580
32 31  0.067  0.586
33 32  0.072  0.548
34 33  0.684  0.620
35 34  0.420  0.654
36 35  0.614  0.683
37 36  0.824  0.848
38 37  0.717  1.112
39 38  0.910  1.309
40 39  1.161  1.581

コメント
  • Twitterでシェアする
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村