多变量数据可视化

首先，我们还是先加载需要的包。

library(ggplot2)

## Warning: 程辑包'ggplot2'是用R版本4.1.3 来建造的

library(cowplot)

## Warning: 程辑包'cowplot'是用R版本4.1.3 来建造的

library(MASS)
library(mvtnorm)
library(mgcv)

## 载入需要的程辑包：nlme

## This is mgcv 1.8-36. For overview type 'help("mgcv-package")'.

library(quantreg)

## Warning: 程辑包'quantreg'是用R版本4.1.3 来建造的

## 载入需要的程辑包：SparseM

## 
## 载入程辑包：'SparseM'

## The following object is masked from 'package:base':
## 
##     backsolve

## 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 updated

library(data.table)

## Warning: 程辑包'data.table'是用R版本4.1.3 来建造的

library(JWileymisc)

## Warning: 程辑包'JWileymisc'是用R版本4.1.3 来建造的

library(magrittr)

分布

评估单变量分布相对比较容易，但是评估多变量分布却比较困难。多个变量的联合分布是由他们的单个变量分布组成的。

多变量正态分布

正态分布由两个参数决定：均值（）和标准差（），分别控制分布的位置和尺度。在这里我们首先使用mvtnorm包里的rmvnorm()函数生成一个多变量正态数据集。然后使用geom_density2d()函数和ggplot2包绘制密度图，如下图所示：

mu <- c(0,0)
sigma <- matrix(c(1,.5,.5,1),2)

set.seed(1234)
d <- as.data.table(rmvnorm(500, mean = mu, sigma = sigma))
setnames(d, names(d), c("x","y"))
ggplot(d, aes(x,y)) + 
  geom_point(colour = "grey60") + 
  geom_density2d(size = 1, colour = "black") + 
  theme_cowplot()

为了对获得的密度与多变量正态分布的期望密度进行比较，我们创建一个由x 和y 值构成的网络。然后使用dmvnorm()函数计算每个位置(x,y)的密度。为了清除起见，删除了原始数据。

testd <- as.data.table(expand.grid(
  x = seq(from = min(d$x), to = max(d$x), length.out = 50),
  y = seq(from = min(d$y), to = max(d$y), length.out = 50)
))
testd[,Density:= dmvnorm(cbind(x,y), mean = colMeans(d), sigma = cov(d))]

ggplot(d,aes(x,y)) + 
  geom_contour(aes(x,y,z = Density), data = testd,
               colour = "blue", size = 1, linetype = 2) + 
  geom_density2d(size = 1, colour = "black") + 
  theme_cowplot()

密度和正态分布密度非常接近，因此认为数据是多变量正态分布是合理的。虽然单个变量是正态分布的，但这不能保证它们的组合是多变量正态分布的。

如下图所示，这两个变量并不是多变量正态分布的。虽然这只是一个极端的例子，但它突出了这样一个事实:我们很容易受到单变量评估的误导。

set.seed(1234)
d2 <- data.table(x = rnorm(500))
d2[, y := ifelse(abs(x)>1, x, -x)]
testd2 <- as.data.table(expand.grid(
  x = seq(from = min(d2$x), to = max(d2$x), length.out = 50),
  y = seq(from = min(d2$y), to = max(d2$y), length.out = 50)
))

testd2[, Density:=dmvnorm(cbind(x,y), mean=colMeans(d2), sigma = cov(d2))]
ggplot(d2,aes(x,y)) + 
  geom_contour(aes(x,y,z=Density), data = testd2,
               colour = "blue", size = 1, linetype = 2) + 
  geom_density2d(size = 1, colour = "black") + 
  theme_cowplot()

当p>2时，直接可视化多变量正态分布是非常困难的，但我们可以利用 Mahalanobis的研究成果。他设计了一种方法，可以计算数据离“中心”的距离，数据和中心都可以是多维的。假设数据是多变量正态分布的，则这些距离可以看成p个自由度的卡方变量(Chi-square variable)。下图1是前面二元正态数据的绘制结果，图2是二元非正态数据的绘制结果。

plot(testDistribution(d, "mvnorm", ncol = 2))

plot(testDistribution(d2, "mvnorm", ncol = 2))

即使变量的数量不断增大，但只要数据通过了多变量正态性检验，它们的图形就是相似的。例如，检验mtcars数据集，结果如下图所示，这说明数据接近于多变量正态分布，但是QQ图微有点偏离标准形态，具有一定程度的非正态性。

plot(
  testDistribution(mtcars, "mvnorm", nco1=2)
)

奇异值

对于多变量数据，某些值可能是某个变量的奇异值，也可能是多变量空间中的奇异值。如果某个值是单变量奇异值，那么它肯定也是多变量奇异值。

