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

Работа с файлами

В прошлый раз мы разобрали, как можно создать базу данных (data.frame) в R и как работать с переменными. Теперь посмотрим, как сохранить созданную базу и как открывать файлы разных форматов.

Возьмем базу, которую мы создали в прошлый раз.

id <- paste('respondent', 1:5, sep = '_')
gender <- c('female', 'male', 'female', 'male', 'female')
gender_bin <- c(0, 1, 0, 1, 1)
age <- c(51, 20, 35, 15, 19)
vote <- c('a', 'b', 'c', NA, 'b')
df <- data.frame(id = id, gender = gender, gender_bin = gender_bin, age = age, votes = vote) 
df
##             id gender gender_bin age votes
## 1 respondent_1 female          0  51     a
## 2 respondent_2   male          1  20     b
## 3 respondent_3 female          0  35     c
## 4 respondent_4   male          1  15  <NA>
## 5 respondent_5 female          1  19     b

Сохранение файлов

Теперь сохраним ее в формате Rda (формат R, от Rdata).

#save(df, file = "C:/Users/AllaT/Desktop/data.Rda")

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

getwd() # wd от working directory
## [1] "/home/oem/Dropbox/ВШЭ - работа/rpubs"

Естественно, ее можно поменять:

#setwd("C:/Users/AllaT/Desktop") # сделать папку Desktop рабочей
save(df, file="data.Rda") # сохранить

Допустим, теперь мы хотим сохранить базу в других форматах (не относящихся к R). Для этого нам потребуется загрузить два пакета.

install.packages("foreign") # Stata, SPSS, SAS
install.packages("xlsx") # Excel

Для того, чтобы использовать функции из пакетов, нужно не только установить сами пакеты, но и обратиться к соответствующей библиотеке. Если не обратиться к библиотеке и не прописать library(), R не будет понимать команды, которые относятся к конкретному пакету, и будет выдавать ошибку: “could not find function …”.

Сохраним нашу базу в формате dta (Stata).

library(foreign) # обращаемся к пакету
write.dta(df, "data.dta")

Аналогично для других форматов. Для создания xlsx файла нужен пакет xlsx.

library(xlsx)
write.xlsx(df, "data.xlsx")

Для того, чтобы сохранить базу в формате csv или txt, дополнительных пакетов не требуется.

write.csv(df, "data.csv")
write.table(df, "data.txt")

Загрузка файлов

Базы R

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

library(MASS)
data()

# конечно, это не всё; многие пакеты имеют свои "вшитые" базы данных

Загружать встроенные базы очень просто. Возьмем, например, базу beaver1 - базу, которая содержит результаты измерения температуры тела бобров в течение двух дней (измерения производились каждые 10 минут).

df <- beaver1
View(df)

Описание переменных из этой базы данных можно посмотреть здесь.

Если говорить о социально-экономических данных, то в R можно подключить пакет wbstats, который позволяет загружать данные Всемирного Банка (World Bank).

install.packages("wbstats")
library(wbstats)

Посмотрим на изменение доли экспорта топлива во всем объеме экспорта в России в 1996-2015 (fuel exports, %). Найдем сам показатель, а затем, узнав его код, выставим остальные параметры - страну и интересующий нас промежуток времени:

wbsearch(pattern = "fuel export", fields="indicator") # ищем точное название (код) показателя
##            indicatorID                               indicator
## 1923 TX.VAL.FUEL.ZS.UN Fuel exports (% of merchandise exports)
fuel <- wb(country = c("RU"), indicator = c("TX.VAL.FUEL.ZS.UN"), startdate = 1996, enddate = 2015) 
View(fuel)

Более подробно про возможности пакета wbstats можно почитать здесь.

Загрузка файлов

Открывать файлы разных форматов можно следующим образом:

STATA

library(foreign)
dta <- read.dta('USArrests.dta') # Stata

Excel

library(xlsx)
xlsx <- read.xlsx('USArrests.xlsx', 1) # Excel

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

Опять, если нам нужно открыть файлы в форматах csv или txt, дополнительных пакетов не требуется.

csv <- read.csv('USArrests.csv') # csv

При работе с txt-файлами необходимо указывать, каким образом столбцы отделены друг от друга (аргумент sep, в данном случае столбцы отделены друг от друга табуляцией), а также учитывать, что представляет собой первая строка: наблюдение или шапку таблицы (имена переменных) – аргумент header. Сравним:

txt <- read.table('USArrests.txt', sep='\t', header = TRUE)  
View(txt) # header=T - первая строка читается как имя переменной

txt2 <- read.table('USArrests.txt', sep='\t', header = FALSE)
View(txt2) # header=F - первая строка читается как наблюдение

Конвертация файлов

