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的處理不若轉成圖形整齊,要如何使用,可依需要取捨。

2018年5月10日 星期四

R 的視覺化之二:3D入門



R的立體3D繪圖,可分為靜態,互動與動態三種。這篇文章介紹靜態的兩種型態,互動式和動態另文介紹。
3D繪圖在財經數據上相當重要,尤其是期間結構商品,如果在接上時間序列,就會更繁複。我們先從簡單的開始。主要用兩個套件:scatterplot3d和plot3D。plot3D有更多可以表現的風格。資料還是用 reshape 套件內建的 tips

I. 使用 scatterplot3d
library(scatterplot3d)
tips=reshape::tips
names(tips)                                     
used=tips[,c("total_bill","size","tip")]

Case 1. 基本3D繪圖
scatterplot3d(used,main="3D Scatter Plot", xlab="Total bill",ylab="table size",zlab="Tip",pch=16, color="steelblue")

Case 2. 分組
   我們依照小費高低,用cut函數分成三群:L=Low tip, M=Medium tip, H=High tip。如下
group=cut(tips$tip,3,labels =c("L","M","H"))
   再給定形狀,就是pch使用的編號
myshapes=c(16,17,18)
   再把編號對應到樣本
myshapes=myshapes[as.numeric(group)] 
   給定三群對應顏色
mycols=c("green","blue", "red")
   再把三個色對應到樣本
mycols=mycols[as.numeric(group)]
scatterplot3d(used,main="3D Scatter Plot", xlab="Total bill",ylab="table size",zlab="Tip",pch=myshapes, color=mycols,box=FALSE,grid=TRUE)


Case 3. 添加 legends
須要上添加物時,必須先把圖物件存起來,下例為my3d
plot.new()
my3d=scatterplot3d(used,main="3D Scatter Plot", xlab="Total bill",ylab="table size",zlab="Tip",pch= 16, color=mycols)
第1個方法是指定座標,使用 xyx.convert()
legend(my3d$xyz.convert(50,5,4),legend=levels(group),col=c("green","blue", "red"),pch=16)
第2個方法是用文字 "bottom","center","right"等等。bty="n" 可取消 legend 框框。
legend("bottomright",legend=levels(group),col=c("green","blue", "red"),pch=16, bty="n")

第3個方法使用 inset指定 legend 和 origin 的距離。下圖同時也指定legend文字水平呈現。
legend("bottom",legend=levels(group),col=c("green","blue", "red"),pch=16:18,inset=-0.15, xpd=TRUE, horiz=TRUE, bty="n")

Case 4. 在不同面添加格線
這功能須要一個小工具graph3D,由git上裝置
devtools::install_github("kassambara/graph3d")
library(graph3d)
plot.new()
s3d=scatterplot3d(used,main="3D Scatter Plot", xlab="Total bill",ylab="table size",zlab="Tip",pch= "", box=FALSE,grid=FALSE)
   指定添加格線平面
s3d_addgrids(used, grid=c("xy","xz","yz"))
s3d$points3d(used,pch=myshapes, col=mycols)
繪圖如下:

上例中的繪圖,風格較少,同時顯示的角度調整不容易(基本上是不能調),plot3D則具備更多的功能。

II. 使用plot3D


library(plot3D)
宣告座標資料。
x=used$total_bill
y=used$size
z=used$tip

Case 1. 標準繪圖,如不需顯示右邊的色棒,colkey=FALSE就可以。
scatter3D(x,y,z,main="3D Scatter Plot", xlab="Total bill", ylab="table size", zlab="Tip", pch=18, clab=c("Tip, US$"), colkey=TRUE)


Case 2. 分組上色
scatter3D(x,y,z,main="3D Scatter Plot", xlab="Total bill",ylab="table size", zlab="Tip", pch=18, col.var=as.integer(group), col=c("#E69F00","#56B4E9", "#B2182B"), clab=c("Tip, US$"))

