在这里,我们先将我们需要的包加载出来。
library(ggplot2)
## Warning: 程辑包'ggplot2'是用R版本4.1.3 来建造的
library(cowplot)
## Warning: 程辑包'cowplot'是用R版本4.1.3 来建造的
library(MASS)
library(data.table)
## Warning: 程辑包'data.table'是用R版本4.1.3 来建造的
library(knitr)
## Warning: 程辑包'knitr'是用R版本4.1.3 来建造的
ggplot2是用来创建图形,cowplot是一款插件,用来使图形更加干净整洁,MASS包提供了很多常用函数,这些函数可以检验各种分布曲线和数据拟合程度。
直方图用于表示数据在某个范围内出现的次数,这个范围由条形的宽度决定。例如,下面的直方图显示的花萼长度的分布,这些数据来自著名的iris数据集。
ggplot(iris,aes(Sepal.Length)) + geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
我们可以发现上图中花萼的长度数据的直方图与正态分布基本相似,尽管并不是完全一样。
如果数据的分布不符合我们所期望的分布,这时我们需要对原始数据进行变换。如下图所示,这里显示的是每年加拿大鲍鱼捕获的分布。
ggplot(data.table(lynx = as.vector(lynx)), aes(lynx)) +
geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
我们可以看出,他是正偏分布,有一个长长的右尾,通常我们也叫做肥尾效应,在金融中经常出现。
对于这样的正偏分布,我们可以用平方根或者对数变换来使其尾部正向偏移,这样才可以使数据更加接近正太分布。下图我们将自然对数变换应用于鲍鱼捕获量上,如下图所示
ggplot(data.table(lynx = as.vector(lynx)), aes(log(lynx))) +
geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
绘制密度图的代码与绘制直方图的代码非常相似。只是吧代码中的geom_histogram() 函数换成geom_density() 函数。下面是相应的代码和图形展示。
ggplot(iris,aes(Sepal.Length)) +
geom_density()
在绘制密度图时,要考虑到分布曲线的光滑度。对于连续的变量,不可能在任意值位置都有观测值,例如,即使在3.281和3.292位置都有观测数据,也不能保证在3.286位置也有观测值。密度图用光滑度表示分布曲线的总体形状。我们可以通过调整光滑度,得到接近于原始数据的曲线或者接近实际分布的曲线。默认情况下,adjust=1。如下图所示,adjust小于1“噪声较大”,图形也不够光滑;如果adjust大于1,则增加光滑度。
ggplot(iris,aes(Sepal.Length)) +
geom_density(adjust=.5)
ggplot(iris,aes(Sepal.Length)) +
geom_density(adjust=5)
对观测数据的分布进行研究很有必要,通常研究分布的目的是判断分布是否满足我们想要使用的统计学分析的模型。
为了评估数据是否拟合或接近某个期望分布,我们可以使用分位数-分位数图,简称Q-Q图。利用Q-Q图可以分析观测数据是否来自某一个分布。正态分布或者高斯分布是我们经常使用的分布,所以ggplot2包在绘制Q-Q图时使用的默认函数分布是正态分布。如下图所示
ggplot(iris,aes(sample = Sepal.Length)) +
geom_qq()
我们手动生成一个观测值数据集。在这里我需要补充的是,R语言内置有许多概率函数,其中有rnorm()可以生成随机数,pnorm()可以生成某个观测值大于或者小于特定值的概率,qnorm()可以计算正太分布的分位数,dnorm()可以得到正态分布在某个值位置的密度。前面分别是r、p、q和d,后面紧跟分布名。
qnorm(p=.1, mean=0,sd=1)
## [1] -1.281552
上述代码可以得到均值为0,标准差为1的正态分布概率为0.1的分位数。这个代码很容易应用于不同均值、不同标准差的正态分布。假设有三个点以0.5为中心均匀分布。
qnorm(p=c(.25,.5,.75), mean=0,sd=1)
## [1] -0.6744898 0.0000000 0.6744898
为了得到等间隔的概率值,可以使用ppoints()函数
ppoints(n=3,a=0)
## [1] 0.25 0.50 0.75
使用ppoints()并不能得到完全等间隔的数据点,可以使用一个默认值作为微小调整量。
等于或者小于10个数据点时,间隔点计算公式是(1:N - 3/8)/(n+1-2* 3/8);大于10个数据点时,使用公式(1:N - 1/2)/(n+1-2* 1/2)。
ppoints(n=3)
## [1] 0.1923077 0.5000000 0.8076923
现在我们使用ggplot2包中的qplot()函数绘制图形。为了更好的视觉效果,我们使用geom_abline()函数来添加一条斜率为1,截距为0的直线。
qplot(
x = qnorm(
p = ppoints(length(iris$Sepal.Length)),
mean = mean(iris$Sepal.Length),
sd = sd(iris$Sepal.Length)),
y = sort(iris$Sepal.Length),
xlab = "正态分布的理论分位数",
ylab = "Sepal Length") +
geom_abline(slope = 1, intercept = 0)
从上图中我们可以看到,图中绘制的点都非常接近于正太直线,而且基本上以这条直线为对称轴,分布在他的两侧,因此我们有理由认为这个数据集属于正态分布。