R语言在经济学中的应用

class: center, middle, inverse, title-slide

# R语言在经济学中的应用
### 南开大学周恩来政府管理学院 吕小康
### 2017-07-13

---

# R简介

### R是一个免费自由且跨平台通用的统计计算与绘图软件。
  
- 它有Windows、Mac、Linux等版本，均可免费下载使用。

### 从[R主页](https://www.r-project.org/)中选择[download R](https://cran.r-project.org/mirrors.html)链接可下载到对应操作系统的R安装程序。

- 打开链接后的网页会提示选择相应的[CRAN](https://cran.r-project.org/mirrors.html)镜像站。目前全球有超过一百个CRAN镜像站 ，用户可选择就近下载。

---

# R与STATA等统计软件的区别

### R为开源免费的软件，其他基本为商业付费软件。

- 如果你有钱，可以只选贵的、不选对的；但如果你没钱……

###  R是一种脚本语言，强调英文命令操作。

- R的学习比较费时、对汉字编码不友好，但掌握之后的自由性更强

### R在数据可视化上的表现更佳，选择更丰富。

- R的统计绘图是它最有标志性的功能，可以制作达到出版的各种图形

---

# R在经济学中的综合应用

R及与之相关的配套开源软件（如RStudio）已构成一个丰富的数据分析网络生态，具有同类软件很难同时满足的多种可能性。

### 用于课程教学

### 用于数据获取与预处理

### 用于数据的计量分析

### 用于数据可视化

### 用于撰写学术报告

……

---

# R 已可方便导入各类型的数据

利用[Hadley](http://hadley.nz/)等人开发的诸多R包，已可方便导入各种类型的数据，这为R成为一种“兼容并蓄”的统计分析软件奠定了重要基础。

### readxl包: xls or xlsx

```r
library(readxl)
read_excel("file.xlsx") 
```

### rvest包：在线文本

### data.table包：导入大数据文件（> 100 G）

---

# R 已可方便导入各类型的数据

### Haven包： SAS, SPSS, STATA

```r
# SAS
read_sas("mtcars.sas7bdat")
write_sas(mtcars, "mtcars.sas7bdat")

# SPSS
read_sav("mtcars.sav")
write_sav(mtcars, "mtcars.sav")

# Stata
read_dta("mtcars.dta")
write_dta(mtcars, "mtcars.dta")
```

利用R，几乎可以分析类型任何类型的数据，而避免在各类统计软件之间相互转化和跟踪。

---

# 作为课堂教学的辅助软件

可以作为两门经济学基础课程的教学辅助软件

- 《概率论与数理统计》
- 《计量经济学》

我本人在清华大学出版社2017年出版的[《R语言统计学基础》](https://www.amazon.cn/%E6%95%B0%E9%87%8F%E7%BB%8F%E6%B5%8E%E5%AD%A6%E7%B3%BB%E5%88%97%E4%B8%9B%E4%B9%A6-R%E8%AF%AD%E8%A8%80%E7%BB%9F%E8%AE%A1%E5%AD%A6%E5%9F%BA%E7%A1%80-%E5%90%95%E5%B0%8F%E5%BA%B7/dp/B06XGR6LJZ/ref=sr_1_1?ie=UTF8&qid=1490843285&sr=8-1&keywords=%E5%90%95%E5%B0%8F%E5%BA%B7)，内容差不多覆盖经济学类入门概率论与数理统计的教学要求，全程使用R作为分析和绘图软件。

---

## 概率论与数理统计的课程教学

这里以抽样分布的教学设计为例进行说明。

#### 抽样分布（sampling distribution）的基础知识

####　标准定义

- 抽样分布是样本统计量的分布。这显然精确而“无用”，即不能帮助人真正理解什么是抽样分布。

- 实质是重复抽样的假想前提下形成的一个统计推论框架，它在现实中是不一定存在的。

---

## 抽样分布的具体解释

- 抽样分布是对同一总体，做出相同样本容量的、重复(无限) 多次的简单随机抽样取遍样本统计量的所有可能值后所体现出来的取值规律性。

- 对这一规律性，采用概率论的基本知识加以描述，即可概括为某一**分布**（distribution），也即 `$F(x) = P(X \leq x)$`
  
- 如果能够找到 `$F(x)$` 的精确数学表达形式，后续的统计推论即可基于这一概念框架而得到概率意义上的精准推导。
  
---

## 抽样分布的教学难点

### 理论框架的“非现实性”
  
- 现实中的研究通常只能有一次抽样，不可能对同一总体进行反复抽样从而得到关于样本统计量的精确分布的直观感知

### 数据并非总是随机抽样获得的

- 通过随机化实验，以及通过普通的观测收集的数据，也需要进行推论统计。但此时很难直接套用基于“重复抽样”获得的抽样分布理论。

### 某些统计量的抽样分布数学推导较为困难

- 况且，能够找到精确数学形式的分布总是少见的。很多统计量本身就是很难找到精确分布，然而推论总是要做下去……

---

## 建立经验感知的方式：模拟

抽样分布的建立需要一对互相联系的概念：总体（population）与样本（sample）。

不妨以这样的思路进行教学:

- 先从假想的理论分布（如正态分布、二项分布、指数分布等）总体进行重复抽样，模拟某一简单样本统计量的分布，并与数学推导的结果进行对比解释；

- 到假想的理论分布总体进行重复抽样，模拟某一很难或无法从数学推导获得精确分布的样本统计量的抽样分布；

- 再到从实际的、不满足特定分布的总体进行重复抽样，模拟对应的样本统计量的分布，验证数学推导的结果是否能够应用于现实，并理解数学推导的局限性与模拟的自由性

---

## 样本均值等常见统计量的抽样分布示例

已知总体 `$X \sim Exp(0.5)$`，即服从某指数分布，注意该分布本身是右偏的，且是一无限总体。现从中抽取样本量为100的样本，重复10000次，每次计算如下样本统计量，再绘制这些样本统计量的的直方图，即可在一定程度上展示该些样本统计量的抽样分布形状，以便形成直观的感知。

- 样本均值（sample mean）
- 样本标准差（samlpe standard deviation）
- 样本方差（sample variance）
- 其他需要的样本统计量

---

左上图为给定总体的概率密度图，右上图为样本均值的抽样分布示意图；左下图为样本标准差的抽样分布示意图，右下图为样本方差的抽样分布示意图。红色虚线为虚线表示总体均值、 总体标准差和总体方差所在位置。

---

## 生成样本均值抽样分布的代码

```r
sample_mean <- numeric(10000)
for (i in 1:10000) {
  set.seed(i)
  samples <- rexp(100, 0.5)
  sample_mean[i] <- mean(samples)
}
hist(
  sample_mean,
  prob = T,
  main = "Sampling Distribution for Sample Mean",
  xlab = "Sample Mean",
  breaks = 30,
  cex.main = 0.8
)
abline(v = 2, col = "red", lwd = 3)
```

其余可交由学生思考复制。

---

### 真实的观测数据：Nile 流量

`Nile`数据是`R`自带的数据，记录了尼罗河在埃及阿斯量(Ashwan) 地
区1871-1970 年这100年间的年流量值。这一数据服从什么特定的精确分布吗？--不清楚。

```r
hist(Nile, breaks = 30)
abline(v = mean(Nile), col = "red", lty = 2, lwd = 3)
```

但若以它为“总体”，再从中进行简单随机抽样，然后观察某些特定统计量（如样本均值、样本方差、样本中位数）的分布，仍可获得经验感知。

---

红色虚线表示“总体”均值所在的位置。这是不是就是中心极限定理（Central Limit Theorem）告诉我们的道理呢？

---

## 部分代码示例

```r
x <-
  data.frame(
  a = numeric(10000),
  b = numeric(10000),
  c = numeric(10000),
  d = numeric(10000)
  )
  par(mfrow = c(2, 2))
  for (i in 1:10000)
  {
  set.seed(i)
  samples <- sample(Nile, 5)
  x$a[i] <- mean(samples)
  }
  
  hist(
  x$a,
  xlab = "Mean flow",
  main = "n=5",
  probability = T,
  xlim = c(600, 1200)
  )
  abline(
  v = mean(Nile),
  lwd = 3,
  lty = 2,
  col = "red"
  )
```

---
### Nile数据，其他样本统计量的抽样分布

<img src="R4Eco201707_files/figure-html/unnamed-chunk-4-1.png" style="display: block; margin: auto;" />
样本最小值、样本最大值、样本中位数、样本第一四分位数的抽样分布示意图（样本容量为30）。红色虚线分别表示“总体”最小值、最大值、中位数和第一四分位数所在的位置。

---

### 关于代码教学的建议

教师只需在课堂讲解说明并提供一个图形的原始代码，即可要求学生仿照此代码，自行绘制其他图形，并作为作业进行布置。如此可加深学生对抽样分布形成过程的直观认识。

实际上关于抽样分布还有一些基于`RStudio`的`shiny`平台搭建的动态呈现模式，例如[Nicole Radziwill](https://radziwill.shinyapps.io/sdm-clt)就制作了相关的简单结果展示。这可以作为课堂教学的参考。

如果教师本人精力允许，可以带领学生自行制作相关网页。这样收获更大。

实际上，有了重复抽样情形下的抽样分布及模拟，教师还可利用`R`进行自助分布（Bootstrap Distribution）、随机化分布（Randomization Distribution）等理论分布的模拟，如此可将推论的背景框架推广至其他观测数据或实验数据的情形。

---

## 计量经济学的课程教学

这里使用一个经常在计量经济学中使用到的数据（`Affairs`）进行示例。这是美国 *Psychology Today* 杂志于1969年采集的关于婚外情的数据。该数据经常用于广义线性模型的示例。

```r
if(!require(AER)) install.packages("AER")
data("Affairs")
head(Affairs)
```

```
##    affairs gender age yearsmarried children religiousness education
## 4        0   male  37        10.00       no             3        18
## 5        0 female  27         4.00       no             4        14
## 11       0 female  32        15.00      yes             1        12
## 16       0   male  57        15.00      yes             5        18
## 23       0   male  22         0.75       no             2        17
## 29       0 female  32         1.50       no             2        17
##    occupation rating
## 4           7      4
## 5           6      4
## 11          1      4
## 16          6      5
## 23          6      3
## 29          5      5
```

---

## OLS回归

```r
fm_ols <- lm(affairs ~ age + yearsmarried + religiousness + occupation + rating,
data = Affairs)
summary(fm_ols)
```

```
## 
## Call:
## lm(formula = affairs ~ age + yearsmarried + religiousness + occupation + 
##     rating, data = Affairs)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0382 -1.7076 -0.7780  0.2086 12.8134 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    5.60816    0.79660   7.040 5.31e-12 ***
## age           -0.05035    0.02211  -2.278   0.0231 *  
## yearsmarried   0.16185    0.03690   4.387 1.36e-05 ***
## religiousness -0.47632    0.11131  -4.279 2.18e-05 ***
## occupation     0.10601    0.07110   1.491   0.1365    
## rating        -0.71224    0.11829  -6.021 3.03e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.087 on 595 degrees of freedom
## Multiple R-squared:  0.1314,	Adjusted R-squared:  0.1241 
## F-statistic:    18 on 5 and 595 DF,  p-value: < 2.2e-16
```

---

## OLS 回归

### 查看模型拟合值

```r
fit <- fitted(fm_ols)
head(fit)
```

```
##          4          5         11         16         23         29 
##  1.8279300  0.7779612  3.2055393 -0.1406526  2.1685688  0.2559938
```

### 查看模型残差

```r
re <- residuals(fm_ols)
head(re)
```

```
##          4          5         11         16         23         29 
## -1.8279300 -0.7779612 -3.2055393  0.1406526 -2.1685688 -0.2559938
```

---

## 查看用于模型诊断的相关图示

```r
opar <- par(no.readonly = T)
par(mfrow = c(2, 2))
plot(fm_ols)
```

```r
par(opar)
```

---

## 关于回归假设诊断的一个常见范例

线性模型最大的优点可能在于数学形式上的简单性。但这仍依赖于许多基本假定（assumptions）。对这些假定进行检验，是经济学教学和实证研究环节中不可忽视的内容。

通过一些“极端”的教学示例，可培养学生检验模型基本假定的良好习惯。

[Frank Anscombe](https://en.wikipedia.org/wiki/Frank_Anscombe)(1918-2001) 是20世纪著名的英国统计学家， 他于1973年发表了一篇具有深远影响的文章，讨论图形在统计检验中的作用。其所构造的一组数据经常被作为演示数据。

![](https://upload.wikimedia.org/wikipedia/en/d/d5/Francis_Anscombe.jpeg)

---

## Anscombe 四重奏(Anscombe’s Quartet)

请观察以下数据。

```r
anscombe
```

```
##    x1 x2 x3 x4    y1   y2    y3    y4
## 1  10 10 10  8  8.04 9.14  7.46  6.58
## 2   8  8  8  8  6.95 8.14  6.77  5.76
## 3  13 13 13  8  7.58 8.74 12.74  7.71
## 4   9  9  9  8  8.81 8.77  7.11  8.84
## 5  11 11 11  8  8.33 9.26  7.81  8.47
## 6  14 14 14  8  9.96 8.10  8.84  7.04
## 7   6  6  6  8  7.24 6.13  6.08  5.25
## 8   4  4  4 19  4.26 3.10  5.39 12.50
## 9  12 12 12  8 10.84 9.13  8.15  5.56
## 10  7  7  7  8  4.82 7.26  6.42  7.91
## 11  5  5  5  8  5.68 4.74  5.73  6.89
```

x1, x2, x3这三列完全相同，x4 与前三例不同；y1, y2, y3, y4各不相同。

---

## Anscombe的“四个”回归方程

```r
fit1 <- lm(y1 ~ x1, data = anscombe)
coefficients(fit1)
```

```
## (Intercept)          x1 
##   3.0000909   0.5000909
```

```r
fit2 <- lm(y2 ~ x2, data = anscombe)
coefficients(fit2)
```

```
## (Intercept)          x2 
##    3.000909    0.500000
```

```r
fit3 <- lm(y3 ~ x3, data = anscombe)
coefficients(fit3)
```

```
## (Intercept)          x3 
##   3.0024545   0.4997273
```

```r
fit4 <- lm(y4 ~ x4, data = anscombe)
coefficients(fit4)
```

```
## (Intercept)          x4 
##   3.0017273   0.4999091
```

---

## Anscombe的“四个”线性相关系数

```r
attach(anscombe)
cor(x1, y1)
```

```
## [1] 0.8164205
```

```r
cor(x2, y2)
```

```
## [1] 0.8162365
```

```r
cor(x3, y3)
```

```
## [1] 0.8162867
```

```r
cor(x4, y4)
```

```
## [1] 0.8165214
```

```r
detach(anscombe)
```

---

## 相同的回归系数，相同的相关系数

从近似角度看，这几乎可以统一为一个回归方程：

$$ \hat{y} = 0.5 x + 3 $$

从近似的角度看，它们还拥有相同的相关系数：0.816。

## 但是，如果绘制各自的散点图……

---

---

## 结论：统计数字会骗人！

###  不能迷信回归系数

###  结合可视化进行模型检验有其优势

###  结合R来做可视化较为便利

---

## 相关代码

```r
attach(anscombe)
opar <- par(no.readonly = T)
par(mfrow = c(2, 2))
plot(x1, y1, col = "orange", pch = 20, cex = 2)
abline(lm(y1~x1), col = "blue", lwd = 2)
plot(x2, y2, col = "orange", pch = 20, cex = 2)
abline(lm(y2~x2), col = "blue", lwd = 2)
plot(x3, y3, col = "orange", pch = 20, cex = 2)
abline(lm(y3~x3), col = "blue", lwd = 2)
plot(x4, y4, col = "orange", pch = 20, cex = 2)
abline(lm(y4~x4), col = "blue", lwd = 2)
par(opar)
detach(anscombe)
```

---

# 广义线性模型

广义线性模型（Generalized Linear Models）的一般形式：

`$$f(\mu_Y)=\beta _0 + \beta _1 X_1 + \beta _2 X_2 + \cdots + \beta _k X_k=\beta _0 + \sum ^k _{j=1} \beta _j X_j$$`

其中
- `$f(\mu_Y)$`表示响应变量的条件均值的某种函数（称为连接函数，link function）。

- 此时对 `$Y$` 不再有服从正态分布的要求，而可以服从任何指数分布族中的某一分布。

- 设定好连接函数与分布类型后，就可以利用极大似然法通过多次迭代推导出各参数值。

---

# 常用的广义线性模型

- Probit/Logistic 回归模型
- Poisson 回归模型
- Negative Binomial 回归模型
- Zero Inflation 回归模型
- Tobit 回归模型
- ……

## 这些都可通过R的相关函数方便求得。

---

# 广义线性模型

## Probit 回归

```r
fm_probit <- glm(I(affairs > 0) ~ age + yearsmarried + religiousness + occupation + rating,
data = Affairs, family = binomial(link = "probit"))
summary(fm_probit)
```

```
## 
## Call:
## glm(formula = I(affairs > 0) ~ age + yearsmarried + religiousness + 
##     occupation + rating, family = binomial(link = "probit"), 
##     data = Affairs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6143  -0.7678  -0.5841  -0.2368   2.4615  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    0.97667    0.36537   2.673 0.007516 ** 
## age           -0.02202    0.01032  -2.134 0.032821 *  
## yearsmarried   0.05990    0.01712   3.499 0.000468 ***
## religiousness -0.18365    0.05172  -3.551 0.000384 ***
## occupation     0.03751    0.03285   1.142 0.253399    
## rating        -0.27298    0.05257  -5.192 2.08e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 675.38  on 600  degrees of freedom
## Residual deviance: 614.59  on 595  degrees of freedom
## AIC: 626.59
## 
## Number of Fisher Scoring iterations: 4
```

注：`I(affairs > 0)`用于生成是否有婚外情的虚拟变量，> 0 则赋值为1，否则为0。

---

## Probit 回归

## 查看模型拟合值

```r
fit <- fitted(fm_probit)
head(fit)
```

```
##          4          5         11         16         23         29 
## 0.26779939 0.16359028 0.47303362 0.07518241 0.33574947 0.11842761
```

---

## Logistic/Logit 回归

```r
fm_logit <- glm(I(affairs > 0) ~ age + yearsmarried + religiousness + occupation + rating,
data = Affairs, family = binomial(link = "logit"))
summary(fm_logit)
```

```
## 
## Call:
## glm(formula = I(affairs > 0) ~ age + yearsmarried + religiousness + 
##     occupation + rating, family = binomial(link = "logit"), data = Affairs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6633  -0.7500  -0.5750  -0.2691   2.4189  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    1.74904    0.62802   2.785 0.005352 ** 
## age           -0.04009    0.01785  -2.245 0.024748 *  
## yearsmarried   0.10543    0.02952   3.572 0.000354 ***
## religiousness -0.32332    0.08950  -3.613 0.000303 ***
## occupation     0.07250    0.05677   1.277 0.201565    
## rating        -0.46842    0.08928  -5.247 1.55e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 675.38  on 600  degrees of freedom
## Residual deviance: 613.70  on 595  degrees of freedom
## AIC: 625.7
## 
## Number of Fisher Scoring iterations: 4
```

---

## Poisson 回归

```r
fm_pois <- glm(affairs ~ age + yearsmarried + religiousness + occupation + rating,
data = Affairs, family = poisson)
summary(fm_pois)
```

```
## 
## Call:
## glm(formula = affairs ~ age + yearsmarried + religiousness + 
##     occupation + rating, family = poisson, data = Affairs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.5968  -1.5728  -1.1627  -0.7067   8.3473  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    2.533905   0.196924  12.867  < 2e-16 ***
## age           -0.032255   0.005851  -5.512 3.54e-08 ***
## yearsmarried   0.115698   0.009908  11.677  < 2e-16 ***
## religiousness -0.354037   0.030892 -11.460  < 2e-16 ***
## occupation     0.079828   0.019449   4.105 4.05e-05 ***
## rating        -0.409443   0.027381 -14.953  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2925.5  on 600  degrees of freedom
## Residual deviance: 2360.1  on 595  degrees of freedom
## AIC: 2866.1
## 
## Number of Fisher Scoring iterations: 7
```

---

## Negative Binomial 回归

```r
if(!require(MASS)) install.packages("MASS")
fm_nb <- glm.nb(affairs ~ age + yearsmarried + religiousness + occupation + rating,
data = Affairs)
summary(fm_nb)
```

```
## 
## Call:
## glm.nb(formula = affairs ~ age + yearsmarried + religiousness + 
##     occupation + rating, data = Affairs, init.theta = 0.142555597, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1139  -0.8067  -0.6943  -0.4533   2.4548  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    2.189666   0.727326   3.011  0.00261 ** 
## age           -0.002624   0.020312  -0.129  0.89722    
## yearsmarried   0.084819   0.034205   2.480  0.01315 *  
## religiousness -0.422227   0.104081  -4.057 4.98e-05 ***
## occupation     0.060443   0.066262   0.912  0.36167    
## rating        -0.431331   0.107449  -4.014 5.96e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.1426) family taken to be 1)
## 
##     Null deviance: 390.89  on 600  degrees of freedom
## Residual deviance: 339.66  on 595  degrees of freedom
## AIC: 1470.5
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.1426 
##           Std. Err.:  0.0160 
## 
##  2 x log-likelihood:  -1456.4880
```

---

## Zero Inflation 回归

```r
if(!require(pscl)) install.packages("pscl")
fm_zero <- zeroinfl(affairs ~ age + yearsmarried + religiousness + occupation + rating | age +
yearsmarried + religiousness + occupation + rating, data = Affairs)
summary(fm_zero)
```

```
## 
## Call:
## zeroinfl(formula = affairs ~ age + yearsmarried + religiousness + 
##     occupation + rating | age + yearsmarried + religiousness + occupation + 
##     rating, data = Affairs)
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -1.4643 -0.5190 -0.3827 -0.2444 14.3993 
## 
## Count model coefficients (poisson with log link):
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    2.201940   0.210148  10.478  < 2e-16 ***
## age           -0.007238   0.006343  -1.141    0.254    
## yearsmarried   0.049224   0.010990   4.479 7.50e-06 ***
## religiousness -0.131668   0.031154  -4.226 2.37e-05 ***
## occupation     0.016029   0.020101   0.797    0.425    
## rating        -0.118672   0.028693  -4.136 3.53e-05 ***
## 
## Zero-inflation model coefficients (binomial with logit link):
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.74075    0.62879  -2.768 0.005633 ** 
## age            0.03977    0.01788   2.224 0.026137 *  
## yearsmarried  -0.10297    0.02961  -3.478 0.000506 ***
## religiousness  0.31696    0.08980   3.530 0.000416 ***
## occupation    -0.07186    0.05689  -1.263 0.206541    
## rating         0.46378    0.08944   5.185 2.16e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Number of iterations in BFGS optimization: 21 
## Log-likelihood: -765.9 on 12 Df
```

---

## Tobit 回归

```r
library(AER)
fm_tobit <- tobit(affairs ~ age + yearsmarried + religiousness + occupation + rating,
data = Affairs)
summary(fm_tobit)
```

```
## 
## Call:
## tobit(formula = affairs ~ age + yearsmarried + religiousness + 
##     occupation + rating, data = Affairs)
## 
## Observations:
##          Total  Left-censored     Uncensored Right-censored 
##            601            451            150              0 
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    8.17420    2.74145   2.982  0.00287 ** 
## age           -0.17933    0.07909  -2.267  0.02337 *  
## yearsmarried   0.55414    0.13452   4.119 3.80e-05 ***
## religiousness -1.68622    0.40375  -4.176 2.96e-05 ***
## occupation     0.32605    0.25442   1.282  0.20001    
## rating        -2.28497    0.40783  -5.603 2.11e-08 ***
## Log(scale)     2.10986    0.06710  31.444  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Scale: 8.247 
## 
## Gaussian distribution
## Number of Newton-Raphson Iterations: 4 
## Log-likelihood: -705.6 on 7 Df
## Wald-statistic: 67.71 on 5 Df, p-value: 3.0718e-13
```

---

# 在数据获取、预处理和可视化中的应用

- tidyverse系列数据处理包
  - dplyr: 数据操纵
  - tidyr: 数据操纵
  - stringr: 文本数据操纵
  - rvest: 在线抓取文本
  - ……
- 可视化系列数据处理包
  - ggplot2
  - ggtheme
  - ggvis
  - shiny
  - wordcloud2
  - ……

---

## 数据处理示例1：一手问卷调查数据

[我们项目组](https://github.com/xkdog/Seminar)目前正在编制《中国医患社会心态调查问卷》，问卷已经基本完成编制并已进行预测试。对初测数据的统计分析工作正在进行。初测问卷使用问卷星填答，要求被调查者使用自身手机或在访问员的手机上完成填答。数据示例见Excel文件。

以下命令可简单地统计被试的地理位置分布。

```r
library(readxl)
library(stringr)
library(tidyverse)
PDSurveyBasic <- read_excel("PDSurveyBasic.xlsx")
ip.location <- str_extract(PDSurveyBasic$ip, "(?<=\$).*(?=\$)") %>%
  str_split("-", n = 2, simplify = TRUE) %>%
  as_tibble %>%
  transmute(province = .[[1]], city = .[[2]]) %>% 
  group_by(province) %>% 
  summarise(n=n()) %>% 
  arrange(desc(n))
```

---
class: center

## 地理位置信息分布结果

```
## # A tibble: 27 x 2
##    province     n
##       <chr> <int>
##  1     天津    97
##  2     云南    71
##  3     辽宁    63
##  4     新疆    57
##  5     河南    47
##  6     山东    44
##  7     北京    36
##  8     四川    34
##  9     重庆    34
## 10     山西    29
## # ... with 17 more rows
```

---

## 数据获取与处理示例2

### 政府工作报告抓取与分析

传统社会科学的量化分析以对数字数据（numeric data）的量化分析为主，对文本数据（text data）的分析较少。这主要是受研究工具的局限所致。

R及Python等开源软件的出现，很大程度上改变这种现状，使得文本分析成为当下社会科学研究的一大潮流。

[中国政府网](http://www.gov.cn/guowuyuan/baogao.htm)提供了自1954年以来所有的政府工作报告全文。这里以中国政府工作报告（2017）为例做一简单的R语言示例（该示例得益于雪晴数据网陈堰平老师的讲座）。

---

## 政府工作报告的抓取与简单分析

[2017政府工作报告](http://www.gov.cn/premier/2017-03/16/content_5177940.htm)

```r
if (!require(rvest)) install.packages('rvest')
if (!require(wordcloud2)) install.packages('wordcloud2')
if (!require(jiebaR)) install.packages('jiebaR')
if (!require(stringr)) install.packages('stringr')
url2017 <-
  "http://www.gov.cn/premier/2017-03/16/content_5177940.htm"
report2017 <- read_html(url2017)
text2017 <- report2017 %>%
  html_nodes("p") %>%
  html_text() %>%
  paste(collapse = "")
writeLines(text2017, "report2017.txt")
library(jiebaR)
cutter <- worker(
  bylines = T,
  user = "./UsrWords.txt",
  stop_word = "./stopWords.txt",
  output = "report2017output.txt"
)
report_seg_file <- cutter["./report2017.txt"]
report_segged <-
  readLines("./report2017output.txt", encoding = "UTF-8")
report <- as.list(report_segged)
doc.list <- strsplit(as.character(report), split = " ")
term.table <- table(unlist(doc.list))
term.table <- sort(term.table, decreasing = TRUE)
del <- term.table < 5 | nchar(names(term.table)) < 2
term.table <- term.table[!del]
vocabDF <- as.data.frame(term.table)
```

---

## 政府工作报告的抓取与简单分析

```r
head(vocabDF, 10)
```

```
##    Var1 Freq
## 1  发展 6125
## 2  改革 3332
## 3  推进 3185
## 4  建设 2646
## 5  经济 2548
## 6  推动 2058
## 7  加快 1960
## 8  政府 1960
## 9  创新 1764
## 10 企业 1715
```

---

## 政府工作报告的抓取与简单分析

```r
library(wordcloud2)
wordcloud2(vocabDF, color = "random-light", backgroundColor = "grey")
```

<div id="htmlwidget-06c1b143915a5eaf3eb3" style="width:504px;height:504px;" class="wordcloud2 html-widget"></div>
<script type="application/json" data-for="htmlwidget-06c1b143915a5eaf3eb3">{"x":{"word":["发展","改革","推进","建设","经济","推动","加快","政府","创新","企业","完善","提高","就业","实施","中国","工作","深化","支持","政策","服务","国家","人民","社会","增长","制度","群众","地区","农村","地方","机制","稳定","治理","保护","保障","持续","合作","市场","重点","落实","试点","投资","消费","综合","产能","加大","教育","结构","领域","维护","城市","改善","国际","国内","生态","我国","降低","力度","农业","强化","提升","引导","增加","制定","质量","标准","创业","国务院","健康","民生","目标","升级","生产","文化","下降","依法","增强","战略","作用","供给","管理","规范","核心","环境","进一步","确保","生活","收费","产业","风险","鼓励","会议","监管","健全","科技","困难","能力","世界","水平","体系","一年","一批","资金","做好","超过","代表","贯彻","行政","减少","金融","经济社会","区域","全国","体制改革","亿元","责任","组织","城镇","党中央","动能","发挥","扶贫","关系","技术","全球","特别","调控","习近平","现代化","协调","需求","优化","预期","制造","中央","专项","资源","部门","城乡","出台","措施","改造","高校","更好","行业","和谐","基础","基础设施","解决","民族","取消","融资","生产总值","实体","特色","脱贫","万人","文明","新型","新增","引领","住房","保险","比例","补短","补助","打造","对外开放","公开","和平","互联网","价格","减税","建立","结构性","精神","开放","领导人","绿色","贸易","面对","能源","农产品","培育","去年","社会主义","涉企","审批","适度","项目","新兴产业","严肃查处","医疗","扎实","支撑","抓好","转型","安排","安置","办法","办好","报告","财政","参与","大国","督查","方式","房地产","覆盖","杠杆","工程","工业","构建","关键","规划","规模","环保","活力","货币政策","基层","基金","家庭","精准","军队","力量","领导","旅游","农民","努力","启动","清洁","情况","融合","实干","实行","市场化","事业","试验区","收入","淘汰","提供","体制","同志","万亿元","效应","新建","压力","研究","一带","一路","医保","优势","优质","着力","整治","政策措施","政务","执法","中心","主体","资本","13","2000","把握","变化","财政赤字","残疾人","成本","成果","传统产业","创造","达标","登记","地位","第十二届","第五次","发生","法律","方向","防治","非公有制","分类","服务业","钢铁","各类","各族人民","公共服务","公平","国防","惠民","货币","激发","既有","坚定不移","监督","建成","降费","接受","金融风险","经济运行","举措","库存","李克强","粮食","民间","内需","排放量","棚户区","贫困地区","平等","平衡","清单","清理","区间","全国人民代表大会","全民","人员","设立","深化改革","深刻","事项","手段","态势","调整","统筹","统一","推广","外贸","外商投资","稳健","问责","污染","相关","小康社会","学生","意见","意识","应对","有效性","责任制","政治","制造业","抓紧","转移","装备","资产","自主","总理","1000","2017","2025","30","北京人民大会堂","比重","毕业生","便利","步伐","财政政策","成就","成效","成长","诚挚","城乡居民","城镇化","充分发挥","出路","传统","大力发展","大众","带动","担当","单位","底线","地下","电信","动力","多层次","多种","遏制","儿童","二是","发展中国家","法治","法治化","繁荣","犯罪","防控","放宽","峰会","负担","覆盖面","感谢","高端","高效","各国","工匠","工具","公共","公共卫生","公里","公立医院","公正","攻坚","股权","关心","广泛开展","国内外","国有","国有企业","过剩","宏观调控","划分","环境保护","回稳","活动","机动车","积极性","激励机制","计划","坚定","交流","缴费","进城","进程","进出口","进展","经费","经营","就医","开发","开幕","考核","科学","科研院所","历史","利益","灵活","迈出","煤电","煤炭","面向","模式","能耗","农民工","排放","贫困人口","品牌","评价","潜力","切实","驱动","权益","全面实施","全面推行","全球化","燃煤","人口","人民币","人文","认真落实","三是","上海","社会治安","胜利","十八","十三","示范区","市场主体","税负","思想","四是","所有制","体现","挑战","条件","铁路","同胞","土壤污染","推开","退出","外资","万户","万亩","网络","稳步","五是","系列","显现","现代农业","香港","享受","销售","效率","效益","新华社","新进展","新设","严格执行","严峻","研发","养老","药品","一道","一是","隐患","优先","有序","债券","这一","支出","支付","职工","职业","中介","重大项目","重组","周年","壮大","准入","资源配置","自贸","自治","总体","作出","1.5","100","12","1300","2030","3.4","4.5","50","5000","600","8000","M2","PM2","澳门","澳门特别行政区","八项","百分点","搬迁","办事","包容","保护主义","北京","并举","补偿","补贴","布局","部署","财富","财力","财政收支","采购","残奥会","产品","产权","产权保护","产业链","产业园","成功","成立","城市群","持之以恒","赤字","充实","处置","传导","创新能力","创造力","创造条件","从严治党","存量","大病","大幅","大局","大气污染","氮氧化物","党和国家","导向","到位","道路","低于","地质","电力","叠加","定向","独特","多边贸易","多渠道","多万人次","多种形式","二氧化硫","二要","发放","法律法规","反对","反腐败","方案","防范","房价","分裂","分流","分配","奋斗","扶持","符合","腐败问题","负责","改革开放","改建","改增","干事","港澳","高度","高水平","革命","格局","各方","根本利益","耕地","更大","工作人员","公民","公平竞争","功能","功能区","共识","共享","贡献","贡献率","构筑","购房","关键环节","关口","管廊","管理制度","贯彻落实","广大干部","广义","规模经营","规则","轨道交通","国防科技","国际贸易","国家主权","国有资本","海外侨胞","海洋","海洋权益","夯实","合法权益","合一","户籍制度","化解","缓解","恢复","辉煌","回顾过去","汇率","混合","货币化","机会","基本面","基础性","基地","基调","激励","集体","记者","加工","坚强","艰巨","监控","兼并","简政放权","建设性","健身","较大","阶段","节能","结构调整","结算","经费支出","经济带","经营性","竞相","竞争力","纠正","救助","居民","居民消费","居住","居住证","局面","决不能","决不允许","决策","均衡","开工","科技进步","科研人员","颗粒物","空间","空气质量","控制","快速","拉动","蓝天","劳动","牢固","利用","联网","两岸","两个","六中全会","落户","率先","矛盾","面积","面临","民主","民主党派","内生","农村公路","农村土地","浓度","盘活","庞兴雷","配套","贫困","品质","平稳","评估","企业家","前提","强国","强军","切实落实","切实做好","侵害","渠道","权力","全面落实","全球性","燃油","人才","人大代表","人民团体","人数","三大","三农","三要","商事","商业银行","上涨","设施","深入开展","审议","生效","失业率","施策","石油","时代","食品药品","始终","示范","市场准入","事后","事业性","事中","守住","树立","数据","数量","数字","水利工程","水污染","税率","税收","四个","随机","台湾同胞","谈判","特大事故","天然气","听取","通信","同步","统计","透明","突出重点","推出","退休","退役军人","托底","拓展","外交","外资企业","万吨","万公里","万套","万元","万众","危房改造","微观","违规","围绕","文化产业","文件","稳妥","稳中求进","污染源","西部","吸纳","系统性","先进","宪法","乡村","相机","协定","协同","新旧","新举措","新一轮","信息","信心","信用","形式","形势","休闲","修订","许可","学习","学校","压减","沿线","养老金","要带头","一般性","一国两制","一律","一要","依法惩治","依法治国","亿吨","义务教育","异地","易地","因地制宜","因素","银行","引向","营造","影响","应有","拥有","勇气","用药","用于","优抚","优惠政策","优秀","玉米","预警","原则","灾害","再创","增收","增值税","债转股","战胜","长江","长期","这是","着眼","整顿","整合","正式","政府职能","政协委员","知识产权","指导","秩序","智能","置换","中国共产党","中华民族","中小企业","中央财政","重点高校","重建","主动","主攻","住院费用","注重","转换","自觉","自贸区","自然资源","自主权","宗教","综合治理","总书记","总体方案","总需求","祖国","最低","1.2","1.55","1.7","1.8","11","11.3","1100","1200","1240","1314","15","16","1600","165","1700","1900","192","2.08","2.3","2.38","2.9","20","2016","2020","21","21.3","220","24","24.5","240","29","300","340","35","3500","380","4.02","4000","420","450","5.6","5076","51.6","550","56.2","5700","6.3","6.5","6.7","6500","6700","70","74","74.4","75","795","8.5","80","8300","8400","85","9.1","9.5","900","95","安康","安全感","安全事故","安全性","安心","安置工作","氨氮","案件","奥运会","澳人治澳","八是","巴黎","板块","办理","办能","办实事","办学","办学条件","伴随","帮扶","包括","薄弱","保护区","保留","保卫战","保险制度","保险资金","保有","保障性","保证","保证率","保值","暴力","北方地区","备案","备战","奔头","本届","本科","本色","迸发","比较完善","必定会","必由之路","边防","边疆","编制","变革","变更","变数","变形","便利化","标本兼治","标志性","标准化","并重","拨付","波动","剥离","博士","不负","不高","不搞","不合理","不堪重负","不靠","不可动摇","不理","不良资产","不容","不容忽视","不辱","不忘","不懈","不懈努力","不易","不折不扣","不正之风","材料","财产权","财产权利","财税","财政性","财政预算","裁量权","采矿业","采取措施","操守","测绘","层层","差别化","差距","产权制度","产业政策","常态","常委会","常住人口","倡导","超出","超低","超额完成","扯皮","沉淀","撑起","撑腰","成蝶","成绩","成熟","成套设备","诚信","承诺","城区","惩处","惩治腐败","持久","充分准备","充满生机","充满希望","充满信心","充裕","充足率","崇尚","筹办","出访","出境","出口","出席","初级阶段","初见成效","初心","除烦","除少数","储备","储蓄率","处于","传承","窗口","闯过","创富","创建","创建活动","创新型","创业板","创意","创造财富","创造性地","垂直","慈善事业","存于","答卷","打好","打通","大胆","大规模","大会","大家庭","大局意识","大陆","大门","大事","大水","大通关","大湾","大型","大业","大义","大灾","大政方针","大中型","带薪休假","单一","当务之急","档案","党的基本","党的领导","党风廉政","党和人民","倒逼","到期","盗窃","得不到","邓小平理论","低保","低估","低迷","低效","底气","砥砺","地方税","地级","地上","地市","地条钢","地震","第二批","第三次","第三方","第十九次","第四代","第五代","电厂","电商","电子商务","顶天立地","定价","定力","定能","东北","东部","东盟","东亚","东中部","冬奥会","冬季","动摇","动用","动员","兜牢","斗争","督察","毒品","短板","对接","对内","对台","对外","对重","多边","多措","多发","多发性","多公里","多国","多难","多年","多万","多万公里","多万户","多亿美元","多亿元","多证","多重","夺取","恶化","发电","发行","发力点","发明创造","发明专利","发债","发展观","发展前景","发展潜力","法定","法规","烦心","繁重","反恐","返贫","返乡","方便群众","方针","防火墙","放开","放心","放在","非金融","非税","非正式","费用","分红","分化","分级","分离","分配制度","分析","分享","分心","奋力拼搏","丰收","丰硕","服务平台","服务体系","服务性","符合规定","幅度","辐射","福祉","腐败分子","付出","负面","负责制","负债","妇联","妇女","附加","复兴","赋予","富民","富强","覆盖率","改变","改进","改饲","改为","感受","干部","干干净净","干字","港人治港","高标准","高度一致","高峰论坛","高级别","高技术","高举","高考","高品质","高速","高速公路","高速铁路","高压","高质量","告别","革命化","个人信息","个体","各界人士","各项","各展","根基","更低","更快","更要","工程质量","工会","工商户","工商联","工资","公布","公道","公共安全","公共设施","公路","公平正义","公司","公司法人","公司制","公务员","公益","公用事业","公园","共建","共进","共青团","共商","共同体","贡献者","供电","供水","供销社","购租","孤寡老人","姑息","鼓劲","关爱","关键性","关停","观光农业","管得","管控","管理体制","贯通","灌溉面积","光缆","光说不练","光纤","广播","广告","广聚","归侨","规划设计","规律","贵在","锅炉","国防建设","国防教育","国际货币基金组织","国际收支","国家级","国六","国内长途","国企","国有资产","国资","过程","过渡","过坎","过快","孩子","海空","海绵","海内外","海上","海外","海域","旱作","捍卫","行动计划","行政村","行政监察","行政长官","杭州","航空","好事","合格","合作项目","和睦相处","和平统一","和衷共济","河长制","核电","黑恶势力","弘扬","红线","宏观政策","宏伟目标","洪涝","洪涝灾害","后备力量","后劲","后勤","厚植","互促","互动","互惠","互联互通","互信","护航","护理","花钱","划定","划转","化学","环境监测","环境污染","缓建","缓中","焕发","荒漠化","黄标车","回报","汇聚","会晤","会越","伙伴关系","货运","机构","机井","机理","积极参与","积极主动","积累","基本保障","基本法","基本方针","基本国情","急需","疾病","集成电路","集群","几个","几千万","挤出","计算","纪念","技术创新","技术改造","继续加强","继续前进","绩效","加成","加大力度","加计","加剧","加快调整","加力","加以解决","佳绩","假冒伪劣","价格上涨","价值","价值观","价值链","坚强意志","坚守","坚信","肩上","艰难险阻","监测","监察","兼容","检验","减半","减不增","减免","减轻","简并","简单","简化","简洁","简陋","简政","建成区","建军","建设者","建议","建筑","建筑业","健体","践行","僵尸","奖补","匠心","降到","降费要","交出","交汇","交通","交通拥堵","交易成本","交易额","较快","较少","教师队伍","教育资源","秸秆","节水","节用裕民","结合实际","结束","竭力","解难","解析","今明两年","金融机构","金融监管","金融市场","金融体制","金砖","紧急","紧密","紧迫","紧日子","尽保","尽责","进口","进社区","禁入","京津冀","经得起","经济渠道","经济体","经济体制","经验","经营机制","精深","精神文明","精益求精","警惕","净化","竞技","竞争","竞争机制","竞争性","敬畏","纠错","九二","旧账","救灾","就地解决","就近","居民收入","居中","举办","举债","巨大成就","具备条件","聚焦","决定性","决心","抉择","绝不允许","崛起","军工","军烈属","军民","军民团结","军委主席","军政","均等化","均衡性","开创性","开发区","开发性","开放型","开局","开启","开拓进取","开越","看齐","考核办法","苛之弊","科技领域","科技人员","科技型","科学决策","科学普及","科学性","科研经费","可靠","可控","可支配","恪尽","客观","空域","恐怖活动","扣除","枯竭","跨界","跨境","快递","快乐","快速增长","框架","困境","困难群众","扩大出口","垃圾","澜沧江","篮子","懒政怠政","劳动力","劳神","劳务输出","老旧","老年","老区","老人","乐观","乐章","累积","离不开","里程","里子","理念","理赔","理政","历史使命","立体化","利率","利润","利息","利用效率","连片","联大","联动","联防","联合国","联合体","联控","廉洁","廉守","练兵","良性","粮食安全","两岸关系","两孩","两项","林场","林区","林权","林田湖","临时","零部件","领土","领土完整","领先","令人","留守","留下","流动","流动性","流通","流向","流域","六是","隆重","隆重庆祝","陆上","路线","轮作","落地","落后","履行","履职","率拟","埋头苦干","迈进","迈向","满意","漫灌","漫游费","慢性病","铆紧","贸易协定","湄公河","煤矿","煤炭行业","每人每年","美好","美丽","梦想","免除","面源","面子","民办教育","民房","民航","民意","民用","民主监督","民主决策","民族大义","民族区域","民族团结","名列","名目繁多","名义","明处","明天","明显提高","明显增加","明显增强","命运","目的","目的地","目录","睦邻","纳入","纳税","难以解决","内地","内蒙古自治区","内外销","内在","年度","年均","年来","年末","年内","年应","凝聚","凝聚力","农村基层","农户","农惠农","农垦","努力完成","爬坡","庞大","培训","培养","配合","配套改革","朋友","棚改","蓬勃","蓬勃生机","碰硬","批评","披荆斩棘","贫困村","贫困家庭","贫困县","频发","品种","平安","平均","平台","平原","迫切","迫切愿望","破产","铺天盖地","普遍推行","普惠","普通高中","七是","期待","期权","欺诈","齐全","其长","企稳","企业财务","气候变化","气象","弃水","汽车","千方百计","签约","前行","前茅","潜能","强刺激","强基","抢劫","抢险救灾","侨眷","侨务","切身感受","锲而不舍","侵犯","勤勉","青年","轻装上阵","倾向","清算","清醒","请予","庆祝","区段","区域性","趋稳","取暖","权利","权威","权威性","权限","全程","全党全军","全方位","全国代表大会","全国人大常委会","全国政协","全力以赴","全力支持","全面推广","全域","确权","群众反映","群众利益","群众体育","热点","热点问题","人才队伍","人才资源","人大","人工智能","人居","人均","人类","人力","人力资源","人民满意","人民政府","人民政协","人命关天","人群","人人有责","人生价值","人士","人事","认定","认可","认真贯彻","任何人","韧性","日电","容错","融资难","三板","三档","三分之一","三个代表","三公","三级","三年","三四","三中","散煤","森林","山水","商谈","上年","上市","上同","上网","上限","上线","少数民族","设备","设市","社保","社会工作","社会科学","社区","涉农","身心健康","深层次","深度","深港","深远","审查","审计","审评","审慎","生产力","生产率","生命财产","生态环境","生物制药","生育","省级","省内","失信","施展才华","施政","十国集团","十九","石漠化","时刻","时限","时有发生","实践","实事","实体店","实效","实质性","使命","始终保持","市场监管","市场调节","市场需求","市场秩序","事故","事关","事件","事权","事务","事业部","事业单位","试验","释放","释放出来","收储","手机","守护","首次","首脑会议","首位","受惠","受益","授权","瘦身","舒心","疏通","属地","属性","树根","数字家庭","双赢","水稻","水价","水运","税费","思潮","思路","四大","四档","四要","四中","诉求","素质","速度","随意","损失","缩小","所得额","所得税","琐事","台独","台海","台湾","抬头","探索","特大","特困","特色产业","特殊教育","梯度","提案","提出","提高质量","提款权","提前完成","提请","提速","体魄","体育","体育产业","体育健儿","天气","天下","田园风光","调查","调动","调减","调优","调整结构","铁腕","停建","停征","通电","通关","通观全局","通用","同级","同心同德","同心协力","同质","同奏","统筹兼顾","痛心","偷排","投产","投入","投向","突发","突破","突破性","土地","团结","团结奋斗","团组织","推进改革","推诿","退耕还林","退化","拖欠","脱实","拓宽","外部环境","外溢","万个","万件","万千瓦","万众一心","网购","为政","违法","违法犯罪","违约","维和","维护者","维稳","伟大旗帜","尾气","委员","卫生","未来","位居","温暖","文化遗产","文明执法","文物","文学","稳定增长","稳是","稳中","稳住","污染物","无党派人士","无穷","无效","五位一体","五要","五中","务实","物流","物质","物质基础","吸引力","希望","下沉","下行","下调","下要","县级","县市","现代文明","现实","现象","限度","线路","相得益彰","相辅相成","相适应","相衔接","香港特别行政区","详查","享誉","向前","向上","向心力","项目管理","消除","消费者","消化","小城镇","小利","小麦","小企业","小时","小微","校际","效能","协作","心头","辛勤","新动力","新高","新格局","新开工","新理念","新能源","新区","新上","新天地","新闻出版","新形势下","新兴","薪酬","信贷","信访工作","信教","信息化","信息网络","信息系统","信用社","形势严峻","形态","兴边","兴起","兴旺发达","幸福","幸福生活","雄厚","休耕","修复","虚假","需氧量","许可证","宣传教育","宣示","选择","选择权","削减","学科建设","学前教育","学杂费","寻常","循环","压倒性","压缩","亚欧","亚太","亚太经合组织","延伸","延续","严而","严格遵守","严控","严厉","严厉打击","严守","严肃","研究生","养老保险","咬定青山不放松","要持","要共担","要活","要实","要素","要稳","要准","业绩","业态","业务","叶茂","一百年","一半","一笔","一场","一二三","一流","一流大学","一人","一视同仁","一体化","一系列","一以贯之","医患","医疗保健","医疗器械","医务人员","医养","医药","依法打击","依法行政","依靠人民","依然","遗留问题","以党","以电代煤","以内","以下几点","以新","亿多","亿人次","艺术创作","议案","议定书","易涝","意义","意愿","因病","因城","引进","饮食","饮水","印记","英才","迎接","迎难而上","营销","赢得","影视","影响力","影子","应保","应急","应用型","硬道理","硬骨头","拥有量","庸政懒","勇往直前","勇于进取","用地","用电","用户","用能","优出","优化结构","优惠","优进","优胜劣汰","优势互补","优异成绩","友好合作","有错必纠","有党","有利于","有所改善","又利","余额","舆论监督","予以","预防","预见性","预决算","预算","预算内","原创性","圆满完成","援助","源头","远海","约定","约法三章","约束","阅读","粤港澳","运动","运营","蕴藏","杂事","再保险","再创新高","在城镇","在思想上","在线","在线教育","早日","造假","增多","增加值","增进","增速","增添","增效","增长极","增值","扎牢","诈骗","债务","占补","战略性","彰显","长城","长江流域","长途电话","长效机制","长远","长征","长征路","长治久安","涨幅","招商引资","招生","招收","找准","哲学","着重","真抓实干","诊疗","阵痛","振兴","征收","整肃","整体","正道","正规化","正确","正确处理","正义","证券化","证照","政策性","政怠政","政府部门","政商","政治经济","之本","之比","之策","之魂","之门","之要","支点","支柱性","知识","织密","直面","职能","职业资格","纸上谈兵","指定","至上","至少","志愿","制度性","制衡","制约","制止","治国","治军","治乱","致力","致贫","智慧","智库","中部","中倒损","中低产田","中国工农红军","中华","中华儿女","中华文化","中起","中外","中西部","中小城市","中性","中医药","种植","众创","众多","重大成果","重大突破","重大意义","重典","重点项目","重度","重视","重托","重心","重要讲话","重在","重中之重","周边国家","周边环境","周密","诸多","诸多矛盾","主板","主办","主产区","主产省","主渠道","主权","主席","主业","助力","助学金","住宅","铸就","筑牢","抓住","专线接入","专业","专业化","专业技能","转化","转为","壮士断腕","状况","追究","准确","卓有成效","资产负债率","资费","资金投入","资源税","资质","资助","自然灾害","自我","自信","自由","自由化","自住","宗教界","宗教事务","综合利用","总量","总体布局","纵容","纵深","纵向","走过","走廊","走势","走向","走样","租房","租赁业","足额","最多","最富","尊崇","做事"],"freq":[6125,3332,3185,2646,2548,2058,1960,1960,1764,1715,1715,1568,1519,1519,1519,1470,1421,1421,1372,1323,1274,1274,1225,1225,1225,1176,1078,980,931,931,931,931,882,882,833,833,833,833,784,784,784,784,784,735,735,735,735,735,735,686,686,686,686,686,686,637,637,637,637,637,637,637,637,637,588,588,588,588,588,588,588,588,588,588,588,588,588,588,539,539,539,539,539,539,539,539,539,490,490,490,490,490,490,490,490,490,490,490,490,490,490,490,490,441,441,441,441,441,441,441,441,441,441,441,441,441,392,392,392,392,392,392,392,392,392,392,392,392,392,392,392,392,392,392,392,392,343,343,343,343,343,343,343,343,343,343,343,343,343,343,343,343,343,343,343,343,343,343,343,343,343,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,294,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,245,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,196,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,147,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,98,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49],"fontFamily":"Segoe UI","fontWeight":"bold","color":"random-light","minSize":0,"weightFactor":0.0293877551020408,"backgroundColor":"grey","gridSize":0,"minRotation":-0.785398163397448,"maxRotation":0.785398163397448,"shuffle":true,"rotateRatio":0.4,"shape":"circle","ellipticity":0.65,"figBase64":null,"hover":null},"evals":[],"jsHooks":[]}</script>

---

如何通过循环来遍历所有年份政府工作报告的链接，留待大家作为思考题。

提示如下：

```r
url <- "http://www.gov.cn/guowuyuan/baogao.htm"
reports <- read_html(url)
links <- reports %>%
  html_nodes(".history_report a") %>%
  html_attr("href") %>%
  str_trim()
head(links)
```

```
## [1] "http://www.gov.cn/premier/2017-03/16/content_5177940.htm"  
## [2] "http://www.gov.cn/premier/2016-03/17/content_5054901.htm"  
## [3] "http://www.gov.cn/guowuyuan/2015-03/16/content_2835101.htm"
## [4] "http://www.gov.cn/guowuyuan/2014-03/14/content_2638989.htm"
## [5] "http://www.gov.cn/premier/2013-03/19/content_2357136.htm"  
## [6] "http://www.gov.cn/premier/2012-03/15/content_2067314.htm"
```

---

## 数据获取、处理与可视化

### 经济学研究的常用数据、世界银行数据可使用两个R包获取：

- [WDI](https://cran.r-project.org/web/packages/WDI/index.html)
- [wbstats](https://cran.r-project.org/web/packages/wbstats/vignettes/Using_the_wbstats_package.html)

### 一个复制Hans Rosling的[Gapminder](https://www.gapminder.org/)软件的动态交互式气泡图

- [Hans Rosling](https://www.ted.com/speakers/hans_rosling)的[TED演讲](https://www.ted.com/talks/hans_rosling_shows_the_best_stats_you_ve_ever_seen)，[中文翻译版](http://open.163.com/movie/2011/12/L/6/M8H1NPQM9_M8H1USJL6.html)
- [R中的复制](https://www.r-bloggers.com/new-r-package-to-access-world-bank-data/)

---

## ggplot系列图形

利用ggplot2及ggthemes、ggsci等包，可便捷产生符合特定杂志风格的图形。

### 常用ggplot系列可视化包

- ggplot2
- ggthemes
- ggsci
- ggcorrplot
- ……

---

### ggplot2 原始风格

```r
library(ggplot2)
ggplot(iris) + 
  geom_boxplot(aes(x = Species, y = Sepal.Length, fill = Species))
```

---
### *The Economist* 风格图形

```r
ggplot(iris) + 
  geom_boxplot(aes(x = Species, y = Sepal.Length, fill = Species)) + 
  ggthemes::theme_economist()
```

---

### *The Wallstreet Journal* 风格图形

```r
ggplot(iris) + 
  geom_boxplot(aes(x = Species, y = Sepal.Length, fill = Species)) + 
  ggthemes::theme_wsj()
```

<img src="R4Eco201707_files/figure-html/wsj style-1.png" style="display: block; margin: auto;" />
---

### Stata风格图形

```r
ggplot(iris) + 
  geom_boxplot(aes(x = Species, y = Sepal.Length, fill = Species)) + 
  ggthemes::theme_stata()
```

---

### *Nature* 风格

```r
library("ggsci")
library("ggplot2")
library("gridExtra")
data("diamonds")
p1 = ggplot(subset(diamonds, carat >= 2.2),
       aes(x = table, y = price, colour = cut)) +
  geom_point(alpha = 0.7) +
  geom_smooth(method = "loess", alpha = 0.05, size = 1, span = 1) +
  theme_bw()
p2 = ggplot(subset(diamonds, carat > 2.2 & depth > 55 & depth < 70),
       aes(x = depth, fill = cut)) +
  geom_histogram(colour = "black", binwidth = 1, position = "dodge") +
  theme_bw()
```

```r
p1_npg = p1 + scale_color_npg()
p2_npg = p2 + scale_fill_npg()
grid.arrange(p1_npg, p2_npg, ncol = 2)
```

---

## Nature 风格

---

### *Lancet* 风格

```r
p1_lancet = p1 + scale_color_lancet()
p2_lancet = p2 + scale_fill_lancet()
grid.arrange(p1_lancet, p2_lancet, ncol = 2)
```

<img src="R4Eco201707_files/figure-html/lancent style-1.png" style="display: block; margin: auto;" />
---
##  更多的R可视化图例

- [RStudio图库](https://shiny.rstudio.com/gallery/)

- [ggplot2图库](http://www.r-graph-gallery.com/portfolio/ggplot2-package/)

- [ggthemes示例](https://cran.r-project.org/web/packages/ggthemes/vignettes/ggthemes.html)

- [ggsci示例](https://cran.r-project.org/web/packages/ggsci/vignettes/ggsci.html)

---

# 用于撰写学术报告

- rmarkdown: html 格式报告
- xaringan: html 格式幻灯片
- rticles: AER 等经济学类顶级刊物LaTeX模板
- stargazer: 生成LaTeX表格

---

# 结论与建议

- R 对于经济学教学与研究来说，是一个“无价”而高效的工具。
   - 无价的两层含义
      - 本身是免费的
      - 作用是巨大的

- 通过RStudio 等IDE（集成开发环境），R已形成一个良好的数据分析生态
  - 可导入SAS、STATA、SPSS等常见格式的数据
  - 可用来生成动态化、交互式报告
  - 几乎可直接用来撰写论文与书籍
  - 有强大的社区支持

- R 自身的学习周期较长，不易上手
 - 这是包括Python等开源软件存在的一种通行问题
 - 可能需要进一步加强基础课程建设来加以解决

---

# 结论与建议

- 尽早让学生接触数据分析的基本流程
  - 在没有概率论支撑的情况下就可引入数据分析的基本理念，培养面向数据的分析思维

- 尽早、全面地引入计算工具，深度参与统计教学
  - 对于科研型院校，可适当强调程序思维和编程操作的基础性地位
  - 要将相关的数据分析的计算机技术当成基础数学三大模块（微积分、线性代数、概率论与数理统计）并列的模块加以重视

- 应更加突出可视化在数据分析中的作用

- 建立统一的、开放的、可编辑的数据展示与教学安全平台，对于节约教师的精力有很大的作用(GitHub可以作为一个有效的平台)。

---

# 常用资源

- 计量经济学中的常用 R 包索引：<https://cran.r-project.org/web/views/Econometrics.html>

- 用R做计量分析网站：<https://econometricswithr.wordpress.com/>

- Using R for Introductory Econometrics(Wooldridge 计量经济学导论配套R语言网站): <http://www.urfie.net/>

- bookdown官方网站：<https://bookdown.org/home/>

- *R for Data Science* 在线版本：<http://r4ds.had.co.nz/>

---

class: center, middle

# 谢谢观看！

### [吕小康](https://xkdog.github.io/)     
### xkdog@126.com

### 南开大学周恩来政府管理学院

本幻灯片由[谢益辉](https://yihui.name/)的 R 包 [**xaringan**](https://github.com/yihui/xaringan) 生成。原始文档可从以下链接下载：

<https://github.com/xkdog/StatsUsingR>

简略版可从以下网址在线观看（图片未能正确显示）：

<https://github.com/xkdog/StatsUsingR/blob/master/R4Eco201707.Rmd>

<http://rpubs.com/xkdog/r4eco2017>