我的学习笔记之数据的描述性分析（R语言）

一、定性数据描述的数值方法

1.某超市为研究不同类型饮料的市场销售情况，随机地调取了50名顾客购买饮料的数据，其性别以及购买饮料的类型如表所示。试计算购买人群（男、女）以及不同种类饮料的频数。

顾客性别及购买的饮料类型

首先将数据存放在drink.data中，然后用scan()函数将数据读出，再使用table()计算出频数.

Lst<-scan("drink.data",what=list(sex="",type=""))   #scan()返回值为列表
Ta<-table(Lst);Ta     # table()返回值为矩阵

#运行结果
sex Fruitjuice greentea others soda water
  F          1        1      6   11     7
  M          2        7      0    9     6

下面计算不同人群（男、女）喝饮料的频数，或者喝不同种类饮料的频数，称为边缘频数。

#margin.table()函数计算边缘频数，其格式为：margin.table(x,margin=NULL)
(sex<-margin.table(Ta,1)) 
(drink<-margin.table(Ta,2))
# Ta为矩阵，1表示计算行边缘频数， 2表示计算列边缘频数
#可以用prop.table()函数计算相对频率，格式为：prop.table(x,margin=NULL)
prop.table(Ta)
prop.table(drink)

#运行结果
> (sex<-margin.table(Ta,1)) 
sex
 F  M 
26 24 
> (drink<-margin.table(Ta,2))
type
Fruitjuice   greentea     others       soda      water 
         3          8          6         20         13 
> prop.table(Ta)
   type
sex Fruitjuice greentea others soda water
  F       0.02     0.02   0.12 0.22  0.14
  M       0.04     0.14   0.00 0.18  0.12
> prop.table(drink)
type
Fruitjuice   greentea     others       soda      water 
      0.06       0.16       0.12       0.40       0.26 

2.为研究新生儿出生时的体重，在某妇产医院随机地选取了100个新生儿，其重量如表（birth.data）. 试对数据进行分组，并计算每组的频数.

#scan()函数读取数据，cut()函数将数据分组table()函数计算频数，
X<-scan("birth.data")/1000   #读取数据并转换为千克
table(cut(X,breaks=10))      #分成10组
b<-c(-Inf,1,1.5,2,2.5,3,3.5,4,4.5,5,Inf)
table(cut(X,breaks=b))    #按照指定的区间长度0.5分成10组

#运行结果
> X<-scan("birth.data")/1000   #读取数据并转换为千克
Read 100 items
> table(cut(X,breaks=10))      #分成10组

(0.956,1.35]  (1.35,1.73]  (1.73,2.12]  (2.12,2.51]   (2.51,2.9]   (2.9,3.28] 
           5            0            5            5           19           14 
 (3.28,3.67]  (3.67,4.06]  (4.06,4.44]  (4.44,4.83] 
          18           16           13            5 
> b<-c(-Inf,1,1.5,2,2.5,3,3.5,4,4.5,5,Inf)
> table(cut(X,breaks=b))    #按照指定的区间长度0.5分成10组

(-Inf,1]  (1,1.5]  (1.5,2]  (2,2.5]  (2.5,3]  (3,3.5]  (3.5,4]  (4,4.5]  (4.5,5] 
       5        0        5        5       25       13       29       13        5 
(5, Inf] 
       0 

3.“Zagat饭店评论”是一个提供世界各地饭店数据的服务机构，如：饭店的质量等级（定性变量）、典型食品的价格（定量变量）。一个样本包含了洛杉矶地区300家饭店的数据（Restaurant.csv）。

rc<-read.csv("Restaurant.csv")  #从CSV文件读取数据
rc$食品价格<-cut(rc$食品价格,breaks=c(9,19,29,39,49))  #定量数据转换成分组数据
(Ta<-table(rc))       #计算交叉列表

#也可用prop.table()计算食品价格关于各种质量等级的频率
round(prop.table(Ta,1)*10000)/100   #Round()是四舍五入

二、描述定性数据的图形法

1.条形图（描述频数）

# 给出相应每一类的频数，长方形的高度（或宽度）与频数成比例. 
# barplot(height,width=1,space=NULL,…)
a<-c(1,2,3,4)
r<-barplot(a)

4.将例1中饮料数据画条形图：

Lst<-scan("drink.data",what=list(sex="",type=""))
Ta<-table(Lst)
r<-barplot(Ta,col=c("lightblue","mistyrose","lightcyan","lavender","cornsilk"),legend.text=c("男","女"),xlab="饮料类型",ylab="频数")

#法2
r<-barplot(Ta,beside=T,col=c("lightblue","mistyrose"),legend.text=c("男","女"),args.legend=list(x=3,y=8),xlab="饮料类型",ylab="频数")

2.饼图

5.将例1中的饮料数据绘制饼图。

