vlambda博客
学习文章列表

十年树木，百年树人！即使你从小就学求平均数，但是现在数据变了，工具变了，所以还要再学怎么求平均数！

---

在上一篇中，我简单地给大家简单介绍了R语言。通过A股的数据，我们介绍了如何利用R语言读取Excel，CVS等格式的文件。在此基础上，我们可以查看数据的特性，例如每个省有多少上市公司。最后还介绍了dplyr，我们可以通过它修改变量名称，删除变量以及合并表格等等。

 

今天我继续为大家介绍如何利用R语言进行描述性统计（Descriptive statistics）。所谓描述性统计，就是利用客观现象数据，通过图表形式，反应出规律性数量特征。通俗地讲，就是通过平均值、最大值、最小值、标准差等参数，描述数据的特征。

 

这个事情听起来简单，谁不会求平均数呢？但是事实上也没有那么简单。因为有些变量为连续型(Continuous)，比方说身高、体重、气温等等。计算这些变量的平均值不难。不过有些变量不是连续型的，而是类别型（Categorical）变量，比方说性别或者上市公司的所在的省份。对类别型变量求平均值就没有意义。所以我们还需要更加适合的方法来描述类别型变量。

1. 连续型变量的描述性统计

我们先从连续型变量开始说起。利用R语言，我们可以很轻松地算出平均值、最大值、最小值等。

 

全部数值型变量的平均值：mean(data)

特定数值变量的平均值：mean(data$variable)

中位数：median(data$variable)

方差：var(data$variable)

标准差：sd(data$variable)

最小值：min(data$variable)

最大值：max(data$variable)

极差：range(data$variable)

百分位：quantile(data$variable)

 

我们还是用A股数据(stock)的例子，来看一下如何求平均值。这个数据集包括3519只股票的6个参数：股票代码（Code）、股票名称（Name）、收盘价（Close_price）、上市公司所在省份（Location）、行业（Industry）和上市日期（List_date），而且没有NA。

比方说，我们想知道这些股票价格的平均值，就可以这样计算。

>mean(stock$Close_price)

[1] 14.16567

 

这就是说A股平均股价在14元左右，所以如果你有1500元，就可以买股票了。

 

股价的标准差

> sd(stock$Close_price)

[1] 27.06492

 

虽然标准差是27， 但是股票价格显然不能为负数。R或者Python这些工具简化了计算，但是我们还要理解数据的含义。不要被自己算出来的数给骗了。

 

当然还有一个简单的方法来做描述性统计，就是用summary()函数：

> summary(stock$Close_price)

  Min.  1st Qu.  Median  Mean  3rd Qu.  Max.

  0.00    5.08    8.58   14.17   15.38  1052.80

 

可见，A股股价最高的是1052元的茅台。有一半股票（也就是1000多只）的价格在8.58元以下，75%的股票价格在15.38元以下。

 

这里还有一个秘籍！我们还可以调用更好的统计库skimr来进行计算。它是一个很轻巧的统计分析包，而且可以和我们之前提到的tidyverse包无缝对接。

library(skimr)

skim(stock)

 

部分输出结果如下：

 

这个方法也可以求出平均值、标准差、以及百分位等等。和我们得出的结果是一致的。对于数值型变量，skimr会画出分布的直方图。最后一列是skimr画出的直方图，但是这个看起来似乎有点问题。

 

为了更好的展示skimr，我们创建两个数组，x和y。x是1到10000，y是均值为1的10000个数，那么我们看看用skim会怎样？

x = seq(10000)

y = rnorm(10000,mean=1)

df = data.frame(x,y)

skim(df)

在后一列，我们可以清晰地看出skim()函数画出了x和y的分布。这就是skimr的特点。

 

上面是整体数据的情况。假如我们想按省份来看股票的平均值、标准差、最大值、最小值，怎么办呢？我们可以用到tapply()函数

 

mean<-tapply(stock$Close_price,stock$Location, mean)

sd<-tapply(stock$Close_price,stock$Location, sd)

med<-tapply(stock$Close_price,stock$Location, median)

min<-tapply(stock$Close_price, stock$Location,min)

max<-tapply(stock$Close_price,stock$Location, max)

cbind(mean, sd, med, min, max)

 

得出的结果如下：

 

我们可以看出，甘肃、宁夏、青海三个省区上市公司的平均股价都比较低，不到7元。贵州上市公司股价的平均值最高，将近46元。这显然是被千元股茅台拉起来的。另外，福建和北京上市公司的股价也比较高。

 

同样，我们可以分行业看上市公司股票的价格。采掘、钢铁、房地产这些行业的平均股票价格都在7元以下。食品饮料行业的股票价格在茅台的带动下，平均值依旧很高。而且食品饮料行业股票价格的中位数在11.55，也不算低。

 

以上是对连续型变量的分析。接下来，我们来看如何对类别型变量进行分析。

2. 类别型变量的描述性统计

 

一般来说，对于类别型变量，我们可以先数数个数，然后看它在行和列中的百分比。对于A股数据集，这个问题就转化为怎么看每一年，每个省有多少公司上市？

