2018年12月6日 星期四

單變量GARCH

再1 個月,我將在CRAN上更新R套件JFE 1.4,我這次將新增單變量 GARCH 的綜合處理。有興趣的公測的同好,可以這樣
2. 解壓縮後,會有幾個檔案,除了資料檔,開啟關鍵 _main.R,設定工作目錄,再執行一行指令 source("JFE.src") 即可
3. 載入的資料必須是 .RData,xts 時間序列格式,原始數據是DJ 30個股收盤價,我們從Transform Price Data選擇目錄和資料,可將之轉換成報酬率。
4. 點選最下方的GARCH,就會出現 如下GARCH pane,上面應該包含所有rugarch套件內所有的選項。
5. 上方左框是mean Eq.的一筆資料,我們選DJ 指數,中間是mean Eq.的其他變數,最右邊是variance Eq.的其他變數。

圖1. GARCH pane
資料選完後,點選您需要執行的模式,如model,
distribution 和 garch order。按OK估計完後,會自動產生兩個物件:
1. 12張圖,如下
圖2. 12 張圖

2. 估計結果將自動儲存在工作目錄內,如 圖1. 主控台的小框框。用 load("apARCH_ged_21.RData") 重新載入後,內存物件名稱都是myFit,可以重複取出所有需要的資訊。
JFE 1.4 即將納入這項功能,JFE 1.5 將再納入multivariate GARCH 的多項功能。歡迎指教

2018年11月18日 星期日

pdR 1.6

     pdR 1.6 is updated and is on CRAN now (install.packages("pdR") works).  In the previous version, pdR contains several econometric routines, for example, panel threshold model and panel unit root (seasonal). Version 1.6 includes the "Hansen B. E. (2000) Sample Splitting and Threshold Estimation. Econometrica, 68, 575-603."  Hansen (2000) is an econometric version of regression tree, but has solid econometric foundation.   
     This version has added three functions:
          SMPLSplit_het() tests the statistical significance of thresholds(sample splitting) by LM statistic and bootstrap p-value and
          SMPLSplit_est() estimates the model. Finally,
          SMPLSplit_example() offers a learning example.

     I modified Hansen's original code and gives a learning function to use Sample Splitting model. SMPLSplit_example() is a function containing learning codes, which executes Hansen's example of Durlauf-Johnson growth data, but has a detailed explanation of  how sample is classified.

Step 1. Prepare the data and declare options
library(pdR)
data("dur_john.rda")
rep <- 200
trim_per <- 0.15
dep=1  # Column number of  dependent  variables
indep=c(2,3,4,5) # Column number of regime-dependent independent  variables
th1=6  #column of the first threshold variable
th2=7  #column of the second threshold variable

Step 2. We then execute the function SMPLSplit_example() below
OUT=SMPLSplit_example(data=dur_john, dep, indep, th1, th2, trim_per, rep, plot=1)
     I gave a output function to collect important information, which can be show below
OUT$TEST
OUT$Hypothesis
OUT$Threshold

stat=matrix(as.numeric(OUT$TEST), byrow = TRUE,8,2)
colnames(stat)=c("F-Stat","P-value")
rownames(stat)=OUT$Hypothesis
stat
In this example, it tests 8 hypothesis and estimates two models, the test can be summarized below:
> stat
                                                                                   F-Stat           P-value
Testing for a First Sample Split, Using GDP60        12.601835   0.085
Testing for a First Sample Split, Using Literacy       10.786273   0.186
Testing for a Second Sample Split, Using GDP60    11.009342   0.161
Testing for a Second Sample Split, Using Literacy  12.091351   0.073
Testing for a Third Sample Split, Using GDP60        7.423359   0.558
Testing for a Third Sample Split, Using Literacy       9.233321   0.186
Testing for a Third Sample Split, Using GDP60        9.415477   0.223
Testing for a Third Sample Split, Using Literacy       9.327215   0.176

     Moreover, the computational details are listed below, which details the implication of sample splitting; for individual need, users can retrieve the source code of SMPLSplit_example() to modify. This code has three levels and I add a conclusion note right below each stage in blue. 

Level 1: Testing for a First Sample Split

<1-1> Testing for a First Sample Split, Using GDP60
Test of Null of No Threshold Against Alternative of Threshold
Allowing Heteroskedastic Errors (White Corrected)
Number of Bootstrap Replications  1000
Trimming Percentage               0.15
Threshold Estimate                833
LM-test for no threshold          12.60184
Bootstrap P-Value                 0.085

