回帰分析の結果から共分散分解
暑すぎてモチベがアレなんで血迷って誰得関数を実装していくシリーズのやつ
Intro. 共分散分解とは
定義
神林龍(2017)『正規の世界, 非正規の世界』*1を読んでいたら共分散分解なるものが出てきた(p.270)。
あんまりなじみのない分解手法だったのだが、やっていることは単純だ。
- 作者: 神林龍
- 出版社/メーカー: 慶應義塾大学出版会
- 発売日: 2017/11/08
- メディア: 単行本
- この商品を含むブログ (4件) を見る
とりあえず出てきた例に即して説明すると、企業jに所属する個人iの賃金をとすると、推定モデルを、
として回帰武分析を得られた推定値から、以下のように「共分散分解」ができる。
神林はこの分解を用いて、賃金分散におけるグループ間のばらつき(第一項)とグループ内のばらつき(第二、三項)のシェアの時系列比較をおこなっている*2。
証明
ちなみに上式が成立することは以下のように証明できる。
簡潔にとすると、右辺はになる。
ここで一般に共分散の定義はであるゆえ、
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)
結論と考察
なんかオチっぽくなってしまうのだが、こうやって算出したがトータルの分散に占める割合は、いわゆる分散説明率と等しい。
じゃあ共分散分解の必要性ってなんなんだ無駄骨じゃんって思ったんだが、単純なOLSではなくランダム効果モデルなどの階層的なモデルにおいて、
Totalの分散を①Lev.1変数による分散、②Lev.2 変数による分散、③ランダム効果による分散、④lev.1残差による分散、みたいに分解する場合には使えるかもしれない。
うーん.....