为了回答这个问题，我们先把上市日期（List_date）转化为上市年份（List_year）。这里需要调用lubridate库。这个库提供了一个简单处理年月的方法，也就是可以通过使用year()，month()函数，直接将日期中的年月提取出来。

 

我们安装并调用这个库

install.packages("lubridate")

library(lubridate)

 

提取年份

year<-year(stock$List_date)

 

我们可以用dplyr库中的mutate()函数，在stock数据集中加入List_year列

stock<-mutate(stock,List_year=year(stock$List_date))

 

这时stock数据集中就出现了List_year这一个列

 

这样就好了，我们来看一下每年每个省有多少公司上市。这里用到table()函数

tbl <- table(stock$Location,stock$List_year)    

 

我们可以看出，90年代初的时候，上市公司主要集中在广东和上海地区。慢慢地向全国扩散。2010年前后，北京、江苏、浙江的上市公司逐渐多了起来。

 

另外的问题是行占比和列占比。也就是说，某省某一年上市公司的数量占这个省总体上市公司数量的比例是多少？在某一年，某省上市公司的数量占全国上市公司数量的百分比是多少？

 

要回答这两个问题，我们用到下面的语句：

prop.table(tbl, 1)

prop.table(tbl, 2)

 

下面的输出是行占比。2014年安徽省上市公司的数量占安徽省总上市公司数量的2.91%，2015年安徽省上市公司的数量占安徽省总上市公司数量的7.77%，依此类推。内蒙古已经好多年没有公司上市了。

 

接下来我们来看一下列占比。2014年，安徽省上市公司的数量，占当年全国上市公司总数量的2.42%。而2014年北京市上市公司的数量，占当年全国上市公司总数量的14.5%，广东省占到了19.3%，浙江省占到了13.7%。

 

这个问题还可以更有意思一点。我们可以看一下每年，每个省，每个行业有多少公司上市。

我们可以这样解决这个问题：

xtabs(~Industry+Location+List_year,data=stock)

 

我们可以看出1990年，广东省有一个公共事业公司上市。

这是哪家公司呢？我们用filter()函数看一下

filter(stock,Location=="广东" & List_year=="1990")

 

这个公司是世纪星源。由于上市时间长，历史上用过的名字都有一大串：深原野A->深星源A->世纪星源->ST星源->*ST星源->ST星源->GST星源->ST星源。 

 

我们再来做一个有趣的事情。为了这个事情，我们先把数据做一些调整。我们把变量Location分为3类。上市公司位于北京的是1，上海的是2，其他的是3。用ifelse()定义新的变量New_location。

stock$New_location<-ifelse(stock$<5,1,

ifelse((stock$log_price>5 & stock$log_price<7),2,3))

 

这样，我们就得到了New_location:

 

下面我们要用到R语言中的“阿森纳”库:arsenal。它主要用于大规模统计计算，是R语言中产生报表的工具。它的主要函数是tableby()。我们就来看看如何使用arsenal和tableby()。

 

当然首先你要安装和调用arsenal

install.packages("arsenal")

library(arsenal)

 

接下来，我们用tableby()函数画出表1

tab1 <- tableby( ~ Close_price +Industry + New_location, data=stock)

summary(tab1, text=T)

 

表1

这个表格计算了股价的平均值和标准差，各个行业上市公司的数量和占比，以及新变量New_location的平均值和极差（虽然毫无意义）。

 

我们如何改变这个表格？我们可以通过tableby.control做一系列操作。

my_controls <- tableby.control(

 total = T,   # 如果不看整体，可以让 total=F                  

 test=F,      # 不显示p值

 numeric.stats = c("meansd", "medianq1q3","range", "Nmiss2"),  #选择连续型变量要看的参数

 cat.stats = c("countpct", "Nmiss2"),   # 选择类别型变量要看的参数

  stats.labels= list(

   meansd = "Mean (SD)",

   medianq1q3 = "Median (Q1, Q3)",

   range = "Min - Max",

   Nmiss2 = "Missing",

   countpct = "N (%)"))

 

为了让显示N (%)，我们把New_location变为因子(factor)

stock_df<-stock %>%

  mutate(New_location=factor(New_location,labels=c("北京","上海","其他")))

 

接下来，我们就可以获得表2。

tab2 <- tableby( ~ Close_price +Industry + New_location, data=stock_df, control=my_controls)

summary(tab2, title = "DescriptiveStatistics: stock Data", text=T)

 

表2

表2和表1差不多，但是仔细看来也有区别。对于连续型变量，我们计算了平均值、标准差、中位数等等参数。对于类别型变量，我们计算了个数和百分比。这相对于表1是一种进步。

 

我们还可以更上一层楼。如果我们按照New_location分类，这个表就会更好一些了。

tab3 <- tableby(New_location ~Close_price + Industry , data=stock_df, control=my_controls)

summary(tab3, title = "DescriptiveStatistics: stock Data", text=T)

表3

 

表3就一目了然了。我们可以看出北京、上海、其他地区以及全国上市公司的股价平均值、标准差等参数。每个行业在不同地区的分布也很清楚了。

好了，先写到这里。欢迎挑毛病！