рекомендации

пятница, 8 мая 2020 г.

Таблицы сопряженности в R

Оригинал: Contingency Tables in R

Как создать таблицу

Во-первых, давайте получим некоторые данные. Пакет MASS содержит данные о 93 автомобилях, поступивших в продажу в США в 1993 году. Они хранятся в объекте Cars93 и включают 27 признаков для каждого автомобиля, некоторые из которых являются категориальными. Итак, давайте загрузим пакет MASS и посмотрим на тип транспортных средств, включенных в cars93:


library(MASS)
Cars93$Type

##  [1] Small   Midsize Compact Midsize Midsize Midsize Large   Large  
##  [9] Midsize Large   Midsize Compact Compact Sporty  Midsize Van    
## [17] Van     Large   Sporty  Large   Compact Large   Small   Small  
## [25] Compact Van     Midsize Sporty  Small   Large   Small   Small  
## [33] Compact Sporty  Sporty  Van     Midsize Large   Small   Sporty
## [41] Sporty  Small   Compact Small   Small   Sporty  Midsize Midsize
## [49] Midsize Midsize Midsize Large   Small   Small   Compact Van    
## [57] Sporty  Compact Midsize Sporty  Midsize Small   Midsize Small  
## [65] Compact Van     Midsize Compact Midsize Van     Large   Sporty
## [73] Small   Compact Sporty  Midsize Large   Compact Small   Small  
## [81] Small   Compact Small   Small   Sporty  Midsize Van     Small  
## [89] Van     Compact Sporty  Compact Midsize
## Levels: Compact Large Midsize Small Sporty Van

У нас там 6 типов автомобилей. Функция table сообщает, сколько у нас их каждого типа:

table(Cars93$Type)
##
## Compact   Large Midsize   Small  Sporty     Van
##      16      11      22      21      14       9

prop.table преобразует их в доли:

prop.table(table(Cars93$Type))

prop.table(table(Cars93$Type))

##
##    Compact      Large    Midsize      Small     Sporty        Van
## 0.17204301 0.11827957 0.23655914 0.22580645 0.15053763 0.09677419

То же самое с происхождением автомобилей:

table(Cars93$Origin)

##
##     USA non-USA
##      48      45

prop.table(table(Cars93$Origin))

##
##      USA  non-USA
## 0.516129 0.483871

Как создать таблицу сопряженности

Отлично, мы увидели, что наш набор данных содержит одинаковое количество автомобилей из США и не из США, и что наиболее распространенными типами являются Midsize и Small. Однако, может быть, США и неамериканцы различаются по типам?

Давайте посмотрим на типы автомобилей с точки зрения их происхождения. Мы можем снова использовать table, но теперь с двумя аргументами. Первый станет переменной строки, а второй станет переменной столбца:

table(Cars93$Type, Cars93$Origin)

##          
##           USA non-USA
##   Compact   7       9
##   Large    11       0
##   Midsize  10      12
##   Small     7      14
##   Sporty    8       6
##   Van       5       4

Теперь мы увидели то, что все знают и так: американцы любят большие машины! В таблице выше показано совместное распределение двух категориальных переменных (Type и Origin). Такие таблицы называются таблицами сопряженности.

Как получить маржинальную форму таблицы сопряженности

Функции rowSums и colSums действительно говорят сами за себя

(tab1<-table(Cars93$Type, Cars93$Origin))

##          
##           USA non-USA
##   Compact   7       9
##   Large    11       0
##   Midsize  10      12
##   Small     7      14
##   Sporty    8       6
##   Van       5       4

rowSums(tab1)

## Compact   Large Midsize   Small  Sporty     Van
##      16      11      22      21      14       9

colSums(tab1)

##     USA non-USA
##      48      45

Как получить проценты из таблицы сопряженности

prop.table, вложенная в table, дает частоты:

prop.table(table(Cars93$Type, Cars93$Origin))

##          
##                  USA    non-USA
##   Compact 0.07526882 0.09677419
##   Large   0.11827957 0.00000000
##   Midsize 0.10752688 0.12903226
##   Small   0.07526882 0.15053763
##   Sporty  0.08602151 0.06451613
##   Van     0.05376344 0.04301075

Преобразование в проценты - это просто умножение на 100:

prop.table(table(Cars93$Type, Cars93$Origin))*100

##          
##                 USA   non-USA
##   Compact  7.526882  9.677419
##   Large   11.827957  0.000000
##   Midsize 10.752688 12.903226
##   Small    7.526882 15.053763
##   Sporty   8.602151  6.451613
##   Van      5.376344  4.301075

Обратите внимание, что это совместное распределение вероятностей, из которого мы видим, например, что около 7,5% автомобилей являются небольшими и имеют американское происхождение.

Чаще всего нас интересует распределение одной переменной внутри групп, созданных другой. Здесь интересным представляется распределение типов автомобилей среди американских и (отдельно) неамериканских автомобилей. Чтобы получить это, мы используем аргумент margin для функции prop.table. Он сообщает, где в строках (margin = 1) или в столбцах (margin = 2) переменная группировки:

prop.table(table(Cars93$Type, Cars93$Origin), margin=2)*100

