Факторный анализ в R

Предположим, у вас есть большой набор утверждений (напр., «человек — это звучит гордо», «все люди — сёстры», «худой мир лучше доброй ссоры» и пр.), своё отношение к которым респонденты оценивали по одинаковому шаблону (напр., «согласен / не знаю / не согласен»). Можно, конечно, в статье дать таблички по каждому пункту, но можно попытаться найти что-то, что объединяет одну часть пунктов в более общую категорию, другую — в ещё одну категорию (безусловно, может оказаться и так, что ваши утверждения ничего не объединяет). Факторный анализ — это один из инструментов, который позволяет найти это общее, если оно там, конечно, есть.

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

Это одновременно сильное и слабое место. Сильное потому, что большой массив данных упрощается и его легче анализировать. А слабое потому, что сильная корреляция, как известно, не указывает на причинность и реальные связи — компьютер покажет вам нечто, но что это значит, насколько находка разумна и правдоподобна, судить только вам. Как написано в одной умной книге «to interpret the factors, which is more like voodoo than science».

Однако перейдём к примеру.

Итак, в 2013 г. Центр социальных экспертиз по заказу ВОО «Гей-Альянс Украины» опрашивал обычных людей (800 чел.) на предмет гомофобии (отчёт). Среди прочего, в опроснике фигурировали и пункты, к гомофобии прямого отношения не имеющие, напр. о доверии к разнообразным политическим и социальным институтам. Вопрос звучал так: «Какой уровень Вашего доверия к следующим социальным институтам? (Дайте один наиболее подходящий ответ по каждой строке)» с вариантами ответов «5. Совсем не доверяю — 4. Скорее не доверяю — 3. Трудно сказать, доверяю или нет — 2. Скорее доверяю — 1. Полностью доверяю». Список институтов, к которым респондент выражал своё отношение, таков:

1. Семье и родственникам
2. Соседям
3. Коллегам
4. Церкви и духовенству
5. Астрологам
6. Средствам массовой информации (телевидение, радио, газеты)
7. Политическим партиям
8. Налоговой инспекции
9. Милиции
10. Прокуратуре
11. Судам
12. Президенту
13. Верховной Раде
14. Правительству
15. Местным органам власти
16. Банкам
17. Страховым компаниям
18. Благотворительным фондам, общественным организациям

Как провести факторный анализ этих данных? (предположим, что таблица с ответами называется dovira)
Присоединяем массив:

>attach(dovira)

Вначале следует убедиться, что в загруженном массиве нет пропусков и ошибок ввода:

>which(is.na(dovira)==T)
integer(0)
>summary(dovira)
p1
Min. :1.000
1st Qu.:2.000
Median :2.000
Mean :2.711
3rd Qu.:4.000
Max. :5.000 ... ... ...


Как видим, всё в порядке (для того, чтобы не загромождать изложение, в выводе оставлен только первый вопрос).
Команда, выполняющая факторный анализ, входит в набор пакетов, устанавливаемых по умолчанию. Она очень проста:

>factanal(dovira,6)
Call:
factanal(x = dovira, factors = 6)

Uniquenesses:

1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18
0.431	0.195	0.379	0.614	0.047	0.672	0.506	0.285	0.174	0.106	0.186	0.215	0.112	0.082	0.464	0.288	0.204	0.533

Loadings:

	Factor1	Factor2	Factor3	Factor4	Factor5	Factor6
1	-0.407	-0.324	0.489	-0.106	-0.213
2			0.879	0.131		-0.112
3			0.784
4		-0.128	0.540	-0.170		0.193
5	0.125	0.171		0.133	0.943
6	0.265	0.122	0.252		0.393	0.139
7	0.522	0.382		0.148	0.151	0.175
8	0.395	0.673	-0.119	0.204	0.182	0.131
9	0.329	0.817		0.181
10	0.297	0.865	-0.113	0.145	0.122
11	0.353	0.769	-0.104	0.277
12	0.805	0.320				0.111
13	0.853	0.318	-0.144	0.151	0.121
14	0.902	0.250			0.125
15	0.582	0.230		0.181		0.325
16	0.196	0.414		0.667	0.139	0.184
17	0.243	0.351		0.694	0.160	0.317
18	0.162	0.109		0.228		0.608

	Factor1	Factor2	Factor3	Factor4	Factor5	Factor6
SSloadings	3.662	3.399	2.079	0.324	1.275	0.765
ProportionVar	0.203	0.189	0.116	0.074	0.071	0.043
CumulativeVar	0.203	0.392	0.508	0.581	0.652	0.695

Test of the hypothesis that 6 factors are sufficient.
The chi square statistic is 257.27 on 60 degrees of freedom.
The p-value is 2.95e-26

Посмотрим на результаты.

Вначале в выводе повторяется отданная машине команда, потом идёт табличка «уникальностей», т. е. долей общей дисперсии, вносимых каждой переменной по отдельности. Следом мы видим таблицу нагрузок, в которой столбцы соответствуют коэффициентам корреляции отдельных переменных с выделенными факторами. Наконец, третья таблица — доля общей дисперсии, объясняемой каждым конкретным фактором и накопление этих дисперсий. Завершает вывод информация о тестировании гипотезы «выбранное число факторов достаточно для описания массива».

Наиболее важными являются таблицы нагрузок и долей объясняемой дисперсии.

Из последней видно, что в сумме 6 выделенных факторов объясняют 70% разброса данных, при этом первый фактор отвечает за пятую часть суммарной дисперсии, второй — 19%, третий — 12% и т. д.
Таблица нагрузок указывает, что в первом факторе объединены 7, 12, 13, 14 и 15 институция (коэффициенты корреляций больше 0.5), во втором — 8, 9, 10, 11, в третьем — 2, 3, 4 и т. д.

Попробуем интерпретировать результаты.

Фактор 1 объединяет доверие к политическим партиям, президенту, Верховной Раде, правительству и к местным органам власти. Иными словами, это доверие к политической сфере в целом.
Фактор 2 объединяет доверие к налоговой инспекции, милиции, прокуратуре и судам. Иными словами, это доверие к фискальным и силовым органам.
Фактор 3 объединяет доверие к соседям, коллегам и, неожиданно, к церкви и духовенству. Эти институции можно обобщить следующим образом — доверие к людям, с которыми респонденты встречаются лицом к лицу. В пользу этого говорит и корреляция с уровнем доверия к родственникам (она лишь ненамного ниже, чем произвольно избранный нами порог коэффициента корреляции 0.5).
Фактор 4 — это доверие к банкам и страховым компаниям, т. е. к финансовым учреждениям.
Фактор 5 стоит особняком — доверие к астрологам (других заметных корреляций нет).
Фактор 6 подобно предыдущему коррелирует только с уровнем доверия только к одной институции — благотворительные фонды и общественные организации.
Лишь одна институция не вошла в эти факторы — средства массовой информации (телевидение, радио, газеты). Доверие к ней приблизительно одинаково «размазано» по выделенным факторам.

Что нам дают эти результаты?

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

Доверие к ...	Средний по фактору балл (шкала от 5 — не доверяю до 1 — доверяю)
людям, с которыми респонденты встречаются лицом к лицу	2.5
благотворительным фондами и общественным организациям	3.1
политической сфере в целом	3.4
астрологам	3.4
финансовым учреждениям	3.6
фискальным и силовым органам	3.6

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

Последний аспект, который не может не вызвать вопросов: откуда мы знаем, что факторов нужно выделить именно 6. Самым, пожалуй, точным ответом будет — ниоткуда. Каждый раз, нужно экспериментировать, опираясь на здравый смысл. Во-первых, количество факторов не может быть большим, чем число переменных. Во-вторых, можно ориентироваться на суммарную объясняемую дисперсию, ибо нет смысла рассуждать о факторах, если они в совокупности не описывают хотя бы её половину (а умные люди рекомендуют добиваться по крайней мере 70%). В-третьих, нужно ориентироваться на возможность подобрать разумное объяснение полученным факторам.

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

Литература

Teetor P. R Cookbook. — O’Reilly, 2011


Читать далее

Даты

Обработка дат в R кажется не вполне тривиальным занятием, но при ближайшем знакомстве выглядит несложной. Итак, есть массив в несколько тысяч записей (назовём его massive), одна переменная (пусть она называется date) которого — даты соответствующих наблюдений в привычном формате день.месяц.год (напр., 12.12.2012). Стоит задача рассортировать все наблюдения по месяцам.

Сначала преобразуем нашу переменную в формат даты:

massive$date.new = as.Date(massive$date,'%d.%m.%Y')

Обратите внимание на параметр формата %d.%m.%Y — мы указываем машине, что наши исходные даты представлены в виде день.месяц (цифрами, а не словами).год (полностью), при этом разделителем выступает точка. Если бы год был в сокращённом виде (двумя цифрами, без столетия), то вместо прописной %Y следовало бы поставить строчную %y.

Полученная переменная date.new пригодна для календарных вычислений. Возможно ли напрямую из неё извлечь месяцы, я пока не знаю, поэтому мы воспользовались кружным путём:

massive$date.sort=strtrim(format(massive$date.new,'%Y-%m-%d'),7)

Иными словами, сначала дату мы преобразуем в строки типа 2012-12-12 командой format(), после чего из каждой такой строки отсекаем первые семь символов (т. е. 2012-12) командой strtrim() и записываем полученное в новую переменную, которая и будет сортировать наши наблюдения.

Читать далее

Снова о логистической регрессии

NB: Этот материал представляет собой сокращённый перевод публикации R Data Analysis Examples: Logit Regression, http://www.ats.ucla.edu/stat/r/dae/logit.htm. Большая часть описанных манипуляций давно автоматизирована в пакете epicalc (см. напр. http://donbas-socproject.blogspot.com/2011/03/blog-post_25.html), однако настоящее описание позволяет взглянуть «под капот» и разобраться в логике построения вывода.

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

Для анализа мы дополнительно воспользуемся пакетами aod и ggplot2, подключаемых командами:

library(aod)
library(ggplot2)

Пример. Исследователь ищет, как переменные GRE (Graduate Record Exam scores — оценки во время обучения в вузе), GPA (grade point average — средний балл) и престиж вуза влияют на поступление в аспирантуру. Т. о. искомая переменная, поступил или не поступил (admit/don’t admit), является бинарной.

В описанном ниже анализе мы воспользовались сгенерированными гипотетическими данными, которые можно получить, используя R:

mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
head(mydata)## посмотрим на первые строки данных

## admit gre gpa rank
## 1 0 380 3.61 3
## 2 1 660 3.67 3
## 3 1 800 4.00 1
## 4 1 640 3.19 4
## 5 0 520 2.93 4
## 6 1 760 3.00 2

Этот масив содержит бинарную зависимую переменную admit, а также три предиктора: gre, gpa and rank. Мы воспользуемся переменными gre и gpa как непрерывными. Переменная rank принимает значения от 1 до 4. У вузов с рангом 1 самый высокий престиж, а у вузов с рангом 4 — самый низкий.

