論理の流刑地

流罪に遭い戸惑う世捨人の雑記

Macで、回帰系の手法の結果を、Excelに貼るだけの状態で出力する関数

俺しか得しない関数シリーズ。要dplyr。

※主な今後の改善点を備忘のために、冒頭に書いておく
✔︎ glmやlmerへの対応を
✔︎ モデル指標のところ、小数点桁を指定可能に
✔︎ 変数をmodelごとに差し引きした時に、表側の変数名とCoefやS.E.の係数推定値の垂直位置がずれるのを直す(今は手動。そんなにめんどくさくないからね)

Motivation

・めっちゃRegression Tableつくらなあかんねんけど、いちいち整形するのめんどい。
・でも私はTex族じゃない→texregとかstargazerとかはぴったりこない。ExcelとかいうMicrosoft様のむっちゃ使いにくい神ソフトで表を整形したい
・でもたくさんあるRegressionのコードごとに関数定義するのはダサい

Source Code

## 便利関数定義 ##
# 引数①:  lm()関数の結果を格納したlist
# 引数②: 数値間の区切り文字(デフォルトはタブ区切り)
# 引数③: 回帰係数をExp()したものを表示するか否か ※対数変換済所得などを従属変数に指定したとき用
# 引数④: fileの出力場所  ※デフォルトはDesktopに"test.txt"という名前で出力する
# 引数⑤: どの関数を使ったか。Defaultはlm() # lm, clm(順序ロジット, ordinalパッケージ), rq(分位点回帰、quantregパッケージ)に対応( As of 17/6/27)

