0. Вводные замечания

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

Был у нас и небольшой опыт анализа данных, включенных в основной модуль (base R) и построения простых графиков.

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

Итак, начнем…

1. Установка необходимых библиотек

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

В качестве небольшого отступления: Кто эти люди, создающие R? Кто эти герои?

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

ТОП-3

  1. Хэдли Уикхэм (Hadley Wickham) - программист из Новой Зеландии, главный научный сотрудник в компании RStudio. Без преувеличения самый известнй, очень трудолюбивый и плодотворный разработчик R, автор почти 100 (!) библиотек, в том числе таких мегапопулярных как ggplot, dplyr, tidyverse, devtools и других. Этими библиотеками воспользовались более 825 тыс. человек.

  1. Дирк Эддельбуэттель (Dirk Eddelbuettel) - канадский ученый в области статистики, программист и сследователь, второй в нашем списке лучших разработчиков библиотек для R. Принял участие в создании более 60 библиотек, наиболее популярными среди которых является Rcpp, позволяющая интегрировать R с другим не менее популярным языком C++, а также RPostgreSQL, предоставляющая интерфейс для работы с системой баз данных PostgreSQL. Работает в Иллинойском университете в Урбане-Шампейне.

  1. Ихуэй Се (Yihui Xie) - еще один плодовитый разрботчик, создавший более 40 библиотек, загруженных более 130 тыс. раз. Его большая заслуга в том, что разработал инструменты для создания интерактивных приложений (knitr, rmarkdown, shiny и htmlwidgets - его детища, и мы тоже с ними познакомимся на следующих занятиях). Yihui Xie также поддерживает систему bookdown, которую можно использовать для написания книг и отчетных документов с помощью R Markdown. Инженер-программист компании RStudio.

Вернемся к нашему разговору.

Прежде, чем загрузить какую-либо библиотеку, ее нужно сначала установить, для этого существует функция install.packages(). Можно также воспользоваться меню: Tools - Install packages.

Нам сегодня понадобится несколько библиотек:

  • haven - позволяет импортировать данные в формате SPSS (что для нас очень важно) и сохранять их в формате .sav (формат файла с данными в SPSS)
  • dplyr - одна из лучших для всевозможных манипуляций над данными - выбора переменных, отбора наблюдений, вычисления различных статистик и создания новых переменных
  • sjmisc - содержит коллекцию функций для трансформации данных, нам полезна для создания таблиц
  • sPlot - коллекция функций для создания таблиц и графиков
  • sjlabelled - нужна для работы с именованными сущностями и очень пригождается при работе с данными, которые созданы в других программах для статистической обработки, таких как SAS или SPSS.
  • likert - создана для анализа и визуального представления шкал Лайкерта (шкал согласия)
  • vcd - создана для визуализации категориальных переменных

Установим их на свой компьютер, введя следующий код: install.packages(c("haven", "dplyr", "sjmisc","sjPlot", "ggpubr", "likert", "sjlabelled", "vcd")).

Далее загрузим эти библиотеки для работы с помощью функции library()

library(haven)
library(dplyr)
library(sjmisc)
library(sjPlot)
library(ggpubr)
library(likert)
library(sjlabelled)
library(vcd)

2. Открываем базу данных SPSS в R

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

Эту базу нужно скачать себе на компьютер в папку. Потом в эту же папку нужно сохранить файл с кодом в R, который Вы сейчас читаете - Lecture 6.Rmd, и установить рабочую директорию в ту папку, где этот код хранится: Session - Set Working Directiory - To Source File Location, чтобы затем не пришлось детально прописывать путь к файлу с базой данных.

Открываем файл с помощью фукции read_sav из библиотеки haven:

df<-read_sav("База_НШ_ВМ_2020_от 21.12.2020.sav", user_na = TRUE)

Объект df появился в вашем рабочем пространстве, и можно увидеть, что он содержит 2538 наблюдений и 221 переменную. Отлично, полдела сделано!

Прежде чем начинать какой-то анализ, необходимо посмотреть, какие переменные (вопросы и отдельные значения) имеются в нашей таблице данных, как эти значения закодированы. Мы, конечно, можем обратиться к документу, содержащему вопросы анкеты, но это значит, что нужно открывать какие-то папки и другие программы, что требует времени. К счастью, в R есть возможность это сделать “не отходя от кассы”:

Должна появиться примерно такая таблица, где можно увидеть, как закодирован вопрос (Name), увидеть метку для этого вопроса (Label), сами значения - коды (Values) и метки значений (Value Labels):

3.Учимся проводить первичный анализ

Наша первая задача - посмотреть одномерные распределения по некоторым интересующим нас вопросам. К примеру, возьмем Q8 - табличный вопрос, исследующий националистические установки и Q9 - “Рискуете ли Вы утратить свою национальную идентичность (раствориться в культуре других народов, забыть национальные традиции, язык)?”. Начнем с Q9, так как в нем только одно утверждение.

Самая простая функция для такого анализа из базового R - table().

Давайте создадим таблицу с частотами по вопросу Q9:

table(df$Q9)
## 
##    1    2    3    4 
##  186  317  744 1276

Ответы есть, но это нам мало что дает, потому что мы привыкли видеть проценты. Что же делать? Для того, чтобы посмотреть относительные частоты, можно немного “поплясать с бубном” и “обернуть” функцию table() в функцию prop.table():

prop.table(table(df$Q9))
## 
##          1          2          3          4 
## 0,07372176 0,12564407 0,29488704 0,50574713

Уже лучше, но это доли единицы, а не проценты. Не беда, нужно просто умножить на 100 и округлить до одной десятой с помощью функции round() - теперь у нас свернуты целых три функции:

round(prop.table(table(df$Q9))*100,1)
## 
##    1    2    3    4 
##  7,4 12,6 29,5 50,6

Фактически это то, что нужно - есть проценты, но нет меток значений.

Зачем так сложно? Почему нельзя это сделать с помощью той же таблицы SPSS, где сразу будет все? Можно, конечно, и в R есть специальные библиотеки, имитирующие стиль работы других программ, не только SPSS. Но на самом деле, базовые функции - очень простые, и иногда они работают быстрее и лучше.Если бы мы не раскладывали по полочкам, по сути, нам понадобилась только одна строка кода.

Давайте попробуем сделать эту же таблицу частот, но уже с помощью функции frq из библиотеки sjmisc:

  frq(df$Q9)
## 
## Рискуете ли Вы утратить свою национальную идентичность (раствориться в культуре других народов, забыть национальные традиции, язык)? (x) <numeric>
## # total N=2538  valid N=2523  mean=3.23  sd=0.93
## 
## Value |                Label |    N | Raw % | Valid % | Cum. %
## --------------------------------------------------------------
##     1 |   Определенно рискую |  186 |  7.33 |    7.37 |   7.37
##     2 |        Скорее рискую |  317 | 12.49 |   12.56 |  19.94
##     3 |     Скорее не рискую |  744 | 29.31 |   29.49 |  49.43
##     4 | Совершенно не рискую | 1276 | 50.28 |   50.57 | 100.00
##  <NA> |                 <NA> |   15 |  0.59 |    <NA> |   <NA>

Тут есть все - метки значений, количество наблюдений и проценты, в том числе учитывающие пропущенные значения (NA). Есть и кумулятивные частоты, позволяющие суммировать проценты по предшествующим строкам, и иногда это может пригоодиться (мы можем, например, сразу же сказать, что 19,94% респондентов рискуют утратить свою национальную идентичность).

Сделаем график для визуализации ответов под этот вопрос:

# установим минимальную тему для графика, без лишних записей и фоновых линий
set_theme(
  base = theme_classic(), 
  axis.linecolor = "white",     # удаляет линии осей
  axis.textcolor.y = "white", # делаем белыми надписи возле осей
  axis.textcolor.x = "black",# делаем белыми надписи возле осей
  axis.tickslen = 0,            # удаляем "засечки"
  axis.title.size = .9,#устанавливаем размер шрифта для названия
  axis.textsize = .9,#устанвливаем размер шрифта для подписей
  legend.size = .7,# размер самой легенды
  geom.label.size = 4 #размер подписей для данных
)
plot_frq(df, Q9, show.n = FALSE)#show.n = FALSE - отключает вывод количества наблюдений

Получился вполне приличный график, который не стыдно вставить в свой отчет или статью.

Что дальше?

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