Методы анализа, которые вы также можете использовать

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

OLS-регрессия, используемая с бинарной зависимой переменной, известна как линейно-вероятностная модель (linear probability model) и может применяться для описания условных вероятностей. Тем не менее, ошибки (т. н. остатки — residuals), даваемые этим методом, не соответствуют предположению о гомоскедастичности (равной дисперсии) и нормальном распределении, заложенном в OLS-регрессии, что приводит к неправильным значениям стандартных ошибок и ложным результатам тестирования гипотез (подробнее см. Long, 1997, p. 38-40).

Дискриминантный анализ двух групп. Многомерный метод для бинарной зависимой переменной.

Hotelling’s T2. Зависимая переменная 0/1 преобразуется в группирующую переменную, а предикторы — в зависимые переменные. Это позволяет выполнить общий тест значимости, но не даёт индивидуальных коэффициентов для каждой переменной, соответственно остаётся неясным вклад каждого такого «предиктора» в эффект прочих «предикторов».

Использование логистической модели

Код, приведённый ниже, выполняет логистическое моделирование с использованием функции glm() function, но сначала мы превратим переменную rank в фактор, чтобы показать, что она должна трактоваться как категориальная переменная:

mydata$rank <- factor(mydata$rank)
mylogit <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial")

Поскольку мы дали нашей модели имя mylogit, R не выдал результаты регресии на экран. Это, однако, легко сделать командой summary():

summary(mylogit)

## Call:
## glm(formula = admit ~ gre + gpa + rank, family = "binomial",
## data = mydata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.627 -0.866 -0.639 1.149 2.079
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.98998 1.13995 -3.50 0.00047 ***
## gre 0.00226 0.00109 2.07 0.03847 *
## gpa 0.80404 0.33182 2.42 0.01539 *
## rank2 -0.67544 0.31649 -2.13 0.03283 *
## rank3 -1.34020 0.34531 -3.88 0.00010 ***
## rank4 -1.55146 0.41783 -3.71 0.00020 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 499.98 on 399 degrees of freedom
## Residual deviance: 458.52 on 394 degrees of freedom
## AIC: 470.5
##
## Number of Fisher Scoring iterations: 4

В приведённых результатах вначале мы видим саму модель (секция «Call») и заданные опции.

Затем видим отклонение остатков (deviance residuals), измеряющий подгонку модели. Ниже мы обсудим, как пользоваться этой строкой, чтобы оценить соответствие модели.

Следующая часть результатов показывает коэффициенты, их стандартные ошибки, z-статистики (называемые иногда z-статистиками Вальда — коэффициенты, делённые на стандартные ошибки) и соответствующие p-значения. Очевидно, что переменные gre и gpa суть статистически значимы, точно также как и три уровня переменной rank. Коэффициенты логистической регресии указывают, насколько изменится вероятность искомого события при увеличении соответствующего предиктора на одну единицу.

Так, возрастание на одну единицу gre увеличивает логарифм шансов поступить в аспирантуру (по сравнению с непоступлением) на 0.002. Возрастание на одну единицу gpa, увеличивает логарифм шансов поступления на 0.804.

Уровни переменной rank несколько отличаются интерпретацией. Напр., у выпускников вуза с рангом 2 (по сравнению с вузом первого ранга), меньше логарифм шансов быть в аспирантуре на 0.675.

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

Мы можем воспользоваться функцией confint(), чтобы получить доверительные интервалы рассчитанных коэффициентов.

confint(mylogit)

## 2.5 % 97.5 %
## (Intercept) -6.2716202 -1.792547
## gre 0.0001376 0.004436
## gpa 0.1602959 1.464143
## rank2 -1.3008888 -0.056746
## rank3 -2.0276713 -0.670372
## rank4 -2.4000265 -0.753543

Мы можем протестировать общее влияние переменной rank функцией wald.test() из библиотеки aod. Порядок, в котором коєффициенты размещены в таблице коэффициентов, повторяет порядок слагаемых в модели. Это важно, поскольку функция wald.test() применяется к коэффициентам согласно этому порядку (опция Terms указывает R, какие именно слагаемые из модели должны быть протестированы, в данном случае это 4, 5 и 6 суть три слагаемых, относящихся к уровням переменной rank).

wald.test(b = coef(mylogit), Sigma = vcov(mylogit), Terms = 4:6)

## Chi-squared test:
## X2 = 20.9, df = 3, P(> X2) = 0.00011

Значение теста хи-квадрат 20.9 с тремя степенями свободы ассоциировано с p-значением 0.00011, показывающим, что общее влияние переменной rank является статистически значимым.

Мы также можем протестировать другие гипотезы о различиях в коэффициентах разных уровней переменной rank. Напр., ниже мы тестируем предположение, что коэффициент для rank=2 равен коэффициенту для rank=3. Первая строка приведённого кода создаёт вектор l, описывающий тот тест, который мы намерены выполнить. В нашем случае мы тестируем разницу слагаемых для rank=2 и rank=3 (т. е., 4-е и 5-е слагаемое в модели). Для контраста мы умножаем один из них на 1, а другой на −1. Прочие слагаемые модели не участвуют в тесте, поэтому мы умножаем их на 0. Вторая линия кода использует равенство L=l для того, чтобы сказать R, что мы решили выполнить тест с использованием вектора l.

l <- cbind(0, 0, 0, 1, -1, 0)
wald.test(b = coef(mylogit), Sigma = vcov(mylogit), L = l)

## Chi-squared test:
## X2 = 5.5, df = 1, P(> X2) = 0.019

Значение теста хи-квадрат 5.5 с одной степенью свободы ассоциировано с p-значением 0.019, показывающим, что разница между коэффициентами для rank=2 и rank=3 статистически значима.

Вы также можете экспоненцировать коэффициенты и интерпретировать их как отношения шансов (OR). R сделает это, если вы укажете ему, что откуда взять коэффициенты — из coef(mylogit). Аналогично вы можете получить и доверительные интервалы отношений шансов. Чтобы вставить это всё в одну таблицу, мы используем команду cbind():

exp(coef(mylogit))

## (Intercept) gre gpa rank2 rank3 rank4
## 0.0185 1.0023 2.2345 0.5089 0.2618 0.2119

exp(cbind(OR = coef(mylogit), confint(mylogit)))

## OR 2.5 % 97.5 %
## (Intercept) 0.0185 0.001889 0.1665
## gre 1.0023 1.000138 1.0044
## gpa 2.2345 1.173858 4.3238
## rank2 0.5089 0.272290 0.9448
## rank3 0.2618 0.131642 0.5115
## rank4 0.2119 0.090716 0.4707

Теперь мы можем сказать, что при возрастании на одну единицу переменной gpa отношение шансов поступить в аспирантуру (в сравнении с непоступлением) возрастает в 2.23 раза.

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

Мы начнём с расчёта презсказанных вероятностей поступления для каждого значения престижности вуза (переменная rank), удерживая переменные gre и gpa на уровне их средних. Создадим и посмотрим на таблицу данных:

newdata1 <- with(mydata, data.frame(gre = mean(gre), gpa = mean(gpa), rank = factor(1:4)))
newdata1

## gre gpa rank
## 1 587.7 3.39 1
## 2 587.7 3.39 2
## 3 587.7 3.39 3
## 4 587.7 3.39 4

Заметьте, что в новой таблице имена переменных должны быть такими же, как и в вашей регрессионной модели (напр., среднее для gre должно быть названо gre). Теперь, имея таблицу, мы можем поручить R создать пресказанные вероятности. Первая линия кода очень компактна — значения rankP должны быть предсказаны [функция predict()] на основе результатов анализа, содержащихся в mylogit и значений предикторов из таблицы newdata1.

newdata1$rankP <- predict(mylogit, newdata = newdata1, type = "response")
newdata1

## gre gpa rank rankP
## 1 587.7 3.39 1 0.5166
## 2 587.7 3.39 2 0.3523
## 3 587.7 3.39 3 0.2186
## 4 587.7 3.39 4 0.1847

Из приведённых результатов следует, что предсказанные вероятности быть принятым в аспирантуру суть 0.52 для студентов из наиболее престижных вузов (rank=1), но 0.18 для студентов из наименее престижных (rank=4) при одинаковых средних показателях успеваемости gre и gpa. Аналогично можно варьировать значения переменных gre и rank. Если мы намерены нарисовать эти зависимости, мы должны создать 100 значений gre в интервале 200 и 800 для каждого значения переменной rank (т. е., 1, 2, 3 и 4):

newdata2 <- with(mydata, data.frame(gre = rep(seq(from = 200, to = 800, length.out = 100), 4), gpa = mean(gpa), rank = factor(rep(1:4, each = 100))))

Код, генерирующий предсказанные вероятности, аналогичен приведённому выше.

newdata2 <- cbind(newdata2, predict(mylogit, newdata = newdata2, type = "link", se = TRUE))
newdata2 <- within(newdata2, { PredictedProb <- plogis(fit) LL <- plogis(fit - (1.96 * se.fit)) UL <- plogis(fit + (1.96 * se.fit)) })
head(newdata2)

## gre gpa rank fit se.fit residual.scale UL LL
## 1 200.0 3.39 1 -0.8115 0.5148 1 0.5492 0.1394
## 2 206.1 3.39 1 -0.7978 0.5091 1 0.5499 0.1424
## 3 212.1 3.39 1 -0.7840 0.5034 1 0.5505 0.1454
## 4 218.2 3.39 1 -0.7703 0.4978 1 0.5512 0.1485
## 5 224.2 3.39 1 -0.7566 0.4922 1 0.5519 0.1517
## 6 230.3 3.39 1 -0.7429 0.4866 1 0.5525 0.1549

## PredictedProb
## 1 0.3076
## 2 0.3105
## 3 0.3134
## 4 0.3164
## 5 0.3194
## 6 0.3224

Рисунок полученной модели делается так:

ggplot(newdata2, aes(x = gre, y = PredictedProb)) + geom_ribbon(aes(ymin = LL,
ymax = UL, fill = rank), alpha = 0.2) + geom_line(aes(colour = rank), size = 1)



Мы можем также захотеть узнать, насколько хорошо наша модель описывает реальные данные. Это может быть особенно полезно, если мы сравниваем конкурирующие модели. Результаты, выводимые на экран командой summary(mylogit), включают показатели подгонки (они выводятся сразу под коэффициентами), в частности разброс остатков нулевой и данной моделей, а также AIC. Один из параметров подгонки — значимость модели в целом. Этот тест показывает, действительно ли модель с предикторами описывает реальные данные значимо лучше, чем модель без предикторов (т. н. нулевая модель). Тестовая статистика — это разница между отклонениями остатков в моделях с предикторами и без них. Эта статистика описывается распределением хи-квадрат с числом степеней свободы равном разнице в степенях свободы тестируемой и нулевой моделью. Чтобы вычислить эту величину, мы можем воспользоваться командой:

