Ясно, что коррелограммы, да и вообще анализ корреляций — не лучший способ анализа данных. Но тем не менее, в качестве разведывательного анализа или как способ отбора переменных для последующего анализа, мне кажется, он может быть полезен. Собственно, в этом посте — информация про то, как построить симпатичные коррелограммы в R с использованием ggplot2. Результат не идеальный, но при желании его вполне можно допилить до того вида, который вам хочется.




В качестве примера будут использованы данные по корреляции между характеристиками текста из недавнего исследования по феномену доверия в блогах. Массив данных состоит из 39 текстов, характеризующихся 5 параметрами.

1
2
3
4
5
6
7
8
> mydata
   word_length words_per_sentence stops_per_text text_length avg_freq
1       5.4538            12.2333         0.0328        2469  10.5943
2       5.4759            11.0312         0.0447        2394   8.1360
3       4.9102             8.4872         0.0423        2082  10.1353
4       4.8058            13.4348         0.0406        1871   9.2711
5       5.0888             8.9412         0.0551        1977  11.2576
...

Для начала нам необходима матрица корреляций. Было бы странно строить коррелограмму без нее.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
> mydata<-as.matrix(mydata)
> mydata
      word_length words_per_sentence stops_per_text text_length avg_freq
 [1,]      5.4538            12.2333         0.0328        2469  10.5943
 [2,]      5.4759            11.0312         0.0447        2394   8.1360
 [3,]      4.9102             8.4872         0.0423        2082  10.1353
 [4,]      4.8058            13.4348         0.0406        1871   9.2711
...
> c <- rcorr(mydata)$r
> c
                   word_length words_per_sentence stops_per_text text_length    avg_freq
word_length          1.0000000          0.2479438    -0.51051790  0.21495645 -0.50763184
words_per_sentence   0.2479438          1.0000000    -0.23731066  0.09888110 -0.06482670
stops_per_text      -0.5105179         -0.2373107     1.00000000 -0.09354671  0.36629570
text_length          0.2149564          0.0988811    -0.09354671  1.00000000 -0.09714153
avg_freq            -0.5076318         -0.0648267     0.36629570 -0.09714153  1.00000000
 
> p <- rcorr(mydata)$P
> p
                    word_length words_per_sentence stops_per_text text_length     avg_freq
word_length                  NA          0.1280325   0.0008983421   0.1887955 0.0009714043
words_per_sentence 0.1280324892                 NA   0.1457634982   0.5492450 0.6949963179
stops_per_text     0.0008983421          0.1457635             NA   0.5710972 0.0218222894
text_length        0.1887955376          0.5492450   0.5710971987          NA 0.5563274677
avg_freq           0.0009714043          0.6949963   0.0218222894   0.5563275           NA

Первая команда делает из dataframe матрицу, поскольку следующая команда, rcorr, требует как раз матриц. Потом мы записываем в переменные c и p, соответственно, матрицу корреляций и матрицу p-уровней этих корреляций.
Теперь, я хочу. чтобы значимые корреляции на графике были отмечены звездочками на основе значений p. Для этого используем команду symnum. Она вместо значений чисел подставляет соответствующие обозначения.

1
2
3
> stars <- as.character(symnum(p, cutpoints=c(0,0.001,0.01,0.05,0.1,1),symbols=c('***', '**', '*', '.','' ), legend=F)) 
> stars
 [1] "?"   ""    "***" ""    "***" ""    "?"   ""    ""    ""    "***" ""    "?"   ""    "*"   ""    ""    ""    "?"   ""    "***" ""    "*"   ""    "?"

Затем мы делаем из корреляций, значений p и звездочек новый dataframe.

1
2
3
4
5
6
7
> molten.data <- cbind(melt(c), stars, melt(p)[,3]) 
> molten.data
                   X1                 X2       value stars melt(p)[, 3]
1         word_length        word_length  1.00000000     ?           NA
2  words_per_sentence        word_length  0.24794379       0.1280324892
3      stops_per_text        word_length -0.51051790   *** 0.0008983421
...

Команда melt переводит dataframe из широкой формы, где отдельные переменные представлены отдельными столбцами, в длинную, где каждая строка содержит три переменных: идентификатор наблюдения, идентификатор переменной и значение переменной, например:

1
2
3
4
5
6
> melt(c)
                   X1                 X2       value
1         word_length        word_length  1.00000000
2  words_per_sentence        word_length  0.24794379
3      stops_per_text        word_length -0.51051790
...

Здесь X1 идентификатор наблюдения (выбран автоматически как первый столбец в c), X2 — идентификатор переменной (названия оставшихся столбцов в c), value — значение на пересечении X1 и X2 в исходных данных.
Команда names записывает вместо невнятных заголовков столбцов в molten.data нормальные.

