2018年6月3日 星期日

如何將統計結果表(Output table)和數據圖合併?一個融資限制的R實做案例

在數據分析時,往往會遇到如何將「統計結果的表格」和「數據圖」合併的問題。例如,做完迴歸時的殘差診斷圖,能否把結果表和圖放在一起。有不少的統計分析圖會有一大塊空白,此時如果能把估計輸出表格嵌入就好了。這篇文的重心在如何將結果表視覺化,以及如何「圖&表」合併。
我們用一個財經融資限制問題 (Financing Constraint) 當範例,做一個含Tobin Q的迴歸分析。

步驟一:載入資料,我們載入套件pdR內建的panel data,再用time average轉成cross-section data。如下:
library(pdR)
data(invest)   # 讀入Tobin Q 數據,美國565家上市公司 15年數據
整理資料。Tobin Q除100是scaling
     myData=data.frame(invest[,c(1,3,4)],Q=invest[,2]/100)
計算投資比率橫斷面平均
     inv=apply(matrix(myData[,1], 565, , byrow=TRUE),1,mean) 
計算現金流量比率橫斷面平均
     cash=apply(matrix(myData[,2], 565, , byrow=TRUE),1,mean)
計算短期債務比率橫斷面平均
     debt=apply(matrix(myData[,3], 565, , byrow=TRUE),1,mean)
計算Tobin Q橫斷面平均
     Q=apply(matrix(myData[,4], 565, ,byrow = TRUE),1,mean) 
產生新的資料集
     myData=data.frame(inv,cash,debt,Q)

接下來執行線性迴歸
     output1=lm(inv~cash+debt+Q,data=myData)

問題一。係數表 Output table 顯示的幾個方法
方法1:由 coef(summary(output1)) 呈現,如Figure 1最上方的顯示
Figure 1
     由 Figure 1 可知Tobin Q對廠商投資率的平均影響,統計上不顯著。

方法2:利用套件papeR,如下

   papeR::prettify(summary(output1),digits=2, confint=TRUE)

上述顯示在Figure 1的中間,這個 Output table 多了信賴區間,且有 rownames,因為papeR只能處理data.frame,去掉 rownames 的方法,就是更換 rownames,如下:

text1=papeR::prettify(summary(output1),digits=2, confint=FALSE)
text11=text1[,-1]
rownames(text11)=text1[,1]
text11

text11 的顯示在Figure 1最下方


問題二。一圖一表:視覺化 Output table 的幾種狀況
     如果要把估計係數視覺化呈現來比較,一個最常用的方法就是使用套件coefplot內的函數coefplot,如下

coefplot::coefplot(output1)


Figure 2. coefplot::coefplot(output1)
      coefplot 將點估計(point estimates)標在中心點,然後左右再增疊兩條95%信賴區間的線。95%的信賴區間跨過0代表統計上相對不顯著。

     第2個方法是利用箭頭方向,來顯示出解釋變數和被解釋變數的關係,可以利用SEM的路徑圖。如下:

library(semPlot)

semPaths(output1, what="est", edge.label.cex=1, rotation=4, nCharNodes=0, shapeMan="ellipse")
Figure 3. 簡單路徑圖
     Figure 3左邊的三角形是截距,右邊圈圈是節點(nodes),正係數是綠色,負係數是紅色。現的粗細反應了數值相對大小。另外,what="std" 則會將參數標準化,從而截距會消失。Figure 3的圖形publication的效果不好,因為左邊上下都有一大塊的空白,而且因為路徑圖沒有統計顯著性。我們可以用legend的技術在左上方(topleft)嵌入Output table。這個技術的關鍵在於利用了capture.output 函數轉換係數表。如下與 Figure 4:

legend("topleft", capture.output(round(coef(summary(output1)),4)), cex=0.8, bty="n", xjust=0, text.col="blue")
Figure 4.

如果用 capture.output 處理text11,就會添增顯著水準的星號。如下與 Figure 5: 
legend("topleft", capture.output(text11), cex=0.8, bty="n", xjust =1)

Figure 5
箭頭顯示圖形,並列Output table有助於整體檢視,尤其是模型稍微龐大一些,如下例包含了交叉效果與Figure 6:
output2=lm(inv~(cash+debt+Q)^3,data=myData)
text2=papeR::prettify(summary(output2),digits=2, confint=FALSE)
text22=text2[,-1]
rownames(text22)=text2[,1]

semPaths(output2, what="est", edge.label.cex=1, rotation=4, nCharNodes=0, shapeMan="ellipse") 
legend("topleft",capture.output(text22), cex=0.8, bty="n", text.col="blue")




Figure 6
     Figure 6 的估計表顯示出單獨的Tobin Q廠商投資率沒有顯著影響,但是,和其他因子結合後,則加大了對投資率影響。當然,這種路徑圖目的在顯示變數結構,事實上,採用 "std" 會將所有的數值標準化,比較也會清楚。



問題三。多圖一表:將 Output 和殘差診斷圖合併

     現在的例子是「多圖一表」,所以和前面「一圖一表」的狀況不同。一圖一表時,表可以用 legend 嵌入,多圖一表時,需要把表轉換成圖,再併圖。

我們用殘差診斷圖為範例。首先,我們先看殘差診斷圖,
par(mfrow=c(2,2))
plot(output1)
par(mfrow=c(1,1))

這四張圖是線性模型很重要的診斷圖,我們來看看如何嵌入結果表
Figure 7
     

Case 1. 選兩張圖,併一個表。
     Step1. 使用layout將圖形頁面分割,然後指定放入位置。

           layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))

     Step 2.將表轉成圖。使用套件 gplots的函數textplot。這個函數比sinkplot好用許多,且不需要宣告一個環境來工作。如下:

     gplots::textplot(text11, halign=c("center"), valign=c("top"), mar=c(0,0,5,0))

     綜合以上,下面四行就產生 Figure 8

layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
gplots::textplot(text11, halign=c("center"), valign=c("top"), mar=c(0,0,5,0))
plot(output1,which=c(1,5))
par(mfrow=c(1,1))


Figure 8
     末行 par(mfrow=c(1,1)) 是 reset 頁面分割,也就是還原。避免之後的分割會受到影響。這務必是一個習慣。

Case 2.  把4 張圖,併一個表
     殘差診斷圖有6個,在不指定情況,會輸出主要的4個圖。如果要這樣呈現,就在於頁面分割法,如下與Figure 9

layout(matrix(c(1,1,2,3,4,5), 3, 2, byrow = TRUE))
gplots::textplot(text11, halign=c("center"), valign=c("top"),mar=c(0,0,5,0))
plot(output1)
par(mfrow=c(1,1))



Figure 9

由上面的介紹,legend的處理不若轉成圖形整齊,要如何使用,可依需要取捨。