Корреляция и взаимосвязь номинальных признаков

[ymr3R9Jge]

Подскажите, пожалуйста, - в какую сторону копать!
Есть таблица данных по метеорологическим наблюдениям такого вида:
num;date;Баштанка;Бобринець;Вел.Олександрівка;Комісарівка;Лозоватка;Лошкарівка
1;01/01/2009/02/00;NW;NA;calm;calm;SW;calm
2;01/01/2009/08/00;WNW;WSW;calm;calm;S;SSW
3;01/01/2009/14/00;WSW;WSW;W;SW;SSW;SW
4;01/01/2009/20/00;WSW;WSW;calm;SSW;NA;SSW
5;01/01/2010/02/00;NW;NNW;calm;WNW;NNW;NW
6;01/01/2010/08/00;NNE;ENE;calm;N;NE;N
7;01/01/2010/14/00;E;ESE;E;E;E;E
8;01/01/2010/20/00;E;E;E;E;E;E
9;01/01/2011/02/00;SW;calm;calm;SSW;SW;S
10;01/01/2011/08/00;SW;WSW;calm;S;SW;SSW
В таблице каждый столбец по метеостанции содержит значения 16-ти румбов, значение "calm" (штиль), и попадаются иногда значения "NA" (пустые). Всего в таблице более 17 тыс. записей. Нужно оценить взаимосвязь направления ветра на разных метеостанциях в одномоментных измерениях и получить некий коэффициент "корреляции" для каждой пары метеостанций. Я предполагаю, что можно воспользоваться таблицами сопряжённости, но как при этом учитывать значения времени???

[gamm]

в основном Хи-квадрат, но мера ненормированная тыц, тыц,тыц

еще можно кросс-энтропию посчитать, тыц

[ymr3R9Jge]

gamm, огромное спасибо за ссылки! То, что было нужно я по ним нашёл!
Теперь по существу: всё решается намного проще. В пакете vcd реализован расчёт как ненормированных (критерий хи-квадрат, отношение правдоподобия) так и нормированных (коэффициент сопряжённости, коэффициент V Крамера). Последовательность следующая:

Код: Выделить всё

tmp.table <- table(my.data.frame$Бобринець, my.data.frame$Лозоватка))#создаём таблицу соответствия для двух переменных
library(vcd)#загружаем пакет vcd
summary(assocstats(tmp.table))#проводим анализ по готовой таблице сопряжённости
#и получаем ответ:
Number of cases in table: 17241 
Number of factors: 2 
Test for independence of all factors:
	Chisq = 25313, df = 256, p-value = 0
                   X^2  df P(> X^2)
Likelihood Ratio 20277 256        0
Pearson          25313 256        0

Phi-Coefficient   : 1.212 
Contingency Coeff.: 0.771 
Cramer's V        : 0.303

Низкое (а по факту - нулевое) p-value критерия Пирсона хи-квадрат указывает, что гипотезу о независимости переменных мы отбрасываем и "завимимость" переменных - статистически значимая. А значение коеффициента V Крамера указывает на наличие слабой связи между переменными.

[gamm]

(1)Низкое (а по факту - нулевое) p-value критерия Пирсона хи-квадрат указывает, что гипотезу о независимости переменных мы отбрасываем и "завимимость" переменных - статистически значимая.
(2)А значение коеффициента V Крамера указывает на наличие слабой связи между переменными.

я бы на реальных данных на асимптотические критерии не полагался, а устроил таки бутстрапчик для очистки совести. Тем более данных много. И значение Хи-кадрат очень большое - на табличку сопряженности глянуть бы ...

И вообще, построил бы модельку, с общей средней вероятностью направлений и случайным эффектом в виде станции, и отклонениями от нее по станциям. Со всеми опциями - средним трендом, сезонностью, автокорреляцией, и т.д. - посмотреть, что дает различие (нет ли аномалий). Для спокойствия и лучшего понимания структуры разницы

Например, вот такая блочно- независимая конструкция

Код: Выделить всё

>> library(vcd)
> my.df<-as.data.frame(cbind(c(sample(1:9+0,8600,replace=T), sample(1:9+9,8600,replace=T)),
+                            c(sample(1:9+0,8600,replace=T), sample(1:9+9,8600,replace=T)) ))
>                           
> names(my.df)<-c("Dir1","Dir2")
> head(my.df)
  Dir1 Dir2