##          
##                 USA   non-USA
##   Compact 14.583333 20.000000
##   Large   22.916667  0.000000
##   Midsize 20.833333 26.666667
##   Small   14.583333 31.111111
##   Sporty  16.666667 13.333333
##   Van     10.416667  8.888889

Теперь мы можем легко увидеть, что маленькие автомобили встречаются за пределами США в два раза чаще, чем в США.

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

(tab2<-prop.table(table(Cars93$Type, Cars93$Origin), margin=2)*100)

##          
##                 USA   non-USA
##   Compact 14.583333 20.000000
##   Large   22.916667  0.000000
##   Midsize 20.833333 26.666667
##   Small   14.583333 31.111111
##   Sporty  16.666667 13.333333
##   Van     10.416667  8.888889

colSums(tab2)

##     USA non-USA
##     100     100

Критерий хи-квадрат

Наиболее распространенный вопрос, который возникает в таблицах сопряженности, заключается в том, являются ли переменные строк и столбцов независимыми. Самый простой способ ответить на него - запустить тест хи-квадрат. Он подробно описан в этом руководстве. Давайте проверим, являются ли Type и Origin независимыми:

chisq.test(Cars93$Type, Cars93$Origin)

## Warning in chisq.test(Cars93$Type, Cars93$Origin): Chi-squared
## approximation may be incorrect

##
##  Pearson's Chi-squared test
##
## data:  Cars93$Type and Cars93$Origin
## X-squared = 14.08, df = 5, p-value = 0.01511

По-видимому, это не так, но мы также получили предупреждение хи-квадрат. Это связано с тем, что статистика хи-квадрат следует распределению хи-квадрат только приблизительно. Чем больше наблюдений мы имеем, тем лучше приближение. Функция chisq.test выдает вышеупомянутое предупреждение всякий раз, когда одно из ожидаемых значений меньше 5 (значение «ожидаемое количество» см. в руководстве, указанном выше).

Точный тест Фишера

Точный критерий Фишера является альтернативой критерию хи-квадрат, используемой в основном, когда аппроксимация хи-квадрат не является удовлетворительной. Давайте запустим его:

fisher.test(Cars93$Type, Cars93$Origin)

##
##  Fisher's Exact Test for Count Data
##
## data:  Cars93$Type and Cars93$Origin
## p-value = 0.007248
## alternative hypothesis: two.sided

Результаты очень похожи на результаты из хи-квадрат, но это не всегда так. Одним из основных недостатков точного теста Фишера является то, что для больших таблиц (или больших выборок) он становится неэффективным в вычислительном отношении.

G-тест

Другой альтернативой является так называемый G-тест. Его статистика также приблизительно распределена по хи-квадрат, но для небольших выборок это приближение лучше, чем при использовании теста хи-квадрат. Для G-теста мы можем использовать функцию GTest из пакета DescTools. Результаты снова очень похожи на два предыдущих теста: Type и Origin не являются независимыми.

library(DescTools)
GTest(Cars93$Type, Cars93$Origin)

##
##  Log likelihood ratio (G-test) test of independence without
##  correction
##
## data:  Cars93$Type and Cars93$Origin
## G = 18.362, X-squared df = 5, p-value = 0.002526

Поправка Йейтса

В таблицах сопряженности 2x2 критерий хи-квадрат может быть улучшен с помощью поправки непрерывности Йейтса. Она просто вычитает 0,5 из каждого члена | Наблюдаемый - Ожидаемый | в статистике хи-квадрат. Снова обратитесь к этому руководству, если вы не понимаете, о чем идет речь. Более того, R применяет поправку Йейтса автоматически, когда это необходимо. Давайте посмотрим на наличие версий с механической коробкой передач для автомобилей США и неамериканских:

(tab3<-table(Cars93$Man.trans.avail, Cars93$Origin))

##      
##       USA non-USA
##   No   26       6
##   Yes  22      39

chisq.test(tab3)

##
##  Pearson's Chi-squared test with Yates' continuity correction
##
## data:  tab3
## X-squared = 15.397, df = 1, p-value = 8.712e-05
Поправка Йейтса была применена автоматически, и R сказал нам об этом.

Таблица размерности 3 (или больше) 

Давайте посмотрим на Man.trans.avail, Origin и Type одновременно. table разбивает набор данных в соответствии с переменной, представленной как третья:

table(Cars93$Man.trans.avail, Cars93$Origin, Cars93$Type)

## , ,  = Compact
##
##      
##       USA non-USA
##   No    2       0
##   Yes   5       9
##
## , ,  = Large
##
##      
##       USA non-USA
##   No   11       0
##   Yes   0       0
##
## , ,  = Midsize
##
##      
##       USA non-USA
##   No    9       4
##   Yes   1       8
##
## , ,  = Small
##
##      
##       USA non-USA
##   No    0       0
##   Yes   7      14
##
## , ,  = Sporty
##
##      
##       USA non-USA
##   No    0       0
##   Yes   8       6
##
## , ,  = Van
##
##      
##       USA non-USA
##   No    4       2
##   Yes   1       2

ftable дает более компактный вид:

