correlation analysis
baorui
Here is correlation analysis:
pairs(~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width,data=iris, main="Simple Scatterplot Matrix")
library(car)
scatterplotMatrix(~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width|Species, data=iris)
# 3d 图形无法展示--略过
# library(rgl)
# scatter3d(iris$Sepal.Length, iris$Petal.Length, iris$Petal.Width) # install.packages('corrgram')
library(corrgram)
#1、设置排序处理
corrgram(mtcars,order=TRUE)
#2、设置上下三角面板形状
corrgram(mtcars,order=TRUE,lower.panel=panel.shade,upper.panel=panel.pie)
#3、只显示下三角部分
corrgram(mtcars,order=TRUE,lower.panel=panel.shade,upper.panel=NULL)
#4、调整面板颜色
corrgram(mtcars,order=TRUE,lower.panel=panel.shade,upper.panel=panel.pie,
col.regions=colorRampPalette(c("darkgoldenrod4","burlywood1","white",
"darkkhaki","darkgreen")))
# install.packages('corrplot')
library(corrplot)
#1、使用不同的method绘制相关矩阵图
# 提取iris的前4个数值列，并进行标准化处理
data0=scale(iris[1:4])
# 计算这4个变量的协方差，由于经过标准化处理，这样得到的也是相关系数
M=cov(data0)
methods <- c("circle","square","ellipse","pie","shade","color")
par(mfrow=c(2,3))
t0 <- mapply(function(x){corrplot(M, method=x, order="AOE")},methods)## Warning in acos(rho): 产生了NaNs
## Warning in acos(rho): 产生了NaNs
par(mfrow=c(1,1))
#2、设置method=color绘制热力矩阵图
corrplot(cor(mtcars), method="color", order = "AOE",tl.col="black",tl.srt=45,
addCoef.col="black",col=colorRampPalette(c("#7F0000","red","#FF7F00",
"yellow","white", "cyan", "#007FFF", "blue","#00007F"))(20))
#3、绘制上下三角及不同色彩的相关矩阵图
library(RColorBrewer)
par(mfrow=c(2,2))
corrplot(cor(mtcars),type="lower")
corrplot(cor(mtcars),type="lower",order="hclust",
col=brewer.pal(n=8,name="RdYlBu"))
corrplot(cor(mtcars),type="upper",order="AOE",
col=c("black","white"),bg="lightblue")
corrplot(cor(mtcars),type="upper",order="FPC",
col=brewer.pal(n=8, name="PuOr"))
par(mfrow=c(1,1)) d<-sqrt(1-cor(mtcars)^2)
hc<-hclust(as.dist(d))
plot(hc)
rect.hclust(hc,k=3)
# install.packages('pvclust')
library(pvclust)
cluster.bootstrap <- pvclust(mtcars, nboot=1000, method.dist="correlation")## Bootstrap (r = 0.5)... Done.## Warning: inappropriate distance matrices are omitted in computation: r =
## 0.5## Bootstrap (r = 0.59)... Done.
## Bootstrap (r = 0.69)... Done.
## Bootstrap (r = 0.78)... Done.
## Bootstrap (r = 0.88)... Done.
## Bootstrap (r = 1.0)... Done.
## Bootstrap (r = 1.09)... Done.
## Bootstrap (r = 1.19)... Done.
## Bootstrap (r = 1.28)... Done.
## Bootstrap (r = 1.38)... Done. plot(cluster.bootstrap)
pvrect(cluster.bootstrap)
#1、提取iris的前4个数值列，并进行标准化处理
data0=scale(iris[1:4])
head(data0)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## [1,] -0.8976739 1.01560199 -1.335752 -1.311052
## [2,] -1.1392005 -0.13153881 -1.335752 -1.311052
## [3,] -1.3807271 0.32731751 -1.392399 -1.311052
## [4,] -1.5014904 0.09788935 -1.279104 -1.311052
## [5,] -1.0184372 1.24503015 -1.335752 -1.311052
## [6,] -0.5353840 1.93331463 -1.165809 -1.048667 #2、计算这4个变量的协方差，由于经过标准化处理，这样得到的也是相关系数
M=cov(data0)
M## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 1.0000000 -0.1175698 0.8717538 0.8179411
## Sepal.Width -0.1175698 1.0000000 -0.4284401 -0.3661259
## Petal.Length 0.8717538 -0.4284401 1.0000000 0.9628654
## Petal.Width 0.8179411 -0.3661259 0.9628654 1.0000000 #3、将M进行分块，1:2两个变量一组，3:4是另外一组，并进行两两组合
X11=M[1:2,1:2]
X12=M[1:2,3:4]
X21=M[3:4,1:2]
X22=M[3:4,3:4]
#4、按公式求解矩阵A和B
A=solve(X11)%*%X12%*%solve(X22)%*%X21
A## Sepal.Length Sepal.Width
## Sepal.Length 0.7308471 -0.3671646
## Sepal.Width -0.3012184 0.1699359 B=solve(X22)%*%X21%*%solve(X11)%*%X12
B## Petal.Length Petal.Width
## Petal.Length 1.3061861 1.1926184
## Petal.Width -0.4554116 -0.4054031 #5、使用eigen函数求解典型相关系数如下
eV=sqrt(eigen(A)$values)
eV## [1] 0.9409690 0.1239369 #6、进行验证
#...比较A与XΛX^(-1)是否相等
round(A-eigen(A)$vectors%*%diag(eigen(A)$values)%*%solve(eigen(A)$vectors),3)## Sepal.Length Sepal.Width
## Sepal.Length 0 0
## Sepal.Width 0 0 ## Sepal.Length Sepal.Width
## Sepal.Length 0 0
## Sepal.Width 0 0
#...比较B与YΛY^(-1)是否相等
round(B-eigen(B)$vectors%*%diag(eigen(B)$values)%*%solve(eigen(B)$vectors),3)## Petal.Length Petal.Width
## Petal.Length 0 0
## Petal.Width 0 0 #...求解A对应的特征向量并计算典型向量C1
C1 <- data0[,1:2]%*%eigen(A)$vectors
head(C1)## [,1] [,2]
## [1,] -1.2214121 0.4937300
## [2,] -0.9989091 -0.6371437
## [3,] -1.3995557 -0.3391746
## [4,] -1.4218362 -0.5984308
## [5,] -1.4217354 0.6427146
## [6,] -1.2435919 1.4756191 #...验证C1对应各变量的标准差是否为1，同时查看均差
apply(C1,2,sd)## [1] 1.041196 0.951045 ## [1] 1.041196 0.951045
apply(C1,2,mean)## [1] -4.880321e-16 -2.759430e-17 ## [1] -4.880321e-16 -2.759430e-17
#...由于均值为0，标准差不为1，这里对特征向量进行伸缩变换
eA <- eigen(A)$vectors%*%diag(1/apply(C1,2,sd))
eA## [,1] [,2]
## [1,] 0.8851871 0.4800626
## [2,] -0.3726619 0.9354889 #...再次验证方差和均值
C1 <- data0[,1:2]%*%eA
apply(C1,2,sd)## [1] 1 1 ## [1] 1 1
apply(C1,2,mean)## [1] -4.667693e-16 -2.745503e-17 ## [1] -4.667693e-16 -2.745503e-17
#...可见，特征向量已经满足要求，同理对B可得
C2 <- data0[,3:4]%*%eigen(B)$vectors
apply(C2,2,sd)## [1] 0.6291236 0.2003530 ## [1] 0.6291236 0.2003530
apply(C2,2,mean)## [1] -1.403572e-17 -9.859870e-18 ## [1] -1.403572e-17 -9.859870e-18
eB <- eigen(B)$vectors%*%diag(1/apply(C2,2,sd))
C2 <- data0[,3:4]%*%eB
apply(C2,2,sd)## [1] 1 1 ## [1] 1 1
apply(C2,2,mean)## [1] -1.598186e-17 5.307097e-17 round(cor(cbind(C1,C2)),3)## [,1] [,2] [,3] [,4]
## [1,] 1.000 0.000 0.941 0.000
## [2,] 0.000 1.000 0.000 0.124
## [3,] 0.941 0.000 1.000 0.000
## [4,] 0.000 0.124 0.000 1.000 x <- as.matrix(iris[,1:2])
y <- as.matrix(iris[,3:4])
cancor(x,y)## $cor
## [1] 0.9409690 0.1239369
##
## $xcoef
## [,1] [,2]
## Sepal.Length -0.08757435 0.04749411
## Sepal.Width 0.07004363 0.17582970
##
## $ycoef
## [,1] [,2]
## Petal.Length -0.06956302 -0.1571867
## Petal.Width 0.05683849 0.3940121
##
## $xcenter
## Sepal.Length Sepal.Width
## 5.843333 3.057333
##
## $ycenter
## Petal.Length Petal.Width
## 3.758000 1.199333