1    1    4
2    3    6
3    5    1
4    6    5
5    9    3
6    2    5
> print(tmp.table <- table(my.df$Dir1,my.df$Dir2))
    
       1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
  1  109  95 125 100 117 114 108 117 120   0   0   0   0   0   0   0   0   0
  2   99 109 102  84 119  94 115 120  92   0   0   0   0   0   0   0   0   0
  3   86 104 115  99 116 106 121 115  99   0   0   0   0   0   0   0   0   0
  4  114 100 109 118 105 113  92 116 106   0   0   0   0   0   0   0   0   0
  5   92 102 112 110 121 109 108 113 101   0   0   0   0   0   0   0   0   0
  6   95 100 102 112 103 102 108  97  95   0   0   0   0   0   0   0   0   0
  7  107  94 117  79 102 109 116  90  93   0   0   0   0   0   0   0   0   0
  8   96 111 114 119  87 109 115 116  96   0   0   0   0   0   0   0   0   0
  9  118 111  86 108 118 105  93 124 112   0   0   0   0   0   0   0   0   0
  10   0   0   0   0   0   0   0   0   0  89 110 103  86 100 114  96 109 107
  11   0   0   0   0   0   0   0   0   0 101 122 110 118 113 109 110 105 113
  12   0   0   0   0   0   0   0   0   0 111 105  96  93  86 108 119 101 108
  13   0   0   0   0   0   0   0   0   0  94 107 135 120 101 106 103 114  93
  14   0   0   0   0   0   0   0   0   0  97 117 114  96 109  96 108 107  93
  15   0   0   0   0   0   0   0   0   0 130 110 105 109 113 105 109  82  95
  16   0   0   0   0   0   0   0   0   0 118  98 109 117 111 114  90 105  91
  17   0   0   0   0   0   0   0   0   0 100 102 116 100 107  92 109 100 114
  18   0   0   0   0   0   0   0   0   0 106 104 133 115 119 103  92 122 103
> 
> summary(assocstats(tmp.table))

Number of cases in table: 17200 
Number of factors: 2 
Test for independence of all factors:
        Chisq = 17463, df = 289, p-value = 0
                   X^2  df P(> X^2)
Likelihood Ratio 23977 289        0
Pearson          17463 289        0

Phi-Coefficient   : 1.008 
Contingency Coeff.: 0.71 
Cramer's V        : 0.244 

дает статистики, близкие к вашим ... хотя в блоках имеем независимость:

Код: Выделить всё

> print(tmp.table2 <- table(my.df$Dir1[1:8600],my.df$Dir2[1:8600]))
   
      1   2   3   4   5   6   7   8   9
  1 109  95 125 100 117 114 108 117 120
  2  99 109 102  84 119  94 115 120  92
  3  86 104 115  99 116 106 121 115  99
  4 114 100 109 118 105 113  92 116 106
  5  92 102 112 110 121 109 108 113 101
  6  95 100 102 112 103 102 108  97  95
  7 107  94 117  79 102 109 116  90  93
  8  96 111 114 119  87 109 115 116  96
  9 118 111  86 108 118 105  93 124 112
> summary(assocstats(tmp.table2))

Number of cases in table: 8600 
Number of factors: 2 
Test for independence of all factors:
        Chisq = 64.24, df = 64, p-value = 0.4681
                    X^2 df P(> X^2)
Likelihood Ratio 65.577 64  0.42187
Pearson          64.238 64  0.46813

Phi-Coefficient   : 0.086 
Contingency Coeff.: 0.086 
Cramer's V        : 0.031 

[ymr3R9Jge]

Вот таблица сопряжённости для двух метеостанций:

Код: Выделить всё

> table(my.data.frame$Бобринець, my.data.frame$Лозоватка)
      
       calm   E ENE ESE   N  NE NNE NNW  NW   S  SE SSE SSW  SW   W WNW WSW
  calm  395 135 190  88 154 289 151 112  88  61  87  51  60  77  62  74  53
  E      27 310 241 150  25 153  59  10   9  20  71  22   4   4   7   4   2
  ENE    22 182 298  46  41 274  89  12   5   7  27  11   6   3   5   3   3
  ESE    23 162  76 134  16  51  26   9   3  11 124  23  10   9   1   5   2
  N      25  36  82  13 255 198 232 157  64  21  11  25  11  15   9  31  14
  NE     20  99 235  44  89 369 224  30   8   5  25  10  12   5   8   6   0
  NNE    16  36  59  14 114 173 167  43  18   5  14   7   6  10   5  11   8
  NNW    47  12  37   6 288 103 174 219 131  10  11   9   8  13  28  73   9
  NW     44   7  25   5 183  58  98 254 245   8   4   7  12  23  51 208  23
  S      31  21   6  37   8  15  11  12   5 198  80 123 117  68   8  11  23
  SE     26 121  43 178  10  25  18   8   6  55 248  94  23  11   5   7   5
  SSE    18  24  19  40   5  11   8   1   1  79 131 116  30  21   3   2  11
  SSW    34  10   6  12  12   7   6   3   6 176  34  78 152 195  25  10  67
  SW     36   8   6  16   9  11   9  10  10  83  20  43 142 273  64  43 184
  W      42   5   6   2  14   7   7  30  54   6   3   6  23  78 226 179 134
  WNW    52  13  17   6  83  26  35 130 219  11   8  11  22  38 175 362  55
  WSW    27  10   5   9   5   2   5   6  10  27   6  11  31 131 172  62 169