with(mylogit, null.deviance - deviance)

## [1] 41.46

А разница степеней свободы может быть получена так:

with(mylogit, df.null - df.residual)

## [1] 5

Наконец, p-значение вычисляется так:

with(mylogit, pchisq(null.deviance - deviance, df.null - df.residual, lower.tail = FALSE))

## [1] 7.578e-08

Хи-квадрат 41.46 с 5 степенями свободы и ассоциированным p-значением, меньшим 0.001, показывают нам, что наша модель, в целом, описывает данные значительно лучше, чем нулевая модель. Этот тест иногда называют «likelihood ratio test».

Примечания

Hosmer, D. & Lemeshow, S. (2000). Applied Logistic Regression (Second Edition). New York: John Wiley & Sons, Inc.
Long, J. Scott (1997). Regression Models for Categorical and Limited Dependent Variables. Thousand Oaks, CA: Sage Publications.

Читать далее

Небольшой алгоритм для импорта данных из SPSS в R

При работе со сторонними данными, набранными в SPSS, у меня возникала проблема: спсс-овские метки переменных не экспортировались в R, что затрудняло работу с получившимся массивом (переменные в этом случае назывались примерно так v1, v2 etc). Решение нашлось в блоге http://strengejacke.wordpress.com/2013/02/22/migrating-from-spss-to-r-rstats/:

my.data=read.spss('data.sav',to.data.frame=FALSE, use.value.labels=T)

my.data.table=as.data.frame(my.data)

my.data.lab=attr(my.data,'variable.labels')

names(my.data.table)=my.data.lab

Чтобы нарисовать картинку с нужным заглавием, можно воспользоваться созданной вспомогательной переменной my.data.lab (напр., для переменной в 11 столбце):

plot(my.data.table[,11], main=my.data.lab[11])


Читать далее

Как показать среднее и погрешности на гистограмме

Представим, что на иллюстрации нужно показать как результаты нескольких экспериментов, так и их среднее с погрешностями. Один из способов в R таков:

data=c(1,2,3,4,3,2,1) # вводим данные
barplot(data) # рисуем нашу иллюстрацию
abline(h=mean(data)) # добавляем линию средней
abline(h=mean(data)+sd(data)) # добавляем линии верхнего и нижнего стандартного отклонений
abline(h=mean(data)-sd(data))

Вуаля:

Вместо функций mean() и sd(), естественно, можно вставить либо конкретные числа, либо другие меры средней тенденции и погрешностей. Например, подключив пакет library(epicalc), можно визуализировать среднее с доверительным интервалом:

barplot(data)
abline(h=mean(data))
abline(h=c(ci(data)$lower95ci,ci(data)$upper95ci))


Читать далее

Упражнения­ по экспоненци­альному моделирова­нию случайных графов на примере флирта в сериале Grey’s Anatomy

Перевод поста http://badhessian.org/lessons-on-exponential-random-graph-modeling-from-greys-anatomy-hook-ups/

Я недавно нашёл прекрасный­ пост Gary Weissman’а Grey’s Anatomy Network of Sexual Relations (перевод на русский тут) и вдохновилс­я им. Для тех, кто ничего не слышал об этом телесериал­е, сообщаю, что он чрезвычайн­о популярен,­ отмечен наградами как лучшая медицинска­я драма, транслируе­мая ABC. Идя навстречу ожиданиям зрителей, шоу регулярно показывает­, как его герои флиртуют друг с другом.

Пытаясь дать студентам-медикам ряд простейших­ концепций анализа социальных­ сетей, Weissman создал набор сетевых данных о сексуальны­х контактах героев этого сериала. Хотя я особо не интересуюс­ь ни этим шоу, а сексуальны­е или воображаем­ые сети лежат вдали от моих исследовательских интересов,­ пост Weissman’а мне кажется отличной демонстрац­ией анализа социальных­ сетей в педагогиче­ских целях.

Я бы хотел, ради развлечения, вернуться к этим данным и рассмотреть их с помощью экспоненциального моделирования случайных графов (ERGM) в аспекте образования связей. Прежде чем я продолжу, я бы хотел сказать, что самые лучшие пособия по ERGM созданы их авторами и доступны на сайте. Моим любимым пособием из-за его прикладной направленности является серия презентаций на февральской конференции Sunbelt Workshop 2011 г. Если вы всерьёз настроены изучить ERGM, то указанные публикации должны стать вашими настольными книгами. Тем не менее, я думаю, что существует та аудитория, которая с бОльшим восторгом наблюдает за приключениями сексуальных врачей на телевидении, чем читает о брачных сетях во Флоренции 16 века.

Экспоненциальные модели случайных графов являются способом понять процесс возникновения сетевой структуры, процесс образования в ней сетей. В нашем примере «образование сетей» и «снять потрахаться» являются синонимами. Эти модели работают, измеряя некий набор известных сетевых характеристик и используя распределение их значений для генерации случайных сетей. Такие случайные сети затем сравниваются с сетью наблюдаемой, чтобы оценить правдоподобие модели. Хорошие ERG-модели продуцируют сети, изоморфно похожие на наблюдаемую сеть, и пользуются небольшим набором характеристик. Плохие ERG-модели продуцируют сети, которые слабо похожи на наблюдаемую сеть, но при этом используется много ненужных характеристик.

Какими характеристикам следует пользоваться? Число рёбер (т. е. «связей», «отношений») является одной из самых важных. Если мы занимаемся образованием связи, то очень важно знать число этих связей. Кроме того, в посте Weissman’а некоторые комментаторы упоминали важность пола персонажей или, к примеру, их должность в больнице. Узлы (т. е., «акторы», «вершины») могут иногда характеризоваться атрибутами, такими как пол, возраст, должность или астрологический знак. При моделировании сети можно использовать атрибуты узлов, чтобы оценить такие ассортативные эффекты как общая гомофилия («Спят ли члены больничных коллективов с сотрудниками, занимающими такие же должности?»), дифференциальная гомофилия («Спят ли штатные врачи с другими такими же? Делают ли они это чаще, чем интерны с другими интернами?») и ассортативное смешивание («Do resident-attending couples occur more frequently than we’d expect vis-à-vis nurse-attending couples?»). Можно анализировать и другие свойства сетей, но это зависит от доступности данных, типа анализируемой сети (направленная, ненаправленная, бимодальная), а также от сходимости модели.

Использованные здесь данные восходят к тем, которые Weissman сделал доступными. Я добавил узлы и рёбра по информации, найденной при просмотре Wikipedia. Я включал лишь те персонажи, которые переспали с одним или несколькими другими персонажами в сериале. Я также добавил такие атрибуты, как пол персонажа, расу, должность в больнице, приблизительный возраст, номер сезона, в котором он появился, и астрологический знак, под которым родился актёр. Я собрал эти атрибуты из Wikipedia и IMDB. Для простоты я кодировал расу каждого персонажа как чёрный (black), белый (white) или другой (other). Хотя я не люблю рассматривать расу именно в этих терминах, однако число персонажей и их разнообразие невелико, чтобы вдаваться в нюансы этого вопроса. Что касается приблизительного возраста персонажей, то я пользовался возрастом актёров, их играющих, а в двух случаях мне пришлось самому сделать предположения. Как у Weissman’а, нынешние данные, вероятно, не содержат какой-то информации об отношениях, узлах, упрощают описание персонажей и их взаимоотношений.

Для начала откройте R-терминал. Установите и загрузите пакеты ergm и RCurl, если вы не сделали этого раньше.

install.packages("RCurl"); install.packages("ergm")

library(RCurl); library(ergm)

Пакет RCurl будет использован для того, чтобы загрузить данные из моего файла на Google Docs. Пакет ergm включает несколько других пакетов, таких как network, поскольку ERGM использует объекты, созданные последним.

Далее, прочтите данные и создайте объект класса network.

#First, read in the sociomatrix

ga.mat<-getURL("https://docs.google.com/spreadsheet/pub?key=0Ai--oOZQWBHSdDE3Ynp2cThMamg1b0VhbEs0al9zV0E&single=true&gid=0&output=txt",ssl.verifypeer = FALSE)

ga.mat<-as.matrix(read.table(textConnection(ga.mat), sep="\t",header=T, row.names=1, quote="\""))

#Second, read in the network attributes

ga.atts<-getURL("https://docs.google.com/spreadsheet/pub?key=0Ai--oOZQWBHSdDE3Ynp2cThMamg1b0VhbEs0al9zV0E&single=true&gid=1&output=txt",ssl.verifypeer=FALSE)

ga.atts<-read.table(textConnection(ga.atts), sep="\t", header=T, quote="\"",stringsAsFactors=F, strip.white=T, as.is=T)

#Third, create a network object using the sociomatrix and its corresponding attributes

ga.net<-network(ga.mat,vertex.attr=ga.atts,vertex.attrnames=colnames(ga.atts),directed=F,hyper=F,loops=F,multiple=F,bipartite=F)

Теперь, когда вы получили данные, мы можем поиграть с ними. Прежде чем перейти к моделированию, давайте посмотрим на схему сети, отметив розовым цветом узлы, означающие женщин, и голубым — мужчин.

plot(ga.net,vertex.col=c("blue","pink")[1+(get.vertex.attribute(ga.net,"sex")=="F")],label=get.vertex.attribute(ga.net,"name"), label.cex=.75)

Расположение узлов на вашем рисунке может несколько отличаться. Здесь то, что получилось у меня:

Следует отметить несколько интересных особенностей этой сети.

Во-1, с точки зрения связности у нас есть большой связный компонент,но в целом граф поделён на четыре отдельных компонента. Если бы шоу продолжалось без появления новых персонажей, я бы думал, что один из самых маленьких компонентов должен бы «соединиться» с самым большим компонентом.

Во-2, большая часть сексуальных отношений в сериале суть гетеросексуальны. Только два ребра (torres-arizona и torres-hahn) были однополыми. По всей видимости, здесь нет нечётно-численных циклов: в сексуальных сетях циклы из трёх, пяти, семи и др. элементов содержат как минимум одну однополую связь внутри цикла.

В-3, в отражение способа построения выборки сеть не содержит сексуально-одиноких изолятов. Минимальная степень (число связей, падающих на узел) равна единице. В-четвёртых, сеть оказывается уязвимой перед табу «бывший бывшей бывшего». В одном из исследований сексуальных сетей учёные обнаружили запрет на создание пары с бывшим партнёром бывшего партнёра бывшего партнёра из-за возможных осложнений. Кажется, что здесь есть как минимум три цикла, отвечающих критериям этого табу. Вместе с тем, также хорошо известно, что нежелательные последствия, возникающие от сексуальных отношений с другими, делают шоу правдоподобным.

