【R】行ごとに引き算を行って列ごとの分散を計算する

自分でBootstrapを行い、各パラメーターの標準誤差を計算するときに、少し行列の演算で引っかかったので備忘録。

最初にまとめておくと

• 行列からベクトルを引くと、自動的に列ベクトルとして計算されてしまう。
• ベクトルを転置してもエラーになる。
• 行列を転置してから列ベクトルを引くと、列ごとに計算してくれる。
• apply関数を使っても同様の結果を得られる。
• colVarsという関数はあるが、不偏分散を計算する。標本分散は計算してくれない。

各列ごとの合計・平均

各列(Column)ごとに合計や平均を計算したければ、baseパッケージにはcolSumsやcolMeansという関数が用意されている。

> # Make a matrix
>   mat = matrix(1:15, nrow = 5, ncol = 3)
>   print(mat)
     [,1] [,2] [,3]
[1,]    1    6   11
[2,]    2    7   12
[3,]    3    8   13
[4,]    4    9   14
[5,]    5   10   15
> 
>   # Column Sum
> 
>   matsum = colSums(mat)
>   print(matsum)
[1] 15 40 65
>   
>   # column mean
>   vec = colMeans(mat)
>   print(vec)
[1]  3  8 13

各列ごとの分散・標準偏差

しかし、列ごとの標準偏差や分散を行う関数は用意されていない。例えばブートストラップで各パラメーターの標準誤差を知りたければ、1000回推定した4つのパラメータを1000 x 4の行列に保存して、列ごとの分散を計算して標準誤差を算出するケースがある。
RfastというパッケージにはcolVarsという関数が用意されているが、これは不偏分散を計算してしまう。(分母がn-1)
各列の平方和は本当は10になり、行数は5なので、標本分散は2になるのだが、5−1の4で割って2.5になっている。
標本分散と不偏分散を使い分けるオプションはcolVarsには用意されていないようだ。

> Rfast::colVars(mat)
[1] 2.5 2.5 2.5

そこで、行ごとに列ごとの平均をひいて、各要素を2乗したものを列ごとに合計し、行数で割る、という手順で分散を計算したい。
ご存知の通り、行列からスカラーを引くと、各要素からスカラーを引いた行列が返される。

> mat - 3
     [,1] [,2] [,3]
[1,]   -2    3    8
[2,]   -1    4    9
[3,]    0    5   10
[4,]    1    6   11
[5,]    2    7   12

しかし、行列からそのままベクトルを引くと、列ごとに平均したベクトルは行ベクトル（横向き）のつもりで、Rでは列ベクトル（縦向き）ベクトルとして扱われるので、左上の要素から、縦向きに引き算を繰り返し行ってしまう。

> mat
     [,1] [,2] [,3]
[1,]    1    6   11
[2,]    2    7   12
[3,]    3    8   13
[4,]    4    9   14
[5,]    5   10   15

> vec
[1]  3  8 13

> mat - vec
     [,1] [,2] [,3]
[1,]   -2   -7    3
[2,]   -6    4   -1
[3,]  -10    0   10
[4,]    1   -4    6
[5,]   -3    7    2

おわかりだろうか？　4行目1列目の下の要素は4だが、ベクトルの1番目の要素の3をひいた1になっていて、5行目1列目が5だったのは、ベクトル2番めの8をひいた−3になっている。次に1行目2列めに移り、ベクトルの3番目の要素である13をひいて−7になり、ベクトルの1番目に戻って繰り返される。

じゃあベクトルを転置して、行ベクトルとして計算すればいいのではないのか、と思うが、うまくいかない。

> # Error: non-conformable
>   mat - t(vec)
Error in mat - t(vec) : non-conformable arrays

元の行列の方を転置して、列ベクトルで列ごとの演算として計算するとうまくいくようだ。

> # Inefficient?
>   t(t(mat) - vec)
     [,1] [,2] [,3]
[1,]   -2   -2   -2
[2,]   -1   -1   -1
[3,]    0    0    0
[4,]    1    1    1
[5,]    2    2    2

Apply 関数を使う方法

Webで見つけたQ&Aでは、apply関数を使う方法が提案されていた。
R help - HELP: How to subtract a vector out of each row of a matrix or array

> apply(mat,1,function(row) row - vec)
     [,1] [,2] [,3] [,4] [,5]
[1,]   -2   -1    0    1    2
[2,]   -2   -1    0    1    2
[3,]   -2   -1    0    1    2

しかし、帰ってくる行列はなぜか転置されているので、結局は転置し直さないといけない。
applyの方を使って計算した各行の分散は以下になる。

> mat2 = t(apply(mat,1,function(row) row - vec))
>   mat22 = mat2^2
>   colMeans(mat22)
[1] 2 2 2