在下面的代码里，我们模拟生成强正相关的多变量正态分布数据，给V3添加两个奇异值，删除奇异值后很难发现还有奇异值。

mu <- c(0, 0, 0)
sigma <- matrix(.7, 3, 3)
diag(sigma) <- 1

set.seed(1234)
d <- as.data.table(rmvnorm(200, mean = mu, sigma = sigma))[order(V1)]
d[c(1,200), V3 := c(2.2, 50)]

plot(
  testDistribution(d$V3, extremevalue = "theoretical", plot = FALSE), 
  labels = c("A")
)

plot(
  testDistribution(d[V3<40]$V3, extremevalue = "theoretical", plot = FALSE),
  labels = c("B")
)

变量之间的关系

对于连续变量，大多数模型都假设变量之间存在某种函数形式，通常是线性关系。

确定两个变量之间函数形式的另一种方法是使用局部带权回归(locally weighted regression，loess)线。局部带权回归线的思想类似于最佳拟合直线，但是能给比较近邻的数据点赋予较大权值。下面的代码会生成一幅图形，它把一张散点图、一条局部带权回归线和一条线性回归线叠加在了一起。局部带权回归线清楚地表明它是一条二次曲线，而线性回归直线无法提供此信息。

set.seed(12345)
d2 <- data.table(x = rnorm(200))
d2[, y := rnorm(200, mean = 2 + x + 2*x^2, sd = 3)]
ggplot(
  d2, aes(x,y)
) + 
  geom_point(colour = "grey50") + 
  stat_smooth(method = "loess", colour = "red") + 
  stat_smooth(method = "lm", colour = "blue", linetype = 2)

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

当大致确定了变量之间的函数关系后，就可以用参数函数来近似表示。在下面的代码中，我们改变线性模型的光滑度，增加一个自定义函数，它表示y回归于x和x,结果得到极大改善，这说明局部带权回线与二次曲线之间相当一致。

ggplot(d2, aes(x,y)) + 
  geom_point(colour = "grey50") + 
  stat_smooth(method="loess", colour ="red") + 
  stat_smooth(method="lm", 
              formula=y~x+I(x^2), 
              colour="blue",linetype=1)

## `geom_smooth()` using formula 'y ~ x'

尽管局部带权回归线有很多优点，但它并不是没有缺点。它需要设置光滑度。这个光滑度可以由span(表示间距)参数控制。这个参数最后需要传递给loess()函数后者才是回归曲线的真正执行者。在下例中，我们绘制了两条loess线，一条是低间距曲线，另一条是高间距曲线，高间距曲线比较光滑。

ggplot(d2, aes(x,y)) + 
  geom_point(colour = "grey50") + 
  stat_smooth(method = "loess", span = .2, 
              colour = "red") +
  stat_smooth(method = "loess", span = 2, 
              colour = "blue", linetype = 2)

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

下面我们从二变量推广到多变量，研究多变量的可视化方法。下面这段代码模拟了一个因变量，它是三个预测变量的函数，其中两个是连续变量，第三个是分类预测量。最初的二变量图形只显示x与y的关系，它说明两者之间存在弱关系，而且可能存在几个奇异点。

set.seed(1234)
d3 <- data.table(
  x = rnorm(500),
  w = rnorm(500), 
  z = rbinom(500, 1, .4)
)
d3[, y := rnorm(500, mean = 3 + 
                  ifelse(z<0&w<0,-2,0) *x +
                  ifelse(x<0, 0, 2) * w * x^2 + 4 * z * w,
                sd = 1)]

ggplot(d3, aes(x,y)) + 
  geom_point(colour = "grey50") + 
  stat_smooth(method = "loess", colour = "blue")

## `geom_smooth()` using formula 'y ~ x'

虽然数据可视化通常只限于完全连续变量的二维空间，但是我们可以使用不同的形状、颜色和子图增加维数。在下面的代码中，再次分析因变量y的预测量。现在要包括所有变量。二值变量z用来表示点和线。我们利用前面的技术把连续值w分割为四组数据。虽然分割w不能完全表示存在的连续关系，但也足以告诉我们正在执行的操作，以及控制w的值是否有必要。

ggplot(d3, aes(x, winsorizor(y, .01), colour = z)) + 
  geom_point() + 
  stat_smooth(method = "loess") + 
  facet_wrap(~ cut(w, quantile(w), include.lowest = TRUE))

## `geom_smooth()` using formula 'y ~ x'

从上图中大概可以看出，当w比较小且x小于0时，y与x之间存在负的线性关系。当w比较大时，就没有这种关系。z值会引起回归曲线的上下相对移动，这具体取决于w的值。

Reference

Wiley, M., & Wiley, J. F. (2019). Advanced R Statistical Programming and Data Models. Berkeley, CA, USA: Apress.