—
单变量数据可视化
在这里，我们先将我们需要的包加载出来。
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 来建造的library(JWileymisc)## Warning in .recacheSubclasses(def@className, def, env): undefined subclass
## "packedMatrix" of class "replValueSp"; definition not updated## Warning in .recacheSubclasses(def@className, def, env): undefined subclass
## "packedMatrix" of class "mMatrix"; definition not updatedggplot2是用来创建图形，cowplot是一款插件，用来使图形更加干净整洁，MASS包提供了很多常用函数，这些函数可以检验各种分布曲线和数据拟合程度。
直方图、密度图、Q-Q图
直方图
直方图用于表示数据在某个范围内出现的次数，这个范围由条形的宽度决定。例如，下面的直方图显示的花萼长度的分布，这些数据来自著名的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图。利用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)
从上图中我们可以看到，图中绘制的点都非常接近于正太直线，而且基本上以这条直线为对称轴，分布在他的两侧，因此我们有理由认为这个数据集属于正态分布。
在实际中，实际数据也可能时其他分布。因此我们使用qlnorm()来显示原始鲍鱼数据与对数正态分布的拟合程度，如下图所示。
ggplot(data.table(lynx = as.vector(lynx)), aes(sample = lynx)) +
geom_qq(distribution = qlnorm)
对于有些不常用的分布，ggplot2也无法给定正确的分布默认参数。这时我们需要通过dparams参数把期望分布所需要的参数传导给ggplot 函数。例如，下图中，我们为了测试鲍鱼数据是否符合泊松分布，可以通过dparams参数把lamda值传递给ggplot函数。
ggplot(data.table(lynx = as.vector(lynx)), aes(sample = lynx)) +
geom_qq(distribution = qpois, dparams = list(lambda = mean(lynx)))
数据密度图与期待分布的密度图
首先，我们用geom_density() 函数绘制数据密度图，然后用stat_function()函数绘制期望分布的密度图。后者是一个泛化方式，可以绘制任何函数。
ggplot(iris, aes(Sepal.Length)) +
geom_density() +
stat_function(
fun = dnorm,
args = list(
mean = mean(iris$Sepal.Length),
sd = sd(iris$Sepal.Length)
), colour = "blue"
)
其他分布拟合
对于正态分布以外的分布，为了得到分位数或者密度，通常需要定义他们的参数。利用MASS包里提供的fitdistr()函数，我们只需要说明分布的名称，就可以拟合许多分布并且让R语言自动估算这些参数。
目前，findistr()函数支持以下分布：
Beta分布
柯西分布
方卡分布
指数分布
伽马分布
几何分布
对数正态分布
Logistic分布
负二项分布
正态分布
泊松分布
t分布
韦伯分布
我们西安根据beta分布生成随机数据，为确保后面的随机数的重复性，我们使用set.seed()函数。
set.seed(1234)
y <- rbeta(150,1,4)
head(y)## [1] 0.13817170 0.03882067 0.11137033 0.09945836 0.37703993 0.38402555fitdistr()函数需要传入一个数据集，一个表达分布名称的字符串和一个表达分布参数的起始值的列表：
y.fit <- fitdistr(y, densfun = "beta",
start = list(shape1 = .5, shape2 = .5))## Warning in densfun(x, parm[1], parm[2], ...): 产生了NaNs
## Warning in densfun(x, parm[1], parm[2], ...): 产生了NaNs
## Warning in densfun(x, parm[1], parm[2], ...): 产生了NaNs
## Warning in densfun(x, parm[1], parm[2], ...): 产生了NaNsy.fit## shape1 shape2
## 1.0830304 4.2790564
## (0.1112984) (0.5182058)y.fit$estimate["shape1"]## shape1
## 1.08303y.fit$estimate["shape2"]## shape2
## 4.279056对数似然值可以告诉我们观测值来自某个分布的可能性有多大。logLik()函数可以得到对数似然值和自由度。
logLik(y.fit)## 'log Lik.' 95.35267 (df=2)y.fit2 <- fitdistr(y, densfun = "normal")
logLik(y.fit2)## 'log Lik.' 67.25708 (df=2)相同自由度下，beta分布拟合的似然值较大，正态分布拟合的你染指较小，以上结果建议我们应该给这个数据集选择beta 分布。
异常值
对于来自正态分布的数据，确定奇异值常用的阈值被定义为z分数+-3，在此范围外的是异常值。假设一个数据集属于正态分布，有0.1%的数据小于z分数-3，大约有0.1%的数据点大于z分数+3。精确的概率值可用pnorm()函数求得。
pnorm(c(-3,3))## [1] 0.001349898 0.998650102通常，比起用定量的准则确定异常值，用可视化方法更容易识别它们。如下有两个图形，这两个图形都在值为5的位置出现三个异常点。然而，图形 A的这三个点显得与其他点格格不入,因为其他点形成了一个几乎连续的一组数据，而它们与异常点有较大的间距。在图形B中，由于其他点之间也存在间距，因此一个点离其他点稍微远一点也不足为奇。在图形B中，数据之间似乎存在一种分离模式，而图形A则没有这种模式。
set.seed(1234)
d <- data.table(
y1 = rnorm(200,0,1),
y2 = rnorm(200,0,.2) + rep(c(-3,-1,1,3), each=50)
)
plot_grid(
qplot(c(d$y1, rep(5,3)), geom = "dotplot", binwidth = .1),
qplot(c(d$y2, rep(5,3)), geom = "dotplot", binwidth = .1),
ncol = 1, lables = c("A", "B")
)## Warning in as_grob.default(plot): Cannot convert object of class character into
## a grob.
难以确定异常值的另一个原因与分布的形状有关。如下图所示，图形A是伽马分布，特征是有一条长长的右尾。虽然有几个点处在极端位置，但是它们与其他数据点并没有明显的“中断”，因为这类分布的特征就是有一条连续的长长的右尾。而且，我们可以很容易看出，这个数据集存在一种频数可压缩的模式，但是几个非常极端的正数除外。相反，图形B是正态分布的，形状比较对称，并且没有出现长尾。在图形B中，添加一两个极端正值，它们就可能成为所谓的异常值。
set.seed(1234)
d2 <- data.table(
y1 = rgamma(200,1,.1),
y2 = rnorm(200,10,10)
)
plot_grid(
qplot(d2$y1, geom = "dotplot", binwidth = 1),
qplot(d2$y2, geom = "dotplot", binwidth = 1),
ncol = 1, lables = c("A","B")
)## Warning in as_grob.default(plot): Cannot convert object of class character into
## a grob.
最后，当奇异值确认后，处理它们时有几种可选方法。如果有可能，最好先检查一下这些值是否准确。一种不同于排除法的方法是温莎法(winsorizing)，这种方法是以Charles Winsor命名的。温莎法接收奇异值，但把它们替换为距离最近的非奇异值。自动实现温莎法的一种办法是计算期望的经验分位数，把任何落在这些经验分位数之外的数值设置为已计算得到的分位数。即使奇异值只存在于分布的尾巴上，这个算法也同样可以作用于上尾和下尾。
温法容易用JWileymisc包的winsorizor函数来实现。唯一需要设置的参数是每端需要缩尾的比例。winsorizor()函数的一个特性就是除了返回被缩尾的变量外，还会添加属性以说明缩尾的界限和百分分位数。
winsorizor(1:10,.1)## [1] 1.9 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 9.1
## attr(,"winsorizedValues")
## low high percentile
## 1 1.9 9.1 0.1Reference
Wiley, M., & Wiley, J. F. (2019). Advanced R Statistical Programming and Data Models. Berkeley, CA, USA: Apress.