Case 3. 背景形式 bty="b2"
scatter3D(x,y,z,main="3D Scatter Plot", xlab="Total bill", ylab="table size", zlab="Tip", bty="b2", colkey=FALSE)
Case 4. 背景形式 bty="g" 是 ggplot2 的風格。
scatter3D(x,y,z,main="3D Scatter Plot", xlab="Total bill",ylab="table size", zlab="Tip", bty="g", cex=1.5, colkey=FALSE)
Case 5. 檢視角度,之前的圖,都是俯視,如果要平行檢視,則 phi=0。如下
scatter3D(x,y,z,main="3D Scatter Plot", xlab="Total bill",ylab="table size", zlab="Tip", bty="g", clab=c("Tip, US$"), pch=20, cex=2, phi=0)
Case 6.  在座標刻度做明確顯示,ticktype="detailed",這個圖也是這筆資料最佳的檢視圖。
scatter3D(x,y,z,main="3D Scatter Plot", xlab="Total bill",ylab="table size", zlab="Tip", bty="g", clab=c("Tip, US$"), pch=20,cex=2,phi=0, ticktype="detailed")


這個圖的分群方式,如果依照 time/day/sex 劃分的話,會呈現出更多的色聚型態,資料科學往往須要檢視型態。這技術就很有用。

資料最佳檢視角度個個不同,所以,互動式圖形就可以利用滑鼠來調整角度,更為簡便。下次的文就來介紹互動式繪圖套件。


2018年5月3日 星期四

R 的視覺化之一:風格美學篇

ggplot2 應該是R社群內最受歡迎的繪圖套件之一。繪圖除了結構清楚,ggplot2 還有就是有很多的風格模版,例如,華爾街日報,經濟學人等等。本文介紹如何使用主題套件ggthemes和其他package套入入風格,程式碼內風格指令反藍,圖形的資料結構反紅,其餘是內建的設定選項。範例資料檔是 tips.csv,這是200筆紀錄小費金額的數據,餐廳發現每張帳單的小費,差距頗大。因此想知道,有哪些因素影響了小費額度。資料如圖


tip =小費金額,美金
total_bill =帳單金額, 美金
sex = 結帳者性別
smoker=結帳者是否吸煙
day=用餐日
time=用餐時段
                         size=帳單顧客人頭數(table size)

require("ggplot2")
tips<-read.csv("tips.csv")
.df <- data.frame(y = tips$tip, x = tips$day, z = tips$sex)

第1個圖是ggplot2內建的盒鬚圖。
.plot <- ggplot(data = .df, aes(x = factor(x), y = y, fill = z)) + 
  stat_boxplot(geom = "errorbar", position = position_dodge(width = 0.9), 
  width = 0.5) +  geom_boxplot(position = position_dodge(width = 0.9)) + 
  xlab("day") + ylab("tip") +  labs(fill = "sex") +
  theme_bw(base_size = 14, base_family = "sans")
print(.plot)
ggplot2 內建風格

第2個圖是經濟學人
.plot <- ggplot(data = .df, aes(x = factor(x), y = y, fill = z)) + 
  stat_boxplot(geom = "errorbar", position = position_dodge(width = 0.9), 
  width = 0.5) +  geom_boxplot(position = position_dodge(width = 0.9)) + 
  xlab("day") + ylab("tip") +  labs(fill = "sex")
  ggthemes::theme_economist(base_size = 14, base_family = "sans")
print(.plot)


經濟學人風格


第3個圖是華爾街日報風格,必須由第三方套件 RcmdrPlugin.KMggplot2 呼叫
.plot <- ggplot(data = .df, aes(x = factor(x), y = y, fill = z)) + 
  stat_boxplot(geom = "errorbar", position = position_dodge(width = 0.9), 
  width = 0.5) +  geom_boxplot(position = position_dodge(width = 0.9)) + 
  xlab("day") + ylab("tip") +  labs(fill = "sex") +
  RcmdrPlugin.KMggplot2::theme_wsj2(base_size = 14, base_family = "sans")
print(.plot)
華爾街日報風格


第4個圖是著名的Few風格,是視覺化大師 Stephen Few的設計風格,Few著有Show me the numbers-- Designing Tables and Graphs to Enlighten一書,超級一流。
.plot <- ggplot(data = .df, aes(x = factor(x), y = y, fill = z)) + 
  stat_boxplot(geom = "errorbar", position = position_dodge(width = 0.9), 
  width = 0.5) +  geom_boxplot(position = position_dodge(width = 0.9)) + 
  xlab("day") + ylab("tip") +  labs(fill = "sex") +
  ggthemes::theme_few(base_size = 14, base_family = "sans")
print(.plot)
Few風格