1
2
3
4
5
6
7
8
> names(molten.data) <- c("M1", "M2", "corr", "stars","p")   
> molten.data
                   M1                 M2        corr stars            p
1         word_length        word_length  1.00000000     ?           NA
2  words_per_sentence        word_length  0.24794379       0.1280324892
3      stops_per_text        word_length -0.51051790   *** 0.0008983421
4         text_length        word_length  0.21495645       0.1887955376
...

Едем дальше. Из всего набора корреляций нам нужен только нижний треугольник матрицы, поскольку матрица корреляций симметрична и отображать ее целиком нет смысла.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
> mi.lower <- subset(molten.data[lower.tri(c),], M1 != M2)
 
> mi.lower
                   M1                 M2        corr stars            p
2  words_per_sentence        word_length  0.24794379       0.1280324892
3      stops_per_text        word_length -0.51051790   *** 0.0008983421
4         text_length        word_length  0.21495645       0.1887955376
5            avg_freq        word_length -0.50763184   *** 0.0009714043
8      stops_per_text words_per_sentence -0.23731066       0.1457634982
9         text_length words_per_sentence  0.09888110       0.5492449543
10           avg_freq words_per_sentence -0.06482670       0.6949963179
14        text_length     stops_per_text -0.09354671       0.5710971987
15           avg_freq     stops_per_text  0.36629570     * 0.0218222894
20           avg_freq        text_length -0.09714153       0.5563274677

Кроме того, для отображения названий в центральной диагонали мы используем еще один кусок molten.data.

1
2
3
4
5
6
7
8
9
> mi.ids <- subset(molten.data, M1 == M2)
 
> mi.ids
                   M1                 M2 corr stars  p
1         word_length        word_length    1     ? NA
7  words_per_sentence words_per_sentence    1     ? NA
13     stops_per_text     stops_per_text    1     ? NA
19        text_length        text_length    1     ? NA
25           avg_freq           avg_freq    1     ? NA

Смысл этого не совсем ясного действа станет понятен чуть позже.
Теперь перейдем собственно к графику.

1
2
3
> library(ggplot2)
 
> plot1<-ggplot(mi.lower, aes(M1, M2, fill=corr))+  scale_fill_gradient(low=rgb(34/255, 94/255, 168/255),high=rgb(235/255, 59/255, 44/255))   + theme_bw() + geom_tile()

Используем нижний треугольник матрицы корреляций (mi.lower),

  1. aes(M1,M2, fill=corr) — в качестве шкалы x используем M1, шкалы y — M2, заливка цветом (fill) задается переменной corr
  2. scale_fill_gradient — заливка цветом меняется от low до high, цвета я задавал через значения rgb, т.е. шкалы «красный, зеленый, синий» со значениями от 0 до 1;
  3. theme_bw() — используем черно-белую тему графиков
  4. geom_tile() — соотношение между x и y задается прямоугольными клетками на пересечении x и y (подробнее про geom_tile)

Пока получается какая-то хрень, это нормально.


Основная проблема — данные не отсортированы, точнее отсортированы, но не так как нам нужно. Решаем ее ручной сортировкой:

1
2
> meas <- as.character(unique(molten.data$M2))
> plot1 <- plot1+scale_x_discrete(limits=meas[length(meas):1]) +  scale_y_discrete(limits=meas)

Уже лучше.

Теперь уберем все лишние подписи и отметки, добавим заголовок.

1
2
3
4
5
6
7
8
> none <-  theme_blank()
> p2<-plot1 + xlab(NULL) + ylab(NULL) + 
     opts(axis.text.x=none) + 
     opts(axis.text.y=none) + 
     opts(axis.ticks=none) + 
     opts(panel.border=none) +
     opts(legend.position='none') +
     opts(title="Показатели удобочитаемости")

Почти готово.

Осталось добавить обозначения:

1
> plot1 <-plot1 + geom_text(data=mi.lower, aes(label=paste(round(corr,2), stars)),size=4, color="white")

paste(round(corr,2), stars) объединяет столбцы corr с округлением до сотых и stars

1
> plot1 <- plot1 + geom_text(data=mi.ids, color="gray40", size=4,aes(label=M2),angle=45,hjust=0.15,vjust=0)

Наконец-то понадобились mi.ids. С их помощью мы на пересечении M1 и M2 добавляем текстовые метки. angle поворачивает надписи на 45 градусов, hjust и vjust слегка модифицируют их положение.

Если вы заметили, после последнего шага изменились цвета. Это связано с тем, что мы добавили к данным mi.ids, где значение corr равно 1, т.е. если раньше градиент заливки был от -0.51 до 0.37, то теперь стал от -0.51 до 1. Можно от этого избавиться, можно так оставить.

1
> plot1+scale_fill_continuous(limits=c(-0.52,0.37))

В итоге получается следующее: