15.2 相关系数和协方差
2018-12-03 本文已影响0人
砖头机的灵感
首稿于2018-12-03
> #install.packages("GGally")
> #chooseCRANmirror()
> #install.packages("RXKCD")
当处理一个以上的变量时，需要测度它们之间的相互关系。两个简单直接的方法是求其相关系数和协方差。先看看ggplot2中的数据集economics。
> require(ggplot2)
> head(economics)
## cor函数计算其相关关系
> cor(economics$pce,economics$psavert)
## 或者根据公式计算
> xPart <- economics$pce - mean(economics$pce)
> yPart <- economics$psavert - mean(economics$psavert)
> nMinusOne <- (nrow(economics) - 1)
> xSD <- sd(economics$pce)
> ySD <- sd(economics$psavert)
> sum(xPart * yPart) / (nMinusOne * xSD* ySD)
##同时比较多个变量，可以对矩阵（必须为数值型变量）用函数cor
> cor(economics[,c(2,4:6)])
## 骚气的用GGally包展示每一个变量和其他变量散点图。
## 这里不是通过require去加载包，而是用 :: 直接调用它包内的函数，这个操作符号可以无须加载包就可以调用其中的函数。
> GGally::ggpairs(economics[,c(2,4:6)])
数据集economics中各对变量之间的散点图并显示其相关系数
## 上图是根据原始数据绘制的，没有直接明了地展现各个变量间的相关系数。
## 为此，绘制相关变量的热图。
## 较高的正相关系数，意味着两变量有密切的正向关系，具有较低的负相关系数，意味着两变量有密切的反向关系，接近于0的相关系数意味着两者没有密切的相关关系。
> require(reshape2)
> require(scales)
> econCor <- cor(economics[,c(2,4:6)])
> econCor
> econMelt <- melt(econCor, varnames=c("x","y"),value.name="Correlation")
> econMelt <- econMelt[order(econMelt$Correlation),]
> econMelt
> ggplot(econMelt, aes(x=x,y=y)) +
geom_tile(aes(fill=Correlation)) +
scale_fill_gradient2(low=muted("red"), mid="white",high="steelblue",guide=guide_colorbar(ticks=FALSE, barheight=10), limits=c(-1,1)) +
theme_minimal() + labs(x=NULL,y=NULL)
数据集economics中变量的相关系数热图。对角线上的图形显示相关系数为1，这是因为每个变量都和自身完全相关。红色代表高度负相关，蓝色代表高度正相关，白色代表没有相关宣传
## 有缺失值时的处理方法
## 函数cor不用设置参数na.rm来除去缺失值，而是应用“all.obs”、“complete.obs”、
## “pairwise.complete.obs”、“everything” 或者 “na.or.complete”其中一种方法来
## 计算多个变量的相关关系。
> m <- c(9,9,NA,3,NA,5,8,1,10,4)
> n <- c(2,NA,1,6,6,4,1,1,6,7)
> p <- c(8,4,3,9,10,NA,3,NA,9,9)
> q <- c(10,10,7,8,4,2,8,5,5,2)
> r <- c(1,9,7,6,5,6,2,7,9,10)
> theMat <- cbind(m,n,p,q,r)
> theMat
cor(theMat, use = "all.obs")
> Error in cor(theMat, use = "all.obs") : cov/cor中有遗漏值
> cor(theMat, use = "na.or.complete")
m n p q r
m 1.0000000 -0.5228840 -0.2893527 0.2974398 -0.3459470
n -0.5228840 1.0000000 0.8090195 -0.7448453 0.9350718
p -0.2893527 0.8090195 1.0000000 -0.3613720 0.6221470
q 0.2974398 -0.7448453 -0.3613720 1.0000000 -0.9059384
r -0.3459470 0.9350718 0.6221470 -0.9059384 1.0000000
> cor(theMat[c(1,4,7,9,10),])
m n p q r
m 1.0000000 -0.5228840 -0.2893527 0.2974398 -0.3459470
n -0.5228840 1.0000000 0.8090195 -0.7448453 0.9350718
p -0.2893527 0.8090195 1.0000000 -0.3613720 0.6221470
q 0.2974398 -0.7448453 -0.3613720 1.0000000 -0.9059384
r -0.3459470 0.9350718 0.6221470 -0.9059384 1.0000000
> identical(cor(theMat, use = "complete.obs"), cor(theMat[c(1,4,7,9,10), ]))
[1] TRUE
> cor(theMat, use = "pairwise.complete.obs")
m n p q r
m 1.00000000 -0.02511812 -0.3965859 0.4622943 -0.2001722
n -0.02511812 1.00000000 0.8717389 -0.5070416 0.5332259
p -0.39658588 0.87173889 1.0000000 -0.5197292 0.1312506
q 0.46229434 -0.50704163 -0.5197292 1.0000000 -0.4242958
r -0.20017222 0.53322585 0.1312506 -0.4242958 1.0000000
> cor(theMat[, c("m","p")], use = "pairwise.complete.obs")
m p
m 1.0000000 -0.3965859
p -0.3965859 1.0000000
> cor(theMat, use = "everything")
> cor(theMat, use = "all.obs")
> cor(theMat, use ="complete.obs")
> cor(theMat, use = "na.or.complete")
> cor(theMat[c(1,4,7,9,10),])
> identical(cor(theMat, use = "complete.obs"), cor(theMat[c(1,4,7,9,10), ]))
> cor(theMat, use = "pairwise.complete.obs")
> cor(theMat[, c("m","p")], use = "pairwise.complete.obs")
> data(tips, package = "reshape2")
> head(tips)
> GGally::ggpairs(tips)
用tips中连续变量和分类变量绘制出的ggpairs图
## 相关性并不意味着存在因果关系。两个变量有相关关系也不一定说明它们是相互影响的。
## 来自动漫深深的鄙视
> require(RXKCD)
> require(RXKCD)
> getXKCD(which = "552")
15-3-RXKCD.png
> cov(economics$pce, economics$psavert)
> cov(economics[, c(2,4:6)])
## cov 和 cor差别在公式上。cor多除了两个标准差
> identical(cov(economics$pce, economics$psavert),
cor(economics$pce, economics$psavert) * sd(economics$pce) * sd(economics$psavert))
[1] TRUE
> savehistory("practice-15.R")
加料：
cor函数中四种方法解读：
1 everthing， 它意味着所有列的元素必须不含缺失值，否则结果是NA。
2 all.obs， 同样要求所有列不含缺失值否者提示错误。
3,4 complete.obs 和 na.or.complete 处理方法类似，它们仅留下不存在缺失值的行。两者的不同在于，如果经过处理后没有一个具有完整数据的行，complete.obs 会发返回错误提示，而 na.or.complete 会反回NA。
5 pairwise.complete.obs 用途更加广泛。它依次比较多对变量，并把两个变量相互之间的缺失列剔除，用余下的数据计算两者的相关系数。这种方法在本质上是与complete.obs在计算每一对变量组合的相关系数是相同的。