<1-2> Testing for a First Sample Split, Using Literacy
Test of Null of No Threshold Against Alternative of Threshold
Allowing Heteroskedastic Errors (White Corrected)
Number of Bootstrap Replications  1000
Trimming Percentage               0.15
Threshold Estimate                10
LM-test for no threshold          10.78627
Bootstrap P-Value                 0.186

=== Estimate First Sample Split, Using GDP60 as Threshold ===

1. Global OLS Estimation
Dependent Variable:      gdpGrowth
Heteroskedasticity Correction Used
Variable       Estimate        St Error
----------------------------------------
Constant       2.8262496      0.72472989
logGDP60      -0.2815736      0.05335405
Inv_GDP       0.4918965      0.10434991
popGrowth      -0.5467001      0.23277830
School       0.2395413      0.06560082

Observations:                       96
Degrees of Freedom:                 91
Sum of Squared Errors:              9.622743
Residual Variance:                  0.1057444
R-squared:                          0.4805041
Heteroskedasticity Test (P-Value):  0.2008698
****************************************************

2. Threshold Estimation
Threshold Variable:                 GDP60
Threshold Estimate:                 863
0.95 Confidence Interval:           [594.0000, 1794.0000]
Sum of Squared Errors:              8.024881
Residual Variance:                  0.09331257
Joint R-squared:                    0.5667667
Heteroskedasticity Test (P-Value):  0.1565995
****************************************************

     Regime 1: GDP60<=863.000000

Parameter Estimates
Variable       Estimate        St Error
----------------------------------------
Constant       4.3120283      1.62679940
logGDP60      -0.6569710      0.21761579
Inv_GDP       0.2277417      0.07160391
popGrowth      -0.2948695      0.33677597
School       0.0180607      0.09685598
----------------------------------------
0.95 Confidence Regions for Parameters
Variable       Low               High
----------------------------------------
Constant       0.68755168       9.5624430
logGDP60      -1.25006934      -0.1464676
Inv_GDP       0.02470878       0.5740113
popGrowth      -1.51316271       0.9224831
School      -0.24700520       0.4397419

Observations:                       18
Degrees of Freedom:                 13
Sum of Squared Errors:              0.6742722
Residual Variance:                  0.05186709
R-squared:                          0.5165107
****************************************************

     Regime 2: GDP60>863.000000

Parameter Estimates
Variable       Estimate        St Error
----------------------------------------
Constant       3.6630685      0.71904747
logGDP60      -0.3233915      0.06144147
Inv_GDP       0.4957500      0.14497428
popGrowth      -0.4876940      0.25532245
School       0.3569407      0.08996972
----------------------------------------
0.95 Confidence Regions for Parameters
Variable       Low               High
----------------------------------------
Constant       1.84478944       5.79544463
logGDP60      -0.52300712      -0.18203224
Inv_GDP       0.18229032       0.95436145
popGrowth      -1.06852144       0.03368612
School      -0.08479933       0.54918938

Observations:                       78
Degrees of Freedom:                 73
Sum of Squared Errors:              7.350609
Residual Variance:                  0.1006933
R-squared:                          0.5492007


####################################################
####################################################
We check output above to determine which way to go.
Because regime ' GDP60 > 863 ' has more obs, Level 2 continues
#####################################################
#####################################################

Level 2
Sub-Sample GDP60  > 863

<2-1> Testing for a Second Sample Split, Using GDP60

Test of Null of No Threshold Against Alternative of Threshold
Allowing Heteroskedastic Errors (White Corrected)
Number of Bootstrap Replications  1000
Trimming Percentage               0.15
Threshold Estimate                1410
LM-test for no threshold          11.00934
Bootstrap P-Value                 0.161


<2-2> Testing for a Second Sample Split, Using Literacy

Test of Null of No Threshold Against Alternative of Threshold
Allowing Heteroskedastic Errors (White Corrected)
Number of Bootstrap Replications  1000
Trimming Percentage               0.15
Threshold Estimate                57
LM-test for no threshold          12.09135
Bootstrap P-Value                 0.073



=== Estimate Second Sample Split, Using Literacy as Threshold ===

1. Global OLS Estimation
Dependent Variable:      gdpGrowth
Heteroskedasticity Correction Used
Variable       Estimate        St Error
----------------------------------------
Constant       3.6630685      0.71904747
logGDP60      -0.3233915      0.06144147
Inv_GDP       0.4957500      0.14497428
popGrowth      -0.4876940      0.25532245
School       0.3569407      0.08996972