reg_to_excel <- function( res_list , blk = "\t" , exp_coef = F , file_name = "/Users/Ronri_Rukeichi/Desktop/test.txt", reg_function = "lm" ){

def_funcs <- c("lm", "rq", "clm")

if(!reg_function %in% def_funcs ){
stop("指定された回帰関数は未対応です, clm or lm or rqのどれかを指定すべきです")
} #if

info_list <- lapply( res_list , function( rslt ){

#実装上の注意 ,,,, (A)回帰係数系(独立変数に付随する情報)と、(B)モデル指標系、は別個にみたほうがよいので、別の情報としてかえす
if( reg_function == "lm" ){
# - - -Case1: lm関数(通常のOLS)- - -#
# summary objectを取得
smr <- summary(rslt) 

## 回帰係数、S.E.を取得
coef_arr <-  coef(smr)[,1] #回帰係数
se_arr <- coef(smr)[, 2] #標準誤差
pval_arr <- coef(smr)[, 4] #p値
var_name_arr <- rownames(coef(smr)) #変数名
coef_df <- data.frame(Var = var_name_arr , Coef = coef_arr, S.E.= se_arr  , Pval = pval_arr )  #DF化

## Model指標を取得 ### 
N <-  length(smr$residuals) # 分析ケース数
nparam <-  slot(logLik((rslt)),"df") #number of parameters
val_logLik <- as.numeric(logLik((rslt))) #対数尤度 
val_AIC <- AIC( rslt) 
val_BIC <-  BIC( rslt)
r_squared <-  summary(rslt )$r.squared
adj_r_squared <- summary(rslt )$adj.r.squared
cri_names <- c("N" , "Number of Parameters" , "log-likelihood" , "AIC" , "BIC" , "R Squared" , "Adj. R Squared")
cri_values <- c( N , nparam,  val_logLik , val_AIC , val_BIC , r_squared , adj_r_squared )

#-- convert to data frame 
cri_df <- data.frame( Var = cri_names , Value = cri_values )

## Tableに記載する情報を返す ##
return(list(  Coef = coef_df ,  Model = cri_df ) )  #

}else if( reg_function == "clm"){
#- - -Case2: clm関数("ordinal" package , ordered logit/probit )  - - -#
# ordered logit modelの場合も同様に実装する..

#--Summary Object --
smr <- summary(rslt) 

# - - Coefに関する情報行列は、lmの場合と同様の実装- -
coef_arr <-  coef(smr)[,1] #回帰係数
se_arr <- coef(smr)[, 2] #標準誤差
pval_arr <- coef(smr)[, 4] #p値
var_name_arr <- rownames(coef(smr)) #変数名
coef_df <- data.frame(Var = var_name_arr , Coef = coef_arr, S.E.= se_arr  , Pval = pval_arr )  #DF化

#- - -Model Information Matirx - - -#
N <-  smr$nobs # 分析ケース数
nparam <-  slot(logLik((rslt)),"df") #number of parameters
val_logLik <- as.numeric(logLik((rslt))) #対数尤度 
val_AIC <- AIC( rslt) 
val_BIC <-  BIC( rslt) 

cri_names <- c("N" , "Number of Parameters" , "log-likelihood" , "AIC" , "BIC" )
cri_values <-  c( N  , nparam , val_logLik , val_AIC , val_BIC ) 

#- - -Convert to Data Frame - - -#
cri_df <- data.frame( Var = cri_names , Value = cri_values ) 

## Tableに記載する情報を返す ##
return(list(  Coef = coef_df ,  Model = cri_df ) )  #

}else if( reg_function == "rq"){
#- - -Case3 : rq関数(分位点回帰モデル , "quantreg" package ) - - -#


##  出力イメージ(他の回帰系との相違点)
# ・OLSやGLMの場合はregression tableに記載するときは, 「Model XX」の下に各情報が並ぶ感じになるけど、
#  Quantile Regressionの場合は、「Model XX」の下に 「YY%点」が併記される形になるので、一階層分深めにdata構造を作っておく

##--- summary object --- ##
smr <- summary( rslt ) 

coef_df_list <- list()
cri_df_list <- list()

for( i in 1: length( smr ) ){
## 回帰係数周りの情報 ##
#print(smr)

coef_arr <- ((smr[[i]])$coefficients)[,1]
se_arr <- (smr[[i]]$coefficients)[,2]
pval_arr <- (smr[[i]]$coefficients)[,4] 
var_name_arr <- rownames((smr[[i]]$coefficients))
coef_df <- data.frame(Var = var_name_arr , Coef = coef_arr, S.E.= se_arr  , Pval = pval_arr )  #DF化

coef_df_list[[i]] <- coef_df

## - - -Model Information --- ##
N <- nrow( rslt$residual)
nparam <- slot( logLik(rslt),"df")
val_logLik <- logLik(rslt)[i]
val_AIC <- AIC(rslt)[i]
val_BIC <-  -2 * val_logLik + log( N ) * nparam
val_tau <- smr[[i]]$tau
cri_names <- c("N" , "Number of Parameters" , "log-likelihood" , "AIC" , "BIC", "Quantile" )
cri_values <-  c( N  , nparam , val_logLik , val_AIC , val_BIC,  val_tau) 

#- - -Convert to Data Frame - - -#
cri_df <- data.frame( Var = cri_names , Value = cri_values ) 
cri_df_list[[i]] <- cri_df

} ##for

## Tableに記載する情報を返す ##
return(list(  Coef = coef_df_list ,  Model = cri_df_list ) )  #

} #if 回帰関数別
} ) #lapply

###====: info_listからテキストに変換:=== ###
## Quantileの場合が難しい → とりあえず、tauの下に各モデルを配置するような感じにしとく, 


# 関数定義: 文字列のvectorをk行空きにする関数
mix_blank <- function( str_arr ,k=1 , blank_str=" " ){
return(as.vector( matrix(c( str_arr , rep( blank_str , length(str_arr) * k  ) ) , nrow= (k + 1) , byrow=T ) ) )
} #function 

# - - 有意確率→Symbolの変換関数
pval_mark <- function( pval){
return(ifelse( pval >= 0.10  , " " , ifelse( pval  >= 0.05 , "†" , ifelse( pval >= 0.01 , "*" , ifelse( pval >= 0.001 , "**" , "***" ) ) ) ))
} #function


#整形用の関数を定義 
#・標準誤差は[ ] で囲んで回帰係数の下に表示. 
#・カスタマイズできるようにしてもいいが、とりあえず後回し。
library(dplyr) #dplyrパッケージを読み込む

coef_to_txt <- function( coef_df, exp_flag = F,  round_digit = 4  ){
#必要情報だけ抜いてくる
info_df <-  dplyr::select( coef_df, Coef, S.E., Pval)

# S.E.を [ ] で囲んでCoefの下にもってきて一列にする
se_str <-   paste("[", formatC(round(info_df$S.E. , round_digit ) , format= "f" , digits = round_digit) , "]" , sep="")
coef_str <- formatC( round( info_df$Coef ,  round_digit ) , format="f" , digits = round_digit )
est_str <-   as.vector( matrix( c( coef_str , se_str ) , byrow=T , nrow = 2 ) )

# 有意確率を表すsymbolの文字列生成 
p_smbl <-  pval_mark( info_df$Pval) 
p_smbl_arr <- mix_blank(p_smbl , k= 1 ) #1行空きにする

# expのCoef乗を返す(Flagに応じて)
if( exp_flag == TRUE ){
exp_beta <-  formatC( round( exp( info_df$Coef ),  round_digit )  , digits = round_digit , format="f" )
exp_beta_str <- mix_blank( exp_beta, k = 1 )  #一行空きにする。

# 行列型にして返す
return( cbind( est_str, p_smbl_arr , exp_beta_str) )

}else{
#行列型にして返す
return( cbind( est_str, p_smbl_arr ) ) #
} #if
} #function( coef_to_txt)

# 簡易関数, Model Criteria version
mdl_to_txt <- function( model_df , exp_flag=F){
cri_values <- formatC( model_df$Value  , format="f" , digits = 3 ) 
## exp(回帰係数)を表示するか否かで、何列空きにするかを変える ##
if( exp_flag ){
ret_mat <- cbind( cri_values , rep( " " , length(cri_values ) ),rep( " " , length(cri_values ) ) )
}else{
ret_mat <- cbind( cri_values , rep( " " , length(cri_values ) ) )
} #if
return(ret_mat ) 
} # function 

# 行数が違う行列をcbindするときに、数が少ない方を空の文字列で埋める関数
cbind_v2 <- function( mat1 , mat2 , blank_str = " "){
n1 <- nrow( mat1 ) 
n2 <- nrow( mat2 ) 

if( n1 > n2 ){
#第1の行列のほうが、第2の行列よりも行数が長い場合..
add_mat <-  matrix( rep(" " ,  ncol(mat2 ) * (n1- n2 )), ncol= ncol(mat2) )
# - -
mat2 <- rbind( mat2 ,add_mat ) 
ret_mat <- ( cbind( mat1 , mat2 )  ) 
}else if( n1 < n2 ) {

add_mat <-  matrix( rep(" " ,  ncol(mat1 ) * (n2- n1 )), ncol= ncol(mat1) )
# - -
mat1 <- rbind( mat1 ,add_mat ) 
ret_mat<- ( cbind( mat1 , mat2 ) )

}else if( n1 == n2 ) {
# なんもしない
ret_mat <- cbind( mat1 , mat2 ) 
}

return( ret_mat )  #統合した行列を返す。
} #function

##- - - 先に外枠を完成させてから、中に埋めるものを完成させる方式- - - ##
#- - とりあえず普通のOLSとかordered logitとかの、シンプルな場合から考えていく - - -#

if(reg_function %in% c( "lm" , "clm")){

##- - -通常パターン- - - ##
## ====表側の情報 ==== ##
# 変数名
var_names <- as.character( info_list[[length(info_list)]]$Coef$Var )  
# モデル指標の名前
model_info_names <- as.character( info_list[[length(info_list)]]$Model$Var )  

#- - 行列形式に変換してk × 1の行列で保存,, 
# 変数名に関しては一行空きにする。 


# 表側の表示テキストを格納した行列
left_side <-  matrix(c(as.vector( matrix(c(var_names , rep(" ",length(var_names))), nrow=2,  byrow=T)), model_info_names) ,ncol=1)

## ==== 表頭の情報 === ##
 if( exp_coef  == T){
title_arr <-  c("β", " ", "Exp(β)" )
}else{
title_arr <-  c("Coef." , " " ) 
}
 #回帰係数のexp()を表記するか(引数③)によって2列or3列が変わる。
model_name_arr <- paste( "Model" ,1:length( info_list )  ,sep = "" ) #
#print( title_arr)
#print(model_name_arr)

# 表頭の表示テキストを格納した行列
upper_side <- matrix( c( mix_blank( model_name_arr , k = ifelse( exp_coef ==T , 2 , 1 ) ) , rep(title_arr,length(res_list )  )) , byrow = T , nrow= 2)

## ==== 表のcontents(回帰係数の推定結果) ==== ##
tbl_content <- NA #初期化
lower_side <- NA #初期化


lapply( info_list ,  function( info_elm  ){
#print(coef_to_txt(info_elm$Coef, exp_flag = exp_coef))
if( is.na( tbl_content  ) ){
#初回の処理
tbl_content <<- coef_to_txt(info_elm$Coef, exp_flag = exp_coef)
lower_side <<- mdl_to_txt( info_elm$Model,exp_flag = exp_coef) 
} else{
# 2回目以降の処理
tbl_content <<- cbind_v2(tbl_content , coef_to_txt(info_elm$Coef , exp_flag = exp_coef)  ) #横につなげていく
lower_side <<- cbind( lower_side ,  mdl_to_txt( info_elm$Model,exp_flag = exp_coef)  )
} #if
}) #lapply

## ここまで生成した文字列行列を統合
ret_mat <-  cbind( rbind( matrix( c(" " , " ") , ncol=1 ) , left_side ) , rbind( upper_side , tbl_content , lower_side ) )

#return(ret_mat) 


}else if(reg_function ==  "rq"){
##- - -Quantile Regressionの場合  - - -##

info_list2 <- list()
#- - 入れ物だけつくっておく
for( i in 1:length(info_list[[1]]$Model) ){

info_list2[[i]] <- list()
} #for

cnt <- 0


lapply(info_list , function(info_elm){
cnt <<- cnt + 1 #インクリメント
#まずCoefの情報を、順番に、格納していく
coef_cnt <- 0 
lapply(info_elm$Coef  , function( coef_info){
coef_cnt <<- coef_cnt + 1 
#print(">>>> <<<<<")
#print(coef_cnt)
info_list2[[coef_cnt]][[cnt]] <<- list()
info_list2[[coef_cnt]][[cnt]][["Coef"]] <<- coef_info


} )#lapply Coef

#print(info_list2)
# 次はModelの情報
model_cnt <- 0 
lapply(info_elm$Model  , function( mdl_info){
model_cnt <<- model_cnt + 1 
info_list2[[model_cnt]][[cnt]][["Model"]] <<- mdl_info
} )#lapply Coef
}) #lapply info_list



ret_mat_list <- lapply( info_list2 ,function(info_list){
## ====表側の情報 ==== ##
# 変数名
var_names <- as.character( info_list[[length(info_list)]]$Coef$Var )  
# モデル指標の名前
model_info_names <- as.character( info_list[[length(info_list)]]$Model$Var )  

#- - 行列形式に変換してk × 1の行列で保存,, 
# 変数名に関しては一行空きにする。 

# 表側の表示テキストを格納した行列
left_side <-  matrix(c(as.vector( matrix(c(var_names , rep(" ",length(var_names))), nrow=2,  byrow=T)), model_info_names) ,ncol=1)

## ==== 表頭の情報 === ##
 if( exp_coef  == T){
title_arr <-  c("β", " ", "Exp(β)" )
}else{
title_arr <-  c("Coef." , " " ) 
}
 #回帰係数のexp()を表記するか(引数③)によって2列or3列が変わる。
model_name_arr <- paste( "Model" ,1:length( info_list )  ,sep = "" ) #
#print( title_arr)
#print(model_name_arr)

# 表頭の表示テキストを格納した行列
upper_side <- matrix( c( mix_blank( model_name_arr , k = ifelse( exp_coef ==T , 2 , 1 ) ) , rep(title_arr,length(res_list )  )) , byrow = T , nrow= 2)

## ==== 表のcontents(回帰係数の推定結果) ==== ##
tbl_content <- NA #初期化
lower_side <- NA #初期化


lapply( info_list ,  function( info_elm  ){
#print(coef_to_txt(info_elm$Coef, exp_flag = exp_coef))
if( is.na( tbl_content  ) ){
#初回の処理
tbl_content <<- coef_to_txt(info_elm$Coef, exp_flag = exp_coef)
lower_side <<- mdl_to_txt( info_elm$Model,exp_flag = exp_coef) 
} else{
# 2回目以降の処理
tbl_content <<- cbind_v2(tbl_content , coef_to_txt(info_elm$Coef , exp_flag = exp_coef)  ) #横につなげていく
lower_side <<- cbind( lower_side ,  mdl_to_txt( info_elm$Model,exp_flag = exp_coef)  )
} #if
}) #lapply

## ここまで生成した文字列行列を統合
ret_mat <-  cbind( rbind( matrix( c(" " , " ") , ncol=1 ) , left_side ) , rbind( upper_side , tbl_content , lower_side ) )

return(ret_mat) 
}) #lapply info_list


#return(ret_mat_list) 


#内容をタテにつなげる
ret_mat <- NULL
lapply( ret_mat_list , function(x){
ret_mat <<- rbind( ret_mat , x ) 
}) #lapply

#return(ret_mat)
} # IF


# 最後にfileに出力する。
write.table(ret_mat , sep=blk , file=file_name , quote=F , row.names = F  , col.names=F )
return( info_list ) 
} #reg_to_excel

How to Use

例えばOLSの場合。

#推定
ols_res1 <- lm( ols_fml , data = test_data  )
ols_res2 <- lm( ols_fml , data = test_data  )
ols_res3 <- lm( ols_fml , data = test_data  )

# 推定結果をlistにまとめる
ols_res_l <- list( ols_res1 , ols_res2,ols_res3) 

# 上記関数により出力
reg_to_excel(ols_res_l , reg_function= "lm", file_name  = "/Users/Ronri_Rukeichi/Desktop/ols_reg.txt" , exp_coef = T )

できあがり

出力したtxtファイルを全コピしてExcelにpasteするだけで
f:id:ronri_rukeichi:20170629112024p:plain

あらふしぎ!


Enjoy!!