Для конвертации файлов необходим пакет rio. Более подробно про него можно прочитать здесь.

install.packages('rio')
library(rio)
convert("USArrests.csv", 
        "USArrests.sav") # на первом месте - 
# что конвертируем, на втором - во что

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

Теперь, когда мы можем открыть практически любую базу, можем приступить к анализу данных. Для начала возьмем базу swiss, содержащую социально-экономические показатели 47 франкоговорящих провинций Швейцарии за 1888 год (встроенная база R).

s <- swiss

Удаление пропущенных значений

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

s <- na.omit(s) # удаляем

Первичный анализ данных

Описательные статистики

Основные описательные статистики

Функция summary() выдает описательные статистики для всех переменных в базе.

summary(s)
##    Fertility      Agriculture     Examination      Education    
##  Min.   :35.00   Min.   : 1.20   Min.   : 3.00   Min.   : 1.00  
##  1st Qu.:64.70   1st Qu.:35.90   1st Qu.:12.00   1st Qu.: 6.00  
##  Median :70.40   Median :54.10   Median :16.00   Median : 8.00  
##  Mean   :70.14   Mean   :50.66   Mean   :16.49   Mean   :10.98  
##  3rd Qu.:78.45   3rd Qu.:67.65   3rd Qu.:22.00   3rd Qu.:12.00  
##  Max.   :92.50   Max.   :89.70   Max.   :37.00   Max.   :53.00  
##     Catholic       Infant.Mortality
##  Min.   :  2.150   Min.   :10.80   
##  1st Qu.:  5.195   1st Qu.:18.15   
##  Median : 15.140   Median :20.00   
##  Mean   : 41.144   Mean   :19.94   
##  3rd Qu.: 93.125   3rd Qu.:21.70   
##  Max.   :100.000   Max.   :26.60

Естественно, можем запросить отдельные статистики для конкретных переменных: максимум, минимум, размах, среднее, медиану, дисперсию и стандартное отклонение.

max(s$Fertility); min(s$Fertility); range(s$Fertility)
## [1] 92.5
## [1] 35
## [1] 35.0 92.5
mean(s$Fertility); median(s$Fertility)
## [1] 70.14255
## [1] 70.4
var(s$Fertility); sd(s$Fertility) # дисперсия и ст. отклонение
## [1] 156.0425
## [1] 12.4917

Кроме того, можно получить набор основных статистик одной командой:

fivenum(s$Fertility) # Tuckey's five numbers: min, 1st quartile, median, 3rd quartile, max
## [1] 35.00 64.70 70.40 78.45 92.50

Как можно заметить, fivenum() выдает те показатели, которые используются для построения “ящика с усами” (разработанном Дж.Тьюки), речь о котором пойдет позже.

Процентили, квартили, децили:

quantile(s$Catholic, 0.65) # 65 процентиль
##   65% 
## 57.54
quantile(s$Catholic, c(0.25, 0.5, 0.75)) # квартили
##    25%    50%    75% 
##  5.195 15.140 93.125
quantile(s$Catholic, seq(0, 1, 0.1)) # децили
##      0%     10%     20%     30%     40%     50%     60%     70%     80% 
##   2.150   2.832   4.610   5.542   9.174  15.140  38.912  90.732  97.094 
##     90%    100% 
##  99.000 100.000

Стандартизация

Для стандартизации переменных в R есть функция scale(). Однако, если помнить смысл стандартизации, всегда можно сделать ее вручную.

Fert_std <- scale(s$Fertility) # стандартизация
Fert_std2 <- (s$Fertility-mean(s$Fertility))/sd(s$Fertility) # то же, но вручную

Сравним результаты (head() - выдает несколько первых значений):

head(Fert_std) 
##           [,1]
## [1,] 0.8051305
## [2,] 1.0372847
## [3,] 1.7897846
## [4,] 1.2534283
## [5,] 0.5409551
## [6,] 0.4769125
head(Fert_std2) # видим, что результаты одинаковы
## [1] 0.8051305 1.0372847 1.7897846 1.2534283 0.5409551 0.4769125

Функции

Часто бывает, что встроенных функций R не хватает для полноценной работы, поэтому приходится дописывать некоторые функции самим. Создадим функцию descriptives, которая будет выдавать набор описательных статистик для какой-либо количественной переменной (как мы убедились, в R такие функции есть, но пусть будет для примера).

При создании функции нужно учитывать три момента: 1) что функция получает на вход (аргументы функции); 2) что функция делает; 3) что функция выдает на выходе (результат).

В данном случае будем считать, что функция descriptives получает на вход некоторую количественную переменную – числовой вектор (вектор x). Затем она считает минимальное значение переменной, ее максимальное значение, размах, среднее и стандартное отклонение. На выходе функция выдает вектор с полученными показателями.