Давайте начнём наш формальный анализ с базовой модели, которая включает только число рёбер и, как второй предиктор, nodematch("sex"), который схватывает гетеросексуальную тенденцию сети.

ga.base<-ergm(ga.net~edges+nodematch("sex")) #Estimate the model summary(ga.base) #Summarize the model

Вот выдача команды summary(ga.base):

==========­==========­======
Summary of model fit
==========­==========­======

Formula: ga.net ~ edges + nodematch(“sex”)

Iterations: 20

Monte Carlo MLE Results:

Estimate Std. Error MCMC % p-value

edges -2.3003 0.1581 NA <1e -04="-04" -3.1399="-3.1399" 0.001="0.001" 0.01="0.01" 0.05="0.05" 0.1="0.1" 0.7260="0.7260" 0="0" 1311.43="1311.43" 1="1" 2="2" 320.47="320.47" 324.47="324.47" 334.17="334.17" 944="944" 946="946" 990.97="990.97" aic:="aic:" bic:="bic:" code="code" codes:="codes:" degrees="degrees" deviance:="deviance:" e-04="e-04" freedom="freedom" na="na" nodematch.sex="nodematch.sex" null="null" of="of" on="on" residual="residual" signif.="signif."

В согласии с ожиданиями число рёбер и гетеросексуальность суть значимые предикторы образования связей, поскольку p-значение существенно меньше, чем условный уровень значимости 0.05. Колонка «MCMC%» слева указана как пропущенная (NA), поскольку эта модель не требует оценки марковских цепей по методу Монте-Карло. Коэффициенты, −2.30 и −3.14, показывают условные логарифмы шансов (log-odds), таким образом логарифмы шансов того, что два персонажа вступили в сексуальные отношения, равны (-2.30)*увеличение числа связей на единицу + (-3.14)*увеличение числа однополых отношений. Логарифм шансов того, что одна гетеросексуальная связь образуется внутри сети, равен −2.30, что даёт её вероятность exp(-2.30)/(1+exp(-2.30)) = 0.09. Аналогично, логарифм шанса образования гомосексуальной связиравен −2.30-3.14 = 5.44, что соответствует вероятности 0.004 = exp(-5.44)/(1+exp(-5.44)).

Давайте посмотрим, какую случайную, расчётную сеть эта модель могла бы создать.

plot(simulate(ga.base),vertex.col=c("blue","pink")[1+(get.vertex.attribute(ga.net,"sex")=="F")])

Ваша картинка может несколько отличаться от того, что получилось у меня:

Эта сеть точно отражает как число рёбер, так и склонность к образованию гетеросексуальных связей, однако здесь не всё в порядке: некоторые из наших узлов суть изоляты. Иными словами, расчёт предсказывает существование в сети сексуально-воздерживающихся индивидов. Поскольку никаких таких людей в нашей реальной сети нет из-за особенностей формирования выборки, построенная модель неудовлетворительна.

Хотя у нас есть значения критерия согласия (goodness of fit), данного функцией правдоподобия, один обычный, систематический способ определения критерия согласия — это сравнение параметров немоделированной, наблюдаемой сети с такими же параметрами рассчитанных ERGM сетей.

ga.base.gof<-gof(ga.base)

summary(ga.base.gof) #Summarize the goodness of fit

par(mfrow=c(3,1)); plot(ga.base.gof) #Plot three windows. It's OK to ignore the warning.

Рисунок воспризводит ту же информацию, которую выдаёт команда summary(). p-Значения в правых колонках суть вероятности того, что рассчитанные сети точно совпадают с сетью наблюдаемой в ряду данных параметров. Мы видим, что в терминах of the edgewise shared partner and minimum geodesic distance measurements represent the data quite well. Однако, если мы посмотрим на такой параметр как центральность, то увидим то, что согласуется с нашими выводами из рисунков — число воздерживающихся изолятов (узлов с центральностью 0) больше, чем нужно, число узлов с центральностью 1 (моногамистов) — меньше, и что нет оценок для узлов с центральностью 3 и 9.

Мы сейчас добавим в модель параметр, контролирующий эффект моногамии (центральность 1).

ga.base.d1<-ergm(ga.net~edges+nodematch("sex")+degree(1)) summary(ga.base.d1) plot(simulate(ga.base.d1),vertex.col=c("blue","pink")[1+(get.vertex.attribute(ga.net,"sex")=="F")]) ga.base.d1.gof<-gof(ga.base.d1); summary(ga.base.d1.gof) par(mfrow=c(3,1)); plot(ga.base.d1.gof)

Если принять значение BIC предыдущей модели 334.17 как исходное, то мы увидим, что модель ga.base.d1 лучше — её BIC равен 316.56. Эффект моногамии также является статистически значимым. Рисунок демонстрирует меньше изолятов, что подтверждает команда gof(). Судя по значениям критерия согласия, модель кажется разумным приближением к реальности, если не принимать во внимание персонаж шоу с девятью сексуальными партнёрами.

В отличие от предыдущей модели, эта использует методику оценки марковских цепей (MCMC). Мы можем оценить результаты MCMC.

mcmc.diagnostics(ga.base.d1)

В выводе следует сначала глянуть на второй пункт, в котором сравниваются параметры расчётной и наблюдаемой моделей.

Are sample statistics significantly different from observed?

edges nodematch.sex degree1 Overall (Chi^2)

diff. -0.0968000 0.0232000 0.0534000 NA
test stat. -0.2942081 0.8328180 0.3379779 1.086016
P-val. 0.7685989 0.4049474 0.7353798 0.780451

Вроде бы всё в порядке: рассчитанные при помощи MCMC графы значимо не отличаются от наблюдаемого. Давайте посмотрим по-другому:

Sample statistics auto-correlation:
Chain 1

edges nodematch.sex degree1

Lag 0 1.0000000 1.00000000 1.0000000
Lag 100 0.8476292 0.47956824 0.7334256
Lag 200 0.7332856 0.26641687 0.6116231
Lag 300 0.6371281 0.17109432 0.5198955
Lag 400 0.5587029 0.12489565 0.4480439
Lag 500 0.4883954 0.08496445 0.3860807

Sample statistics burn-in diagnostic (Geweke):
Chain 1

Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5

edges nodematch.sex degree1

-1.590 -1.068 1.372

P-values (lower = worse):

edges nodematch.sex degree1

0.1117691 0.2856841 0.1701817

В целом, нет ничего подозрительного в результатах диагностики MCMC. Параметр автокорелляции между лагами несколько больше, чем предпочитаемый мною (лаги должны быть ближе к нулю и дальше от единицы, за исключением нулевого лага). Кроме того, статистика Geweke, приблизительно соответствующая z-статистике, когда её p-значения незначимы, должна быть значительно ближе к нулю. К счастью, мы можем так изменить некоторые параметры по умолчанию в MCMC для того, чтобы получить более приемлемый результат.

ga.base.d1<-ergm(ga.net~edges+nodematch("sex")+degree(1),control=control.ergm(MCMC.burnin=50000, MCMC.interval=5000))

При увеличении MCMC.burnin с 10,000 (по умолчанию) до 50,000 должен уменьшаться параметр Geweke. При увеличении MCMC.interval со 100 (по умолчанию) до 5000 должна уменьшаться автокорреляция между лагами. Если бы параметры нашего расчётного образца действительно значимо отличались от наблюдаемых, мы бы увеличили размер выборки в MCMC с 10,000 (по умолчанию) до 20,000 или даже 50,000.

ga.base.d1<-ergm(ga.net~edges+nodematch("sex")+degree(1),control=control.ergm(MCMC.burnin=50000,MCMC.interval=5000, MCMC.samplesize=50000))

Нам следовало бы поупражняться в увеличении этих параметров MCMC, но это может резко увеличить время расчётов.

summary(ga.base.d1) #The figures are mostly the same

mcmc.diagnostics(ga.base.d1)

Ура! MCMC диагностика заметно улучшилась, что указывает на большую жёсткость нашей последней модели.

Теперь, когда мы понимаем чуть больше, как оценивать качество модели, мы можем вернуться к интригующему вопросу секса! Только что мы моделировали сеть по числу рёбер, склонности к гетеросексуальности и по моногамии, но могут ли другие аспекты привнести смысл в анализ? Я бы предложил начать с возраста, поскольку большинство любовных связей возникают между сверстниками.

ga.base.d1.age<-ergm(ga.net~edges+nodematch("sex")+degree(1)+absdiff("birthyear"),control=control.ergm(MCMC.burnin=50000, MCMC.interval=5000))

summary(ga.base.d1.age)

mcmc.diagnostics(ga.base.d1.age)

summary(gof(ga.base.d1.age))

Кажется, модель ощутимо улучшилась! Значение BIC снизилось с 316.56 до 297.51, MCMC-диагностика выглядит хорошо, значение критерия согласия выглядит, в целом, разумным, а разница в возрасте значима и отрицательна, что указывает на то, что чем больше разница в возрасе акторов, тем меньше шансы у них флиртовать друг с другом.

Наконец, давайте рассмотрим ассортативность с точки зрения группировки по расам. Предыдущие исследования юношеских сексуальных и любовных сетей открыли такую ассортативность. Создательница сериала, Shonda Rhimes, сознательно построила отбор актёров так, чтобы представить расовое разнообразие. Поскольку шоу включает ряд афроамериканских (наряду с обилием других не-белых) персонажей, можно попытаться установить, привело ли такое разнообразие при отборе персонажей к расовому разнообразию в сексуальных отношениях. Для начала, давайте попробуем построить модель общей гомофилии по расе.

ga.base.d1.age.race<-ergm(ga.net~edges+nodematch("sex")+degree(1)+absdiff("birthyear")+nodematch("race"),control=control.ergm(MCMC.burnin=50000,MCMC.interval=5000))

summary(ga.base.d1.age.race)

mcmc.diagnostics(ga.base.d1.age.race)

summary(gof(ga.base.d1.age.race))

Мы в самом деле нашли эффект расового отбора, предполагающий, что образование связи более вероятно между персонажами одной и той же расы. В этом шоу сексуальные отношения являются скорее внутри-, нежели межрасовыми по сравнению с тем, что было бы, если бы персонажи разных рас комбинировались совершенно случайно. MCMC-диагностика выглядит хорошо, а измерения критерия согласия в пределах разумного. Хотя полученная модель выглядит хорошо по многим индикаторам, но заметно небольшое увеличение BIC (298.21) по сравнению с предыдущей моделью (297.51). Если бы нам были интересны только параметры сети, то предыдущая модель даёт лучшие значения.

Но возможно, что белые персонажи проявляют большую гомофилию, а чёрные — меньшую? Или наоборот? Мы можем посмотреть на это как на дифференциальную гомофилию. В модель добавляем параметр diff=T к слагаемому nodematch("races").

ga.base.d1.age.racediff<-ergm(ga.net ~ edges + nodematch("sex") + degree(1) +absdiff("birthyear") + nodematch("race", diff=T), control=control.ergm(MCMC.burnin=50000,MCMC.interval=5000))

summary(ga.base.d1.age.racediff)

Хотя мы можем моделировать дифференциальную гомофилию, ergm() не смог интерпретировать категорию «other», поскольку в наблюдаемой сети нет сексуальных отношений между латиноамериканцами и азиатоамериканцами. ergm() указывает на эту сложность следующей фразой:

Observed statistic(s) nodematch.race.Other are at their smallest attainable values. Their coefficients will be fixed at -Inf.

Аналогично, the AIC, BIC, отклонение и эффекты nodematch.race.Other помечены как NA в summary. Давайте попробуем создать новую модель, исключив гомофилию между латиноамериканцами и азиатоамериканцами. Чтобы это сделать, мы включаем параметр keep=c(1,3) в слагаемое nodematch("race"), поскольку «Black» и «White» были первой и третьей категориями, перечисленными в summary(ga.base.d1.age.racediff).

ga.base.d1.age.racediff<-ergm(ga.net~edges+nodematch("sex")+degree(1)+absdiff("birthyear")+nodematch("race", diff=T, keep=c(1,3)),control=control.ergm(MCMC.burnin=50000, MCMC.interval=5000))

summary(ga.base.d1.age.racediff)

mcmc.diagnostics(ga.base.d1.age.racediff)

summary(gof(ga.base.d1.age.racediff))

Как и в предыдущей модели, мы видим значимую ассортативность по расовым группам в сериале. Как чёрные, так и белые персонажи проявляют склонность к расовой гомофилии бОльшую, чем можно было бы ждать, если бы выбор партнёра был случайным. Эта модель также указывает, что эффекты гомофилии сильнее среди афроамериканских персонажей (на это указывает более высокий коэффициент и меньшее p-значение при nodematch.race.Black по сравнению с nodematch.race.White term). Проверка MCMC-диагностики и статистик критерия согласия говорит, что полученная модель является хорошим отражением оригинальной сети. Вместе с тем, BIC ещё больше сдвигается к 300.73, указывая что модель менее экономна по сравнению с моделью общей гомофилии и модели с возрастом.

Давайте посмотрим на картинку.

plot(simulate(ga.base.d1.age.racediff),vertex.col=c("blue","pink")[1+(get.vertex.attribute(ga.net,"sex")=="F")])

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

Хотя число съёмов в сериале Grey’s Anatomy может казаться удивительным, в аспекте образования диад видно скорее традиционное восприятие сексуальных отношений: доминируют нормы моногамии, а бОльшая часть связей суть гетеросексуальны и соединяют пары близкого возраста и одной расы. Следует, однако отметить, что описанный анализ обошёл ряд альтернативных версий, которые могли бы лучше отразить свойства сети. Я приглашаю читателей самостоятельно поиграть с данными и с другими моделями. Несколько дополнительных примеров:

#Perhaps men and women have different tendencies to form sexual partnerships?

ga.base.d1.age.sex<-ergm(ga.net~edges+nodematch("sex")+degree(1) +absdiff("birthyear")+nodefactor("sex"), control=control.ergm(MCMC.burnin=100000, MCMC.interval=5000)) #Maybe less "traditional" than previously concluded?

#Perhaps you're interested in the assortative mixture among roles in the hospital?
#Resident-resident sexual contacts (28) are the reference group,
#and unobserved pairings between positions (3-5, 9, 12-15, 17-21, 24, 27) are omitted.

ga.base.d1.age.rolemix<-ergm(ga.net~edges+nodematch("sex")+degree(1)+absdiff("birthyear") +nodemix("position", base=c(3:5, 9, 12:15, 17:21, 24, 27, 28)), control=control.ergm(MCMC.burnin=50000, MCMC.interval=5000))

list.vertex.attributes(ga.net) #List the different vertex attributes in the network object

get.vertex.attribute(ga.net, "sign") #Provides each actor's astrological sign.

?ergm.terms #This command will provide a list of other terms that ergm() can model

?control.ergm #This command will present different options on the estimation methods.


Читать далее

Быстро построить линейную модель

Быстро построить картинку линейной зависимости двух непрерывных переменных можно так:

plot(Y~X)
abline(lm(Y~X))

Ну а параметры этой модели выводятся на экран как обычно: summary(lm(Y~X))
Читать далее

Диаграмма Парето и недоступность тестирования

Популярное правило «20% усилий определяет 80% результата» было установлено эмпирически экономистом и социологом Вильфредом ПАРЕТО в 1897 г. Оно часто применяется при оценке качества деятельности, её эффективности и поиске путей её оптимизации. С помощью R можно легко выполнить анализ Парето и визуализировать результат (справедливости ради отмечу, что не менее легко это делается в любом табличном редакторе, поскольку сам по себе анализ прост до неприличия).

Анализ Парето заключается в классификации проблем на немногочисленные, но существенно важные, и многочисленные, но несущественные. Он позволяет распределить усилия и установить основные факторы, с которых нужно начинать действовать, чтобы преодолеть возникающие проблемы.

Рассмотрим анализ на примере данных Мониторинга сексуального поведения МСМ и их знаний в отношении ВИЧ за 2009 г. (краткие результаты тут, полный отчёт тут).

В опроснике было два вопроса, касающихся доступности тестирования на ВИЧ: «Является ли тестирование на ВИЧ для Вас доступным?» (да/нет) и «Почему для Вас лично тестирование недоступно?» (вопрос задавался тем, кто на предыдущий ответил «нет», в качестве вариантов респондентам предложено 10 формулировок, в том числе «другое», при этом каждый опрошенный мог выбрать сразу несколько вариантов).

В результате получено такое распределение ответов:

1) «не знаю, к кому обратиться» — 44%
2) «в нашем населённом пункте нет учреждения, где тестируют» — 2%
3) «не знаю, где находится пункт тестирования» — 24%
4) «нет денег на тестирование» — 22%
5) «неудобный график работы учреждения или пункта тестирования» — 3%
6) «неудобное расположение пункта тестирования» — 1%
7) «не устраивает отношение персонала» — 1%
8) «боюсь разглашения своего статуса» — 18%
9) «другое» — 3%
10) «затрудняюсь ответить» — 21%

При анализе доли полученных ответов сортируются по мере убывания, а кумулятивная кривая показывает тот момент, когда сумма самых больших долей достигает 80%.

Проделаем это в R.

1. Введём наши данные в вектор:
> nedostup = c(44,2,24,22,3,1,1,18,3,21)

2. Для удобства припишем данным названия — номера ответов:
> names(nedostup)=c(1:10)

3. Установим в R нужный пакет и подключим его:
> install.packages('qcc')
> library(qcc)

4. Нарисуем, собственно, диаграмму и получим табличку:
> pareto.chart(nedostup)

Pareto chart analysis for nedostup
Frequency Cum.Freq. Percentage Cum.Percent.
1 44.0000000 44.0000000 31.6546763 31.6546763
3 24.0000000 68.0000000 17.2661871 48.9208633
4 22.0000000 90.0000000 15.8273381 64.7482014
10 21.0000000 111.0000000 15.1079137 79.8561151
8 18.0000000 129.0000000 12.9496403 92.8057554
5 3.0000000 132.0000000 2.1582734 94.9640288
9 3.0000000 135.0000000 2.1582734 97.1223022
2 2.0000000 137.0000000 1.4388489 98.5611511
6 1.0000000 138.0000000 0.7194245 99.2805755
7 1.0000000 139.0000000 0.7194245 100.0000000

5. Собственно, всё — техническая часть закончена. Теперь смотрим на содержание. Из колонки «Cum.Percent.» в табличке видно, что самыми важными (80% всех случаев недоступности) являются такие причины как «не знаю, к кому обратиться», «не знаю, где находится пункт тестирования», «нет денег на тестирование» и «затрудняюсь ответить».

Поскольку вариант «нет денег на тестирование» означает, по сути, то же незнание (тестирование на ВИЧ в Украине бесплатно), постольку напрашивается вывод, что основным препятствием при тестировании для опрошенных МСМ является недостаток информации о том, где и как пройти тест (отметим также, что вариант «затрудняюсь ответить» в целом согласуется с предыдущими — у таких респондентов нет вообще представления о том, что у них спрашивают, т. е. нет информации о тестировании).

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

Примечание:

В заметке использован материал http://sixsigmaonline.ru/load/22-1-0-255
Читать далее

Нелинейная регрессия

В своей работе химика-кинетика я часто встречался с задачей описать набор данных заданной нелинейной функцией. В R это делается достаточно просто.

1. Загружаем наши данные (предположим, что они содержатся в файле 1.csv):
> test=read.csv('~/1.csv',header=T)

2. На всякий случай проверяем, что получилось:
> test
X Y
1 0.5 0.00
2 1.0 0.01
3 2.0 0.02
4 3.0 0.05
5 4.0 0.13
6 5.0 0.18
7 6.0 0.24
8 7.0 0.29
9 8.0 0.35
10 9.0 0.37
11 10.0 0.41
12 11.0 0.44
13 12.0 0.47
14 15.0 0.57
15 20.0 0.71
16 25.0 0.80
17 30.0 0.86
18 35.0 0.90
19 43.0 0.95
20 50.0 0.97
21 60.0 0.99
22 70.0 0.99
23 80.0 1.00

Переменная X — время в минутах, Y — объём поглощённого кислорода в мл.

3. Попробуем построить график:
> attach(test)
> plot(Y~X)

4. Данные должны описываться кинетическим уравнением первого порядка Y = a*(1-exp(-k*X)), пробуем:
> test.mod=nls(Y~a*(1-exp(-k*X)),data=test,start=list(a=1,k=0.05))

Здесь nls() — команда для выполнения нелинейной регрессии, атрибутами которой являются:
а) модель, в нашем случае это Y~a*(1-exp(-k*X))
б) экспериментальные данные, у нас это test
в) стартовые значения параметров модели, которые подбираются исследователем эмпирически.

Если всё сделано правильно и параметры подобраны удачно, то в test.mod будет содержаться информация о полученной модели:
> test.mod
Nonlinear regression model
model: Y ~ a * (1 - exp(-k * X))
data: test
a k
1.05594 0.04967
residual sum-of-squares: 0.03621

Number of iterations to convergence: 4
Achieved convergence tolerance: 4.59e-06

5. Попробуем визуализировать нашу модель:
> test.predict=predict(test.mod,newdata=data.frame(X=0:100))

Здесь мы на промежутке от 0 минут до 100 вычисляем Y

Рисуем:
> plot(Y~X)
> lines(0:100,test.predict)

Ура! Всё получилось!

Читать далее

03. Просто R — Использование R для ознакомительной статистики

Продолжение серии моих переводов из учебного пособия «Просто R — Использование R для ознакомительной статистики» американского автора John’a VERZANI.
http://www.math.csi.cuny.edu/Statistics/R/simpleR

Раздел 15: Анализ различий (дисперсионный анализ, ANOVA)

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

Разберём вначале пример однофакторного анализа.

Классификация стипендий

Предположим, что школа должна рассмотреть 300 заявок на стипендии. Эта работа слишком велика для одного человека, поэтому школа нанимает шестерых классификаторов. Комитет по рассмотрению стипендий хотел бы убедиться, что все классификаторы пользуются одной и той же шкалой при оценивании заявок, поскольку в ином случае заявители не будут находиться в равных условиях. Одним из способов сделать это является распределить в случайном порядке между классификаторами задания (каждому из шестерых достанется по 50) и затем сравнить оценки, выставленные каждым из шести, зная, что различия между классификаторами будут минимальны (т. е. будут в пределах случайной ошибки) в том случае, если все пользуются одинаковой шкалой оценок.

Для того, чтобы проиллюстрировать работу метода, положим, что в нашем распоряжении есть 27 тестов и три классификатора (а не 300 и 6 — для простоты). Кроме того, пусть классификаторы пользуются пятибальной шкалой. Выставленные ими оценки выглядят так:

grader 1: 4, 3, 4, 5, 2, 3, 4, 5
grader 2: 4, 4, 5, 5, 4, 5, 4, 4
grader 3: 3, 4, 2, 4, 5, 5, 4, 4

Введём эти данные в R и объединим затем в одну таблицу:
> x = c(4,3,4,5,2,3,4,5)
> y = c(4,4,5,5,4,5,4,4)
> z = c(3,4,2,4,5,5,4,4)
> scores = data.frame(x,y,z)
> boxplot(scores)

Вначале мы построили диаграммы размахов, позволяющие зрительно сравнить три распределения. Из рисунка следует, что второй классификатор отличается от двух других.

Дисперсионный анализ позволяет нам установить, действительно ли средние оценки всех классификаторов одинаковы в пределах ошибки. Функция в R, позволяющая проводить однофакторный дисперсионный анализ (oneway.test) требует данных в ином формате, а именно — данные должны составлять один вектор, а во второй должен быть фактором, описывающим классификаторов или категорию (иными словами в одном столбце должны быть все оценки классификаторов, а во втором — указание, какая оценка какому классификатору принадлежит). Команда stack() сделает это преобразование:

> scores = stack(scores)
> names(scores)
[1] "values" "ind"

Посмотрев на имена векторов (names), мы узнали, что значения собраны в переменной «values», а категории — в переменной «ind». Чтобы вызвать oneway.test(), мы должны воспользоваться такой же записью, какая применяется в формулах:

> oneway.test(values ~ ind, data=scores, var.equal=T)
One-way analysis of means
data: values and ind
F = 1.1308, num df = 2, denom df = 21, p-value = 0.3417

Мы видим, что p = 0.34, а это значит, что мы принимаем нулевую гипотезу о равенстве средних. Более детальную информацию об анализе можно получить, воспользовавшись функциями anova() и aov().

Дополнение из книги Paul TEETOR. R Cookbook.

Итак, вы узнали, что есть значимые отличия в средних, относящихся к разным группам. Однако сам по себе тест не показывает, какие именно группы отличаются.

Функция TukeyHSD() позволяет посчитать эти отличия и указывает вам на самые большие из них. В ней реализован т. н. метод «по-честному значимых отличий» («honest significant differences»), предложенный Джоном ТЬЮКИ.

Выполнение дисперсионного анализа функцией aov() даёт результаты, к которым применяют функцию TukeyHSD():
> TukeyHSD(aov(x ~ f))

Здесь x — это ваши данные, а f — группирующий их фактор. Вы можете визуализировать результаты для наглядности:
> plot(TukeyHSD(m))

Тест Краскела–Уоллиса

Тест Краскела–Уоллиса является непараметрическим тестом, используемым вместо однофакторного дисперсионного анализа в том случае, если данные не распределены нормально. Его применяют подобно тому, как тест Вилкоксона вместо t-теста.

Тест Краскела–Уоллиса используется в R подобно команде oneway.test:
> kruskal.test(values ~ ind, data=scores)
Kruskal-Wallis rank sum test
data: values by ind
Kruskal-Wallis chi-squared = 1.9387, df = 2, p-value = 0.3793



Читать далее

Преобразование таблицы в массив для логлинейного анализа

Описанная в предыдущем материале команда loglm() работает со специфической формой исходных данных — массивом (на языке R — array), отличающимся от обычной таблицы. Задача преобразовать обычную таблицу (на языке R — data.frame) в array оказалась не вполне тривиальной.

1) импортируем таблицу в R так, чтобы стоящие в ней NA были значимой категорией:

> lxx=read.csv('путь к файлу',header=T,na.strings=' ')

2) делаем, как указано в http://r-statistics.livejournal.com/5467.html, но с соответствующими смысловыми заменами:

> names <- list(age=c('15-19','20-29','30-39','40-49','50.and.older','NA'), family.status.legal=c('single','married','divorsed','NA'), family.status.actual=c('alone','parent','with.man','with.woman','NA'), clubA=c(0,1,'NA'), partiesA=c(0,1,'NA'), cruisingA=c(0,1,'NA'), internet=c(0,1,'NA'), newspaper=c(0,1,'NA'), teletext=c(0,1,'NA'), on.street=c(0,1,'NA'))

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

> lx=array(c(age=lxx$age,
family.status.legal=lxx$family.status.legal,
family.status.actual=lxx$family.status.actual,
clubA=lxx$clubA,
partiesA=lxx$partiesA,
cruisingA=lxx$cruisingA,
internet=lxx$internet,
newspaper=lxx$newspaper,
teletext=lxx$teletext,
on.street=lxx$on.street),dim=c(6,4,5,3,3,3,3,3,3,3),dimnames=names)
Читать далее

Логлинейный анализ

Перевод материала William B. KING, http://ww2.coastal.edu/kingw/statistics/R-tutorials/loglin.html

Вступительные замечания

Мы будем упражняться на встроенном в R наборе данных «Titanic». Он описывает последствия катастрофы лайнера «Титаник» в 1912 г. и представляет из себя массив с четырьмя столбцами...

> data(Titanic)
> dimnames(Titanic)
$Class
[1] "1st" "2nd" "3rd" "Crew"

$Sex
[1] "Male" "Female"

$Age
[1] "Child" "Adult"

$Survived
[1] "No" "Yes"


Как видим, здесь есть четыре категориальных переменных, а именно «Class» описывает класс каюты путешествовавших пассажиров (значения «1st», «2nd», «3rd», «Crew»), «Sex» описывает пол пассажиров (значения «Male», «Female»), «Age» описывает их возраст по стратам («Child», «Adult») и, наконец, выжил ли пассажир или погиб в этой катастрофе, указано в переменной «Survived» («No», «Yes»). В этой таблице...

> margin.table(Titanic)
[1] 2201

...записей, каждая из которых соответствует человеку. Эти данные были собраны British Board of Trade при расследовании катастрофы.

Логлинейный анализ (Log linear analysis) позволяет взглянуть на взаимоотношения между переменными в многомерной таблице сопряжённости. Подобно своему старшему брату, критерию хи-квадрат для двумерной таблицы, логлинейный анализ основывается на допущении, что все представленные в таблице наблюдения являются независимыми (возможно, это не совсем так в случае наших данных), а ожидаемые частоты в таблицах 2 на 2 будут достаточно высоки, как правило 5 и больше. До тех пор пока ожидаемые частоты составляют 1 или больше (при этом не более 20% из них меньше пяти), точность процедуры будет приемлемой, однако её мощность будет очень мала! (собственно, это такие же допущения, как и у пирсоновского критерия хи-квадрат).

Посмотрим на данные, собранные в переменных «Class» и «Age»...

> margin.table(Titanic, c(2,4)) # эта команда показывает столбцы 2 и 4
Survived
Sex No Yes
Male 1364 367
Female 126 344

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

> (344/126) / (367/1364)
[1] 10.14697

... величину большую 10 (т. е. десять к одному, 10:1). Отношение правдоподобия (relative likelihood или likelihood ratio, LR) соответственно...

> (344/(344+126)) / (367/(367+1364))
[1] 3.452165

... (отношение правдоподобия — это доля выживших женщин, делённая на долю выживших мужчин). Итак, женщины более чем в 3.5 раз чаще выживали. В этом и состоит язык логлинейного анализа. Пирсоновский критерий хи-квадрат для этой таблицы показал бы, без сомнения, высокую значимость связи (в этом контексте мы называем её «взаимодействием» (interaction) между рассмотренными двумя факторами).

Немного теории

Хи-квадрат является основным кирпичиком логлинейного анализа, однако в несколько иной форме, называемой хи-квадрат отношения правдоподобия (likelihood ratio chi square). Её преимуществом является то, что хи-квадраты отношения правдоподобия обладают свойством аддитивности. Это значит, что хи-квадраты, полученные из отдельных действий, складываются вместе и дают значения сложного эффекта (или «модели»). R (насколько мне известно) не содержит теста, позволяющего вычислить хи-квадрат отношения правдоподобия для таблицы 2 на 2, однако нам ничего не мешает пользоваться этим тестом для такой простейшей таблицы, поэтому ниже я привожу скрипт:

### начало текста скрипта
> likelihood.test = function(x) {
nrows = dim(x)[1] # no. of rows in contingency table
ncols = dim(x)[2] # no. of cols in contingency table
chi.out = chisq.test(x,correct=F) # do a Pearson chi square test
table = chi.out[[6]] # get the OFs
ratios = chi.out[[6]]/chi.out[[7]] # calculate OF/EF ratios
sum = 0 # storage for the test statistic
for (i in 1:nrows) {
for (j in 1:ncols) {
sum = sum + table[i,j]*log(ratios[i,j])
}
}
sum = 2 * sum # the likelihood ratio chi square
df = chi.out[[2]] # degrees of freedom
p = 1 - pchisq(sum,df) # p-value
out = c(sum, df, p, chi.out[[1]]) # the output vector
names(out) = c("LR-chisq","df","p-value","Pears-chisq")
round(out,4) # done!
}
### конец скрипта

Выделите, скопируйте и вставьте текст в окно терминала R. При необходимости нажмите клавишу Enter в конце. Это создаст функцию в вашем рабочем пространстве с именем likelihood.test(). Посмотрим, как она применяется:

> sex.survived = margin.table(Titanic, c(2,4)) ### создаём таблицу 2 на 2
likelihood.test(sex.survived)
LR-chisq df p-value Pears-chisq
434.4688 1.0000 0.0000 456.8742

Результат выполнения команды состоит из значения отношения правдоподобия, числа степеней свободы, значения p-критерия, и, для сравнения, значения пирсоновского критерия хи-квадрат.

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