Observations:                       78
Degrees of Freedom:                 73
Sum of Squared Errors:              7.350609
Residual Variance:                  0.1006933
R-squared:                          0.5492007
Heteroskedasticity Test (P-Value):  0.362994
****************************************************

2. Threshold Estimation
Threshold Variable:                 Literacy
Threshold Estimate:                 45
0.95 Confidence Interval:           [19.0000, 57.0000]
Sum of Squared Errors:              6.198249
Residual Variance:                  0.09115072
Joint R-squared:                    0.6198728
Heteroskedasticity Test (P-Value):  0.551557
****************************************************

     Regime 1: Literacy<=45.000000

Parameter Estimates
Variable       Estimate        St Error
----------------------------------------
Constant       2.0922814      1.8701913
logGDP60      -0.1163168      0.1648712
Inv_GDP       0.1726423      0.2161463
popGrowth      -0.3902176      0.5161286
School       0.4525333      0.1168289
----------------------------------------
0.95 Confidence Regions for Parameters
Variable       Low               High
----------------------------------------
Constant      -1.8538485      6.1610409
logGDP60      -0.4708739      0.3447183
Inv_GDP      -0.3181809      0.6554835
popGrowth      -1.5932356      0.8683475
School       0.1953806      0.7459717

Observations:                       30
Degrees of Freedom:                 25
Sum of Squared Errors:              2.609994
Residual Variance:                  0.1043998
R-squared:                          0.5780366
****************************************************

     Regime 2: Literacy>45.000000

Parameter Estimates
Variable       Estimate        St Error
----------------------------------------
Constant       4.31048461      0.96515976
logGDP60      -0.39502686      0.06102754
Inv_GDP       0.83363986      0.13936990
popGrowth      -0.41800939      0.26956672
School       0.09457796      0.13489238
----------------------------------------
0.95 Confidence Regions for Parameters
Variable       Low               High
----------------------------------------
Constant       1.5818922       6.2291568
logGDP60      -0.5354340      -0.2495195
Inv_GDP       0.4271144       1.1318791
popGrowth      -1.1133443       0.1145814
School      -0.3443145       0.4020814

Observations:                       48
Degrees of Freedom:                 43
Sum of Squared Errors:              3.588255
Residual Variance:                  0.08344778
R-squared:                          0.5821175


####################################################
####################################################
We check outputs above to determine which way to go.
Because both sub-regimes by Literacy have similar obs, Level 3 continues
############################################################
############################################################

Level 3
3A: Given GDP60 > 863 , Sub-Sample Literacy <= 45

<3A-1> Testing for a Third Sample Split, Using GDP60
Test of Null of No Threshold Against Alternative of Threshold
Allowing Heteroskedastic Errors (White Corrected)
Number of Bootstrap Replications  1000
Trimming Percentage               0.15
Threshold Estimate                1618
LM-test for no threshold          7.423359
Bootstrap P-Value                 0.558


<3A-2> Testing for a Third Sample Split, Using Literacy
Test of Null of No Threshold Against Alternative of Threshold
Allowing Heteroskedastic Errors (White Corrected)
Number of Bootstrap Replications  1000
Trimming Percentage               0.15
Threshold Estimate                26
LM-test for no threshold          9.233321
Bootstrap P-Value                 0.186


3B: Given GDP60 > 863 , Sub-Sample Literacy > 45

<3B-1> Testing for a Third Sample Split, Using GDP60
Test of Null of No Threshold Against Alternative of Threshold
Allowing Heteroskedastic Errors (White Corrected)
Number of Bootstrap Replications  1000
Trimming Percentage               0.15
Threshold Estimate                3493
LM-test for no threshold          9.415477
Bootstrap P-Value                 0.223


<3B-2> Testing for a Third Sample Split, Using Literacy
Test of Null of No Threshold Against Alternative of Threshold
Allowing Heteroskedastic Errors (White Corrected)
Number of Bootstrap Replications  1000
Trimming Percentage               0.15
Threshold Estimate                83
LM-test for no threshold          9.327215
Bootstrap P-Value                 0.176

####################################################
####################################################
Because, by the Bootstrap P-Value, the LM tests of both sub-regimes do not have
statistical significance, we then stop here without estimating the third level estimation.
#####################################################################
#####################################################################

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