# descriptives - имя функции, x - аргумент (их может быть несколько)

descriptives <- function(x){
  
  # в фигурных скобках перечисляем, что должна делать функция
  
  minimum <- min(x)
  maximum <- max(x)
  range <- maximum - minimum
  m <- mean(x)
  s <- sd(x) 
  allstat <- c(minimum, maximum, range, m, s) 
  
  # возвращает вектор с рассчитанными значениями
  return (allstat)
}

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

descriptives(s$Agriculture) 
## [1]  1.20000 89.70000 88.50000 50.65957 22.71122

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

Создадим матрицу, содержащую доходы пяти различных людей за три месяца (5 строк, 3 столбца). Немного странная матрица, но для примера подойдет:

income_1 <- c(25000, 15000, 8000, 12000, 4500)
income_2 <- c(24000, 32000, 17000, 46000, 53000)
income_3 <- c(36000, 80000, 12000, 27000, 6500)
M <- cbind(income_1, income_2, income_3) # получили матрицу из 3 столбцов
M
##      income_1 income_2 income_3
## [1,]    25000    24000    36000
## [2,]    15000    32000    80000
## [3,]     8000    17000    12000
## [4,]    12000    46000    27000
## [5,]     4500    53000     6500

Теперь попробуем применитьapply(): сначала указываем, к чему применяем, потом - что.

apply(M, 1, mean) # считаем среднее для каждой строки (1 - отвечает за строки)
## [1] 28333.33 42333.33 12333.33 28333.33 21333.33

Получили средний доход каждого человека.

apply(M, 2, mean) # считаем среднее для каждого столбца (2 - отвечает за столбцы)
## income_1 income_2 income_3 
##    12900    34400    32300

Получили среднемесячный доход.

apply(M, c(1,2), log) # логарифмируем значения в каждой ячейке (с(1,2) - отвечает за ячейки)
##       income_1  income_2  income_3
## [1,] 10.126631 10.085809 10.491274
## [2,]  9.615805 10.373491 11.289782
## [3,]  8.987197  9.740969  9.392662
## [4,]  9.392662 10.736397 10.203592
## [5,]  8.411833 10.878047  8.779557

Можем осуществлять сразу несколько действий. Например, логарифмировать значения в каждой ячейке и округлять результат до второго знака после запятой:

round(apply(M, c(1,2), log),2) 
##      income_1 income_2 income_3
## [1,]    10.13    10.09    10.49
## [2,]     9.62    10.37    11.29
## [3,]     8.99     9.74     9.39
## [4,]     9.39    10.74    10.20
## [5,]     8.41    10.88     8.78

То же можно применять к базам данных.

Проверка нормальности

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

  1. Нормальная вероятностная бумага
qqnorm(s$Catholic)
qqline(s$Catholic) # добавляем reference line, чтобы было видно, с чем сравнивать

Как можно заметить, наблюдения располагаются почти на S-образной кривой. Видно невооруженным глазом, что распределение показателя “доля католиков в кантоне” не является нормальным. Тем не менее, посмотрим на более строгие критерии.

  1. Критерий Шапиро-Уилка
shapiro.test(s$Catholic)
## 
##  Shapiro-Wilk normality test
## 
## data:  s$Catholic
## W = 0.7463, p-value = 1.205e-07
  1. Критерий Колмогорова-Смирнова
ks.test(s$Catholic, "pnorm")
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  s$Catholic
## D = 0.98422, p-value < 2.2e-16
## alternative hypothesis: two-sided

Как показывают критерии, действительно, есть основания отвергнуть гипотезу о нормальности распределения доли католиков среди населения кантонов Швейцарии.

Стоит обратить внимание на то, что при расчете критерия Колмогорова-Смирнова R выдал предупреждение “ties should not be present for the Kolmogorov-Smirnov test”. Он сообщает о том, что для точного расчета p-value в выборке не должны присутствовать повторяющиеся значения (ties).

Сравнение средних и распределений

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

Подготовим данные. Будем сравнивать процент мужчин, занятых в сельском хозяйстве, в кантонах, где доля католиков выше медианного значения (50-я процентиль) и в кантонах, где эта доля не превышает медианное значение.

q <- quantile(s$Catholic, 0.5) # 50 процентиль (медиана)
Gr1 <- subset(s, s$Catholic > q) # две маленькие базы
Gr2 <- subset(s, s$Catholic <= q)

t-test

Применение:

  1. 2 выборки
  2. Нормальное распределение данных
  3. Есть 2 варианта теста: для независимых выборок и для связных выборок

t-test для 2 независимых выборок

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

