Корреляционные матрицы

Пример функции, позволяющей создавать корреляционные матрицы и оценивать уровни статистической значимости коэффициентов корреляции. Предусмотрено вычисление объёма совокупности при попарном исключении вариант.

corrsign <- function(x){ 
  require(Hmisc) 
  x <- as.matrix(x)
  r_coeff <- rcorr(x,type=c("pearson"))$r #коэффициенты корреляции Пирсона
  p_lev <- rcorr(x, type=c("pearson"))$P  #уровни статистической значимости
  n <- rcorr(x, type=c("pearson"))$n    #n при попарном исключении
  sr<-sqrt((1-r_coeff^2)/(n-2)) #стандартная ошибка для r
  sr <- format(round(cbind(rep(-1.111, ncol(x)), sr), 3))[,-1]  #округление ошибки
  stars <- ifelse(p_lev <= .001,"***", (ifelse (p_lev <= .01, "**", 
                                                ifelse(p_lev <= .05, "*", "")))) 
  r_coeff <- format(round(cbind(rep(-1.111, ncol(x)), r_coeff), 3))[,-1]    #округление для r
  r_coeff_new <- ifelse (is.infinite(r_coeff),,
                         matrix(paste(r_coeff,"±",sr,stars,"(",n,")", sep=""), ncol=ncol(x)))
  diag(r_coeff_new) <- paste(diag(r_coeff),"", sep="-")
  rownames(r_coeff_new) <- colnames(x)
  colnames(r_coeff_new) <- paste(colnames(x),"", sep="")
  r_coeff_new <- as.data.frame(r_coeff_new)
  return(r_coeff_new)   #На выходе - таблица корреляций с уровнями значимости
}

#Пирсоновские корреляции по данным вида "Setosa" базы данных "iris"
a<-subset(iris,Species=="setosa")
corrsign(a[1:3])

## Loading required package: Hmisc
## Loading required package: grid
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## Loading required package: ggplot2
## 
## Attaching package: 'Hmisc'
## 
## Следующие объекты скрыты от 'package:base':
## 
##     format.pval, round.POSIXt, trunc.POSIXt, units

##                      Sepal.Length          Sepal.Width      Petal.Length
## Sepal.Length               1.000-  0.743± 0.097***(50)  0.267± 0.139(50)
## Sepal.Width   0.743± 0.097***(50)               1.000-  0.178± 0.142(50)
## Petal.Length     0.267± 0.139(50)     0.178± 0.142(50)            1.000-

#Пирсоновские корреляции по данным вида "Versicolor" базы данных "iris"
a<-subset(iris,Species=="versicolor")
corrsign(a[1:3])

##                      Sepal.Length          Sepal.Width      Petal.Length
## Sepal.Length               1.000-  0.526± 0.123***(50)  0.754± 0.095***(50)
## Sepal.Width   0.526± 0.123***(50)               1.000-  0.561± 0.120***(50)
## Petal.Length  0.754± 0.095***(50)  0.561± 0.120***(50)               1.000-

#Пирсоновские корреляции по данным вида "Versicolor" базы данных "iris"
a<-subset(iris,Species=="virginica")
corrsign(a[1:3])

##                      Sepal.Length          Sepal.Width      Petal.Length
## Sepal.Length               1.000-  0.457± 0.128***(50)  0.864± 0.073***(50)
## Sepal.Width   0.457± 0.128***(50)               1.000-  0.401± 0.132**(50)
## Petal.Length  0.864± 0.073***(50)   0.401± 0.132**(50)              1.000-