回帰分析の結果から共分散分解

暑すぎてモチベがアレなんで血迷って誰得関数を実装していくシリーズのやつ

Intro. 共分散分解とは

定義

神林龍（2017）『正規の世界, 非正規の世界』*1を読んでいたら共分散分解なるものが出てきた（p.270）。
あんまりなじみのない分解手法だったのだが、やっていることは単純だ。

作者: 神林龍
出版社/メーカー: 慶應義塾大学出版会
発売日: 2017/11/08
メディア: 単行本
この商品を含むブログ (4件) を見る

とりあえず出てきた例に即して説明すると、企業jに所属する個人iの賃金を $w_{ij}$ とすると、推定モデルを、
$w_{ij} = X_{ij}\beta + Y_j \gamma + \epsilon_{ij}$
として回帰武分析を得られた推定値から、以下のように「共分散分解」ができる。

$Var(w_{ij}) = Cov(w_{ij}, X_{ij}\beta ) + Cov( w_{ij} , Y_j \gamma) +Cov( w_{ij} , \epsilon_{ij})$

神林はこの分解を用いて、賃金分散におけるグループ間のばらつき（第一項）とグループ内のばらつき（第二、三項）のシェアの時系列比較をおこなっている*2。

証明

ちなみに上式が成立することは以下のように証明できる。

簡潔に $w = a + b+ c$ とすると、右辺は $Cov( w,a) + Cov( w, b) +Cov( w, c)$ になる。
ここで一般に共分散の定義は $Cov( X, Y) = E[XY] - E[X]E[Y]$ であるゆえ、
$Cov( w,a) + Cov( w, b) +Cov( w, c)$
$= E[wa] + E[wb] + E[wc] - E[w]( E[a] + E[b] + E[c])$
$= E[w (a + b+ c) ] - E[w]E[w]= E[w^2] - ( E[w])^2 = Var(w)$

Implementation

実装してみる。

covar_dcp <- function(lm_result, grp_idx_list=NULL){

#---- 内部関数定義 -----#
  # 共分散を得る関数
  calc_covar <- function(x , y){
    e_x <- mean(x, na.rm=T)
    e_y <- mean(y, na.rm=T)
    e_xy <- mean(x*y , na.rm=T)
    return( e_xy - e_x*e_y)
  }
  # 分散を得る関数
    calc_v <- function(x){
    n <- length(x)
    return( ( var(x) * ( (n-1) /n )) )
  } #内部関数

#---- 前処理 -----#
  Ys <- model.frame(lm_result )[,1] #従属変数取得
  Xs <- model.matrix(lm_result)[, -1] #独立変数取得
  betas <- coef(lm_result)[-1] #係数取得
  resids <- residuals(lm_result)
  
  if(is.null(grp_idx_list)){
    grp_idx_list <- list(1:(ncol(Xs)))
  }else if( !(base::setequal( 1:(ncol(Xs)), unlist( grp_idx_list))) ){
    stop("Error, ただしく変数が分割されません")
  } #if grp_idx_list 
  
#---- 分解の実行 -----#

  #Total variance
  y_var <- calc_v(Ys)
  
  x_beta <-  Xs %*% diag(betas) #係数を独立変数にかける
  xy_cov_df  <- apply( x_beta, 2, calc_covar , y=Ys)
  residy_cov <- calc_covar( Ys , resids)
  
  xy_cov <- lapply( grp_idx_list, function(idx){
    sum( xy_cov_df[idx])
  })#lapply
  xy_cov <- as.data.frame(xy_cov)
  xy_cov_clm <- unlist( lapply( grp_idx_list, function( idx){
   if( max(idx)==min(idx)){
     return(paste0("Var", max(idx)))
   }else{
     return(paste0("Var", min(idx), "-", max(idx)))
   } #fun
  })) #unlist / lapply
 colnames(xy_cov) <- xy_cov_clm  
  
 
return( list(Total= y_var , Cov_xy = xy_cov , Cov_residy = residy_cov))
}# function

単純な処理なんでとくに書くこともないんだけど、文系的には独立変数の行列と係数のベクトルをどうやって掛け合わせるのかってときに

 x_beta <-  Xs %*% diag(betas)

で対角行列を右からかけると、各列をスカラー倍できるのを利用できたのはえらいのでは。と、自分をほめる。

結論と考察

なんかオチっぽくなってしまうのだが、こうやって算出した $Cov(w_{ij}, X_{ij}\beta )$ がトータルの分散に占める割合は、いわゆる分散説明率 $R^2$ と等しい。
じゃあ共分散分解の必要性ってなんなんだ無駄骨じゃんって思ったんだが、単純なOLSではなくランダム効果モデルなどの階層的なモデルにおいて、
Totalの分散を①Lev.1変数による分散、②Lev.2 変数による分散、③ランダム効果による分散、④lev.1残差による分散、みたいに分解する場合には使えるかもしれない。