В традиционном критерии хи-квадрат таблицы 2 на 2 более простая модель обычно не имеет смысла, но для того, чтобы показать вам саму идею, допустим, что предлагается следующая модель (нулевая гипотеза): все частоты в клетках таблицы равны. Иными словами, нет не только каких-либо взаимодействий между двумя переменными, но также нет и действия какого-либо фактора. Ничего не мешает нам протестировать нулевую гипотезу с использованием критерия хи-квадрат и мы сделаем его, приняв ожидаемые частоты равными 2201/4=550.25. На языке логлинейного анализа эта модель, в которой все частоты равны, называется нулевой (null model).

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

Логлинейный анализ

Наша цель — объяснить полученные частоты в массиве «Titanic» возможно более простой моделью.

Эффекты, что они означают?

Class в некоторых классах пассажиров больше, чем в других
Sex пассажиров одного пола больше, чем пассажиров другого
Age пассажиров одной возрастной группы больше, чем пассажиров другой
Survived численности выживших и погибших в катастрофе не одинаковы
Class × Sex переменные «Class» и «Sex» не независимы (напр., мужчины зарабатывают больше женщин, поэтому могут позволить себе путешествовать более дорогим классом)
Class × Age переменные «Class» и «Age» не независимы (напр., члены команды корабля [значение «Crew»] ехали без детей)
Class × Survived переменные «Class» и «Survived» не независимы (напр., члены команды корабля способны дольше продержаться на воде после катастрофы в силу профессиональной подготовки, пассажиры верхних палуб, т. е. более дорогого класса, находились в момент катастрофы ближе к шлюпкам и спасательным кругам)
Sex × Age переменные «Sex» и «Age» не независимы (напр., среди путешествовавших взрослых больше мужчин)
Sex × Survived переменные «Sex» и «Survived» не независимы (напр., в силу культурных традиций в шлюпки сажали прежде всего женщин или, наоборот, мужчины, будучи сильнее, могли первыми захватить спасательные средства)
Age × Survived переменные «Age» и «Survived» не независимы (напр., взрослые выживают чаще, так как они сильнее и спасательные средства рассчитаны на них)
Class × Sex × Age эти три переменные связаны друг с другом
Class × Sex × Survived эти три переменные связаны друг с другом
Class × Age × Survived эти три переменные связаны друг с другом
Sex × Age × Survived эти три переменные связаны друг с другом
Class × Sex × Age × Survived все четыре переменные связаны друг с другом

Заметьте, мы не рассматриваем факторы как зависимые или независимые переменные. Мы это могли бы сделать, но не в этот раз. Важно, чтобы вы об этом помнили. Вполне объяснима склонность трактовать «Survived» как зависимую переменную (отклик), на которую влияют те или иные факторы. Но мы сейчас НЕ будем этого делать. Мы попытаемся смоделировать (объяснить, предсказать, оценить) частоты в таблице сопряжённости.

Начнём, пожалуй, с такого:

> summary(Titanic)
Number of cases in table: 2201
Number of factors: 4
Test for independence of all factors:
Chisq = 1637.4, df = 25, p-value = 0
Chi-squared approximation may be incorrect

Вы только что сделали логлинейный анализ модели, которая включает в себя лишь главные эффекты, т. е. эффекты каждого из факторов без каких-либо взаимодействий между ними (первые четыре строчки в подразделе «Эффекты, что они означают»). Таким образом, эту модель (нулевую гипотезу) мы отбрасываем. В некоторых случаях существуют такие взаимодействия, которые мы должны найти, чтобы адекватно объяснить частоты в ячейках.

В R можно найти несколько функций для построения логлинейной модели. Среди них loglin() в библиотеке "stats" (она загружается по умолчанию), loglm() в библиотеке "MASS", а также glm() в библиотеке "stats". Давайте немного поэкспериментируем...

> library("MASS")
> loglm( ~ Sex + Survived, data=sex.survived)
Call:
loglm(formula = ~ Sex + Survived, data = sex.survived)

Statistics:
X^2 df P(> X^2)
Likelihood Ratio 434.4688 1 0
Pearson 456.8742 1 0

Мы видим, что в таблице 2 на 2 «Sex» от «Survived», когда мы предполагаем только чистые эффекты индивидуальных факторов, мы получаем те же самые результаты, какие у нас получились при использовании моей команды likelihood.test(). Иными словами, мы воспользовались логлинейным моделированием, чтобы применить критерий хи-квадрат на независимость переменных в этой таблице сопряжённости. Модель, содержащая лишь чистые эффекты, оказалась неприемлемой (нулевая гипотеза отброшена, p очень близко к 0). Следовательно, тут есть связь между двумя переменными (мы это предполагали).

Проделаем то же самое с полной таблицей...

> loglm(~ Class + Sex + Age + Survived, data=Titanic)
Call:
loglm(formula = ~Class + Sex + Age + Survived, data = Titanic)

Statistics:
X^2 df P(> X^2)
Likelihood Ratio 1243.663 25 0
Pearson 1637.445 25 0

Результат такой же, как и при использовании команды summary(Titanic). Существенно, что мы применили критерий хи-квадрат на независимость четырёх переменных и отбросили нулевую гипотезу. Какие-то из этих факторов взаимодействуют друг с другоми, давая наблюдаемые частоты в ячейках. Насыщенная модель, с другой стороны, не даёт ничего полезного, поскольку она полностью объясняет (предсказывает, моделирует) наблюдаемые частоты:

> loglm(~ Class * Sex * Age * Survived, data=Titanic)
Call:
loglm(formula = ~Class * Sex * Age * Survived, data = Titanic)

Statistics:
X^2 df P(> X^2)
Likelihood Ratio 0 0 1
Pearson NaN 0 1

Нулевая гипотеза не отбрасывается (p очень близко к 1). Истина лежит где-то посредине между этими крайностями.

Попытаемся обойтись без четырёхпеременного взаимодействия:

> loglm(~ Class * Sex * Age * Survived - Class:Sex:Age:Survived, data=Titanic)
Call:
loglm(formula = ~Class * Sex * Age * Survived - Class:Sex:Age:Survived,
data = Titanic)

Statistics:
X^2 df P(> X^2)
Likelihood Ratio 0.0002728865 3 0.9999988
Pearson NaN 3 NaN

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

С другой стороны...

> loglm(~ Class + Sex + Age + Survived + Age:Survived, data=Titanic)
Call:
loglm(formula = ~Class + Sex + Age + Survived + Age:Survived,
data = Titanic)

Statistics:
X^2 df P(> X^2)
Likelihood Ratio 1224.103 24 0
Pearson 1596.846 24 0

... модель со всеми факторами плюс взаимодействием Age × Survived предсказывает частоты, которые значимо отличаются от реальных. Нужно добавить больше вариантов взаимодействий, чтобы получить «истинную правду» (т. е. чтобы создать адекватную статистическую модель того, что произошло на Титанике).

Тестирование конкретных гипотез

Предположим, сначала нам следует проверить гипотезу о том, что пол связан с выживаемостью в катастрофе на Титанике. Из критерия хи-квадрат для двух переменных, который мы сделали выше, следует, что это похоже на правду. Однако, переменная «Sex» связана также с переменными «Class» (p близко к нулю) и «Age» (p близко к нулю), и обе они значимо связаны с «Survived». Следовательно, взаимодействие переменных «Class» и «Age» может испытывать влияние взаимодействия переменных «Sex» и «Survived». Подобно множественной регрессии (multiple regression) для числовых переменных логлинейная регрессия позволит нам рассмотреть отдельный вклад этих эффектов.

Если мы удалим все обозначения взаимодействий, включающие как «Sex», так и «Survived», то модель по-прежнему будет адекватно описывать наблюдаемые частоты, следовательно можно сделать вывод, что пол и выживаемость не являются связанными...

> sat.model = loglm(~ Class * Sex * Age * Survived, data=Titanic)
sat.model
Call:
loglm(formula = ~Class * Sex * Age * Survived, data = Titanic)

Statistics:
X^2 df P(> X^2)
Likelihood Ratio 0 0 1
Pearson NaN 0 1
> ### already knew this...
> model2 = update(sat.model, ~.-(Class:Sex:Age:Survived+Sex:Age:Survived+
+ Class:Sex:Survived+Sex:Survived))
> model2
Call:
loglm(formula = ~Class + Sex + Age + Survived + Class:Sex + Class:Age +
Sex:Age + Class:Survived + Age:Survived + Class:Sex:Age +
Class:Age:Survived, data = Titanic)

Statistics:
X^2 df P(> X^2)
Likelihood Ratio 436.2715 8 0
Pearson NaN 8 NaN

Модель плохо описывает наблюдаемые частоты, поскольку в ней нет взаимодействия между «Sex» и «Survived». Заметьте, что использование команды update() даёт новую модель. Конечно, мы могли бы просто перенабрать формулу модели, но команда update() является более удобным средством. Использованный нами в этой команде синтаксис означает «обновить sat.model с использованием всех тех же предикторов (точка = „те же самые“), за исключением (знак минуса) следующих за знаком минуса». Давайте сделаем её по-другому, убрав только два взаимодействия между тремя переменными...

> model3 = update(sat.model, ~.-(Class:Age:Sex:Survived+Sex:Age:Survived+
+ Class:Sex:Survived))
> model3
Call:
loglm(formula = ~Class + Sex + Age + Survived + Class:Sex + Class:Age +
Sex:Age + Class:Survived + Sex:Survived + Age:Survived +
Class:Sex:Age + Class:Age:Survived, data = Titanic)

Statistics:
X^2 df P(> X^2)
Likelihood Ratio 76.90406 7 5.884182e-14
Pearson NaN 7 NaN

Неа! И эту модель отбрасываем. Очевидно, взаимодействие переменных «Sex» и «Survived» обуславливается (зависит от, изменяется с) ещё одним фактором.

Как вы можете представить, такой путь поиска логлинейной модели является очень утомительным и чреватым ошибками. R предоставляет нам способ сделать это автоматически.

Функция step()

Хотя я, в целом, не одобряю автоматизированное построение статистических моделей, в данном случае я склолнен сделать исключение...

> step(sat.model, direction="backward")
Start: AIC=64
~Class * Sex * Age * Survived

Df AIC
- Class:Sex:Age:Survived 3 58
64

Step: AIC=58
~Class + Sex + Age + Survived + Class:Sex + Class:Age + Sex:Age +
Class:Survived + Sex:Survived + Age:Survived + Class:Sex:Age +
Class:Sex:Survived + Class:Age:Survived + Sex:Age:Survived

Df AIC
- Sex:Age:Survived 1 57.685
58.000
- Class:Sex:Age 3 61.783
- Class:Age:Survived 3 89.263
- Class:Sex:Survived 3 117.013

Step: AIC=57.69
~Class + Sex + Age + Survived + Class:Sex + Class:Age + Sex:Age +
Class:Survived + Sex:Survived + Age:Survived + Class:Sex:Age +
Class:Sex:Survived + Class:Age:Survived