# pie(x,labels=names(x),edges=200,radius=0.8,…)
r<-pie(Ta,radius=0.9)  #绘制简单饼图，下边语句可以绘制详细饼图
Name<-paste0(names(Ta),'',prop.table(Ta)*100,'%')
pie(Ta,radius=0.9,labels=Name,col=rainbow(length(Ta)),font=2,cex=1.2)
Prop.table()*100 #计算每种类型数据的百分比
Font=2 #表示使用黑体
Cex=2 #表示字体大小似乎原来的2倍

3. Pareto图

按条形的高度顺序排列的条形图. Pareto图通常在过程和质量控制中使用，容易辨别哪些因素最重要.。
R没有绘制Pareto图的函数，但是我们可以调用自己编写的函数。
6.绘制例1的饮料数据的Pareto图。

Lst<-scan("drink.data",what=list(sex="",type=""))
Ta<-table(Lst)
Ta
pareto_chart<-function(height,col=NULL,xlab=NULL){
  height<-sort(height,decreasing = TURE)
  s<-sum(height)
  par(mai=c(0.9,0.9,0.9,0.9))
  r<-barplot(height,col=col,xlab=xlab,ylab="数量",ylim=c(0,s+s/20))
  text(r,height/2,format(height))
  cum<-cumsum(height)
  lines(r,cum,type="o",lwd=2,col='red')
  box()
  axis(1,at=r,labels=FALSE)
  at<-seq(0,s,by=s/5)
  label<-paste0(seq(0,100,bu=20),'%')
  axis(4,at=at,label=label)
  abline(h=at,lty=2)
  mtext("累积百分比",4,2.5)
}
pareto_chart(Ta,col=c("lightblue","mistyrose","lightcyan","lavender","cornsilk"),xlab="饮料类型")

三、描述定量数据的图形法

直方图

直方图是展示连续数据分布的常用工具，用来估计数据的分布. 用hist()函数.
7.绘出例2中新生儿体重数据的直方图。

# hist(x,break=“Sturges”,freq=NULL,probability=…)
X<-scan("birth.data")/1000
hist(X,col="lightblue",border="red",labels=TRUE,ylim=c(0,34))   
#频数直方图，增加了直方图和外框的颜色，以及相应的频数
hist(X,freq=FALSE)
lines(density(X),col="blue",lwd=2)
x<-seq(from=.5,to=5,by=.1)
lines(x,dnorm(x,mean(X),sd(X)),col="red",lwd=2)
#密度直方图，并用lines()函数绘出数据的密度曲线和正态分布密度曲线

四、常用的数值描述性度量

1.集中趋势的度量

众数：通过计算频数的最大值，可用max()和which.max()，分别返回数值向量的最大值，以及哪一个分量是最大值.
8.计算例1中饮料数据的众数.

Lst<-scan("drink.data",what=list(sex="",type=""))
Ta<-table(Lst)
(drink<-margin.table(Ta,2))
max(drink)
which.max(drink)

中位数
在R中用median()计算中位数，格式为:median(x,na.rm=FALSE)
9.计算例2中新生儿体重数据的中位数.

X<-scan("birth.data") 
median(X)

算术平均

10.计算例2中新生儿体重数据的均值.

#mean(x,trim=0,na.rm=FALSE,…)
X<-scan("birth.data") 
mean(X)

加权平均数

对分组数据计算平均数不合适，需要加权平均，在计算时需要先将每组变量分别乘以各自的频数，然后再除以总次数（或总频数）.
11.某幼儿园共有学童458名，其中3至6岁儿童分别是90、130、120、118名. 计算该幼儿园儿童的平均年龄.

#weighted.mean(x,w,na.rm=FALSE,…)
x<-3:6; f<-c(90,130,120,118)
weighted.mean(x,w=f)

几何平均数

R没有计算几何平均数的函数，但可以通过对算术平均取对数和加权平均取对数.
12.某市2010-2014年GDP的增长率分别是12%、8%、14%、16%和13%. 计算该市5年的平均增长率.

x<-c(12,8,14,16,13)/100
m<-mean(log(1+x))
x_bar<-exp(m);x_bar

分组数据
13.某市2010-2014年GDP的增长率分别是12%、8%、14%、16%和13%. 计算该市5年的平均增长率.

x<-c(12,8,14,16,13)/100
m<-mean(log(1+x))
x_bar<-exp(m);x_bar

结果为1.027127，即平均增长率为2.71%

调和平均数

是各观测数据的倒数的平均数的倒数.
对于分组数据，是个分组数据倒数的加权平均数的倒数.
与几何平均数一样，可借助于算术平均数或加权平均数，计算调和平均数.
14. 某菜市场3种蔬菜的价格分别是5元/500g、4元/500克和4.5元/500克，如果消费者每种蔬菜各买3元，计算该消费者购买蔬菜的平均价格.
计算平均价格的本质，实际上是计算这3组数据的调和平均数.

x<-c(5,4,4.5); (x_bar<-1/mean(1/x))

15.某水果市场某日的苹果、梨和橘子的价格分别为3.6元/500克、4.2元/500克和2.8元/500克，成交金额分别是2.88万、5.04万和4.20万元，计算这3种水果的平均价格.

x<-c(3.6,4.2,2.8) 
f<-c(2.88,5.04,4.20)*10000
(x_bar<-1/weighted.mean(1/x,w=f))

