文章来源于:炼数成金;摘自《数据分析:R语言实战》
第5章 数据的描述性分析
通过前面两章的学习,我们知道,数据收集是取得统计数据的过程,数据预处理是将数据中的问题清理干净,那么接下来的步骤就是统计分析了。数据分析是通过统计方法研究数据的过程,所用的方法分为描述性统计和统计推断两部分。描述性统计用编制图表、计算统计量等形式对数据进行加工处理和显示,进而综合、概括和分析,得出反映所研究数据的一般性特征——数字特征、分布状态等。
描述性统计分析是统计分析的第一步,做好这一步是下面进行正确统计推断的先决条件。
本章将分别从统计量和图形的角度描述样本的分布特征,重点介绍如何运用R中的函数完成。
5.1 R内置的分布
分布是描述一个样本数据最核心、最重要的方式。R内嵌了很多常用的统计分布,提供了四类函数:概率密度函数(density)、累积分布函数(probability)、分位数(quantile)和伪随机数(random)。在R中分别用d、p、q、r表示这4个项目,后面接分布的英文名称或缩写。
例如正态分布的4个函数分别表示为dnorm、pnorm、qnorm及rnorm。表5.1列出了R中的分布函数,每个分布函数都有各自的参数,有些有默认值(如正态分布默认为N(0,1)),有些则没有。
最常用的函数一般都来自于统计分析程序包stats,它不需要手动安装,打开R的同时其就已经自动加载。另外,有些分布较为复杂但在统计上应用较广,例如Wishart分布和Pareto分布,使用它们需要加载新的程序包。
表5.1 R内置的分布及对应函数
分 布 | R函数 | 参数及默认值 | 所属程序包 |
贝塔Beta | _beta | shape1, shape2, ncp=0 | stats |
二项Binomial | _binom | size, prob | |
柯西Cauchy | _cauchy | location=0, scale=1 | |
卡方Chi-square(x2) | _chisq | df, ncp | |
指数Exponential | _exp | rate | |
F分布Fisher-Snedecor | _f | df1, df2, ncp | |
伽马Gamma | _gamma | shape, scale=1 | |
几何Geometric | _geom | prob | |
超几何Hypergeometric | _hyper | m, n, k | |
对数正态Lognormal | _lnorm | meanlog=0, sdlog=1 | |
逻辑斯蒂Logistic | _logis | location=0, scale=1 | |
负二项Negative binomial | _nbinom | size, prob | |
多项式Multinomial | _multinom | size, prob | |
正态Normal | _norm | mean=0, sd=1 | |
泊松Poisson | _pois | lambda | |
学生Student’s t | _t | df | |
均匀Uniform | _unif | min=0, max=1 | |
威布尔Weibull | _weibull | shape, scale=1 | |
威尔考克森Wilcoxon | _wilcox | m, n | |
帕累托Pareto | _pareto | shape, scale | actuar (含有许多逆函数) |
布尔Burr | _burr | shape1, shape2, rate=1 (scale=1/rate) | |
逆指数Inverse Exponential | _invexp | rate | |
狄利克雷Dirichlet | _dirichlet | alpha | MCMCpack (蒙特卡罗,马氏链) |
威沙特Wishart | _wish | v, S | |
逆威沙特Inverse Wishart | _iwish | v, S | |
广义极值 Generalized Extreme Value | _gev | xi, mu, sigma | evir (极值统计分析) |
广义帕累托 Generalized Pareto | _gpd | xi=1, mu=0, sigma=1 | |
多元正态 Multivariate Normal | _mvnorm | mean, sigma | mvtnorm |
多元t分布Multivariate-t | _mvt | sigma=diag(2), df=1 | mvtnorm |
5.2 集中趋势的分析
5.2.1 集中趋势的测度
样本观测值来自于总体,其中含有总体各方面的信息,乍看这些信息有时显得杂乱无章,需要对样本进行加工,通过探索性分析、计算统计量来提炼有用的信息。
描述数据分布特征的统计量可分为两类:一类表示数量的中心位置;另一类表示数量的变异程度(或称离散程度)。两者相互补充,共同反映数据的全貌。
描述统计分布集中趋势的指标主要是平均数、中位数、众数,也称为“平均指标”。这些指标的主要作用包括:
反映总体各单位变量分布的集中趋势和一般水平;
便于比较同类现象在不同单位之间的水平;
便于比较同类现象在不同时期的发展变化趋势或规律;
用于分析现象之间的依存关系。
(1)均值(平均数)
均值是最常用于度量数据平均水平的统计量,记为 ,定义为
(2)众数
一组数据中出现次数最多的观测值叫做众数,用 表示。众数测度数据的集中性趋势,一般在数据量较大的情况下,众数比较有意义。众数不受极端值的影响,一组数据分布的最高位置对应的数值即众数。
如果数据没有明显的集中趋势,那么众数可能不存在;如果有两个最高峰点,那么这组数据就有两个众数。
(3)中位数
中位数简单来讲就是数据排序位于中间位置的值,记为 ,即
中位数描述数据中心位置,对于对称分布的数据,均值接近中位数;偏态分布是指频数分布不对称,集中位置偏向一侧,数据的均值则与中位数不同。
中位数是描述分析中重要的统计量,与众数类似,它的显著特点在于不受异常值的影响,具有稳健性。
平均数、众数和中位数三者一起可以对数据的分布形式做简单描述。
如果数据具有单一众数且分布对称,那么均值、众数和中位数相等,即 。
若集中位置偏向数值小的一侧,称为正偏态分布,也称右偏分布(可以简单地记为分布的“尾巴”在右边)。此时数据存在极大值,必然拉动均值向极大值一方靠拢,因此 ,如图5.1(b)所示。
如果集中位置偏向数值大的一侧,称为负偏态分布,也叫左偏分布,即 ,如图5.1(a)所示。
(4)分位数
中位数是将数据等分后中间位置的点,与之类似的还有四分位点、十分位点和百分位点。常用的是数据排序后处于25%和75%位置上的值,即四分位点。
5.2.2 R语言实现
R中的函数mean()用于计算样本均值,其调用格式如下:
mean(x, trim = 0, na.rm = FALSE, ...)
内部的参数x是计算对象,可以为向量、矩阵、数据或数据框;trim用于设置计算前去掉与均值差别较大数据的比例,默认为0,表示计算中包括全部数据;na.rm默认为FALSE,表示不允许数据有缺失。
我们使用2007.01-2012.09时间段中国人寿股票价格的样本作为实例,具体说明用R如何做描述统计。这一数据集读入后包含两列,第一列为交易日期,第二列为当日收盘价:
> data=read.csv("d:/data/中国人寿股价.csv")
>attach(data)
The following object is masked from data (position 3):
Clsprc, Trddt
>price=Clsprc[!is.na(Clsprc)] #获得价格向量,同时去掉缺失数据
>mean(price)
[1] 27.83326
有时计算均值时,我们希望赋予各变量的权重是不一样的,也就是要计算加权平均数,这时用到的指令是weighted.mean(x,w,…),其中x是样本数据,w是与x长度相等的数值型向量,用于指出赋予样本各数据的权重值。
中位数通过函数median()计算,该函数与mean()的使用方法相同。由计算结果发现,均值略大于中位数,可以初步判断,中国人寿股价样本呈右偏分布。
>median(price)
[1] 24.02
分位数通过quantile()计算,x是数值向量,probs给出相应的百分位数,默认值是(0,0.25,0.5,0.75,1),na.rm=FALSE表示不允许包含缺失数据。
默认条件下,可以直接计算出“五数”:最小值、25%的四分位数、中位数、75%的四分位数和最大值。函数fivenum()用来计算五数。
quantile(x, probs = seq(0, 1, 0.25), na.rm = FALSE,names = TRUE, type = 7, ...)
>quantile(price)
0% 25% 50% 75% 100%
14.740 19.285 24.020 31.160 75.080
>fivenum(price)
[1] 14.74 19.27 24.02 31.16 75.08
>min(price);max(price)
[1] 14.74
[1] 75.08
这样一个一个地输入函数非常不方便,R提供了总体描述的函数。函数summary()可以计算出一组数据的五数和均值。
>summary(price)
Min.1st Qu.Median Mean3rd Qu. Max.
14.74 19.28 24.02 27.83 31.16 75.08
另外,在R中没有直接的函数用来计算众数。对于离散型样本数据,可以简单地手动计算,使用指令which.max(table(x))。但像本例中这样的连续型样本,众数的定义是其分布的密度函数峰值对应的取值,这就需要我们自己写一个函数(function)来计算,然而出于描述数据分布形态的目的,我们大可以直接画出直方图来观察众数的大概位置。
5.3 离散趋势的分析
集中趋势只是数据分布的一个概括性度量,它反映的是各变量值向中心值聚集的程度。而变量的分散程度如何度量呢?分散程度,即各观测值远离中心值的程度,也称为“离中趋势”。离散程度越小,数据的代表性就越好。
5.3.1 离散趋势的测度
数据分布的离散程度主要靠极差、四分差、平均差、方差、标准差等统计指标来度量。在实际分析中,离散程度分析主要有以下作用:
-
衡量平均指标的代表性;
-
反映社会经济活动的均衡性;
-
研究总体标志值分布偏离正态的情况;
-
抽样推断等统计分析的一个基本指标。
(1)极差
首先从最简单的统计量入手,极差是描述数据离散程度最简单的测度值,是样本数据中两个极端值之差,也称为全距。数据越分散,极差越大。
极差计算简单又容易理解,但它只是利用了数据两端的信息,容易受极端值的影响,并且没有充分利用数列的信息
(2)四分位差
数据的四分位差是两个四分位点之差,反映了中间50%数据的离散程度,其数值越小,说明中间的数据越集中。
虽然四分位差避免了极端数据的影响,但“掐头去尾”后仍然丢失很多信息。
(3)方差与标准差
描述离散程度,最常用的指标是方差和标准差,它们利用了样本的全部信息去描述数据取值分散性。方差是各样本相对均值的偏差平方和的平均,即
方差的平方根,称为标准差。与方差不同的是,标准差是具有量纲的,它与变量值的计量单位相同,因此具有较强的实际意义,在实际分析中会比较常用。
(4)变异系数
一组数据的标准差与相应平均数之比,称为变异系数,也叫做离散系数。它是刻画数据相对分散性的一种度量,记为CV。既然是“相对”,说明它去除了单位的影响,因此是个无量纲的统计量,用百分数表示。在实际应用中可以消除由于不同计量单位、不同平均水平所产生的影响。
5.3.2 R语言实现
可以通过函数range()计算极差。给出最小值和最大值两个点,再相减得到。
> m=range(price);m[2]-m[1]
[1] 60.34
>max(price)-min(price)
[1] 60.34
四分位差同样需要手动计算,比较便捷的方法是直接使用函数fivenum()。
> q=fivenum(price);q[4]-q[2]
[1] 11.89
R中的方差函数和标准差函数分别是var()和sd(),非常简单。还可以手动输入公式,计算验证一下。
var(x, y = NULL, na.rm = FALSE, use)
sd(x, na.rm = FALSE)
>var(price);sd(price)
[1] 138.7829
[1] 11.78061
>sqrt((sum(price^2)-(sum(price))^2/length(price))/(length(price)-1))
[1] 11.78061
除了上面几个最常用的离散趋势统计量外,R还有一个比较特殊的函数,即离差mad(),它用于计算中位数绝对偏差,具有渐近正态的一致性。
mad(x, center = median(x), constant = 1.4826, na.rm = FALSE, low = FALSE, high = FALSE)
其中的参数center可选,默认为中位数,即计算样本偏离中位数的程度;constant是比例因子,默认为一个常数;参数low/high默认取值为FALSE,若是TRUE,表示对于偶数个样本数据,中位数不取两个中间值的平均,而是采用较小/大的那一个。在默认状态下,R中离差的计算公式相当于
其中系数1.4826约等于1/qnorm(3/4),这样离差mad()在估计方差时才具有渐近正态的一致性。
>mad(price)
[1] 7.813302
5.4 数据的分布分析
5.4.1 分布情况的测度
仅通过集中趋势和离散趋势两个指标,仍无法判断数据的分布形态。数据比较集中的区域在分布图上就会形成一个“峰”,这个“峰”可能偏左,也可能偏右;峰值可能较高,也可能较低,这时就需要用到偏态(Skewness)和峰度(Kurtosis)两个统计方法对数据分布特征作进一步描述。
(1)偏度
前面提到,通过均值和平均数、众数之间的关系可以大致判断出数据点分布状态是对称、左偏还是右偏。显然,判断偏斜方向并不难,但要判断偏斜的程度却需要更精确的测度,即偏度系数。样本的偏度系数记为SK,计算公式为
其中,s是标准差, 是样本的3阶中心矩。
具体的判断方法是:如果样本的偏度系数SK=0,则说明数据关于均值对称;当SK<0时,样本是左偏分布(Negative Skew);当SK>0时,表示正偏离程度较大,说明样本是右偏分布(Positive Skew),如图5.2所示。
(2)峰度
偏度系数判断尖峰在水平方向上往哪边偏移,可以说它是从横向的角度描述分布特点,那么从纵向的角度,描述数据分布的平坦或尖峰程度(称为峰态),就要通过峰度系数来衡量。记为K。
其中, 是样本4阶中心矩, 是样本标准差的四次方。对于分组数据,利用组中值 计算峰度系数, 代表组内的频数。
峰态通常是与标准正态分布相比较而言的,如果一组数据服从标准正态分布,则峰度系数K=0;当分布比正态更加尖峰时,K>0;当分布比较扁平时,K<0。如图5.3所示,当K>0时,分布为尖峰形状,因此分布两侧的尾部数据较少(Thin tails);而K<0时,分布的峰值较低,造成尾部较厚(Fat tails)。
5.4.2 R语言实现
在程序包timeDate中(或直接加载fBasics程序包),有直接计算偏度和峰度系数的函数,为skewness()和kurtosis()。
也可以自己写一个函数function()进行验证,由于计算方法上的一些差别,两种方法计算的结果有出入,但不影响对样本分布形态的判断。
>library(timeDate) #加载程序包
>skewness(price) #计算偏度系数
[1] 1.683351
attr(,"method")
[1] "moment"
> SK=function(x){
+ n=length(x);m=mean(x)
+ C=n/((n-1)*(n-2))*sum((x-m)^3)/(sd(x))^3;C
+ }
>SK(price) #手动计算偏度系数
[1] 1.686996
>kurtosis(price)
[1] 2.578657
attr(,"method")
[1] "excess"
> K=function(x){
+ n=length(x);m=mean(x);s=sd(x)
+ g2=((n*(n+1))/((n-1)*(n-2)*(n-3))*sum((x-m)^4)/s^4-(3*(n-1)^2)/((n-2)*(n-3)))
+ g2
+ }
>K(price)
[1] 2.600382
偏度系数SK=1.687>0,说明样本是右偏分布,并且右偏程度比较严重。峰度系数K= 2.60>0,说明分布比正态更尖峰。
5.5 图形分析及R实现
电影《达芬奇密码》中有这样一句话:“A picture says more than a thousand words!”意思是,利用图形描述数据规律一目了然,胜过千言万语。图形显示的结果可以直观地概括上述指标所反映的分布特点,对于一组数据的描述,主要通过直方图、茎叶图等。
5.5.1 直方图和密度函数图
获得了一组样本数据,我们最直接的想法是通过频率分布直方图观察其样本的分布,常与直方图配套出现的是核密度估计函数density(),它可以根据已知样本,去估计样本的密度函数曲线。结合直方图和密度估计曲线,可以得到和统计指标所反映的相似结论:样本呈现严重的右偏分布,右尾拖得很长;数据有一较高的峰值。
>hist(price,breaks=50,prob=T) #参数breaks设置直方图的组距,prob=T规定绘制密度直方图
>lines(density(price),col="blue") #利用核密度估计函数density(),绘制密度曲线图
绘制结果如图5.4所示。
5.5.2 QQ图
QQ图用于直观验证一组数据是否来自某个分布,或者验证某两组数据是否来自同一族的分布。在教学和软件中常用QQ散点图来检验数据是否来自于正态分布。假定总体为正态分布 ,对于样本 ,其顺序统计量为 。 是标准正态分布 的分布函数,而 是其反函数,QQ图即是由以下的点构成的散点图:
R软件中的函数qqnorm()和qqline()可以分别用于绘制正态QQ散点图和相应直线,其使用方法为:
qqnorm(y, ylim, main="Normal Q-Q Plot", xlab="Theoretical Quantiles",
ylab="Sample Quantiles", plot.it = TRUE, datax = FALSE, ...)
qqline(y, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75), qtype = 7, ...)
qqplot(x, y, plot.it = TRUE, xlab = deparse(substitute(x)), ylab = deparse(substitute(y)), ...)
其中x,y为数据样本;xlab,ylab和main是图标签,分别用于设置X轴和Y轴的标签、图像标题。其它参数用法参见第四章plot()的参数设置。
对股票价格样本绘制qq散点图和曲线图,判断它是否来自于正态总体。
> qqnorm(price)
> qqline(price)
绘制结果如图5.5所示。
从图5.5可以看出,QQ散点图并不均匀地分布在直线两侧,因此股价样本不服从正态分布。接下来举一个正态分布样本的示例。
> a=c(21.0,21.4,22.5,21.6,19.6,20.6,20.3,19.2,20.2,21.3)
> qqnorm(a);qqline(a)
绘制结果如图5.6所示。
图 5.6 正态QQ图
从图5.6可以看出,样本数据是来自于正态总体。
5.5.3 茎叶图
茎叶图也是考察样本分布的重要方法,由统计学家约翰·托奇设计。它是将数组中的数按位数进行比较,将数的大小基本不变或变化不大的位作为一个“茎”,将位变化大的数作为“叶”,列在茎的后面,这样就可以清楚地看到每个茎后面有多少样本数,以及具体数值是多少。
茎叶图更适用于描述样本量较小的分布,在本例中,股票价格样本有逾千个,茎叶图不能完全显示,所以抽取其中的50个样本作图,R中用函数stem()绘制茎叶图。
stem(x, scale = 1, width = 80, atom = 1e-08)
其中,x是数据向量,scale控制茎叶图的长度,width控制绘图的宽度,atom是容差。
>set.seed(111) #设置抽样种子
>s=sample(price,50) #从数据集price中抽取50个样本
> stem(s)
结果如下所示:
The decimal point is 1 digit(s) to the right of the |
1 |
1 | 5888999999999
2 | 01112223334444
2 | 555556788889
3 | 022
3 | 788
4 | 1
4 | 8
5 | 3
5 | 69
在上面的结果中,左边的“茎”是股价的十位数字,右边的“叶”是个位数字。茎叶图从外观上来讲,比较像是横放的直方图,我们所抽取的样本也是一个右偏分布。“叶”的每一个数字都代表一个样本,显示频数的个数,所以从上向下数第25、26个数字的平均可以体现样本中位数的值。
5.5.4 箱线图
箱线图在描述统计中也是非常重要的图形,可以说它的重要程度不亚于直方图。箱线图是“五数”的一个图形概括,可以用来对数据分布的形状进行大致判断,是直观展现数据分布主要特征的方式。在上一章已经介绍过,R中使用函数boxplot()绘制箱线图。
函数boxplot()有三种使用方法,第一种最为直观和简单,因此也最常用:
boxplot(x, ...)
其中,x是数值型向量、列表或数据框。第二种形式为:
boxplot(formula, data = NULL, ..., subset, na.action = NULL)
formula是公式,形如y~grp,这里y是数值型向量,而grp通常为因子,用来表示数据的分组,这种形式的绘图结果是分类的箱线图,便于比较。第三种形式比较复杂:
boxplot(x, ..., range = 1.5, width = NULL, varwidth = FALSE, notch = FALSE, outline = TRUE,
names, plot = TRUE, border = par("fg"), col = NULL, log = "",
pars = list(boxwex = 0.8, staplewex = 0.5, outwex = 0.5),
horizontal = FALSE, add = FALSE, at = NULL)
其中,x是数值型变量、列表或数据框;range是“触须”的范围,默认值为1.5;notch是逻辑变量,默认值FALSE表示绘制的箱线图不带有切口;outline也是逻辑变量,默认值为TRUE,即标明异常值点;col用于设置不同颜色(根据不同的数值);horizontal是逻辑变量,默认值为FALSE,设置成TRUE时将箱线图绘制成水平状。
要绘制出本章股票价格样本的箱线图,在R中输入如下指令:
> boxplot(price,main="Boxplot of price")
绘制结果如图5.7所示。
在图5.7的箱线图中,75%和25%四分位数分别为中间箱体的上下两条线,箱体中间的粗线段是中位数 所在的位置。箱体外的两条线是数据的最大值和最小值,但不得超过1.5倍四分位数间距,若是超出了则标记为异常值点,用“o”表示。
在本例中,股价分布的右尾很长,因此就出现了较多的异常值点,这是在实际数据分析中很有可能出现的情况。
5.5.5 经验分布图
直方图主要适用于总体为连续分布的样本。对于更一般的总体分布,如果要估计其分布函数,可以使用经验分布函数(edf)。在R中函数ecdf()给出样本的经验分布,通过plot()绘制。
ecdf(x)
plot(x, ..., ylab="Fn(x)", verticals = FALSE,col.01line = "gray70", pch = 19)
通过ecdf()生成经验分布函数 在各样本点的值,是一个数值型向量;plot()中的verticals是逻辑变量,若为TRUE表示画出竖线,默认是不画竖线。
>plot(ecdf(price),main="empirical cdf",xlab="price",ylab="Fn(price)")
> x=min(price):max(price)
> lines(x,pnorm(x,mean(price),sd(price)),col="red") #正态分布函数曲线
> legend(60,0.4,legend=c("样本分布","正态分布"),lty=1,col=1:2)
绘制结果如图5.8所示。
5.6 多组数据分析及R实现5.6.1 多组数据的统计分析当数据分为多个组别时,其描述统计与单组数据有相通之处,但在体现变量关系方面还是有一定差别的。上文中使用的例子是一组数据,假设按照年份将之分组,重新读入2009、2010和2011这三年的完整数据,在前面的理论基础上,对三组数据进行描述性统计分析。由于每年的交易日不完全相同,为简单处理,把有缺失值的样本删去。 > group=read.csv("d:/data/09-11股价.csv") > group=na.omit(group) #忽略缺失样本 首先,与单组情形类似,直接使用summary()可以得到各组数据的均值和“五数”。通过观察横向数据,我们可以对比不同年份股票价格的基本统计量。以均值为例,它反映出股票价格的平均水平是逐年下降的。 >summary(group) Price09 Price10 Price11 Min. :18.67 Min. :20.90 Min. :14.74 1st Qu. :23.08 1st Qu. :23.25 1st Qu. :16.98 Median :27.45 Median :24.70 Median :18.53 Mean :26.51 Mean :25.20 Mean :18.72 3rd Qu. :29.96 3rd Qu. :27.40 3rd Qu. :21.05 Max. :34.01 Max. :31.42 Max. :22.33 函数var()应用在多组数据上,得到的计算结果是一个协方差阵,其每个元素是各个向量之间的协方差。使用指令cor(group)也得到相同结果。 > options(digits=3) #设置显示格式:数字只显示3位 >var(group) Price09 Price10 Price11 Price09 16.55 -7.18 -6.75 Price10 -7.18 5.99 4.16 Price11 -6.75 4.16 4.84 协方差的大小在一定程度上反映了变量之间的相互关系,但它还受变量本身度量单位的影响,因此我们还要计算相关系数来度量变量之间的线性相关程度。在R中使用函数cor()计算相关系数矩阵。 cor(x, y = NULL, use = "everything",method = c("pearson", "kendall", "spearman")) 其中,x,y是计算的对象,当x是一个数据框或列表时y可以省略;use指定如何处理缺失样本;method给出计算哪一种相关系数:默认的皮尔逊(Pearson)系数度量线性相关性,如果数据呈现的不是线性关系,而是单调的,则可以用肯德尔(Kendall)或斯皮尔曼(Spearman)相关系数,它们描述的是秩相关性。 本例中没有理论基础可以说明两年的股价会呈线性关系,所以计算Spearman系数。 >cor(group,method="spearman") Price09 Price10 Price11 Price09 1.000 -0.687 -0.710 Price10 -0.687 1.000 0.752 Price11 -0.710 0.752 1.000 在上面的矩阵中,对角线都是1,因为变量与自身是完全相关的。而2009年股价与2010年股价的相关系数是-0.687,说明这两年的股票价格呈负相关的关系。 5.6.2 多组数据的图形分析 常用的图形表示方法大多用于一维数据的描述,但实际问题中的数据集往往是多变量样本。接下来的这一节,将介绍在实践中常用的多组数据绘图。 (1)二维散点图 首先,两组数据的关系用散点图可以清楚地描述。在散点图中加入一条拟合曲线有助于更好地把握变量关系。 R中的函数lowess()可通过加权多项式回归对散点图进行平滑,拟合一条非线性的曲线,但其只能适用于二维情况。与之类似的loess()用于处理多维情况。 lowess(x, y = NULL, f = 2/3, iter = 3, delta = 0.01 * diff(range(x))) x,y指定两个向量;f是平滑的跨度,值越大,曲线的平滑程度越高;iter控制应执行的迭代数,值越高平滑越精确,但使用较小的值会使程序跑得比较快。 > attach(group)。 >plot(Price10~Price09,xlab="Price of 2009",ylab="Price of 2010") > lines(lowess(Price09,Price10),col="red",lwd=2) #拟合曲线 绘制结果如图5.9所示。
从散点图中可以看出,2009年和2010年的股票价格并不是完全线性的关系,这与由Spearman系数得出的结论相同,两年的股价具有负的相关性,价格变动趋势相反。 (2)等高线图 有时候数据量很大,散点图上的数据点就会非常集中,不容易看出变量的关系或趋势,这就需要借助二维等高线图来描述。首先利用程序包MASS中的函数kde2d()来估计出二维数据的密度函数,再利用函数contour()画出密度的等高线图。如果不想画出图上的数据标签,可以将参数drawlabels=FALSE去掉。函数kde2d()的使用方法如下: 其中x,y分别为横轴和纵轴的数据;n指定每个方向上的网格点数量,可以是标量或长度为2的一个正数向量;参数lims表示横纵轴的范围。 绘制结果如图5.10所示。
与地理上的等高线地形图相同,即等高线越靠近内部说明密度越高,散点分布越集中。 (3)矩阵散点图 或 绘制结果如图5.11所示。
|
(4)矩阵图
在处理多组数据时,常将各组数据放在一起进行比较,matplot()可将各变量的散点图放在同一个绘图区域中。
上文计算的相关系数矩阵反映出2009年和2010年的股价成负相关,而2010年2011年的股价则变动趋势相同,这些结论从图5.10中也可以得到,黑色曲线表示2009年股价呈上升趋势,而2010和2011年股票价格则整体下跌,在2011年更是创出了历史新低。
>matplot(group,type="l",main="Matplot") #type=“l”表示绘制曲线而不是散点
> legend(0,35,legend=2009:2011,pch="——",cex=0.6,col=1:3) #cex指定字体大小
绘制结果如图5.12所示。
(5)箱线图
在5.5.4节中已经详细介绍过箱线图,函数boxplot()也可以用于多组数据的绘图,在同一区域中画出各变量的箱线图,便于对比。箱线图的好处在于可以同时对比最值、均值、中位数和整体分布情况,尤其适用于时间序列的对比,从图5.13中的几条“线”可以看出2011年的股票价格从最大值、75%分位点到中位数、最小值等都低于2009年和2010年,因此2011年的股价水平整体低于前两年。
>boxplot(group,cex.axis=0.6)
(6)星图(雷达图)
假设变量是p维数据,有n个观测值,故第k次的观测值写作向量形式为
星图作为一种多元数据的图示方法,可以对p维数据进行比较。由于星图(如图5.14所示)看起来很像雷达图像,因此也称为雷达图;星图也像一个蜘蛛网,所以有时也称为蜘蛛图。其绘制原理是:
(1)作一个圆,并将圆周p等分。
(2)连接圆心和圆周上的p个分点,把这p条半径依次视为变量坐标轴,并标上适当的刻度。
(3)对给定的每一次观测值,在上述p条坐标轴上分别标出p个变量的取值,然后将它们连接成一个p边形,即形成了一个“星”。
(4)对n次观测值都重复(3)的步骤,就形成了n个星图。
R中的函数stars()可以绘制星图,其使用方法是:
stars(x, full = TRUE, scale = TRUE, radius = TRUE, labels = dimnames(x)[[1]],
locations = NULL, nrow = NULL, ncol = NULL, len = 1,
key.loc = NULL, key.labels = dimnames(x)[[2]], key.xpd = TRUE,
xlim = NULL, ylim = NULL, flip.labels = NULL, draw.segments = FALSE,
col.segments = 1:n.seg, col.stars = NA, col.lines = NA, axes = FALSE,
frame.plot = axes, main = NULL, sub = NULL, xlab = "", ylab = "",
cex = 0.8, lwd = 0.25, lty = par("lty"), xpd = FALSE,
mar = pmin(par("mar"), 1.1+ c(2*axes+ (xlab != ""), 2*axes+ (ylab != ""), 1, 0)),
add = FALSE, plot = TRUE, ...)
stars()函数中包括很多的参数设置,这里着重介绍其中比较重要的参数,其他请参考help文档。其中x是数据矩阵或数据框;full是逻辑变量,默认值为TRUE表示星图绘制成整圆的,FALSE表示星图绘制成上半圆的图形;scale也是逻辑变量,默认值TRUE表示数据矩阵的每一列是独立的,且每列的最大值为1,最小值为0,若值为FALSE则所有的星图会叠在一起;radius是逻辑变量,默认值为TRUE,表示画出构成星图半径的连线,FALSE则画出的星图没有半径构成的连线;len是半径刻度因子,设置星图的比例,默认值为1;key.loc是由x和y坐标值构成的向量(默认值为1),设置星的位置;draw.segments同样是逻辑变量,默认值为FALSE,即用直线连接p个半径上的点,若为TRUE则画出的星图是一段一段的弧。
但是星图比较适用于维度较低的数据,方便比较。例如学校随机抽取6名学生的5门考试成绩,如表5.2所示,从D盘读入该数据集,我们将画出这6名学生成绩的星图。
表5.2 学生考试成绩
科目1 | 科目2 | 科目3 | 科目4 | 科目5 |
89 | 90 | 92 | 95 | 99 |
90 | 77 | 74 | 86 | 89 |
79 | 88 | 94 | 98 | 97 |
100 | 90 | 85 | 100 | 81 |
83 | 87 | 70 | 77 | 73 |
100 | 94 | 97 | 90 | 89 |
在R软件中输入如下指令:
> score=read.table("d:/data/score.txt",header=T) #读入数据
> stars(score) #绘制星图,所有参数均设置为默认值
绘制结果如图5.14所示。
星图的半径反映数值的大小,因此从图5.14中可以看出学生编号1和6的成绩较好,而第5名学生成绩最差,第2名其次。
设置函数stars()中的参数,可以将上面的星图绘制成图5.15所示的另一种形式。在R中输入指令:
> stars(score,full=FALSE,draw.segments=TRUE,key.loc=c(5,0.5),mar=c(2,0,0,0)) #full=F
# 表示绘制成半圆的图形,draw.segments=T表示画出弧线
(7)折线图
和星图相比,折线图是一种更加直观地展示多维图形的方法,但同样地,折线图比较适用于数据量较小的数据样本。
折线图由以下步骤绘制而成:
(1)在直角坐标系的横轴上取p个点,表示p个变量;
(2)对每一次观测值,p个点的纵坐标即表示变量的取值;
(3)连接p个点可以得到一条折线,即为该次观测值的折线图;
(4)对于n次观测值,均重复上述步骤,即可得到n次观测值的折线图。
在R中并没有现成的绘制折线图的函数,因此需要我们自行编写,这里命名为函数outline()。
outline=function(x){
if(is.data.frame(x)==TRUE)
x=as.matrix(x) #若x为数据框,则先转换为矩阵形式
m=nrow(x);n=ncol(x) #提取x的行列数
plot(c(1,n),c(min(x),max(x)),type="n",
main="The outline graph of data",xlab="Number",ylab="Value")
for(i in 1:m){
lines(x[i,],col=i)
}
}
其中的参数x是数据矩阵或数据框,函数的运行结果是n次观测值的折线图。例如对上面学生成绩的例子绘制折线图,则应在R中输入如下指令:
> outline(score)
绘制结果如图5.16所示。
从图5.16中可以看出,越是靠近图形上方的学生成绩越好。本节仅给出一个编写函数的示例,若样本量较多,编写函数时应在绘图中加入样本编号。这种折线图可以用于聚类分析的初步判断。
(8)调和曲线图
调和曲线图是D.F.Andrews在1972年提出的三角多项式作图方法,所以又称为三角多项式图,其思想是把高维空间中的一个样本点对应于二维平面上的一条曲线。
设p维数据 ,对应的曲线是
上式中,当t在区间[-π,π] 上变化时,其轨迹是一条曲线。
在多项式的图表示中,当各变量的数值太悬殊时,最好先标准化后再作图。这种图对聚类分析帮助很大,如果选择聚类统计量为距离的话,则同类的曲线非常靠近拧在一起,不同类的曲线相互分开,非常直观。
调和曲线图有两点优良数学性质:一是保持线性关系,二是与一般的欧式距离之间的关系。第一点从上面的公式中就可以直观得到结论,第二个性质可通过计算得到。按照一般常用的范式,定义样本x和y之间的距离为如下的平方积分形式
根据三角函数积分,通过运算最终可以得到这一距离与欧氏距离之间具有如下关系:
按照上式的描述,我们可以在R中自己编写调和曲线图函数unison()。
unison=function(x){
if (is.data.frame(x)==TRUE)
x=as.matrix(x) #若x为数据框,则先转换为矩阵形式
t=seq(-pi, pi, pi/30) #设置t的变化范围
m=nrow(x); n=ncol(x) #提取x的行列数
f=array(0, c(m,length(t))) #f赋值为一个数组
for(i in 1:m){
f[i,]=x[i,1]/sqrt(2)
for( j in 2:n){
if (j%%2==0)
f[i,]=f[i,]+x[i,j]*sin(j/2*t)
else
f[i,]=f[i,]+x[i,j]*cos(j%/%2*t)
}
}
plot(c(-pi,pi),c(min(f),max(f)),type="n",main ="The Unison graph of Data",xlab="t",ylab="f(t)")
for(i in 1:m) lines(t,f[i,],col=i)
}
其中,x是数据矩阵或数据框。函数的运行结果是调和曲线图(如图5.17所示)。
图5.17虽然不能用于判断学生成绩的好坏,但这种调和曲线图对聚类分析有很大帮助,可以用于初步判断,即同类的曲线会聚集在一起,而不同类的曲线则会拧成不同的曲线束,比较直观。