俺しか得しない関数シリーズ。要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 )