[gamm]

соображение такое - номинальная табличка для анализа не годится, поскольку данные ординальные. И нужно не просто сравнивать равно/не равно, а похоже/не похоже. Например, штиль - он на всех похож, я его пока выкинул.

На остальные посмотрел глазками - загнал я вашу табличку в файл foo.dat

Код: Выделить всё

Name calm   E ENE ESE   N  NE NNE NNW  NW   S  SE SSE SSW  SW   W WNW WSW
calm  395 135 190  88 154 289 151 112  88  61  87  51  60  77  62  74  53
E      27 310 241 150  25 153  59  10   9  20  71  22   4   4   7   4   2
ENE    22 182 298  46  41 274  89  12   5   7  27  11   6   3   5   3   3
ESE    23 162  76 134  16  51  26   9   3  11 124  23  10   9   1   5   2
N      25  36  82  13 255 198 232 157  64  21  11  25  11  15   9  31  14
NE     20  99 235  44  89 369 224  30   8   5  25  10  12   5   8   6   0
NNE    16  36  59  14 114 173 167  43  18   5  14   7   6  10   5  11   8
NNW    47  12  37   6 288 103 174 219 131  10  11   9   8  13  28  73   9
NW     44   7  25   5 183  58  98 254 245   8   4   7  12  23  51 208  23
S      31  21   6  37   8  15  11  12   5 198  80 123 117  68   8  11  23
SE     26 121  43 178  10  25  18   8   6  55 248  94  23  11   5   7   5
SSE    18  24  19  40   5  11   8   1   1  79 131 116  30  21   3   2  11
SSW    34  10   6  12  12   7   6   3   6 176  34  78 152 195  25  10  67
SW     36   8   6  16   9  11   9  10  10  83  20  43 142 273  64  43 184
W      42   5   6   2  14   7   7  30  54   6   3   6  23  78 226 179 134
WNW    52  13  17   6  83  26  35 130 219  11   8  11  22  38 175 362  55
WSW    27  10   5   9   5   2   5   6  10  27   6  11  31 131 172  62 169

Код: Выделить всё

p<-read.table("W:/tmp/00/foo.dat",header=T,as.is=T)
dir.names<-c("N","NNE","NE","ENE","E","ESE","SE","SSE","S","SSW","SW","WSW","W","WNW","NW","NNW")
p.tab<-as.matrix(p[-1,c(-1,-2)])
rownames(p.tab)<-colnames(p.tab)
pos<-match(dir.names,colnames(p.tab))
p.tab<-p.tab[pos,pos]
image(p.tab)

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

Главное понять, насколько изменчивым может быть значение в данных - что означает ENE в ячейке исходной таблицы - какое ему соответствует распределение вероятностей для всех значений, его же наверняка усредняли, и брали одно значение за временной интервал в 6 часов, а за это время направление наверняка менялось. И вот это распределение нужно сравнивать с тем, которое получилось (замкнув его на тор). Грубо говоря, можно корреляцию оценить (через R2). И сравнивать нужно с усредненным распределением, тут Хи-квадрат не годится ни разу, тут нужно по диагонали распределение построить, усредненное (с учетом изменчивости за срок наблюдения), и с ним все остальные станции сравнивать. Я маргиналы посмотрел, они у этих станций практически одинаковые, так что гипотеза H0: среднее распределение станций с собой (с учетом вариабельности за срок), и между станциями совпадают. Т.е. проверять нужно не независимость (чего быть не может), а одинаковость (следование одному шаблону изменения во времени, с учетом изменчивости самих данных).

[gamm]

попробую сформулировать, о чем речь.