t.test(Gr1$Agriculture, Gr2$Agriculture) # просто две переменных (два вектора)
## 
##  Welch Two Sample t-test
## 
## data:  Gr1$Agriculture and Gr2$Agriculture
## t = 1.453, df = 43.18, p-value = 0.1534
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.701877 22.796442
## sample estimates:
## mean of x mean of y 
##  55.53478  45.98750

По значению t и по p-value видно, что гипотезу о равенстве средних не следует отвергать.

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

s$group <- ifelse(s$Catholic > q, 1, 0) # группирующая переменная
View(s)
t.test(s$Agriculture ~ s$group) # через ~ указываем, по какой переменной группируем
## 
##  Welch Two Sample t-test
## 
## data:  s$Agriculture by s$group
## t = -1.453, df = 43.18, p-value = 0.1534
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -22.796442   3.701877
## sample estimates:
## mean in group 0 mean in group 1 
##        45.98750        55.53478

По умолчанию проверяется двусторонняя гипотеза (two-sided). Если нужна односторонняя гипотеза, задаем значение аргумента alternative.

t.test(s$Agriculture~s$group, alternative = "less") # левосторонняя альтернатива (H0: mean1 - mean2 > 0)
## 
##  Welch Two Sample t-test
## 
## data:  s$Agriculture by s$group
## t = -1.453, df = 43.18, p-value = 0.07672
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##      -Inf 1.497246
## sample estimates:
## mean in group 0 mean in group 1 
##        45.98750        55.53478
t.test(s$Agriculture~s$group, alternative = "greater") # правосторонняя альтернатива (H0: mean1 - mean2 < 0)
## 
##  Welch Two Sample t-test
## 
## data:  s$Agriculture by s$group
## t = -1.453, df = 43.18, p-value = 0.9233
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  -20.59181       Inf
## sample estimates:
## mean in group 0 mean in group 1 
##        45.98750        55.53478

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

t-test для 2 связных выборок

Представим теперь, что у нас появились данные за 1889 год (напомню, в базе s - данные за 1888 год). Измерения в конце каждого года производились одинаково. Можем считать выборку значений Agriculture в 1888 и выборку значений Agriculture в 1889 связными выборками. Сгенерируем нашу вымышленную переменную Agriculture89 - возьмем ее значения из нормального распределения с таким же средним Agriculture, как и в 1888 году:

s$Agriculture89 <- rnorm(n = length(s$Agriculture), mean = mean(s$Agriculture)) # о функции rnorm() мы поговорим более подробно, но позже
# paired = TRUE - отвечает за то, что выборки являются связными

t.test(s$Agriculture, s$Agriculture89, paired = TRUE)
## 
##  Paired t-test
## 
## data:  s$Agriculture and s$Agriculture89
## t = 0.082077, df = 46, p-value = 0.9349
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -6.423852  6.969992
## sample estimates:
## mean of the differences 
##               0.2730697

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

Критерий Уилкоксона (Манна-Уитни)

Применение:

  1. 2 выборки
  2. Нет нормального распределения данных
  3. Выборки независимы

Нулевая гипотеза: выборки взяты из одного и того же распределения (иногда это представляется как равенство медиан двух распределений).

wilcox.test(s$Agriculture ~ s$group)
## 
##  Wilcoxon rank sum test
## 
## data:  s$Agriculture by s$group
## W = 202, p-value = 0.1185
## alternative hypothesis: true location shift is not equal to 0

ANOVA

Применение:

  1. Больше 2 выборок
  2. Нормальное распределение данных

Воспользуемся встроенной в R базой с данными по эффективности спреев от насекомых InsectSprays - отойдем от политологии:)

example <- InsectSprays
View(example)

Будем смотреть на разницу в эффективности 6 спреев (A–F). Описание базы данных можно найти, набрав help(InsectSprays). В качестве показателя эффективности возьмем переменную count, которая, видимо, означает число насекомых, исчезающих при использовании спрея той или иной марки.

Нулевая гипотеза: средние по группам равны.

estimates <- aov(example$count ~ example$spray) # ANOVA, выдает сумму квадратов
summary(estimates) # все показатели
##               Df Sum Sq Mean Sq F value Pr(>F)    
## example$spray  5   2669   533.8    34.7 <2e-16 ***
## Residuals     66   1015    15.4                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

Критерий Краскала-Уоллиса

Применение:

  1. Больше 2 выборок
  2. Нет нормального распределения данных
  3. Выборки независимы

Нулевая гипотеза: выборки взяты из одного и того же распределения (иногда это представляется как равенство медиан распределений).

kruskal.test(example$count ~ example$spray)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  example$count by example$spray
## Kruskal-Wallis chi-squared = 54.691, df = 5, p-value = 1.511e-10

И опять нулевую гипотезу следует отвергнуть.