第5個圖是著名的網站538風格,fivethirtyeight是預測高手Nate Silver所創立,Nate著有Noices and Signals一書,中譯「精準預測」
.plot <- ggplot(data = .df, aes(x = factor(x), y = y, fill = z)) + 
  stat_boxplot(geom = "errorbar", position = position_dodge(width = 0.9), 
  width = 0.5) +  geom_boxplot(position = position_dodge(width = 0.9)) + 
  xlab("day") + ylab("tip") +  labs(fill = "sex") +
  ggthemes::theme_fivethirtyeight(base_size = 14, base_family = "sans")  
print(.plot)  
538風格

最後一個是利用 igray 風格呈現新的資料結構,把前面的圖,再依照星期和吸煙與否做成2x2=4格。
.df <- data.frame(y = tips$tip, x = tips$day,z = tips$sex, s = tips$time, t = tips$smoker)
.plot <- ggplot(data = .df, aes(x = factor(x), y = y, fill = z)) + 
  stat_boxplot(geom = "errorbar", position = position_dodge(width = 0.9), 
  width = 0.5) + 
  geom_boxplot(position = position_dodge(width = 0.9)) + 
  facet_grid(s ~ t) + 
  xlab("day") +  ylab("tip") + labs(fill = "sex") + 
  ggthemes::theme_igray(base_size = 14, base_family = "sans") + 
  theme(panel.spacing = unit(0.3, "lines"))
print(.plot)
 igray 風格

2018年4月18日 星期三

mice::md.pattern() 缺值型態解讀

數據規模龐大時,對於缺值的認識就很重要。R有多個套件檢視缺值狀態,如VIM和mice。VIM將缺值型態視覺化,mice的函數 md.pattern()可以呈現一個比較清楚的表,md.pattern()也是GUI rattle 所採用,如下圖解釋如下:




1.      裡面的0=缺值,1=有值
2.      最右邊的直欄從0-4,則是缺值個數。例如,0就是8個變數同時都無缺值的型態,1就是8個變數只缺1個的例子。但是,缺哪一個則不詳。
3.      最左邊的直欄,代表上述型態有多少列。例如,1575代表有1575列的資料是圓滿的,88代表雙缺Occupation Employment這兩個,有88(Rows)。詳細見下面第4點括弧。
4.      內列(row) {0,1} 代表了缺值出現在哪個:例如,第1列無缺,故8個變數都是1。第2列代表單缺Age(承上,有39);第3列代表單缺Employment(承上,有40);第9列代表單缺Deduction(承上,有34) ;第11列代表雙缺Occupation Employment (承上,有88)
5.      末列代表行(Column)變數總缺值。例如,圈圈42=Age這個變數的紀錄中,有42個缺值。
6.      右下角560代表整體資料的區缺值個數。

VIM的視覺化,則可以呈現對比的比較
更多: http://rstudio-pubs-static.s3.amazonaws.com/4625_fa990d611f024ea69e7e2b10dd228fe7.html

2018年4月11日 星期三

使用 sparklyr 的經驗分享之一: 連接失敗多半是 Java_Home 沒指定

R處理大規模數據時,ETL往往需要克服一些問題。以讀取資料為例,需要一些技巧來提高讀資料的速度。原始資料讀取後,轉換存取用 feather 或 pigz 都是高效能的好方法,或以記憶體管理為方向的 bigmemory 都是好方法。本文主要解說啟動 sparklyr 連接失敗的問題如何克服,多半是Java_home沒有適當的指定,或版本須要更新。

以 1.5G 的範例檔 exampleFile.csv (1千萬列,10個行變數),硬體 64位元 PC Intel Core(TM) i7-4700 CPU 3.40 GHz,32 GB 記憶體。我們先看看三個原始資料的子載入速度:

read.csv 讀取需 8分鐘。
> system.time(csvData <- read.csv("exampleFile.csv"))
   user  system elapsed 
 473.68    4.11  494.14

data.table::fread 需 55秒:
> system.time(freadData <- fread("exampleFile.csv")) 
   user  system elapsed 
  53.96    0.41   55.31

版本2.2.1的 Spark,R 套件 sparklyr 則只20秒
> library(sparklyr)
> sc <- spark_connect(master="local")
> system.time(sparkDF <- spark_read_csv(sc, name= "spDF", path= "exampleFile.csv"))
   user  system elapsed 
   0.10    0.00   19.22