Для этого нам нужно перекодировать исходную переменную так, чтобы первые два варианта ответа образовывали одну категорию - “Рискую”, а вторые два - другую категорию “Не рискую”.

Воспользуемся библиотекой dplyr и функцией mutate(), которая позволяет создать новые переменные (подробнее об этой функции с примерами: https://dplyr.tidyverse.org/reference/mutate.html).

! Обратите внимание на знак %>% в коде ниже - это так называемый pipe operator - оператор пайп (от англ. pipe - трубка), который позволяет “нанизывать” друг на друга разные функции, что существенно сокращает объем кода и увеличивает скорость его написания.

df<-df %>% mutate(Q9n = case_when(Q9==1 | Q9==2   ~ "Рискую", Q9==3 | Q9==4   ~ "Не рискую")) #Знак | обозначает логический оператор "или"
frq(df, Q9n)
## 
## Q9n <character>
## # total N=2538  valid N=2523  mean=1.20  sd=0.40
## 
## Value     |    N | Raw % | Valid % | Cum. %
## -------------------------------------------
## Не рискую | 2020 | 79.59 |   80.06 |  80.06
## Рискую    |  503 | 19.82 |   19.94 | 100.00
## <NA>      |   15 |  0.59 |    <NA> |   <NA>

Видим, что 80% респондентов не испытывают опасений по поводу утраты своей идентичности и 20% - испытывают.

Делать какой-то график две цифры не очень целесообразно, но ради практики сделаем вариацию круговой диаграммы - кольцевую диаграмму - donut chart (“пончик”):

#сначала создадим таблицу частот и сохраним ее как отдельный объект, удалив данные о пропущенных значениях
table<-as.data.frame(frq(df$Q9n, show.na = FALSE))

Объект table появился в рабочем пространстве, если кликнуть на него, то можно увидеть, что это таблица частот:

Переменные этой таблицы следующие:

val - варианты ответа

label - какие-то метки (здесь их нет)

frq - частоты - количество наблюдений

raw.prc - проценты по строке

cum.prc кумулятивные проценты

Теперь на основе этой таблицы создадим график:

ggdonutchart(table, "raw.prc",
                     label = "raw.prc", color="val",
                     fill = "val",lab.pos = "in", borders ='n', 
                    palette = c("dodgerblue3", "darkorange2"),legend.title = "")

Давайте разберем код, который был использован для создания графика:

ggdonutchart - функция, которая создает кольцевую диаграмму (пончик)

Что в скобках?

table - источник данных, наша таблица частот

“raw.prc” - указываем данные, на основе которых строится диаграмма

label = “raw.prc” - просим отобразить проценты на дольках нашего “пончика”

color=“val” - переменная, на основе которой программа определяет количество цветов, у нас это варианты ответа, их два, значит, требуется два цвета

fill = “val” - аргумент, указывающий, каким цветом заполнит “куски” нашего “пончика” - то есть отдельные области диаграммы, если этот аргумент убрать - внутри будет пустота, просто белый цвет

lab.pos = “in” - где будут располагаться надписи, in - внутри каждой области

borders =‘n’ - убираем границы, если убрать этот аргумент, то у диаграммы будет черная “обводка”

palette = c(“dodgerblue3”, “darkorange2”)- задаем цвета для элементов диаграммы - синий и оранжевый

Разбираем второй вопрос - Q8

Второй вопрос - интереснее, это три утверждения с разной степенью согласия. Что это за вопросы и какие у них ответы?

Воспользуемся функцией get_label из библиотеки sjlabelled:

df %>% 
  get_label(starts_with("Q8_SQ"))
##                                                                                                                                       Q8_SQ001 
## "[В любой стране власть должна в основном находиться в руках представителей коренной национальности] Насколько Вы согласны или не согласны <U+FFFD>" 
##                                                                                                                                       Q8_SQ002 
## "[Для человека естественно и правильно думать, что его нация лучше остальных] Насколько Вы согласны или не согласны с каждым из следующих утв" 
##                                                                                                                                       Q8_SQ003 
##  "[Национальные меньшинства имеют слишком много власти и влияния в России] Насколько Вы согласны или не согласны с каждым из следующих утверж"
get_labels(df$Q8_SQ001)
## [1] "Совсем не согласен" "Скорее не согласен" "Скорее согласен"   
## [4] "Полностью согласен"

Как свести их воедино и сделать красивую визуализацию так, чтобы сразу все стало понятно?

Воспользуемся библиотекой likert.

Предварительно создадим мини-таблицу, где будет только этот вопрос, вернее три вопроса, которые оцениваются одинаково:

Q8<-df %>% select(Q8_SQ001,Q8_SQ002, Q8_SQ003) # функция select() - отбирает нужные нам переменные
Q8<-as.data.frame(Q8)# переведем в формат таблицы данных

Теперь, проанализируем вопрос с помощью функции likert() из библиотеки likert (да, фантазии иногда у разработчиков маловато, с именами просто беда:

#Переведем наши переменные в факторные с четырьмя уровнями
Q8[,1:3]<- lapply(Q8[,1:3], factor, levels=1:4, labels =c("Совсем не согласен", "Скорее не согласен", "Скорее согласен", "Полностью согласен"))
colnames(Q8)<-c("В любой стране власть должна в основном находиться в руках представителей коренной национальности", "Для человека естественно и правильно думать, что его нация лучше остальных", "Национальные меньшинства имеют слишком много власти и влияния в России")#присвоим столбцам имена меток
# Обработаем с помощью функции likert (она позволяет сделать таблицу сразу по всем вопросам)
Q8_likert<-likert(Q8, nlevels = 4)
Q8_likert
##                                                                                                Item
## 1 В любой стране власть должна в основном находиться в руках представителей коренной национальности
## 2                        Для человека естественно и правильно думать, что его нация лучше остальных
## 3                            Национальные меньшинства имеют слишком много власти и влияния в России
##   Совсем не согласен Скорее не согласен Скорее согласен Полностью согласен
## 1           15,91270           28,45238        36,07143          19,563492
## 2           40,75547           33,95626        19,52286           5,765408
## 3           26,02631           43,56317        21,68194           8,728577

Описывать можно и таблицу, но мы сделаем еще и визуализацию:

set_theme(base=theme_classic())#вернем тему на классическую
plot(Q8_likert, centered=T, wrap=30)

Полезные ссылки о библиотеке `likert’:

https://github.com/jbryer/likert/blob/master/demo/likert.R

4. Сравниваем два вопроса вместе

Мы посмотрели ответы на оба вопроса, а теперь хотим понять, связаны ли как-то оценки об утрате идентичности с националистическими убеждениями. Для этого давайте перекодируем первый подвопрос вопроса Q8 в бинарный вид (Согласен - Не согласен) и таблицу сопряженности, а также оценим, насколько значимы различия в двух группах (на основе статистики хи-квадрат Пирсона).

df<-df %>% mutate(Q8_s1 = case_when(Q8_SQ001==1 | Q8_SQ001==2   ~ "Не согласен", Q8_SQ001==3 | Q8_SQ001==4   ~ "Согласен"))
sjt.xtab(df$Q8_s1, df$Q9n, show.col.prc=T,  var.labels = c("Мнения о власти","Потеря идентичности"), encoding="UTF-8",  show.exp = F, show.obs=F)

Мы видим из таблицы, что среди рискующих потерять свою идентичность 61,2% согласны с тем, что власть должна быть в руках какой-то одной национальности, очевидно, их собственной (то есть демонстрируют негативное отношение), тогда как среди тех, кто не рисует - доля согласившихся составляет 54,2% - и это статистически достоверное различие по статистике хи-квадрат, которая благодаря функции sjt.xtab выводится непосредственно под таблицей (так гораздо удобнее), p < 0,006.

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

Сделаем мы такую диаграмму с помощью библиотеки vcd - аббревиатура которой расшифровывается как visualization of categorical data. Как раз то, что нам нужно!

mosaic(~Q9n+Q8_s1, data = df, direction = c("v", "h"),   gp = shading_max, labeling_args = list(set_varnames = c(Q9n="Риск утраты идентичности", Q8_s1="Мнение о власти")))

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

P.S.: Мы сделали только несколько самых первых шагов, чтобы проанализировать данные в R. Для того, чтобы полностью погрузиться в анализ, требуется время и, главное, постоянная практика.

Вашим домашним заданием будет воспроизведение кода урока и построение дополнительных таблиц и мозаичных диаграмм по 2 и 3 подвопросу вопроса 8.