2.离散程度的度量

极差（全距）

16.计算例2中新生儿体重数据的极差.

#可用max()和min()函数计算极差，或者用range()函数。
X<-scan("birth.data")
max(X)-min(X)
R<-range(X); R[2]-R[1]

分位数与四分位差
17.计算例2中心生儿体重数据的四分位差，以及0%、20%、40%、60%、80%和100%的百分位数.

# quantile()函数计算分位数（或百分位数），格式为：
#quantile(x,probs=seq(0,1,0.25),na.rm=FALSE,names=TRUE,type=7,…)
X<-scan("birth.data")
Q<-quantile(X,p=c(0.25,0.75),names=F)
IQR<-Q[2]-Q[1]; IQR     #四分位差
quantile(X,p=seq(0,1,by=0.2))  # 百分位数

18.计算例2中新生儿体重数据的方差和标准差.
方差与标准差

# var()函数计算方差，sd()函数计算标准差，
# var(x,y=NULL,na.rm=FALSE,use)
# sd(x,na.rm=FALSE)
X<-scan("birth.data")
var(X); sd(X)

3. 分布形态的度量

偏度系数：描述数据的对称性。
关于均值对称的数据其偏度系数为0，右（左）侧更分散的数据偏度系数为正（负）。
R中没有计算样本偏度系数的基本函数，可编写程序skew.R

skew<-function(x){
  mu<-mean(x);n<-length(x)
  m2<-sum((x-mu)^2)/n
  m3<-sum((x-mu)^3)/n
  m3/sqrt(m2^3)
}

19.计算例2中新生儿体重数据的偏度系数.

X<-scan("birth.data") 
source("skew.R") 
skew(X)

结果为-0.6558841，说明数据左侧更分散.

峰度系数：正态分布的偏离程度.
峰度系数是反映分布峰的尖峭程度的重要指标.
当峰度系数大于0时，两侧极端数据较多；
当峰度系数小于0时，两侧极端数据较少.
R中没有计算样本峰度系数的基本函数，可编写程序kurt.R

kurt<-function(x){
  mu<-mean(x);n<-length(x)
  m2<-sum((x-mu)^2)/n;m4<-sum((x-mu)^4)/n
  m4/m2^2-3
}

20.计算例2中新生儿体重数据的峰度系数.

X<-scan("birth.data");source("kurt.R");kurt(X)

五、检验异常值的方法

异常值通常归咎于以下原因：
1）观测、记录或录入计算机时不正确；
2）测量值来自不同的总体；
3）观测值是正确的，但代表一个稀有（或偶然）的事件.
介绍两种检测异常值的方法—标准分数法和图形法.

1.标准分数法

变量值与其平均数的差除以标准差的值称为标准分数，或称为Z得分。
当Z_i的绝对值大于某个数值时，可以将第i个样本看成异常值.
判别数据分布的两个法则：
<1>经验法则：若数据集近似于丘形对称分布，则
(1)大约有68%的测量值位于平均值的1个标准差的范围内；
(2)大约有95%的测量值位于平均值的2个标准差的范围内；
(3)几乎所有的测量值位于平均值的3个标准差的范围内.
<2>Chebyshev(切比雪夫)法则：对任意数据集，无论数据的频数分布是什么形状的，则
(1)可能有很少的测量值位于平均值的1个标准差的范围内；
(2)至少有75%的测量值位于平均值的2个标准差的范围内；
(3)至少有8/9的测量值位于平均值的3个标准差的范围内；
(4)对于任意大于1的数k，至少有1-1/k^2的测量值落在平均值的k个标准差的范围内.
通过z得分和这两个法则，可以判断哪些样本是异常的.
21.考虑新生儿体重数据中哪些新生儿的体重属于异常值.

X<-scan("birth.data");hist(X) #数据不对称，用切比雪夫法则
X<-scan("birth.data")
names(X)<-1:length(X)
Xbar<-mean(X)
S<-sd(X)
Z<-(X-Xbar)/S
X[abs(Z)>3]

统计学上一般认为Z得分绝对值大于3，是异常值.

#运行结果
named numeric(0)

2.箱线图法（箱形图）

从四分位数的角度来描述数据的分布，它由上边缘线、上四分位数、中位数 、下四分位数和下边缘线组成。盒子的上限边缘线为上、下四分位数，盒子中间的线是中位数。
22.使用箱线图的方法考查例3.2的新生儿体重数据的异常值.

#boxplot()绘制箱线图.
X<-scan("birth.data")
boxplot(X)

还可以用boxplot()函数画出多组数据的箱线图.
23.已知由两种方法得到的数据如表，试画出两组数据的箱线图.

A<-c(79.98, 80.04, 80.02, 80.04, 80.03, 80.03, 80.04, 79.97, 80.05, 80.03, 80.02, 80.00, 80.02)
B<-c(80.02, 79.94, 79.98, 79.97, 79.97, 80.03, 79.95, 79.97 )
boxplot(A,B,names=c('A','B'),col=c(2,3))