ftable(Cars93$Man.trans.avail, Cars93$Origin, Cars93$Type)

##              Compact Large Midsize Small Sporty Van
##                                                    
## No  USA            2    11       9     0      0   4
##     non-USA        0     0       4     0      0   2
## Yes USA            5     0       1     7      8   1
##     non-USA        9     0       8    14      6   2

Тест Кокрана-Мантеля-Хензеля

Из приведенных выше результатов должно быть ясно, что зависимость между Origin и Man.trans.avail отличается от Type. Например, для небольших и спортивных автомобилей ассоциация отсутствует: каждый американский, а также любой неамериканский маленький или спортивный автомобиль выпускается с ручной коробкой передач. С другой стороны, большинство автомобилей среднего размера в США не имеют такой версии, в то время как большинство неамериканских автомобилей этого типа имеют. Чтобы учесть возможные различия, можно использовать критерий Кохрана-Мантеля-Хензеля:

mantelhaen.test(Cars93$Man.trans.avail, Cars93$Origin, Cars93$Type)

##
##  Mantel-Haenszel chi-squared test with continuity correction
##
## data:  Cars93$Man.trans.avail and Cars93$Origin and Cars93$Type
## Mantel-Haenszel X-squared = 8.0153, df = 1, p-value = 0.004638
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##   2.226531 76.891307
## sample estimates:
## common odds ratio
##          13.08438

Третий аргумент, переданный mantelhaen.test, определяет страты. Сравните приведенные выше результаты с результатами без стратификации (пример поправки Йейтса выше). Ассоциация все еще присутствует, но доказательства для этого гораздо слабее.

Меры ассоциации

Как только мы обнаружили некоторые связи между переменными, пришло время измерить их силу. Существует множество различных типов возможных мер. Многие из них описаны здесь. Сейчас я остановлюсь на двух наиболее широко используемых.

V-коэффициент Крамера

V основан на статистике хи-квадрат:

где:

N - общий итог таблицы сопряженности (сумма всех ее ячеек),
C - количество столбцов
R - количество строк.
V [0; 1]. Чем больше V, тем сильнее связь между переменными. V = 0 можно интерпретировать как независимость (поскольку V = 0 тогда и только тогда, когда χ2 = 0). Основным недостатком V является отсутствие точной интерпретации. Является ли V = 0,6 сильной, умеренной или слабой ассоциацией?

Функция CramerV от DescTools может рассчитать это для нас:

CramerV(Cars93$Type, Cars93$Origin)

## [1] 0.3890967

Еще раз: это сильная, умеренная или слабая ассоциация?

Лямбда Гудмена — Краскэла

Лямбда Гудмена — Краскэла - пример меры, основанной на пропорциональном уменьшении вариаций. Эти виды мер направлены на имитацию R2 - коэффициента детерминации линейной регрессии:

они принимают значения из [0, 1]
они выражают долю вариаций, объясняемых независимой переменной.

Давайте использовать столбец-переменную как независимую. Тогда формула лямбда Гудмена — Краскэла:

где:

Lj - сумма немодальных частот в j-м столбце,
L - сумма немодальных частот в столбце «total»

Это лучше всего объяснить на примере. Давайте посмотрим на Type и Origin, последняя будет независимой:

table(Cars93$Type, Cars93$Origin)

##          
##           USA non-USA
##   Compact   7       9
##   Large    11       0
##   Midsize  10      12
##   Small     7      14
##   Sporty    8       6
##   Van       5       4

Мода столбца US - 11, поэтому L1 = 7 + 10 + 7 + 8 + 5 = 37. Мода столбца non-US - 14, поэтому L2 = 9 + 0 + 12 + 6 + 4 = 31.

Столбец «Total»

rowSums(table(Cars93$Type, Cars93$Origin))

## Compact   Large Midsize   Small  Sporty     Van
##      16      11      22      21      14       9
Мода равна 22, поэтому L = 16 + 11 + 21 + 14 + 9 = 71 и

λ=(71−(37+31))/71=0.042

Вместо этого мы можем использовать функцию Lambda из пакета DescTools:

Lambda(Cars93$Type, Cars93$Origin, direction='row')

## [1] 0.04225352

Параметр direction определяет, где (в строке или в столбце) находится зависимая переменная.

Origin объясняет только около 4% изменений в Type. Обратите внимание на формулу, что лямбда определяет вариацию как дихотомию между принадлежностью и не принадлежностью к самой большой группе.

Стоит отметить, что лямбда равна нулю, когда модальная категория каждого столбца одинакова. Посмотрите на эту таблицу:

(lambdaTab<-cbind(c(0,100), c(49,51)))

##      [,1] [,2]
## [1,]    0   49
## [2,]  100   51

chisq.test(lambdaTab)

##
##  Pearson's Chi-squared test with Yates' continuity correction
##
## data:  lambdaTab
## X-squared = 62.279, df = 1, p-value = 2.981e-15

Модальная категория одинакова для каждого столбца (это вторая строка). Таким образом, лямбда должна быть нулевой, несмотря на значительную и видимую связь:

Lambda(lambdaTab, direction='row')

## [1] 0

Комментариев нет:

Отправка комментария