sparkDF 是 Spark data frame, sparklyr 無法存成正規的 Spark 資料庫,所以 sparklyr 只能以 data.frame 存在工作空間。但是, feather pigz 都可以有效管理Spark data frame的存取。只要是 spark data frame, Spark 的機器學習工具都可以有效使用,sparklyr 的環境也可以把 R data frame轉成 spark data frame.

  sparklyr 有一些優點,但是啟動 spark_connect(master="local")會遇到一些錯誤訊息,這些錯誤訊息,關鍵在於Java_Home路徑沒有清楚的在電腦中被指定。如果使用時出現錯誤,最主要的問題都是Java_Home不存在。檢查方法如下:在 cmd 環境執行 echo %Java_Home% 查一下Java_Home路徑是否存在。

  如果沒有裝好 Java,則請先安裝Java。
  如果你還沒裝置套件,則:
install.packages("sparklyr") #裝套件

sparklyr::spark_install(version = "2.2.1") #安裝Spark

  再由控制台->系統,如下圖 


  將正確的 Java_Home 路徑書寫進去,退出後,重新啟動R,目前為止,都可以成功執行 sc <- spark_connect(master="local")的連結。一旦連結成功,建立 Spark data frame的方法很簡單:
                 RDF_tbl <- copy_to(sc, R_data.frame)

  只要是Spark data frame,所有 Spark 的機器學習工具都可以使用。 15 GB 的資料(美國上市公司財報資料庫),很快就可以完成演算,不會當機。詳細的說明,下一則再來繼續。









2018年3月31日 星期六

JFE 1.2: Portfolio Backtesting

  1. Portfolio Backtesting here is basically based on R package fPortfolio. The method is designed for portfolio optimization, JFE provides more covariance estimators and GMVP strategy for backtesting. JFE offers a comprehensive computation(Backtesting All in One) for 6 covariance estimators combined with 2 strategies, which is a little bit time-consuming, 3-min for DJ30 dataset.
  2. To use this function, you must have a multivariate time series dataset with R format, xts is most encourgaed; and the file is saved in .RData or .rda.
  3. If the loaded data is price, then you have to pull down the menu and choose Transform Price Data, else, Load Returns Data.
  4. The Next-Month Advice is the output bottom is the assets weights suggestion computed by backtesting for the next period from the end of data. The rolling length is 1 month and estimation is 1 year, which are not allowed to change so far.
  5. Download the dataset DJ30.RData to practice, which is close price of Dow Jones 30 component Stocks.
While data inputting is done, click Backtesting button, and access the pane below. First of all, you have to "Pick 1 bench asset", if you choose "None", JFE will compute cross-section average as Bench. Secondly, you may pick asset(s) you do not need. The following 4 blocks are its specification, we offer 3 multivariate covariance estimators: sample covariance, Ledoit-Wolf Baysian Shrinkage and the one based on Student t distribution.


Clicking OK, we give out an output in R console, "Next-Month Advice" is the diversified portfolio suggested by backtesting.  For our dataset, it suggests more weights on DJI (Dow Jones Index) components. 


Whether the backtesting strategy is good enough to follow, we can check the graph below: the middle right one shows our optimized portfolio (red line) outperforms the bench one where we pick Russia (RTS), although drawdown is not satisfactory.

Clicking "Backtesting All in One" pop out the pane below. Users have only to pick assets and portfolio risk, including risk free rate and smooth lambda, it will give a comprehensive computation including two strategies (Tangency & GMVP) and 6 multivariate covariance estimators.  If dataset is not large, it may take roughly 3 minutes.  If the dataset is somewhat large, for example, 50-year daily with 100 assets, it is time-consuming.  In our experiment with parallel computation(Clusters), it takes 10 minutes.




JFE 1.2: Assets Selection

Assets Selection

  1. Assets Selection here is basically designed for International Assets Selection, which is based on, inequality of specific performance index, for example, Sharpe ratio.  Interested readers please refer to Bekaert and Hodrick(2009, International Financial Management, PP.466-468).  We provide more performance indices based on the R package PerformanceAnalytics.
  2. To use this function, you must have a multivariate time series dataset with R format, xts is most encourgaed; and the file is saved in .RData or .rda.
  3. If the loaded data is price, then you have to pull down the menu and choose Transform Price Data, else, Load Returns Data
  4. Download the dataset world20.RData to practice, which is 20 country price indices of MSCI world index.