置信椭圆与R画法
2020-08-14 本文已影响0人
小潤澤
置信椭圆
在二维或者更高维的空间里，数据的聚类往往需要添加一个“置信区间”。仿照一维空间的数据，置信区间往往相对于点估计而来的，在统计学中，一个概率样本的置信区间（Confidence interval）是对这个样本的某个总体参数的区间估计。置信区间展现的是这个参数的真实值有一定概率落在测量结果的周围的程度，其给出的是被测量参数的测量值的可信程度，一般我们用95%置信区间来表示。
那么二维空间的数据也是如此，二维空间的置信区间往往利用置信椭圆来描述，展现真实值的可信程度，一般我们用95%置信椭圆来表示。
为什么是置信椭圆呢，往往是因为我们的而数据分布形似椭圆
置信椭圆绘制原理
本文若无特殊说明，椭圆的中心都在左边原点，与图片略微不符
(1).数据中心在原点且轴对齐型
我们先看最简单的，轴对齐型的，比方说椭圆长轴平行于x轴
平行x轴
假设我有两列随机变量 x 和 y
对于该二维数据，x和y，我们不妨计算下x和y的协方差矩阵：
二维协方差矩阵计算公式：
摘自百度
由于我们的数据，形似椭圆，并且椭圆长轴平行于x轴，那么协方差矩阵中x轴方向的方差为8.4213，y轴方向的方差为0.9387。由于x方向与y方向正交，所以x和y这两个随机变量的协方差为0，也可以理解为相关性为0。
由于我们的数据点分布形似椭圆，所以我们定义置信区间就把它定义为椭圆，我们知道椭圆的标准方程为：
那么我们不妨将椭圆的半长轴a定义为σx，半短轴b定义为σy，那么有：
显然，该椭圆的数据中心为(0,0)。
其中s定义椭圆的规模，可以是任意的数（例如，s=1）。现在的问题是如何选择s，使得所得到的椭圆规模代表我们所选择的置信水平（例如，95％的置信水平对应于s =5.991）。
这个s=5.991可以通过卡方分布（随机变量平方和）进行计算（假设随机变量x和y服从高斯分布），并且自由度为2.
联想到卡方分布，上式左侧可以看作为随机变量的平方和，通过查询卡方表，有：
我们知道当 s = 5.991 时，P(s < 5.991) = 0.95，也就是说，我们要寻找一个在椭圆上的点，使得 s = 5.991，那么将 s = 5.991 带入椭圆方程：
这个等式表明数据点小于5.991的概率为95%（在椭圆内表示数据点小于5.991）。
那么，s=5.991即为95%置信区间，表示有一组数据，其x方向方差为σx，y方向方差为σy时，如果s=5.991，那么95%的数据点在该椭圆内。
那么置信椭圆的长轴2a：
置信椭圆的短轴2b：
(2).数据中心在远点且任意方向的置信椭圆
我们首先计算随机变量x和y的协方差矩阵，根据协方差矩阵计算该矩阵的特征向量（上图绿色为椭圆短轴方向的特征向量；粉红色为椭圆长轴方向的特征向量），我们定义椭圆长轴方向的特征向量为v1，椭圆短轴方向的特征向量为v2；而特征向量大小为特征值，我们定义椭圆长轴方向的特征值为λ1，椭圆短轴方向的特征值为λ2。
此时椭圆的长轴为：
椭圆的短轴为：
并且定义椭圆长轴与x轴正方向夹角为α：
那么，我们就可以根据角度来对置信椭圆进行“旋转”。
(3).任意位置置信椭圆
有了上面的介绍，那么任意位置的椭圆无非是满足数据中心不在原点，且椭圆长轴与x轴正方向夹角为任意角α。并以此来建立椭圆方程，画出置信椭圆（可见数据中心化可以简便很多计算）
R绘制置信椭圆
如果是利用ggplot画置信椭圆，例如对PCA或者PCoA降维后的二维结果进行画图
我们利用stat_ellipse()画置信椭圆
pcoa_plot <- ggplot(pcoa, aes(PCoA1, PCoA2, group = group, color = group)) +
geom_point() +
stat_ellipse(level = 0.95, show.legend = F) +
annotate('text', label = 'A', x = -0.3, y = 0.06, size = 5, colour = '#f8766d') +
annotate('text', label = 'B', x = 0.1, y = 0.35, size = 5, colour = '#00ba38') +
annotate('text', label = 'C', x = 0.3, y = -0.13, size = 5, colour = '#619cff')
或者利用一些R包“car”来画
pdf(" ")
dataEllipse(data[,1],data[,2],levels=0.95,col="black",xlim=c(-1,1),ylim=c(-4,2),xlab="DNMT3A",ylab="PGLYRP2",lwd=1,
grid=F,pch=20,lty=2)
dev.off()