Между станциями Y1,Y2 заведомо существует зависимость, индуцированная общими факторами, т.е. это не зависимость Y1 от Y2 (или наоборот), а зависимость их от общего фактора X (направление ветра в целом на больших территориях, которое определяется глобальными атмосферными явлениями). Поэтому обычно стремятся убрать влияние фактора X, и сравнивать Y1\X с Y2\X - в простейшем случае строится регрессия на X (или на его суррогат - в данном случае время, с учетом сроков и сезонности), и сравнивают остатки. Регрессия строится общая (одновременно для всех станций), временные ряды на станциях рассматриваются как реализации одного процесса, коррелированного во времени. Учитывая особенности данных (они ординальные), придется использовать что-то типа VGAM, приделав туда корреляцию между разными направлениями, или использовать распределение Дирихле - но нужно знать, какое реально распределение скрывается под записанными румбами.

для начала можно просто вычислить главные компоненты временных рядов, нарезав их соответствующим образом (как минимум, объединив одинаковые сроки в один ряд), и сравнить проекции на компоненты, отличные от первой - первая компонента и есть общий тренд.

[nickleb]

.. с учетом сроков и сезонности...

ох уж эта сезонность... как её грамотно вычленить... иногда вихрь синоптического масштаба напорочь съедает эту сезонную... в атмосфере м.б. попроще - сеть мощная-равномерная наблюдений, ряды длинные, репрезентативные... а в океане?: мониторинг съёмок недёшев? дрейфующие буи - с ними да (ARGO-эра) революция.. но они дрейфуют (буйковые заякоренные - дорогое удовольствие)... с буями типа ARGO - да - вал данных - ктр. мало кто понимает как обрабатывать: "которая" в них изменчивость по времени, а "которая" по месту??? почитаешь Вас, gamm, и хочется опять дрель с утра в руки взять... спасибо, что свёрла правильные помогаете выбирать... голову-то не поменять, а свёрла-то можно...

[gamm]

ох уж эта сезонность... как её грамотно вычленить...

все дело в волшебных пузырьках (с)

главная проблема - это правильно задать распределение, т.е. структуру модели. А потом Байесовские методы способны творить чудеса. Все эти ваши вихри нужно просто включить в модель.

что касается сезонности, то можно пытаться вычленить не календарную, а феноменологическую сезонность, опираясь на какие-то ключевые события. Или совмещать ряды, считая корреляцию (опять на основе модели). Иногда получается.

[nickleb]

что касается сезонности, то можно пытаться вычленить не календарную, а феноменологическую сезонность, опираясь на какие-то ключевые события

иногда такого рода сезон - это сезон удобного времени наблюдений - светлое время года и ещё лёд крепкий - не стал подтаивать и "гулять" - март-май - так для высоких широт - когда на "самолётках-вертолётках" летали и съёмки делали - вот это за "зиму" и считали... а глухой-то настоящей ночью-зимой (декабрь-февраль) - что там делалось?

[gamm]

иногда такого рода сезон - это сезон удобного времени наблюдений - светлое время года и ещё лёд крепкий - не стал подтаивать и "гулять" - март-май - так для высоких широт - когда на "самолётках-вертолётках" летали и съёмки делали - вот это за "зиму" и считали... а глухой-то настоящей ночью-зимой (декабрь-февраль) - что там делалось?

а снимки Модиса на что нам дадены добрым американским правительством, тепловой канал и ночью работает - наверняка, что-то можно вытянуть А если доброе канадское правительство Радарсат даст - так и вообще щастье наступит

[nickleb]

иногда такого рода сезон - это сезон удобного времени наблюдений - светлое время года и ещё лёд крепкий - не стал подтаивать и "гулять" - март-май - так для высоких широт - когда на "самолётках-вертолётках" летали и съёмки делали - вот это за "зиму" и считали... а глухой-то настоящей ночью-зимой (декабрь-февраль) - что там делалось?

а снимки Модиса на что нам дадены добрым американским правительством, тепловой канал и ночью работает - наверняка, что-то можно вытянуть А если доброе канадское правительство Радарсат даст - так и вообще щастье наступит

"знам" про то... да сейчас ещё и вмороженные в лёд ITP-буи дрейфуют - наши "кораблики" - тоже их ставили... правда самописцы у них бегают не до нижней границы тёплых атлантических вод в Канадском секторе... но это очень здорово... надо будет- анимашку их дрейфа вышлю... а так они для всего сообщества доступны и расположены на стр. http://www.whoi.edu/website/itp/overview