Df AIC
57.685
- Class:Sex:Age 3 71.953
- Class:Age:Survived 3 95.899
- Class:Sex:Survived 3 126.904
Call:
loglm(formula = ~Class + Sex + Age + Survived + Class:Sex + Class:Age +
Sex:Age + Class:Survived + Sex:Survived + Age:Survived +
Class:Sex:Age + Class:Sex:Survived + Class:Age:Survived,
data = Titanic, evaluate = FALSE)

Statistics:
X^2 df P(> X^2)
Likelihood Ratio 1.685479 4 0.7933536
Pearson NaN 4 NaN

Эта процедура применяется для насыщенной модели ("sat.model"), подвергая её тому, что называется «регрессия методом исключения» (backward regression). Функция step() применяется с опцией "direction=backward" (другие опции суть "forward" и "both"). Она упрощает модель, используя критерий Akaike's Information Criterion (AIC). Чем меньше значение AIC, тем лучше. В нашем случае, R приходит к модели, в которой только два взаимодействия удалены, а именно — четырёхпеременное и Sex:Age:Survived. Это значит, что отношение между переменными «Sex» и «Survived» обуславливаются переменной «Class». Давайте это проверим...

> margin.table(Titanic, c(2,4,1))
, , Class = 1st

Survived
Sex No Yes
Male 118 62
Female 4 141

, , Class = 2nd

Survived
Sex No Yes
Male 154 25
Female 13 93

, , Class = 3rd

Survived
Sex No Yes
Male 422 88
Female 106 90

, , Class = Crew

Survived
Sex No Yes
Male 670 192
Female 3 20

Отношения шансов очень легко рассчитать вручную, исходя из данных этих таблиц...

> ### odds ratio 1st class
> (141/4) / (62/118)
[1] 67.08871
> ### odds ratio 2nd class
> (93/13) / (25/154)
[1] 44.06769
> ### odds ratio 3rd class
> (90/106) / (88/422)
[1] 4.071612
> ### odds ratio crew
> (20/3) / (192/670)
[1] 23.26389

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

Получение большего объёма информации из модели

Допустим, что мы хотим работать с моделью, созданной командой step() как с моделью случившегося на Титанике. Вначале, я помещу модель в data object (я попросту скопирую и вставлю команду, продуцирующую эту модель из вывода команды step())...

> loglm(formula = ~Class + Sex + Age + Survived + Class:Sex + Class:Age +
+ Sex:Age + Class:Survived + Sex:Survived + Age:Survived +
+ Class:Sex:Age + Class:Sex:Survived + Class:Age:Survived,
+ data = Titanic) -> step.model

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

> print(step.model) # то же самое, как если бы мы просто напечатали step.model
Call:
loglm(formula = ~Class + Sex + Age + Survived + Class:Sex + Class:Age +
Sex:Age + Class:Survived + Sex:Survived + Age:Survived +
Class:Sex:Age + Class:Sex:Survived + Class:Age:Survived,
data = Titanic)

Statistics:
X^2 df P(> X^2)
Likelihood Ratio 1.685479 4 0.7933536
Pearson NaN 4 NaN
>
> fitted(step.model) # ожидаемые частоты, вычисленные по нашей модели
Re-fitting to get fitted values
, , Age = Child, Survived = No

Sex
Class Male Female
1st 0.00000 0.00000
2nd 0.00000 0.00000
3rd 37.43281 14.56719
Crew 0.00000 0.00000

, , Age = Adult, Survived = No

Sex
Class Male Female
1st 118.0000 4.0000
2nd 154.0000 13.0000
3rd 384.5672 91.4328
Crew 670.0000 3.0000

, , Age = Child, Survived = Yes

Sex
Class Male Female
1st 5.00000 1.00000
2nd 10.98493 13.01507
3rd 10.56718 16.43282
Crew 0.00000 0.00000

, , Age = Adult, Survived = Yes

Sex
Class Male Female
1st 57.00000 140.00000
2nd 14.02291 79.97709
3rd 77.43281 73.56719
Crew 192.00000 20.00000

> resid(step.model) # стандартные остатки (standardized residuals)
Re-fitting to get frequencies and fitted values
, , Age = Child, Survived = No

Sex
Class Male Female
1st 0.000000e+00 0.000000e+00
2nd 0.000000e+00 0.000000e+00
3rd -4.020602e-01 6.208006e-01
Crew 0.000000e+00 0.000000e+00

, , Age = Adult, Survived = No

Sex
Class Male Female
1st 0.000000e+00 0.000000e+00
2nd 0.000000e+00 0.000000e+00
3rd 1.239264e-01 -2.555637e-01
Crew 0.000000e+00 0.000000e+00

, , Age = Child, Survived = Yes

Sex
Class Male Female
1st 2.142148e-08 -4.552313e-08
2nd 4.546268e-03 -4.178434e-03
3rd 7.221252e-01 -6.159455e-01
Crew 0.000000e+00 0.000000e+00

, , Age = Adult, Survived = Yes

Sex
Class Male Female
1st -6.825888e-08 6.182584e-08
2nd -6.118921e-03 2.561368e-03
3rd -2.779358e-01 2.820973e-01
Crew 0.000000e+00 0.000000e+00

Излишне отмечать, что в наших ожидаемых частотах есть те, которые равны нулю. Именно поэтому, существует опасность неточности нашей модели.

Использование glm() для логлинейного моделирования

Мы также можем выполнить логлинейный анализ, воспользовавшись командой glm(), которая используется для общих линейных моделей. Для этого нам понадобится таблица с данными (data frame), а не таблица сопряжённости, однако нет ничего сложного в том, чтобы её получить:

> ti = as.data.frame(Titanic)
> ti
Class Sex Age Survived Freq
1 1st Male Child No 0
2 2nd Male Child No 0
3 3rd Male Child No 35
4 Crew Male Child No 0
5 1st Female Child No 0
6 2nd Female Child No 0
7 3rd Female Child No 17
8 Crew Female Child No 0
9 1st Male Adult No 118
10 2nd Male Adult No 154
11 3rd Male Adult No 387
12 Crew Male Adult No 670
13 1st Female Adult No 4
14 2nd Female Adult No 13
15 3rd Female Adult No 89
16 Crew Female Adult No 3
17 1st Male Child Yes 5
18 2nd Male Child Yes 11
19 3rd Male Child Yes 13
20 Crew Male Child Yes 0
21 1st Female Child Yes 1
22 2nd Female Child Yes 13
23 3rd Female Child Yes 14
24 Crew Female Child Yes 0
25 1st Male Adult Yes 57
26 2nd Male Adult Yes 14
27 3rd Male Adult Yes 75
28 Crew Male Adult Yes 192
29 1st Female Adult Yes 140
30 2nd Female Adult Yes 80
31 3rd Female Adult Yes 76
32 Crew Female Adult Yes 20

После чего логлиненый анализ выполняется так:

> glm.model = glm(Freq ~ Class * Age * Sex * Survived, data=ti, family=poisson)

Это — насыщенная модель. Команда summary(glm.model) вывалила бы нам огромный и запутанный кусок информации, но я предпочитаю пользоваться...

> anova(glm.model, test="Chisq")
Analysis of Deviance Table

Model: poisson, link: log

Response: Freq

Terms added sequentially (first to last)

Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL 31 4953.1
Class 3 475.8 28 4477.3 8.326e-103
Age 1 2183.6 27 2293.8 0.0
Sex 1 768.3 26 1525.4 4.169e-169
Survived 1 281.8 25 1243.7 3.078e-63
Class:Age 3 148.3 22 1095.3 6.048e-32
Class:Sex 3 412.6 19 682.7 4.126e-89
Age:Sex 1 6.1 18 676.6 1.363e-02
Class:Survived 3 180.9 15 495.7 5.634e-39
Age:Survived 1 25.6 14 470.2 4.237e-07
Sex:Survived 1 353.6 13 116.6 7.053e-79
Class:Age:Sex 3 4.0 10 112.6 0.3
Class:Age:Survived 3 35.7 7 76.9 8.825e-08
Class:Sex:Survived 3 75.2 4 1.7 3.253e-16
Age:Sex:Survived 1 1.7 3 4.237e-10 0.2
Class:Age:Sex:Survived 3 0.0 0 4.463e-10 1.0

Девиация (deviance) является терминологическим вариантом критерия хи-квадрат для отношения правдоподобия. Нулевая модель, в которой все частоты одинаковы, находится вверху. Очевидно, что нулевая модель не подходит для описания этого набора данных, поскольку добавление следующих эффектов значительно уменьшает девиацию. При переходе от верхних к нижним строкам в таблице добавляется другие предикторы и это уменьшает значение девиации, указанное в столбце Deviance. Значимость p-value в последнем столбце основывается на изменениях во втором столбце (Deviance) и Df (первый столбец), происходящих при добавлении новых эффектов в модель. Здесь мы видим значимое уменьшение девиаций, свидетельствующее о том, что каждый новый эффект помогает предсказывать наблюдаемые частоты.

Эта таблица даёт возможность предположить, что целесообразно исследовать исключение трёх сочетаний переменных: четырёхпеременное взаимодействие и два трёхпеременных. Вспомним, что команда step() сохранила взаимодействие Class:Age:Sex ...

> anova(update(glm.model,.~.-(Class:Age:Sex:Survived+Age:Sex:Survived
+ +Class:Age:Sex)),test="Chisq")
Analysis of Deviance Table

Model: poisson, link: log

Response: Freq

Terms added sequentially (first to last)

Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL 31 4953.1
Class 3 475.8 28 4477.3 8.326e-103
Age 1 2183.6 27 2293.8 0.0
Sex 1 768.3 26 1525.4 4.169e-169
Survived 1 281.8 25 1243.7 3.078e-63
Class:Age 3 148.3 22 1095.3 6.048e-32
Class:Sex 3 412.6 19 682.7 4.126e-89
Age:Sex 1 6.1 18 676.6 1.363e-02
Class:Survived 3 180.9 15 495.7 5.634e-39
Age:Survived 1 25.6 14 470.2 4.237e-07
Sex:Survived 1 353.6 13 116.6 7.053e-79
Class:Age:Survived 3 29.2 10 87.4 2.024e-06
Class:Sex:Survived 3 65.4 7 22.0 4.066e-14
> 1-pchisq(22,df=7)
[1] 0.002540414

В этой модели последняя строка содержит остаточное отклонение (residual deviance), равное 22 при 7 степенях свободы. Судя по значению p = .0025, это неподходящая модель. Поэтому один или больше удалённых эффектов должны быть возвращены назад.

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

Последнее замечание

Логлинейные модели рассчитываются с использованием итеративного алгоритма. Соответственно, этот метод можно рассматривать как «машинно-нагруженный». Различные программные пакеты по разному реализуют этот алгоритм или даже используют совсем другие алгоритмы. Поэтому результаты, полученные в разных программах, могут не совпадать с абсолютой точностью. К счастью, они всё-таки будут похожи!


Читать далее