Файлы для работы можно скачать здесь.

Kачественные данные

Возьмем базу данных, в которой содержатся результаты опроса, проводившегося в 2006 году в рамках исследования национализма и антисемитизма в современной Германии (более полные данные обсуждались на курсе Structural Equation Modelling, ECPR Summer School, 2014).

df <- read.csv('C:/Users/AllaT/Desktop/Rclass/survey06.csv')
View(df)

Большинство вопросов – вопросы вида “В какой степени вы согласны с утверждением…”. Соответственно, есть 7 градаций – степеней согласия (1 - полностью не согласен, 7 - полностью согласен). Для удобства разделим все эти степени на 2 группы: от 1 до 3 - не согласен, от 4 до 7 - согласен). Деление, конечно, очень условное. Для примера возьмем переменную j_shame - “чувство стыда за нацистское прошлое страны”.

df$j_shame_bin <- rep(NA, length(df$j_shame)) # как всегда, создадим переменную из NA

# и ее заполним

for (i in 1:length(df$j_shame)){
if (df$j_shame[i]<4)  df$j_shame_bin[i] <- 'disagree'
else df$j_shame_bin[i] <- 'agree'
}

Tаблицы сопряженности и хи-квадрат

Построим таблицы сопряженности (contigency tables) для двух признаков: пол и чувство стыда за нацистское прошлое.

table(df$gender, df$j_shame_bin) # в скобках - два признака
##    
##     agree disagree
##   0    62      181
##   1    61      196

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

chisq.test(table(df$gender, df$j_shame_bin))
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(df$gender, df$j_shame_bin)
## X-squared = 0.128, df = 1, p-value = 0.7205

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

Посмотрим, что еще можем вытащить из результатов теста:

summary(chisq.test(table(df$gender, df$j_shame_bin)))
##           Length Class  Mode     
## statistic 1      -none- numeric  
## parameter 1      -none- numeric  
## p.value   1      -none- numeric  
## method    1      -none- character
## data.name 1      -none- character
## observed  4      table  numeric  
## expected  4      -none- numeric  
## residuals 4      table  numeric  
## stdres    4      table  numeric
chisq.test(table(df$gender, df$j_shame_bin))$observed # наблюдаемые частоты
##    
##     agree disagree
##   0    62      181
##   1    61      196
chisq.test(table(df$gender, df$j_shame_bin))$expected # ожидаемые частоты
##    
##      agree disagree
##   0 59.778  183.222
##   1 63.222  193.778
chisq.test(table(df$gender, df$j_shame_bin))$stdres # стандартизированные остатки
##    
##          agree   disagree
##   0  0.4616429 -0.4616429
##   1 -0.4616429  0.4616429

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

Графики для качественных данных: bar plot и pie chart

Столбчатая диаграмма (bar plot)

Для начала построим обычную столбчатую диаграмму (bar plot) для возрастных групп респондентов. Сначала надо сделать таблицу с разбиением на группы:

groups <- table(df$age_group) 
groups
## 
##   1   2   3   4   5   6 
##  62 135 141 127  33   2
barplot(groups) # строится для групп (table), а не для просто переменной

Теперь начинаем приводить график в приличный вид. С помощью ? узнаем, что можно изменить.

?barplot # help
## starting httpd help server ...
##  done

Для начала поменяем цвет графика.

barplot(groups, col = 'navy')

А теперь зададим несколько цветов (хотя бы два):

colors = c('navy', 'red', 'red', 'red', 'navy', 'navy') # надо прописать вектор
barplot(groups, col = colors)

Кроме того, можем задать цвет границ столбиков:

barplot(groups, col = colors, border = 'white') # border

Добавим названия осей (x и y):

barplot(groups, col = colors, xlab = 'groups', ylab = 'number of respondents')

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

labs = c('18-29', '30-44', '45-59', '60-74', '75-89', 'more 90') # названия групп

# names.arg - подписи для каждого столбика 

barplot(groups, col = colors, xlab = 'age, years', ylab = 'number of respondents', names.arg = labs)

Осталось дать графику название.

# main - основной заголовок, название графика

barplot(groups, col = colors, xlab = 'age, years', ylab = 'number of respondents',names.arg = labs, main = 'Age of respondents')

Вроде минимальные требования к хорошему графику мы учли. Что бы еще сделать?

Можем что-нибудь попереворачивать.

barplot(groups, col = colors, xlab = 'age, years', ylab = 'number of respondents',names.arg = labs,
        main = 'Age of respondents', las = 2) 

# las = 2 - означает, что названия групп  (под каждым столбиком) должны идти вертикально. 
# Если ничего не пишем или пишем las = 1, то подписи по умолчанию остаются горизонтально.

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

А теперь повернем сам график.

barplot(groups, col = colors, xlab = 'age, years', ylab = 'number of respondents',names.arg = labs, main='Age of respondents', horiz = TRUE) 

