理学療法士がまとめたノート

統計学備忘録(R言語のメモ)

since2016 ときどきTEXのメモ

分散共分散行列 相関のあるサンプル作成

データセットから分散共分散行列を求めてみます

x1 <- c(151, 164, 146, 158)
x2 <- c(48, 53, 45, 61)
x3 <- c(8, 11, 8, 9)
data <- data.frame(x1,x2,x3)

#分散共分散行列
var(data) 
         x1        x2        x3
x1 62.25000 38.250000 10.333333
x2 38.25000 48.916667  4.333333
x3 10.33333  4.333333  2.000000

#相関行列
cor(data) 
          x1        x2        x3
x1 1.0000000 0.6931597 0.9260955
x2 0.6931597 1.0000000 0.4381055
x3 0.9260955 0.4381055 1.0000000

定義に従って相関係数を求めてみます. 例)cor(data)[2]を算出してみます

cor(data)[2]
var(data)[2]/sqrt(var(data)[1]*var(data)[5])

シミュレーション

平均や分散を指定した乱数を作成するのは簡単です.例えば正規分布であればrnorm(10,5,3)で、平均5分散9の乱数が10個生成されます.
しかし繰り返しデータの場合に、それぞれの回数の相関を設定した乱数を作成するのはちょっと難しくなります.
5回繰り返しデータ(X1~X5)を想定します.また、それぞれの平均値が(0,3,5,-2,-1)、分散8、共分散3となるように乱数を発生させます.

#まず、次のような共分散行列を作成します
   x1 x2 x3 x4 x5
x1  8  3  3  3  3
x2  3  8  3  3  3
x3  3  3  8  3  3
x4  3  3  3  8  3
x5  3  3  3  3  8
  
sigma <- matrix(rep(6,25),ncol = 5)   # 3を25個作成
sigma  

     [,1] [,2] [,3] [,4] [,5]
[1,]    6    6    6    6    6
[2,]    6    6    6    6    6
[3,]    6    6    6    6    6
[4,]    6    6    6    6    6
[5,]    6    6    6    6    6

#次に分散8を挿入します
diag(sigma)=rep(8,5)       #斜めに挿入
sigma

     [,1] [,2] [,3] [,4] [,5]
[1,]    8    6    6    6    6
[2,]    6    8    6    6    6
[3,]    6    6    8    6    6
[4,]    6    6    6    8    6
[5,]    6    6    6    6    8

#X1~X5の平均値を設定します
m <- c(0, 3, 5, -2, -1) 

#パッケージをインストールします
install.packages("mvtnorm")
library(mvtnorm)

#平均m, 分散8, 共分散3になるようなX1~X5を10組の乱数を発生させます.
y <- rmvnorm(10, m, sigma, method = "chol")      #乱数
colnames(y) <- name
#イメージのために書き出しておきます.乱数なので算出毎に数値は変化します.
 y
               x1         x2         x3         x4         x5
 [1,]  3.52189310  6.3816131  8.5906099  0.3448154  2.4660227
 [2,]  0.92752113  4.9547725  7.9595828  0.3717669  1.3875238
 [3,] -5.19449966 -0.3019742 -1.4509597 -4.3760477 -5.0233539
 [4,] -1.01591408  4.2919216  5.4850598 -1.0672082 -1.4707111
 [5,] -3.80629633 -0.6012603  4.7040327 -6.4885712 -3.7804199
 [6,] -3.10433704  5.7128433  2.2531806 -1.1011307  0.8055750
 [7,]  2.05904069  4.8775650  6.2192196  0.3335801  0.4234229
 [8,] -0.62623689  4.2532797  3.7670864 -1.5564523 -1.8985524
 [9,] -0.09027824  3.3205490  5.0517354 -1.2558948 -1.0637946
[10,] -1.30173038 -3.0114726  0.8776159 -8.0027438 -3.9780232

#各列に適当に色をつけて、プロットしてくれるmatplotを使って、作図します.
#デフォルトでは col の順番なのでx1=黒 ,x2=赤, x3=緑, x4=青, x5=水色
matplot(y,type="l") 

f:id:yoshida931:20180423183232p:plain:w500

互いに相関のある群(相関係数0.75)を作成

#分散共分散行列よりそれぞれの相関係数は、6/√8*√8=0.75 が読み取れます
   x1 x2 x3 x4 x5
x1  8  3  3  3  3
x2  3  8  3  3  3
x3  3  3  8  3  3
x4  3  3  3  8  3
x5  3  3  3  3  8

しかしグラフからは、X1~X5の相関が理解できません.
ここからは、上記yをコピーしたy2を使用して進めていきます.
相関行列を確認してみたら6/√8*√8=0.75に近い相関になっています.

cor(y2)
          x1        x2        x3        x4        x5
x1 1.0000000 0.6014473 0.8170968 0.6413491 0.7693525
x2 0.6014473 1.0000000 0.6684236 0.9623296 0.8923660
x3 0.8170968 0.6684236 1.0000000 0.6635989 0.7835190
x4 0.6413491 0.9623296 0.6635989 1.0000000 0.8511389
x5 0.7693525 0.8923660 0.7835190 0.8511389 1.0000000

nのサイズを大きくすれば、相関係数はさらに6/√8*√8=0.75に近づきます