# horiz = TRUE - столбцы будут идти горизонтально
# в этом случае поворот тоже не очень удачный

Круговая диаграмма (pie chart)

Построим для тех же данных круговую диаграмму.

pie(groups) # тоже строим для таблицы частот, не просто для переменной

Пока выглядит грустно. Начинаем исправлять.

colors <- c('snow', 'royalblue', 'slategray', 'thistle', 'slateblue', 'palevioletred')

# да, я специально подобрала нетипичные названия цветов

pie(groups, col = colors)

Кстати, для любителей эстетики, ссылка на таблицу цветов в R.

Аналогично предыдущему графику, добавим необходимые подписи:

labs = c('18-29', '30-44', '45-59', '60-74', '75-89', 'more 90')

# в отличие от barplot, для надписей используем аргумент labels

pie(groups, col = colors, main = 'Age of respondents', labels = labs)

Сильно красивее не стало.

Поэтому сделаем по-другому: на самом графике подпишем проценты (% респондентов в каждой возрастной группе), а остальное вынесем в легенду.

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

groups
## 
##   1   2   3   4   5   6 
##  62 135 141 127  33   2

Поделим каждый элемент в groups на общее число респондентов (sum) и умножим на 100 - получим проценты. Потом для удобства округлим до первого знака после запятой:

perc <- round(groups/sum(groups)*100, 1)
perc
## 
##    1    2    3    4    5    6 
## 12.4 27.0 28.2 25.4  6.6  0.4

Теперь добавим знаки % - создадим вектор надписей.

labs <- paste(perc, '%', sep = '') # paste - присоединим % ко всем элементам perc и сделаем так, чтобы между ними ничего не было - пустой разделитель sep=''
labs
## [1] "12.4%" "27%"   "28.2%" "25.4%" "6.6%"  "0.4%"

Наконец, возвращаемся к самому графику. Добавляем легенду: позиция на графике, список названий, cex - размер, fill - параметр для цветов, должен совпадать с вектором цветов.

pie(groups, col = colors, main = 'Age of respondents', labels = labs)

legend("topright", c('18-29','30-44','45-59','60-74','75-89','more 90'), cex = 0.8, fill = colors)

Положение легенды на графике можно задавать двумя способами: числами (см.?legend) или выбрать одну из стандартных позиций:

topright - верхний правый угол

topleft - верхний левый угол

bottomright - нижний правый угол

bottomleft - нижний левый угол

Как мы уже убедились, смотреть на графики, используя окно Zoom не всегда удобно. Более того, выгружать из него графики в таком виде тоже печально. Поэтому просто сохраним график в формате png:

# сначала указываем формат, потом - имя файла (путь к файлу и его имя)

dev.copy(png,'C:/Users/AllaT/Desktop/pie.png')
## png 
##   3
# перечисляем, что строим
pie(groups, col = colors, main = 'Age of respondents', labels = labs)
legend("topright", c('18-29','30-44','45-59','60-74','75-89','more 90'), cex = 0.8, fill = colors)

dev.off() # закрываем файл
## png 
##   2

А можно изобразить еще такую симпатичную (хотя для кого как) круговую 3D-диаграмму:

# для этого потребуется пакет plotrix

install.packages('plotrix')
library(plotrix)

pie3D(groups,col = colors, main = 'Age of respondents', labels = labs)

Теперь вернемся к количественным данным. Но для начала несколько базовых вещей. Как вообще строятся обычные графики (функций, плотности распределения и другие)?

Графики функций

Математические функции

Построим, например, график \(y=cos(x)\):

curve(cos, xlim = c(-pi, 3*pi), col ="navy") 

# cos - сама функция
# xlim - границы по оси Ox
# ylim - границы по оси Oy
# col - цвет

Со встроенными функциями (cos, sin, log и др.) все просто. Если нужно построить график более сложной, “своей”, функции, ее надо сначала задать (как это делать, мы уже обсуждали).

f <- function(x){x^2-3*x+1} # задаем функцию
curve(f, xlim = c(-10, 10), ylim = c(-3, 3), col = "navy") # строим график

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

Пример:

y <- function(x){sin(1/x)} 
curve(y, xlim = c(-1, 1), ylim = c(-1, 1),col = "navy") 
## Warning in sin(1/x): NaNs produced

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

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

# n - число точек, в которых R будет рассчитывать значение функции

curve(y, xlim = c(-1, 1), ylim = c(-1,1), col  = "navy", n = 2000)

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

Графики функции плотности распределения и функции распределения

Построим графики плотностей часто используемых распределений: нормальное распределение, распределение Стьюдента (t-распределение) и хи-квадрат распределение.

Note: английские эквиваленты: функция плотности распределения (вероятности) - probability density function (pdf), функция распределения - cumulative density function (cdf).