y3 <- rmvnorm(50, m, sigma, method = "chol")
cor(y3)
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 1.0000000 0.8068844 0.7566464 0.7541858 0.7390697
[2,] 0.8068844 1.0000000 0.7345947 0.6725027 0.6586652
[3,] 0.7566464 0.7345947 1.0000000 0.7810799 0.7455624
[4,] 0.7541858 0.6725027 0.7810799 1.0000000 0.7226452
[5,] 0.7390697 0.6586652 0.7455624 0.7226452 1.0000000

y3をグラフにしてみます

matplot(y3,type="l") 

f:id:yoshida931:20180423223431p:plain:w700

x軸をX1~X5に移行することで相関が可視化できます

matplot(t(y3),type="l")    

f:id:yoshida931:20180423223335p:plain:w700
x1~x5まで、ほぼ同じ範囲に入っており相関の強さが視覚的に確認できました

復習のためにy3を残しておきます

             [,1]         [,2]       [,3]       [,4]       [,5]
 [1,]  0.10151443  2.499704243  5.9057675 -3.4943595 -1.3900978
 [2,]  1.50663613  5.280020390  6.0871375 -1.6427891  4.1792270
 [3,]  5.15807547  9.240978804 10.1524971  1.6396755  4.2384627
 [4,]  0.52636740  7.362142992  8.2131410 -2.7272272 -0.1657287
 [5,] -1.82768487  2.078275786  4.8978490 -1.3363880 -2.6741024
 [6,]  4.43524347  8.060044313  8.8144178  2.4425139  2.0623879
 [7,]  0.37328244  5.317456919  3.7818383 -1.0436226 -0.7600355
 [8,] -4.91121977 -0.627691232  3.4766397 -4.2966815 -1.9693309
 [9,]  3.08701302  4.490715808  8.0062255  1.3102477  3.4494656
[10,] -1.07464046  4.124389139  4.1597881 -3.2331504 -1.5478108
[11,] -4.98187707 -2.960175387 -1.1852294 -8.4042762 -6.7175592
[12,]  0.55455335  4.811030651  4.4954359 -2.1924783 -1.6559525
[13,] -3.25343072  0.216598109  2.5139412 -4.0381462 -2.4717920
[14,]  0.71190379  0.664704194  7.4396890  1.2862593 -0.7421359
[15,]  2.05860252  2.432888164  3.5347403 -1.8087834 -1.5556714
[16,] -4.35657559  0.316841111  1.6688628 -5.3447686 -4.5653558
[17,]  0.52655249  0.326658121  3.3274462 -5.0033717 -1.4797982
[18,] -0.51997740  1.522057422  3.9754390 -2.2568758 -3.4186446
[19,] -2.48896458  0.609983346  4.4815961 -4.4651604  0.5507060
[20,] -1.27166463  0.860673117  4.5316817 -1.6695152 -1.5116476
[21,]  2.36579566  8.115692097 10.1842162  1.5469158  1.2745295
[22,] -2.61885567  2.264994539  1.8073273 -4.6059818 -5.2668700
[23,] -2.87195778 -0.354779514  4.4273497 -2.1824915 -2.3847656
[24,] -0.31689104  4.490830114  7.2127744  0.8737950  3.2108026
[25,]  0.86603870  3.574589801  4.4276743 -1.6189355 -0.6536883
[26,] -1.08641447  1.650209214  0.2495224 -3.6113148 -3.8762204
[27,] -1.06718456  1.333869467  2.7389823 -6.0151694 -1.0285024
[28,]  2.92878284  3.572365454  7.5983883  4.4834977  1.1370509
[29,]  1.25628266  3.077046228  4.4049753 -3.9389680  0.8752170
[30,] -2.25266309  1.031792708  2.0312006 -5.2208742 -1.1417857
[31,] -1.83386252  1.464612026  5.4984764 -2.5442451 -2.5926546
[32,] -3.74237017  0.948552374  1.6345401 -5.3518601 -4.3701054
[33,]  1.57641371  2.949825751  7.6793848  2.8385187  2.5930688
[34,] -1.70707918  4.328429267  4.8674954 -2.3957196 -0.6474211
[35,] -0.68729957  0.186676655  6.3561554 -4.9844379 -2.6635800
[36,] -5.70176744  0.002329385  3.3561491 -3.6218192 -2.5791208
[37,]  0.73014712  2.200125852  3.0212866 -1.4821494 -3.4750934
[38,]  2.51743524  4.256512695  3.2907951  1.0304869  3.3052088
[39,] -3.71975203 -0.761618136  3.5496253 -2.0532424 -1.8430078
[40,]  4.00292393  6.456935876  9.6564018  1.9699281  0.9976009
[41,] -0.76466675  5.343766540  4.4333739 -0.7760534 -1.2543057
[42,] -5.23139631 -2.144231142  0.3253171 -7.3690128 -5.2907999
[43,]  0.22047353  0.598735545  5.6964360 -3.8599177  0.1664410
[44,]  0.25324554  4.186434455  3.3374806 -2.2359334 -2.5759674
[45,] -0.33142225  3.259191383  4.9409355 -4.1641440 -3.3400309
[46,]  4.28859451  7.946808493  9.4197735  1.7908695  1.2385851
[47,]  0.91987986  4.440995369  4.2562333 -4.3308892 -1.1173378
[48,] -3.44169133  0.376280288  2.3650166 -3.1425327 -3.2670228
[49,]  0.07415149  2.969230453  5.1548999 -1.9221235  2.5954580
[50,]  2.51582740  4.186680476  6.5627895  1.7226084  3.3340985