# dnorm, d - от density

curve(dnorm(x), xlim = c(-5, 5)) # стандартное нормальное распределение

curve(dnorm(x, mean = 2, sd = 0.5), xlim = c(-5, 5)) # нормальное распределение, среднее равно 2, ст.отклонение равно 0.5

curve(dt(x, df = 2), xlim = c(-5, 5)) # распределение Стьюдента с 2 степенями свободы (df=2)

curve(dchisq(x, df = 1), xlim = c(0, 5)) # хи-квадрат распределение с одной степенью свободы (df=1)

Теперь мы хотим показать три плотности распределения на одном графике:

curve(dnorm(x), xlim = c(-5, 5), col = "green")

# add=TRUE показывает, что мы хотим изобразить еще одну кривую на том же графике

curve(dt(x, df = 2), xlim = c(-5, 5), col = "blue", add = TRUE) 
curve(dchisq(x, df = 3), xlim = c(0, 5), col = "red", add = TRUE)

Хотя график самой функции распределения нужен довольно редко, пусть будет:

curve(pnorm(x), xlim = c(-5, 5)) # стандартное нормальное распределение

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

dt<-read.csv('C:/Users/AllaT/Desktop/Rclass/elections_Leningrad.csv')
View(dt)

Гистограммы (histograms)

Для начала построим самую простую гистограмму.

hist(dt$turnout)

Начинаем приводить ее в порядок:

hist(dt$turnout, breaks = 20, freq = TRUE, col = "lightgreen", 
     main = "Явка", xlab = "явка", ylab = "частоты")

# freq=TRUE (по оси Оy - частоты, а не вероятности)
# col - цвет
# main - главный заголовок (название графика)
# xlab - название оси Ox
# ylab - название оси Oy

Теперь хотим сделать то же самое, но чтобы вместо частот по оси Oy были указаны вероятности:

hist(dt$turnout, breaks = 20, freq = FALSE, col = "lightgreen", 
     main = "Явка", xlab="явка", ylab = "вероятности")

А еще мы можем добавить график плотности нормального распределения, чтобы посмотреть, распределены ли данные нормально:

hist(dt$turnout, breaks = 20, freq = FALSE, col = "lightgreen", 
     main = "Явка", xlab="явка", ylab = "вероятности")

# в качестве параметров распределения мы берем среднее и ст.отклонение переменной, которую изучаем

curve(dnorm(x, mean = mean(dt$turnout), sd = sd(dt$turnout)), add = TRUE, col = "navy", lwd = 2)

# lwd - толщина линии

Кроме того, мы можем просто построить график плотности распределения переменной (аналог kdensity в Stata). Для этого нам потребуется пакет car.

install.packages('car')
library(car)
density(dt$turnout) # функция распределения переменной turnout
## 
## Call:
##  density.default(x = dt$turnout)
## 
## Data: dt$turnout (963 obs.); Bandwidth 'bw' = 4.298
## 
##        x                  y            
##  Min.   : -0.5133   Min.   :1.779e-06  
##  1st Qu.: 27.8388   1st Qu.:1.197e-03  
##  Median : 56.1910   Median :7.839e-03  
##  Mean   : 56.1910   Mean   :8.809e-03  
##  3rd Qu.: 84.5432   3rd Qu.:1.585e-02  
##  Max.   :112.8953   Max.   :2.102e-02
plot(density(dt$turnout), main = 'Density') # график плотности распределения

Ящики с усами (boxplots)

“Ящик с усами” для процента голосов за кандидата (R допускает русскоязычные подписи к графикам, но иногда с их отображением могут возникать проблемы).

# sub - подзаголовок

boxplot(dt$Drozdenko, col = 'greenyellow', main = "Процент голосов",
        sub = 'Дрозденко ("Единая Россия")')

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

# разобьем переменную turnout (явка) по значениям группирующей переменной OIK (округ)

boxplot(split(dt$turnout, dt$OIK), col='yellow')

Возникли проблемы с русскоязычными надписями. Исправим:

labels <- c(1, 2, 3, 5, 6, 7, 11, 12, 13, 15, 16, 17, 20, 21, 24, 25, 26, 27) 

# names - подписи под каждым boxplot'ом

boxplot(split(dt$turnout, dt$OIK), names = labels, col = "yellow")

Добавим красоты к нашим “yellow submarines”:

# в последней строчке cex.lab=0.75 означает, 
# что мы хотим, чтобы шрифт у подписей к осям составлял 75% стандартного размера
# аналогично для заголовка и подписей к отдельным boxplot'ам

boxplot(split(dt$turnout,dt$OIK), names = labels, col = 'yellow',
main = 'Выборы губернатора Ленинградской области',
xlab = 'номер окружной избирательной комиссии',
ylab = 'явка, %',
cex.lab = 0.75, cex.title = 0.55, cex.names = 0.75)