Что такое R?

Начнем с краткого введения в R.

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

Чтобы скачать R, перейдите по ссылке на сайт CRAN и выберите версию R, которая совместима с вашим программным обеспечением: Linux, Mac или Windows.

Мы выбираем версию для Windows:

После скачивания файла с программой запустите его и нажимайте кнопку “Next” во всех открывающихся окнах интерфейса, а затем кнопку “Finish” в конце. Теперь все готово к работе!

Откройте R и наберите в командной строке код print("Hello!"). Затем нажмите кнопку “Enter”, как показано ниже:

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

Чтобы создать новый скрипт, нажмите “File” и выберите в выпадающем меню “New script”. Теперь в открывшемся окне выведите фразу “Hello!”. Для этого наберите команду print("Hello!"), выделите ее, кликните правой кнопкой мыши и нажмите “Run line or selection” или используйте сочетание клавиш “Control + R” на ПК или “Command + Enter” на Mac. Результат отобразится в окне консоли красным цветом.

С языком программирования R удобнее работать в программе RStudio. Ее можно скачать с сайта RStudio. Скачайте бесплатную версию программы - ее вполне достаточно. Наш обучающий модуль построен исключительно с использованием бесплатной версии программы RStudio. Выберите версию программы, которая совместима с вашей операционной системой. Мы выбираем версию для Windows.

Ниже показано, как открыть новый скрипт в RStudio и выполнить простой код. Результат выполнения кода отображается в консоли (внизу слева).

Прогнать код можно также выделив его мышкой и нажав “Control + Enter” на ПК или “Command + Enter” на Mac.

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

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

Чтобы загрузить файл с данными, нажмите на кнопку “Import dataset” в верхнем правом углу окна и выберите тип файла, который требуется загрузить. Откроется окно, в котором можно выбрать папку и файл с данными:

Удобнее, однако, настроить в RStudio рабочую папку и загружать файлы с данными из нее. Этот вопрос мы рассмотрим подробнее в Разделе 2.

Теперь же давайте поработаем с кодом. Чтобы попрактиковаться, просто копируйте код в командное окно RStudio и нажимайте ""Run" или "Ctrl+Enter"!

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

Помимо тестов, мы рекомендует самостоятельно прогнать код, представленный в разделах модуля, в RStudio на вашем компьютере. Для этого вы можете просто скопировать код в командное окно RStudio. При этом не забудьте указать правильный абсолютный или относительный путь к папкам на вашем компьютере. Чтобы не запутаться, мы рекомендуем прогонять весь код одного раздела в ходе одной сессии работы в RStudio. Если на изучение одного раздела вам понадобятся две отдельные сессии, закрывая программу RStudio обязательно сохраняйте скрипт кода данного раздела и рабочую среду (workspace). В каждом разделе мы используем новые наборы данных. Чтобы пройти весь модуль, может потребоваться от 2 недель до 1 месяца, в зависимости от того, сколько часов в день вы готовы тратить на обучение, а также от того, насколько вы уже знакомы с языком R или с программированием.

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

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

1.1 Создание переменных и арифметические операции в R

Начнем мы с простых арифметических операций: сложение, вычитание, умножение и деление.

Сначала создадим переменные “a” и “j” со значениями 34 и 88 соответственно. Для этого введите следующие команды в R:

a <- 34
j <- 88

(Напомним: чтобы выполнить код, можно либо набрать код в консоли и нажать "Enter", либо набрать код в окне скрипта, выделить его и нажать сочетание клавиш “Ctrl + Enter” на ПК или “Command + Enter” на Mac.)

Обратите внимание, что мы используем функцию присвоения <-, чтобы задать значение 34 для переменной "а".

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

a <- 34
a = 34 

Мы можем выполнять арифметические действия с нашими переменными. Например, ниже представлен код для выполнения операций сложения и деления:

a + j

## [1] 122

a/j

## [1] 0.3863636

Обратите внимание, что результат операции отображается под строчкой с кодом после ## [1]. После первых двух строчек кода не отображается никакого результата, так как в данном случае код просто создает две переменные. Далее, когда мы будем оценивать статистические модели или проверять наборы данных без создания новых объектов R, в консоли будет отображаться результат выполнения команд. Это станет понятнее по мере нашего продвижения через модуль.

Внимание! Если вы попробуете выполнить приведенный ниже код в R, программа перейдет в режим ожидания кода, так как предполагает, что после знака плюс должно быть что-то еще. Чтобы выйти из этого режима, нужно нажать кнопку “Escape”/“Esc”.

a + 

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

a

## [1] 34

Мы можем рассчитать значение логарифма при помощи функции log:

log(a)

## [1] 3.526361

По умолчанию функция "log" вычисляет натуральный логарифм.

#Ответ А неверный. Такой код только вычисляет значение натурального логарифма числа 1000. Но, согласно заданию, нам необходимо создать новую переменную "export", равную значению натурального логарифма числа 1000. Соответственно, правильным ответом является ответ С с кодом "export <- log(1000)." 

#Ответ B неверный. Новая переменная указана рядом с кодом "log(1000)", но ей не присваивается значение, равное log(1000). Чтобы присвоить переменной "export" значение, равное log(1000), необходимо добавить код "<-". Соответственно, правильным ответом является ответ С с кодом "export <- log(1000)." 

#Ответ C верный. Молодец! 

1.2 Типы данных в R

Используемые в R данные могут быть различных типов. Они могут быть числовыми, например: "1011", и текстовыми, например: "abc". Тип данных в R называется классом. Класс данных можно проверить с помощью функции class. Различные классы данных в R имеют различные свойства.

Переменная a является числом, и поэтому относится к классу числовых данных:

a <- 34
class(a)

## [1] "numeric"

Добавив к имени функции “?”, можно просмотреть справочную информацию о ней и о значениях ее аргументов:

?Classes_Details

Вы можете познакомиться с различными классами данных в R, выполнив код ниже и затем прокрутив страницу Справки до раздела “Basic classes”:

Справка отображается в нижней правой части окна RStudio.

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

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

j <- "China, PRC"
class(j)

## [1] "character"

Теперь, если мы попытаемся сложить две переменные a + j, в консоли отобразится сообщение об ошибке:

a + j

Error in a + j : non-numeric argument to binary operator

Сообщения об ошибках могут быть очень полезными, так как они сообщают, почему не выполняется код, что позволяет понять, как его исправить. В примере выше в сообщении указано "non-numeric argument to binary operator". Это значит, что арифметическая операция сложения невозможна с "нечисловым объектом" (в нашем случае j <- "China, PRC"). Если текст ошибки непонятен, введите его в поисковую систему Google и добавьте ключевое слово "stackoverflow".

1.3 Генерирование распределений

Мы можем сгенерировать набор целых чисел в диапазоне от числа А до числа B, как показано ниже:

1:10

##  [1]  1  2  3  4  5  6  7  8  9 10

Также можно сгенерировать случайную совокупность нормально распределенных чисел с помощью функции rnorm. Мы можем открыть справочную информацию об этой функции с помощью кода ниже:

?rnorm

Обязательно просматривайте справочную информацию о каждой новой функции в R. В данном модуле мы даем краткое описание для каждой новой функции. Но более детальную информацию можно получить только в справочном файле.

Мы можем использовать функцию rnorm, чтобы сгенерировать 100 случайных нормально распределенных чисел со средним значением 1 и стандартным отклонением 5.

y <- rnorm(100,1,5)

Теперь мы можем построить гистограмму с помощью функции hist.

hist(y)

Теперь сгенерируем еще одну переменную “x”, которая является результатом умножения значений переменной “y” на 2 и сложения результатов со 100 случайными числами со средним значением 0 и стандартным отклонением 1.

x <- y * 2 + rnorm(100,0,1)

Теперь мы можем графически отобразить отношение между двумя переменными “x” и “y”, создав график при помощи функции plot. Аргумент main позволяет добавить название графика, а с помощью аргумента col можно настроить цвет точек на графике. 2 - код для красного цвета.

plot(x,y, main="My first R plot", col=2)

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

plot(x,y, main="My first R plot", col=1:10)

Ниже приведены примеры создания графиков с точками различных цветов. Код “#000000” соответствует черному цвету.

plot(x,y, col="red")

plot(x,y, col="#000000") #черный

`"#черный" является комментарием. Использование комментариев позволяет добавлять в скрипт кода пояснения или вопросы. При выполнении кода комментарии игнорируются. Чтобы добавить комментарий, перед текстом комментария необходимо ввести “#”.

# Ответ верный. Супер!

#Ответ B неверный. Код не указывает, чему должны быть равны переменные x и y. Код для создания графика написан верно. Соответственно, правильный ответ А.

#Ответ C неверный, так как в коде для построения графика не указаны аргументы для названия графика, а в аргументе для настройки цвета точек графика указан черный цвет. Поэтому, правильный ответ А.

2.1 Загрузка данных

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

Мы можем загрузить файл .csv с интернет-сайта с помощью кода ниже:

trade_data <- read.csv("https://r.tiid.org/R_tute_data2007-2017.csv")

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

Вы также можете скопировать путь к папке из “File explorer”. Для пользователей ПК: при копировании пути к рабочей папке в окно RStudio необходимо заменить \ на /, иначе при попытке выполнить код в консоли отобразится сообщение об ошибке.

setwd("C:/path/to_your/working_directory")

Чтобы проверить путь к рабочей папке, используйте функцию getwd.

getwd()

## [1] "C:/path/to_your/working_directory"

После настройки рабочей папки вы можете загрузить файл .csv в R с помощью функции read.csv, указав в аргументе функции только название файла .csv. (Вы можете выполнить код ?read.csv, чтобы просмотреть все аргументы этой функции.)

trade_data <- read.csv("R_tute_data2007-2017.csv")

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

Error in file(file, “rt”) : cannot open the connection
In addition: Warning message:
In file(file, “rt”) :
cannot open file ‘C:/path/to_your/working_directory/r_tute_data2007-2017.csv’: No such file or directory

Файл trade_data, который мы загрузили в R, является электронной таблицей. В R такой объект называется “data.frame” или "таблица данных".

class(trade_data)

## [1] "data.frame"

С помощью функции names можно проверить названия переменных величин (колонок) в загруженном наборе данных.

names(trade_data)

## [1] "reporter" "flow"     "partner"  "year"     "value"

Просмотреть таблицу данных можно также кликнув на ее имя. Как было упомянуто во Введении, все объекты с переменными и таблицами данных отображаются в верхнем правом углу окна RStudio.

С помощью функции names мы увидели, что таблица данных trade_data содержит столбец reporter, в котором указаны страны, сообщившие об объемах своих торговых потоков, стоимостные значения которых указаны в столбце value. В столбце partner указаны страны-торговые партнеры. В столбце flow указано, является ли торговый поток импортом или экспортом товаров по отношению к стране в столбце "reporter". В столбце year указан год, в котором было зарегистрировано значение торгового потока, указанное в столбце value.

Теперь займемся созданием выборок данных!

2.2 Основы создания выборок данных

Мы можем просмотреть первые несколько строк таблицы данных с помощью функции head:

head(trade_data)

##   reporter flow partner year      value
## 1      213    I     186 2008  156130202
## 2      213    E     174 2008  117661208
## 3      213    E     172 2008   31986252
## 4      213    E     134 2008 1507966544
## 5      213    I     181 2008    2407260
## 6      213    I     182 2008   80414681

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

Если нужно просмотреть конкретное количество строк таблицы, можно задать это число в соответствующем аргументе функции, как показано ниже:

head(trade_data,10)

##    reporter flow partner year       value
## 1       213    I     186 2008   156130202
## 2       213    E     174 2008   117661208
## 3       213    E     172 2008    31986252
## 4       213    E     134 2008  1507966544
## 5       213    I     181 2008     2407260
## 6       213    I     182 2008    80414681
## 7       213    I     676 2008      991884
## 8       213    E     258 2008   107246580
## 9       614    I     110 2008 15302709034
## 10      213    I     311 2008

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

first10 <- head(trade_data,10)

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

first10[1,]

##   reporter flow partner year     value
## 1      213    I     186 2008 156130202

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

first10(1,)

Error in first10(1, ) : could not find function “first10”.

Круглые скобки () используются для указания параметров функций, а квадратные скобки [] используются для указания параметров таблиц данных.

Просмотреть первые 5 строк таблицы можно с помощью обоих вариантов кода, приведенных ниже: Во втором варианте кода 1:5 указывает на строчки с 1 по 5 включительно.

head(first10, 5)

##   reporter flow partner year      value
## 1      213    I     186 2008  156130202
## 2      213    E     174 2008  117661208
## 3      213    E     172 2008   31986252
## 4      213    E     134 2008 1507966544
## 5      213    I     181 2008    2407260

first10[1:5,]

##   reporter flow partner year      value
## 1      213    I     186 2008  156130202
## 2      213    E     174 2008  117661208
## 3      213    E     172 2008   31986252
## 4      213    E     134 2008 1507966544
## 5      213    I     181 2008    2407260

Конкретные строки можно указать с помощью функции c(), создающей вектор значений. Например, код ниже позволит просмотреть только первую и пятую строки:

first10[c(1,5),]

##   reporter flow partner year     value
## 1      213    I     186 2008 156130202
## 5      213    I     181 2008   2407260

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

Например, код ниже выдаст содержимое столбца 1:

first10[,1]

##  [1] 213 213 213 213 213 213 213 213 614 213

Указав соответствующие индексы, мы можем просмотреть содержимое первого и третьего столбцов (столбцы "reporter" и "partner"):

first10[,c(1,3)]

##    reporter partner
## 1       213     186
## 2       213     174
## 3       213     172
## 4       213     134
## 5       213     181
## 6       213     182
## 7       213     676
## 8       213     258
## 9       614     110
## 10      213     311

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

first10[,c("reporter","partner")]

##    reporter partner
## 1       213     186
## 2       213     174
## 3       213     172
## 4       213     134
## 5       213     181
## 6       213     182
## 7       213     676
## 8       213     258
## 9       614     110
## 10      213     311

Обратиться к определенному столбцу таблицы данных можно также указав название соответствующего объекта с данными и знак $ перед названием нужного столбца. В примере ниже мы указали название нашего объекта first10 и название столбца “flow”:

first10$flow

##  [1] I E E E I I I E I I
## Levels: E I

Ниже приведены все три способа обратиться к столбцу "reporter":

first10[,1]

##  [1] 213 213 213 213 213 213 213 213 614 213

first10[,"reporter"]

##  [1] 213 213 213 213 213 213 213 213 614 213

first10$reporter

##  [1] 213 213 213 213 213 213 213 213 614 213

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

Если в названии столбца/переменной есть пробел, его обязательно необходимо указывать в кавычках:

first10$'reporter code' 

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

first10$'flow'

##  [1] I E E E I I I E I I
## Levels: E I

first10$"flow"

##  [1] I E E E I I I E I I
## Levels: E I

С помощью кода можно проверить класс данных в столбце:

class(first10$flow)

## [1] "factor"

В этом случае переменная “flow” относится к классу факторов. Фактором в R является переменная с ограниченным количеством возможных значений, называемых в R "уровнями" ("levels"). Понятие фактора в терминологии R используется для обозначения категориальной (качественной) переменной.

Теперь используем код для одновременного выбора строк и столбцов. Мы можем выбрать только пять первых стран из столбца "reporter" (страны, которые подают данные по статистике своей внешней торговли со странами-партнерами, указанными в таблице в столбце "partner"):

first10[1:5,"reporter"]

## [1] 213 213 213 213 213

Код ниже выводит на экран первые пять стран-reporters и стран-partners.

first10[1:5,c(1,3)]

##   reporter partner
## 1      213     186
## 2      213     174
## 3      213     172
## 4      213     134
## 5      213     181

first10[1:5,c("reporter", "partner")]

##   reporter partner
## 1      213     186
## 2      213     174
## 3      213     172
## 4      213     134
## 5      213     181

#Ответ A неверный, так как он обращается к ячейке во второй строке и пятом столбце: `mydata[5,2]`, в то время как нам нужна ячейка в пятой строке и втором столбце `mydata[2,5]`. Правильный ответ B. 

#Ответ B верный! 

#Ответ C неверный, так как в аргументе функции "head" не указано, что нужно выбрать первые 20 строк объекта. Если в аргументе к функции "head" не указано, что нужно выбрать 20 строк, по умолчанию отображается только шесть строк.  Правильный ответ B. 

2.3 Логические операторы

Логические операторы используются в программировании в логических операциях, то есть выражениях, которые проверяют данные на соответствие определенным критериям и которые возвращают значения TRUE/FALSE (правда/ложь). Эти выражения и операторы применяются в том числе для выборки данных из таблиц с использованием более сложных критериев, которые мы рассмотрим в следующем разделе.

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

a == 34

## [1] TRUE

Если выполнить код ниже, то отобразится результат FALSE.

a == 35 

## [1] FALSE

Оператор != означает "не равно".

a != 34

## [1] FALSE

Если мы используем этот оператор с числом 34, то получим результат FALSE, так как переменная a равна 34. Если мы используем этот оператор с числом 35, то получим результат TRUE, так как наша переменная не равна 35.

a != 35

## [1] TRUE

Другие логические операторы включают: >: больше, чем, <: меньше, чем, >=: больше, либо равно и <=: меньше, либо равно.

a < 34

## [1] FALSE

a <= 34

## [1] TRUE

Логические операторы можно использовать с векторами и таблицами данных.

one2five <- 1:5
three2seven <- 3:7

Теперь можно проверить, содержат ли новые векторы число 1.

one2five == 1

## [1]  TRUE FALSE FALSE FALSE FALSE

three2seven == 1

## [1] FALSE FALSE FALSE FALSE FALSE

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

one2five==three2seven

## [1] FALSE FALSE FALSE FALSE FALSE

К счастью, есть еще один способ проверить содержатся ли какие-либо элементы вектора one2five среди элементов вектора three2seven. Для этого можно воспользоваться оператором %in%.

one2five %in% three2seven

## [1] FALSE FALSE  TRUE  TRUE  TRUE

Логические операторы оказываются очень полезными при создании выборок из таблиц данных. Например, давайте выберем только те строки объекта first10, в столбце "reporter" которых содержится значение "614". Мы указываем название нашего объекта и затем определяем параметры выборки в квадратных скобках: first10[]. Сначала мы указываем критерий выборки для строк и ставим запятую. Если мы не указываем никаких критериев для столбцов - мы выбираем все столбцы таблицы данных. Если мы хотим обратиться к данным из конкретных столбцов, мы указываем критерии для них после запятой.

first10[first10$reporter==614,]

##   reporter flow partner year       value
## 9      614    I     110 2008 15302709034

Сам по себе код для критерия выборки дает логическую оценку соответствия данных в указанном столбце этому критерию. В данном случае критерий указан для столбца "reporter".

first10$reporter==614

##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE

Код ниже позволяет выбрать все строки таблицы, в которых указаны данные по экспорту (flow==“E”).

first10[first10$flow=="E",]

##   reporter flow partner year      value
## 2      213    E     174 2008  117661208
## 3      213    E     172 2008   31986252
## 4      213    E     134 2008 1507966544
## 8      213    E     258 2008  107246580

Теперь мы можем применить критерий выше, одновременно выбирая только те строки, для которых значение объема экспорта >100,000,000, применив оператор &:

first10[first10$flow=="E"&
          first10$value>100000000,]

##      reporter flow partner year value
## NA         NA <NA>      NA   NA  <NA>
## NA.1       NA <NA>      NA   NA  <NA>
## NA.2       NA <NA>      NA   NA  <NA>
## NA.3       NA <NA>      NA   NA  <NA>

Оператор & значить "и".

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

Код выше сейчас выдаст нам сообщение об ошибке, так как данные в столбце “value” относятся к классу факторов, а не чисел.

class(first10$value)

## [1] "factor"

Мы можем поменять класс данных в столбце с фактора на числа:

first10$value <- as.numeric(as.character(first10$value))

Чтобы конвертировать факторы в числа, нужно сначала конвертировать факторы в текст, а только затем в числа, иначе код обратится к набору вариантов факторов (уровням), а не к самим числам. Конвертация в текст осуществляется с помощью функции as.character. Для конвертации в числа мы используем функцию as.numeric.

Теперь в столбце числовые данные:

class(first10$value)

## [1] "numeric"

И теперь мы можем использовать логический оператор > для создания выборки данных:

first10[first10$flow=="E"&
          first10$value>100000000,]

##   reporter flow partner year      value
## 2      213    E     174 2008  117661208
## 4      213    E     134 2008 1507966544
## 8      213    E     258 2008  107246580

Можно создать выборку из тех строк, в которых значение объёма экспорта партнеру под номером 174 выше 100 миллионов:

first10[first10$flow=="E"&
          first10$value>100000000&
          first10$partner==174,]

##   reporter flow partner year     value
## 2      213    E     174 2008 117661208

Ниже мы создаем такую же выборку, как выше, но теперь включаем партнеров по экспорту под номерами 174 и 134:

first10[first10$flow=="E"&
          first10$value>100000000&
          first10$partner %in% c(174,134),]

##   reporter flow partner year      value
## 2      213    E     174 2008  117661208
## 4      213    E     134 2008 1507966544

Выберем все строки для экспорта и импорта и для стран-партнеров 174 и 34 с пользованием оператора “|”, который значит “или”:

first10[first10$partner==174|first10$partner==134,]

##   reporter flow partner year      value
## 2      213    E     174 2008  117661208
## 4      213    E     134 2008 1507966544

Используя оператор “или” (|), мы выбираем все ряды соответствующие, одному критерию, а также все ряды, соответствующие другом критерию.

Два варианта кода, приведенные ниже, эквивалентны:

first10[first10$partner %in% c(174,134),]

##   reporter flow partner year      value
## 2      213    E     174 2008  117661208
## 4      213    E     134 2008 1507966544

first10[first10$partner==174|first10$partner==134,]

##   reporter flow partner year      value
## 2      213    E     174 2008  117661208
## 4      213    E     134 2008 1507966544

#Ответ А верный. Молодец! 

#Ответ B неверный. В нем правильно указаны критерии выбора партнеров (коды 186 и 181), но не указан критерий выбора только тех строк, которые содержат данные импорта ("I") и в которых стоимостный объем импорта выше 1000. А значит, ответ А верный. 

rm(first10)

3.1 Выборка данных (продолжение)

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

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

trade_data <- read.csv("R_tute_data2007-2017.csv")

Убедитесь, что вы работаете в той же рабочей папке, что и в предыдущем разделе! Проверить текущую рабочую папку можно с помощью функции getwd.

Теперь проверим количество строк в таблице данных с помощью функции nrow.

nrow(trade_data)

## [1] 1244330

В нашей таблице 1244330 строк!

Важное примечание! Когда вы работаете с файлом .csv, в котором такое большое количество данных, не открывайте его в программе Excel, так как Excel по умолчанию ограничен 1 048 576 строками, а значит, открыв файл в Excel, часть строк вы потеряете. В этом случае файл, загруженный в R, будет содержать только 1 048 576 строк.

И снова просмотрим содержащиеся в этом наборе данных переменные:

names(trade_data)

## [1] "reporter" "flow"     "partner"  "year"     "value"

Мы можем просмотреть значения только одной переменной, содержащиеся в первых нескольких строках:

head(trade_data$year)

## [1] 2008 2008 2008 2008 2008 2008

Использовав функцию unique, мы также можем проверить, данные для каких уникальных лет содержатся в нашей таблице данных:

unique(trade_data$year)

##  [1] 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017

В столбце "partner" кодом “1” обозначен мир, а кодом “924” - Китай.

Теперь давайте найдем строку, в которой отображен весь экспорт Китая во все страны мира в 2016 году:

trade_data[trade_data$flow=="E"&    #выбираем только экспорт
trade_data$reporter==924&        #выбираем только Китай в столбце "reporter"
trade_data$year==2016&           #выбираем только 2016 год
trade_data$partner==1, ]         #выбираем только "World" (мир) в столбце ""partner"

##         reporter flow partner year       value
## 1019644      924    E       1 2016 2.13659E+12

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

options(scipen = 999)

С помощью кода ниже мы можем конвертировать данные в столбце "value" в числовой формат:

trade_data$value <- as.numeric(as.character(trade_data$value))

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

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

trade_data <- trade_data[!is.na(trade_data$value),]

Обратите внимание, что сначала мы указываем название объекта с таблицей данных, так как мы хотим ее изменить. Затем мы используем функцию is.na по отношению к значениям интересующей нас переменной. Функция is.na проверяют столбец на наличие отсутствующих значений. Перед is.na мы включаем оператор !, который исключает строки с отсутствующими значениями переменной "value".

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

trade_data[trade_data$flow=="E"&
trade_data$reporter==924&
trade_data$year==2016&
trade_data$partner==1,"value"
]/1000000000

## [1] 2136.59

Мы также можем суммировать весь экспорт Китая во все страны (а не миру) и проверить, совпадет ли цифра с полученным выше значением (2136,59). Для этого нужно выбрать все строки, где партнер Китая по экспорту не мир. В критериях выборки мы можем указать, что в столбце "partner" должен быть указан не мир, использовав оператор !=.

head(trade_data[trade_data$flow=="E"&
trade_data$reporter==924&
trade_data$year==2016&
trade_data$partner!=1,
])

##         reporter flow partner year        value
## 1016163      924    E     293 2016   6081495000
## 1016164      924    E     672 2016   1236859000
## 1016166      924    E     744 2016   1320129000
## 1016184      924    E     399 2016    217025000
## 1016193      924    E     156 2016  27854748000
## 1016255      924    E     405 2016 102783000000

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

CN_X_2016<-trade_data[trade_data$flow=="E"&
trade_data$reporter==924&
trade_data$year==2016&
trade_data$partner!=1,
"value" # добавим название столбца "value" , так как нас интересуют величины только этой переменной
]
sum(CN_X_2016)/1000000000

## [1] 8012.502

Полученное значение оказалось выше, чем то, что у нас получилось ранее (2136,59). Это связано с тем, что эта цифра содержит данные об объемах экспорта Китая в различные региональные экономические группировки стран, а значит некоторые страны-партнеры по экспорту были учтены два или более раз. Чтобы узнать, какой код используется для отдельных стран, а какой для региональных экономических группировок, необходимо использовать словарь данных МВФ. Скачайте его отсюда.

codes <- read.csv("imf_codes.csv")
head(codes)

##   code iso2                  country isCountry                 region
## 1  110 XR29 Advanced Economies (IMF)         0
## 2  512   AF              Afghanistan         1           Asia Pacific
## 3  605   F1                   Africa         0
## 4  799  F19     Africa not specified         1 Africa and Middle East
## 5  914   AL                  Albania         1                 Europe
## 6  612   DZ                  Algeria         1 Africa and Middle East
##          income_level                 AP_region LLDC LDC SIDS
## 1                                                 NA  NA   NA
## 2          low_income South and South-West Asia    1   1    0
## 3                                                 NA  NA   NA
## 4                                                  0   0    0
## 5 upper_middle_income                              0   0    0
## 6 upper_middle_income                              0   0    0

Из этой таблицы данных нам нужны столбцы: code, country, isCountry, region.

Перезапишем объект с таблицей данных "codes", сохранив только эти столбцы.

codes <- codes[, c("country","code", "isCountry", "region" )]

Мы можем изменить порядок столбцов:

codes <- codes[, c("code","country", "isCountry", "region" )]

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

Как это сделать? В объекте "codes" есть столбец isCountry. Если в его ячейке указано значение 1, то относящийся к ней трехзначный код описывает отдельную страну. В объекте "trade_data" мы также хотим сохранить данные об объемах торговли с миром. Как установить оба этих критерия?

1. Сначала мы используем код codes$isCountry==1, где isCountry==1 является фильтром, который выбирает только строки, содержащие значение 1 в столбце "isCountry".
   
2. Затем мы используем оператор "или" - ( | ), вводящий критерий, согласно которому мы также выбираем строки, содержащие значение 1 в столбце code.
   
3. После запятой мы указываем название столбца code, так как только он нам нужен из таблицы данных "codes". (Если мы не укажем название столбца после запятой, то мы получим данные всех столбцов, и тех строк, которые соответствуют нашим критериям. Можете сами проверить! codes[codes$isCountry==1|codes$code==1,])

Используем этот код, чтобы создать новый объект из объекта codes:

cntry_w <- codes[codes$isCountry==1|codes$code==1,"code"]

Таким образом, объект cntry_w содержит список кодов МВФ, которые соответствуют отдельным странам (а не группировкам стран), а также код 1 для мира.

Мы можем использовать новый объект , чтобы отфильтровать данные объекта CN_X_2016 так, чтобы он содержал только данные об объемах экспорта Китая в отдельные страны и в мир в 2016 году. Для этого мы используем код trade_data$partner %in% cntry_w, указывая на то, что мы хотим сохранить данные только для тех партнеров по экспорту, коды для которых содержатся в объекте cntry_w.

CN_X_2016<-trade_data[trade_data$flow=="E"&   #выбираем только экспорт
trade_data$reporter==924& # из Китая
trade_data$year==2016& # в 2016 году
trade_data$partner!=1& #партнерам, не отмеченным кодом 1
trade_data$partner %in% cntry_w, #но отмеченным другими кодами из объекта "cntry_w"
"value" #сохраняем только данные столбца "value"
]
sum(CN_X_2016)/1000000000

## [1] 2135.245

Теперь мы получили значение 2135,25, которое гораздо ближе к значению 2136,59, чем 8012,50. Небольшое расхождение связано с тем, что некоторые регионы мира не отмечены в качестве стран.

#Ответ A неверный, так как он не содержит критерия выбора только тех строк, которые содержат данные по экспорту (`trade_data$flow="E"&`). Такой критерий установлен в ответе D, а значит он верный. 

#Ответ B неверный, так как код содержит "=" вместо "==". Логический оператор равенства "==" используется для установки критерия проверки значений переменных, а функция "=" используется для присвоения значений переменным (например, код "trade = 100" присваивает значение 100 переменной "trade", а код "trade == 1000" проверяет, равно ли значение уже существующей переменной "trade" 1000.) Таким образом, ответ D верный. 

#Ответ C неверный. Переменные, содержащиеся в столбцах таблиц данных, не являются объектами в R, а значит R не сможет их найти при обращение к переменным только по их имени. Код должен содержать название объекта "trade_data", содержащего таблицу данных с этой переменной: "trade_data$flow". Ответ D верный.

#Ответ D верный. Молодец! 

3.2 Создание графиков

Попробуем построить график значений экспорта Китая в разные годы.

Сначала создадим таблицу данных, содержащую значения валового экспорта Китая в мир (1) во все годы, для которых эти данные доступны:

CN_X <- trade_data[trade_data$flow=="E"&
trade_data$reporter==924&
trade_data$partner==1,
]
head(CN_X)

##        reporter flow partner year         value
## 24180       924    E       1 2008 1428870000000
## 148613      924    E       1 2009 1202050000000
## 273046      924    E       1 2010 1578440000000
## 397479      924    E       1 2011 1899280000000
## 521912      924    E       1 2012 2050090000000
## 646345      924    E       1 2013 2210660000000

Мы можем воспользоваться функцией plot для построения графика, показывающего экспорт Китая в мир в разные годы:

plot(CN_X[,c("year", "value")])

Выразим эти значения в триллионах?

CN_X$value <- CN_X$value/1000000000000

Так график выглядит лучше:

plot(CN_X[,c("year", "value")])

Добавим название графика и подпишем ось y:

plot(CN_X[,c("year", "value")],
main="China's exports",
xlab="",
ylab="USD, trillions")

Обратите внимание, что мы используем аргумент main, чтобы добавить название графика. С помощью аргументов xlab и ylab можно добавить описание осей x и y соответственно.

Точки не очень понятны. Попробуем использовать линию:

plot(CN_X[,c("year", "value")],
main="China's exports",
xlab="",
ylab="USD, trillions",
type="l") # указываем тип линии

Чтобы выбрать линию, используем аргумент type = l (где буква “l”, а не цифра 1, указывает на линию). По умолчанию этот аргумент настроен на p, что указывает на использование точек. (Используйте код ?plot, чтобы посмотреть информацию об этой функции.)

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

CN_M <- trade_data[trade_data$flow=="I"&
trade_data$reporter==924&
trade_data$partner==1,
]
CN_M$value <- CN_M$value/1000000000000


plot(CN_X[,c("year", "value")],
main="China's exports and imports",
xlab="",
ylab="USD, trillions",
type="l",
ylim=c(1,2.5), #добавим ограничения по оси y
col=2, #выберем красный цвет линии
lty=1 # выберем сплошную линию; более подробную информацию см. выполнив код ?par
)
lines(CN_M[,c("year", "value")],col=4, lty=2) #выберем синюю пунктирную линию

legend("bottomleft",
legend = c("Exports", "Imports"),
col=c(2,4),
lty=c(1,2)
)

Сначала мы создали объект с таблицей данных CN_M. Для этого мы применили фильтр к объекту trade_data, выбрав только данные по объемам импорта: trade_data$flow=="I", в Китай: trade_data$reporter==924 из мира: trade_data$partner==1. Мы выразили числовые значения в триллионах: CN_M$value <- CN_M$value/1000000000000.

Мы объединили график экспорта CN_X и график импорта CN_M в один, сначала создав график для первого набора данных, а затем использовав функцию lines() для добавления в этот же график второй линии. Для первого графика на основе объекта CN_X мы выбираем переменные год и объем экспорта: CN_X[,c("year", "value")]. Мы указываем название графика: main="China's exports and imports", и затем подписываем ось: ylab="USD, trillions". Указываем, что хотим создать график в виде линии: type="l". С помощью аргументов xlim и ylim можно настроить координатные оси x и y соответственно. Цвет графика настраиваем с помощью аргумента col, а сплошную линию выбираем с помощью аргумента lty=1. Выполнив код ?par, можно просмотреть все настраиваемые параметры графиков (с помощью кода ?plot также можно просмотреть некоторые настраиваемые параметры, но полный список параметров находится в par).

Чтобы добавить второй график, нужно указать переменные второго объекта с данными: CN_M[,c("year", "value")]. Пусть линия будет синего цвета: col=4 и пунктирной: lty=2.

С помощью функции legend мы добавляем легенду к графику. Указываем, где на нашем графике мы хотим поместить легенду (внизу слева в нашем случае), даем наименования элементам легенды: legend = c("Exports", "Imports"), а также указываем параметры для цвета и типов линии обоих графиков: col=c(2,4) и lty=c(1,2) соответственно. Обратите внимание, что параметры одного типа мы указывает с использованием функции c, если их нужно указать вместе. В Разделе 4 мы рассмотрим вопросы построения графиков более подробно.

Теперь давайте построим график, на котором отметим 10 основных партнеров Китая по экспорту в 2016 году.

Создадим объект с таблицей данных по экспорту Китая в отдельные страны в 2016 году:

CN_X_countries <-
trade_data[trade_data$flow=="E"&
trade_data$reporter==924&
trade_data$year==2016&
trade_data$partner!=1,]

head(CN_X_countries)

##         reporter flow partner year        value
## 1016163      924    E     293 2016   6081495000
## 1016164      924    E     672 2016   1236859000
## 1016166      924    E     744 2016   1320129000
## 1016184      924    E     399 2016    217025000
## 1016193      924    E     156 2016  27854748000
## 1016255      924    E     405 2016 102783000000

Мы хотим указать названия стран на графике. Для этого мы используем функцию merge, чтобы объединить объекты CN_x_countries и codes.

CNX_merged <- merge(CN_X_countries,
codes, by.x="partner",
by.y="code")
CNX_merged <- CNX_merged[, c("country", "value", "isCountry","region")]

head(CNX_merged)

##                    country         value isCountry        region
## 1    Export earnings: fuel  167185000000         0
## 2 Export earnings: nonfuel 1969410000000         0
## 3 Advanced Economies (IMF) 1403500000000         0
## 4            United States  389714000000         1 North America
## 5           United Kingdom   56645733000         1        Europe
## 6                  Austria    2245591000         1        Europe

При использовании функции merge сначала необходимо указать объекты с таблицами данных, которые мы хотим объединить. В нашем случае: CN_X_countries и codes. Затем нужно указать названия столбцов, по которым мы объединяем две таблицы данных. Названия этих столбцов разные, но содержат они значения одной и той же категориальной переменной (факторы). С помощью аргумента by.x указываем название столбца в первом объекте, а с помощью аргумента by.y - название столбца во втором объекте. (Если столбцы, по которым производится объединение двух таблиц данных, имеют одно название, то указывать значения аргументов by, by.x или by.y не требуется.)

Затем мы сохраняем только столбцы с переменными country, value, isCountry, region. Столбец isCountry мы используем, чтобы сохранить строки с данными только для отдельных стран, а столбец region сохраняем, чтобы в графике отметить страны одного региона одним цветом.

Сохраняем строки только с данными для отдельных стран:

CNX_merged <- CNX_merged[CNX_merged$isCountry==1,]

Сортируем данные с помощью функции order. Код ниже сортирует строки в порядке возрастания значения в столбце "value ":

CNX_merged <- CNX_merged[order(CNX_merged$value),]
head(CNX_merged)

##                 country   value isCountry                      region
## 65           Montserrat  110000         1 Latin America and Carribean
## 60            Greenland  993000         1               North America
## 11           San Marino 1506000         1                      Europe
## 178               Nauru 1644000         1                Asia Pacific
## 174       Faroe Islands 1737000         1                      Europe
## 68  St. Kitts and Nevis 4850000         1 Latin America and Carribean

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

CNX_merged <- CNX_merged[order(-CNX_merged$value),]
head(CNX_merged)

##                country        value isCountry        region
## 4        United States 389714000000         1 North America
## 101   Hong Kong, China 293997000000         1  Asia Pacific
## 19               Japan 129617000000         1  Asia Pacific
## 105 Korea, Republic of  95815599000         1  Asia Pacific
## 10             Germany  66044117000         1        Europe
## 116           Viet Nam  62040928000         1  Asia Pacific

Создаем новый объект CNX_merged_top10, который содержит выборку из первых 10 стран. Затем строим круговую диаграмму с помощью функции pie.

CNX_merged_top10 <- head(CNX_merged,10)

pie(CNX_merged_top10$value,
CNX_merged_top10$country)

Также можно построить столбчатую диаграмму с помощью функции barplot. Укажем, что страны одного региона должны быть отображены на графике в одном цвете: col=CNX_merged_top10$region.

barplot(CNX_merged_top10$value,
names=CNX_merged_top10$country,
col=CNX_merged_top10$region)

Рассмотрим еще один пример перед тем, как мы перейдем в подраздел о циклах в R. Построим аналогичные графики на основе данных о торговле Филиппин (566) в 2016.

Сначала извлечем данные об объемах экспорта из Филиппин в отдельные страны в 2016 году.

PH_X_countries <-
trade_data[trade_data$flow=="E"&
trade_data$reporter==566&
trade_data$year==2016&
trade_data$partner!=1,]

head(PH_X_countries)

##         reporter flow partner year       value
## 1081095      566    E     112 2016   475800328
## 1081100      566    E     158 2016 11674107669
## 1081102      566    E     542 2016  2095041077
## 1081103      566    E     544 2016      684932
## 1081108      566    E     343 2016      429046
## 1081109      566    E     433 2016    14130931

Объединим новый объект c объектом "codes":

PHX_merged <- merge(PH_X_countries,
codes, by.x="partner",
by.y="code")
PHX_merged <- PHX_merged[PHX_merged$isCountry==1, c("country", "value", "isCountry","region")]

head(PHX_merged)

##          country      value isCountry        region
## 4  United States 8670654131         1 North America
## 5 United Kingdom  475800328         1        Europe
## 6        Austria   67115695         1        Europe
## 7        Belgium  399641671         1        Europe
## 8        Denmark   35623543         1        Europe
## 9         France  726782604         1        Europe

Отсортируем полученную таблицу данных в порядке убывания значений переменной "value":

PHX_merged <- PHX_merged[order(-PHX_merged$value),]

head(PHX_merged)

##              country       value isCountry        region
## 18             Japan 11674107669         1  Asia Pacific
## 4      United States  8670654131         1 North America
## 97  Hong Kong, China  6582766881         1  Asia Pacific
## 190            China  6192432905         1  Asia Pacific
## 109        Singapore  3700590527         1  Asia Pacific
## 10           Germany  2293130299         1        Europe

Выберем 10 первых стран-партнеров по экспорту и построим столбчатую диаграмму.

PHX_merged_top10 <- head(PHX_merged,10)

barplot(PHX_merged_top10$value,
names=PHX_merged_top10$country,
col=PHX_merged_top10$region,
main="PHL")

Из объекта "codes", используя код 566, можем извлечь название страны "Philippines":

codes[codes$code==566,"country"]

## [1] Philippines
## 245 Levels: Advanced Economies (IMF) Afghanistan Africa ... Zimbabwe

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

barplot(PHX_merged_top10$value,
names=PHX_merged_top10$country,
col=PHX_merged_top10$region,
main=codes[codes$code==566,"country"],
las=2) # повернем текст перпендикулярно к оси

Аргумент las = 2 указывает, как должны быть расположены подписи к осям. Значение 2 значит “всегда перпендикулярно к оси”. Дополнительную информацию о других аргументах и параметрах можно прочитать, выполнив код ?par.

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

par(mar=c(11,8,1,1))
barplot(PHX_merged_top10$value,
names=PHX_merged_top10$country,
col=PHX_merged_top10$region,
main=codes[codes$code==566,"country"],
las=2)

Теперь мы можем сохранить график в формате .png с помощью функции png.

png(filename="exports.png") #укажем имя файла
par(mar=c(11,8,1,1)) #настроим поля
barplot(PHX_merged_top10$value,
names=PHX_merged_top10$country,
col=PHX_merged_top10$region,
main=codes[codes$code==566,"country"],
las=2) # повернем текст перпендикулярно к оси

dev.off()

## png
##   2

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

#Сначала найдите код вашей страны в объекте "codes". Условно обозначим, что название вашей страны "ABC", а ее трехзначный код "XYZ". Чтобы воспользоваться кодом ниже, замените "ABC" на название вашей страны, а вместо "XYZ" укажите трехзначный код.
ABC_X_countries <- trade_data[trade_data$flow=="E"&
trade_data$reporter==XYZ&
trade_data$year==2017&
trade_data$partner!=1,]

# Объедините объект "trade_data" с объектом "codes", чтобы добавить к данным о стоимостных объемах экспорта соответствующие названия экономик, регионов и переменную, указывающую на то, является ли экономика страной.
ABC_X_merged <- merge(ABC_X_countries,
codes, by.x="partner",
by.y="code")
ABC_X_merged <- ABC_X_merged[ABC_X_merged$isCountry==1, c("country", "value", "isCountry","region")]

# Отсортируйте строки полученной таблицы данных в порядке убывания значений стоимостных объемов экспорта. Не забудьте добавить знак "-" (минус), иначе сортировка будет по умолчанию осуществлена в порядке возрастания значений.
ABC_X_merged <- ABC_X_merged[order(-ABC_X_merged$value),]

# Выберите 10 стран, являющихся основными импортерами продукции вашей страны
ABC_X_merged_top10 <- head(ABC_X_merged,10)

# Постройте столбчатую диаграмму. Не забудьте заменить "ABC" на название вашей страны.
par(mar=c(11,8,1,1))
barplot(ABC_X_merged_top10$value,
names=ABC_X_merged_top10$country,
col=ABC_X_merged_top10$region,
main="ABC", las=2) 

3.3 Циклы со странами

Это цикл:

for (c in 1:3){
print(c*5)
}

## [1] 5
## [1] 10
## [1] 15

Мы указали, что переменная c принимает значения от 1 до 3 включительно. Затем мы последовательно умножаем каждое из этих значений на 5. Цикл выполняет отдельные итерации данного действия для каждого значения, двигаясь от одного к трем. Отдельный цикл для каждого значения переменной "с" называется итерацией. Для c=1, код умножает значение c на 5, печатает результат (в нашем случае "5"), и затем переходит к следующей итерации, где c=2, и совершает с данным значением переменной эти же действия. Цикл продолжается, пока для последней итерации, где c=3, значение "с" не будет умножено на 5 и результат не будет напечатан.

Изучите подробно синтаксис цикла и использование различных скобок. Сразу после этого выполните тест для практики!

#Ответ A неверный, так как цикл печатает все значения переменной "с" без умножения их на 2, как указано в задании теста. Ответ B верный. 

#Ответ B верный! Супер! 

#Ответ C неверный, так как пропущена фигурная скобка в конце цикла! Обратите внимание на синтаксис цикла "for": for (условие) {выполняемый код}. Ответ В правильный.

#Ответ D неверный. Ошибка в синтаксисе цикла, как и в ответе C, но в данном случае условие цикла не заключено в круглые скобки. Ответ В верный, так как условие цикла заключено в круглые скобки. 

3.4 Циклы (продолжение)

Мы можем использовать цикл, чтобы построить графики сразу для всех стран, аналогичные тем, что мы создали выше! Это очень удобно и экономит время и количество кода, который нужно написать! Циклы очень популярны.

Давайте создадим цикл, который построит столбчатые графики сразу для разных стран. Воспользуемся тем же кодом, который мы использовали для Китая и Филиппин. Ниже приведен пример с кодом для Филиппин. Просмотрите комментарии к коду, чтобы вспомнить, как мы создавали графики!

# создадим выборку с данными по экспорту Филиппин в отдельные страны в 2016 году
PH_X_countries <-
trade_data[trade_data$flow=="E"&
trade_data$reporter==566&
trade_data$year==2016&
trade_data$partner!=1,]

# зададим имя файла и сохраним график в формате .png
png(filename="exports.png") #укажем имя файла
par(mar=c(11,8,1,1)) #настроим поля
# укажем аргументы для настройки параметров графика
barplot(PH_X_countries$value,
names=PH_X_countries$partner,
main=codes[codes$code==566,"country"],
las=2) # повернем текст перпендикулярно к оси

dev.off()

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

Разберемся с каждым этапом отдельно. Начнем с первого. Как должен выглядеть код для цикла, который в каждой своей итерации создает новую выборку данных для каждой отдельной страны-экспортера? Рассмотрим пример для Филиппин ниже:

PH_X_countries <-
trade_data[trade_data$flow=="E"&
trade_data$reporter==566&
trade_data$year==2016&
trade_data$partner!=1,]

Строчка кода, содержащая код страны-экспортера (trade_data$reporter==566), дает R указание создать выборку данных для этой страны. Значит именно эту строчку мы должны отредактировать так, чтобы в каждой итерации формировалась новая выборка.

Для этих целей нам необходимо создать список всех уникальных кодов отдельных стран. С помощью кода ниже мы создаем такой список кодов стран на основе объекта codes.

unique_ctry_codes <- codes[codes$isCountry==1,"code"]
head(unique_ctry_codes)

## [1] 512 799 914 612 859 614

Теперь напишем код цикла, который в каждой своей итерации создает разные выборки данных!

Подробно рассмотрим каждую строчку кода ниже. Цикл будет выполнен для первых пяти кодов стран: for(i in 1:5). Итерации цикла будут выполняться для каждого из значений i в заданном диапазоне. Таким образом, в объекте CN_X_countries будут поочередно содержаться пять разных набора данных. Для создания выборки данных в CN_X_countries укажем нужный нам критерий - только экспортные потоки: trade_data[trade_data$flow=="E". Используем возможные значения переменной i для настройки итераций цикла: trade_data$reporter==unique_ctry_codes[i]. Согласно этому коду, объект с выборкой данных CN_X_countries создается для тех стран в столбце "reporter" объекта trade_data, трехзначный код которых соответствует коду i-го элемента объекта unique_ctry_codes. Мы поместили переменную "i", содержащую индекс элемента из списка unique_ctry_codes, в квадратные скобки, чтобы извлечь значение нужного элемента в зависимости от итерации цикла. И наконец, с помощью элементов кода trade_data$year==2016 и trade_data$partner!=1мы указываем, что нас интересуют данные только за 2016 год и все торговые партнеры, кроме мира.

for(i in 1:5){
CN_X_countries <-
trade_data[trade_data$flow=="E"&
trade_data$reporter==unique_ctry_codes[i]&
trade_data$year==2016&
trade_data$partner!=1,]
}

Мы разобрались, как использовать цикл для создания отдельных выборок данных. Теперь давайте посмотрим, как использовать цикл для создания отдельных графиков!

Ниже приведем код для создания файла с графиком для одной страны.

png(filename="exports.png")
par(mar=c(11,8,1,1))
barplot(PH_X_countries$value,
names=PH_X_countries$partner,
main=codes[codes$code==566,"country"],
las=2)

dev.off()

Как нужно изменить код, чтобы он генерировал отдельный файл с графиком для создаваемой в каждой итерации цикла выборки данных?

Необходимо изменить строку filename="exports.png", так как мы хотим сохранять графики для разных стран под разными именами. Также мы хотим создавать разные заголовки графиков в каждой итерации цикла (main=codes[codes$code==566,"country"]). Остальные параметры зависят от конкретной выборки данных, формируемой в каждой отдельной итерации цикла. (Поэтому нам не нужно ничего менять в этих строках кода.)

Начнем со строки, определяющей заголовок графика. Мы можем отредактировать ее аналогично тому, как мы поменяли код trade_data$reporter==566 на trade_data$reporter==unique_ctry_codes[i] выше. В этом случае мы меняем main=codes[codes$code==566,"country"] на main=codes[codes$code==unique_ctry_codes[i],"country"].

Все немного сложнее в случае с именем файла, так как нам нужно не только изменить название страны (codes[codes$code==unique_ctry_codes[i],"country"]), но и добавить расширение файла .png.

Для этого мы можем использовать функцию paste0, чтобы совместить название страны с расширением “.png” и конвертировать полученное имя файла в текст. Соответственно наш код filename="exports.png" нужно поменять на filename=paste0( codes[codes$code==unique_ctry_codes[i],"country"], ".png" ).

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

for(i in 1:5){
CN_X_countries <-
trade_data[trade_data$flow=="E"&
trade_data$reporter==unique_ctry_codes[i]&
trade_data$year==2016&
trade_data$partner!=1,]

png(filename=paste0(
codes[codes$code==unique_ctry_codes[i],"country"], ".png"
))
par(mar=c(11,8,1,1))
barplot(CN_X_countries$value,
names=CN_X_countries$partner,
main=codes[codes$code==unique_ctry_codes[i],"country"],
las=2)

dev.off()
}

Error in plot.window(xlim, ylim, log = log, …) :
need finite ‘xlim’ values
In addition: Warning messages:
1: In min(w.l) : no non-missing arguments to min; returning Inf
2: In max(w.r) : no non-missing arguments to max; returning -Inf
3: In min(x) : no non-missing arguments to min; returning Inf
4: In max(x) : no non-missing arguments to max; returning -Inf

В результате получаем сообщение об ошибке, приведенное выше, что скорее всего говорит о том, что одна из итераций цикла генерирует пустую выборку данных, а значит график построен быть не может. Но не надо отчаиваться! Это можно исправить. Давайте подробно рассмотрим создаваемые циклом выборки данных. Как проверить, является ли выборка данных пустой? Мы можем проверить количество строк в выборке с помощью функции nrow и напечатать результат с помощью функции print:

for(i in 1:5){
CN_X_countries <-
trade_data[trade_data$flow=="E"&
trade_data$reporter==unique_ctry_codes[i]&
trade_data$year==2016&
trade_data$partner!=1,]

print(nrow(CN_X_countries))
}

## [1] 135
## [1] 0
## [1] 130
## [1] 137
## [1] 98

Как видно из результатов, вторая выборка данных содержит ноль строк. Рассмотрев код страны подробнее, можно обнаружить, что он относится к значению “Africa not specified”. Это не страна, и данные для этого кода не содержатся в объекте с таблицей данных trade_data. Если трехзначный код не относится к стране, итерация создает объект CN_X_countries пустым (результат выполнения команды nrow(CN_X_countries) - 0). Для того, чтобы пропустить этот код, в цикл можно включить условный оператор if, который проверит, является ли созданная выборка пустой, с помощью кода: if (nrow(CN_X_countries)==0){...}. Если это условие истинно, код переходит к следующей переменной в последовательности: if (nrow(CN_X_countries)==0){next}. Код next в R прекращает выполнение программы для текущей итерации цикла и переходит к следующей итерации.

Код цикла выше можно отредактировать следующим образом:

for(i in 1:5){
CN_X_countries <-
trade_data[trade_data$flow=="E"&
trade_data$reporter==unique_ctry_codes[i]&
trade_data$year==2016&
trade_data$partner!=1,]

if (nrow(CN_X_countries)==0){next}

print(nrow(CN_X_countries))
}

## [1] 135
## [1] 130
## [1] 137
## [1] 98

Теперь добавим условный оператор if в код, который создает выборки данных и файлы с диаграммами:

for(i in 1:5){
CN_X_countries <-
trade_data[trade_data$flow=="E"&
trade_data$reporter==unique_ctry_codes[i]&
trade_data$year==2016&
trade_data$partner!=1,]

if (nrow(CN_X_countries)==0){next}

png(filename=paste0(
codes[codes$code==unique_ctry_codes[i],"country"], "2.png"
))
par(mar=c(11,8,1,1))
barplot(CN_X_countries$value,
names=CN_X_countries$partner,
main=codes[codes$code==unique_ctry_codes[i],"country"],
las=2)

dev.off()
}

А если мы хотим, чтобы диаграммы включали только 10 самых крупных импортеров продукции рассматриваемых нами стран? Логика программного кода будет такой же, но нам нужно добавить еще два действия! Нам нужно отсортировать страны по значению объема их импорта в убывающем порядке, а также объединить объекты с выборками данных с объектом "codes", чтобы получить названия стран-импортеров (столбец "partner") и регионов, к которым они относятся.

Сделаем это только для первых пяти элементов списка в объекте unique_ctry_codes. Код цикла начинается с указания элементов последовательности и кода для создания выборки данных CN_X_countries в каждой итерации цикла. Затем мы объединяем выборку данных с объектом codes, чтобы извлечь названия стран-импортеров и названия регионов, к которым они относятся. Далее мы сортируем объединенную таблицу данных и выбираем только 10 самых крупных стран-импортеров. И, наконец, создаем столбчатую диаграмму.

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

f <- 5
for(i in 1:f){

CN_X_countries <-
trade_data[trade_data$flow=="E"&
trade_data$reporter==unique_ctry_codes[i]&
trade_data$year==2016&
trade_data$partner!=1,]
if (nrow(CN_X_countries)==0){next}

#объединяем таблицу, содержащую данные об объемах экспорта, с объектом "codes"

CNX_merged <- merge(CN_X_countries,
codes, by.x="partner",
by.y="code")
CNX_merged <- CNX_merged[CNX_merged$isCountry==1, c("country", "value", "isCountry","region")]

# сортируем строки в порядке убывания

CNX_merged <- CNX_merged[order(-CNX_merged$value),]

# выбираем 10 крупнейших стран-импортеров

CNX_merged_top10 <- head(CNX_merged,10)

png(filename=paste0(
codes[codes$code==unique_ctry_codes[i],"country"], ".png"
))
par(mar=c(11,8,1,1))
barplot(CNX_merged_top10$value,
names=CNX_merged_top10$country,
col=CNX_merged_top10$region,
main=codes[codes$code==unique_ctry_codes[i],"country"],
las=2)
dev.off()
}

Снова подробно рассмотрим код выше.

Цикл выполняется для элементов последовательности от 1 до 5 (f <- 5 for(i in 1:f){...}), но в результате мы получаем только 4 диаграммы. Это связано с тем, что один из элементов последовательности не является страной. Для таких элементов формируется пустая выборка данных. Чтобы учесть формирование пустых выборок, мы используем условный оператор if, который переводит цикл к выполнению следующей итерации, если сформированная в текущей итерации выборка оказывается пустой: if (nrow(CN_X_countries)==0){next}.

Далее мы используем уже знакомый нам код. Мы объединяем таблицу, содержащую данные об экспорте обрабатываемой в текущей итерации страны в другие отдельные страны, с объектом "codes", чтобы извлечь названия стран-импортеров: CNX_merged <- merge(CN_X_countries, codes, by.x="partner", by.y="code"). Далее мы выбираем переменные, которые хотим сохранить: CNX_merged <- CNX_merged[CNX_merged$isCountry==1, c("country", "value", "isCountry","region")]. Сортируем данные: CNX_merged <- CNX_merged[order(-CNX_merged$value),]. Выбираем 10 крупнейших импортеров: CNX_merged_top10 <- head(CNX_merged,10).

Если Вы хотите сгенерировать подобные диаграммы для всех стран-экспортеров, используйте первую строчку кода, чтобы задать значение индекса последнего элемента последовательности цикла: f <- nrow(unique_ctry_codes). (Использование отдельной задаваемой переменной индекса для последнего элемента последовательности, позволяет легко его менять.)

Супер! В Разделе 4 мы научимся строить еще более крутые графики!

#Ответ A неверный, так как мы не изменили аргументы функции "barplot". Ответ C такой же, как ответ A, но в нем указаны верные значения аргументов для функции "barplot", а значит - ответ С верный.  

#Ответ B неверный, так как не задано нужное значение для переменной f! В ответе C нужное значение для данной переменной задано, а значит он верный. 

#Ответ C верный. Супер!

4.1 Базовые графики

Для работы с этим разделом скачайте этот файл. Это ссылка на оригинальную базу данных.

Загрузим файл с данными в RStudio. Так как это файл .csv, используем функцию read.csv. Указываем аргумент stringsAsFactors = FALSE, так как мы не хотим, чтобы текстовые данные были интерпретированы как факторы. (Помните, что вы можете выполнить код ?read.csv, чтобы подробно изучить все возможные аргументы этой функции?)

df <- read.csv("sample_data.csv",
              stringsAsFactors = FALSE) #добавим данный аргумент, чтобы текстовые данные не были конвертированы в факторы

Посмотрим наименования столбцов и класс данных в каждом из них.

names(df)

## [1] "iso3"           "year"           "country"        "gdp_per_capita"
## [5] "pop"            "fertility"

class(df$country)

## [1] "character"

head(df)

##   iso3 year   country gdp_per_capita      pop fertility
## 1  ARG 1968 Argentina       6434.954 23261278     3.049
## 2  ARG 1969 Argentina       6954.764 23605987     3.056
## 3  ARG 1970 Argentina       7056.848 23973058     3.073
## 4  ARG 1971 Argentina       7335.759 24366439     3.104
## 5  ARG 1972 Argentina       7329.921 24782949     3.148
## 6  ARG 1973 Argentina       7407.367 25213388     3.203

Как мы видим, этот набор данных содержит такие переменные как ВВП на душу населения (gdp_per_capita), размер населения (pop) и коэффициент рождаемости (fertility) для каждой страны (country) в разные годы (year).

Построим диаграмму рассеяния для зависимой переменной "ВВП на душу населения" и независимой переменной "размер населения". При использовании функции plot(x,y), первой мы указываем переменную для оси x.

plot(df$pop,df$gdp_per_capita)

#С помощью кода ниже найдите последний год, для которого имеются данные в таблице:
max(df$year)

#В результате получаем 2016 год. Сначала сделаем выборку данных из объекта "df", выбрав только те строки, которые содержат данные за 2016 год:
df2016 <- df[df$year==2016,]

#Чтобы построить график, используем следующий код:
plot(df2016$pop,df2016$gdp_per_capita)

4.2 Графики (продолжение)

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

df2016 <- df[df$year==2016,]
plot(log(df2016$pop, base=10),df2016$gdp_per_capita)

Наиболее удобным пакетом функций для построения диаграмм в R является пакет ggplot2:

#install.packages("ggplot2")
library("ggplot2")

Первая линия кода, оформленная как комментарий (#install.packages("ggplot2")), используется для установки пакета функций ggplot2. Если этот пакет не установлен, удалите # и выполните код install.packages("ggplot2")для установки пакетаggplot2. Несколько основных пакетов функций уже установлены в R, но не все.

Вторая строчка кода library("ggplot2") загружает пакет функций в RStudio. Для каждой новой сессии работы в RStudio необходимо заново загружать нужные вам пакеты функций. Если нужный пакет функций не установлен, то при попытке обратится к функции из этого пакета отобразится следующее сообщение об ошибке:

Error in ggplot(XXX) : could not find function “ggplot”

Теперь попробуем построить график с помощью функции ggplot.

В аргументах функции ggplot сначала необходимо указать объект c данными, которые будут использоваться для построения графика. В нашем случае это объект df. Следующий аргумент aes используется для обращения к переменным и для настройки других параметров, связанных с ними. Параметр geom_point используется для построения диаграмм рассеяния. Ниже приведен пример построения графика с переменными "размер населения" и "ВВП на душу населения" с помощью функции ggplot:

ggplot(df,aes(pop,gdp_per_capita)) + # укажем источник данных и нужные переменные
  geom_point() #укажем, что нужно сделать с данными и переменными

То же самое на основе данных только за 2016 год:

ggplot(df[df$year==2016,],aes(pop,gdp_per_capita)) +
  geom_point()

Сравните оба варианта кода и найдите различие. (Ответ: df[df$year==2016,])

Мы можем построить график зависимости между переменными "ВВП на душу населения" и "коэффициент рождаемости" на основе данных за 2016 год:

ggplot(df[df$year==2016,],aes(gdp_per_capita,fertility)) +
geom_point() 

Мы можем добавить график линейной аппроксимации с помощью функции geom_smooth и указать метод построения графика линейной регрессии. В данном случае мы укажем модель линейной регрессии lm, которая соответствует формуле “y = x + c” или отношению y~x:

ggplot(df[df$year==2016,],aes(gdp_per_capita,fertility)) +
  geom_point() +
geom_smooth(method='lm',formula=y~x) 

Используя данные за 2016 год, построим график зависимости двух переменных: ВВП на душу населения, выраженный через логарифм для учета асимметрии распределения, и коэффициент рождаемости:

ggplot(df[df$year==2016,],aes(log(gdp_per_capita),fertility)) +
  geom_point()

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

Но полученные единицы измерения по оси x сложно интерпретировать. Это можно исправить, выразив масштаб оси x через логарифм по основанию 10:

ggplot(df[df$year==2016,],aes(gdp_per_capita,fertility)) +
  geom_point() +
  scale_x_log10()

Мы можем устранить экспоненциальное представление чисел с помощью кода options(scipen = 999):

options(scipen = 999)

ggplot(df[df$year==2016,],aes(gdp_per_capita,fertility)) +
  geom_point() +
  scale_x_log10()

(Если ранее в ходе текущей сессии работы в RStudio вы уже использовали код options(scipen = 999), после выполнения этого кода сейчас вы не увидите никаких изменений, так как R уже получил соответствующую инструкцию.)

Мы можем отрегулировать размер точек на графике так, чтобы они соответствовали значениям переменной "размер населения". Для этого используем аргумент size=pop:

ggplot(df[df$year==2016,],
       aes(gdp_per_capita,fertility,size=pop)) +
  geom_point() +
  scale_x_log10()

Мы можем настроить разный цвет точек для разных стран с помощью аргумента color=iso3:

ggplot(df[df$year==2016,],
       aes(gdp_per_capita,fertility,size=pop,colour=iso3)) +
  geom_point() +
  scale_x_log10()

#Объект "codes" содержит код стран и соответствующие названия регионов. Добавим в таблицу объекта "df" столбец, который будет содержать названия регионов:
df_merged <- merge(df, codes)

#Затем выполним код с функцией "ggplot" ниже:
ggplot(df_merged[df_merged$year==2016,],
       aes(gdp_per_capita,fertility,size=pop,colour=region)) +
  geom_point() +
  scale_x_log10()

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

Продолжим! Теперь удалим легенду и добавим подписи к осям с помощью аргумента labs:

ggplot(df[df$year==2016,],
       aes(gdp_per_capita,fertility,size=pop,colour=iso3)) +
  geom_point(show.legend = FALSE) + # убираем легенду
  scale_x_log10() +
  labs(x="GDP per capita (constant 2010 USD)",       # добавляем подпись к оси х
       y="Fertility rate, total (births per woman)") # добавляем подпись к оси у

Мы можем отразить на графике временной аспект, указав в качестве источника данных объект с оригинальной таблицей данных ‘df’ (не извлекая из нее данные только для 2016 года):

ggplot(df,
       aes(gdp_per_capita,fertility,size=pop,colour=iso3)) +
  geom_point(show.legend = FALSE) +
  scale_x_log10() +
  labs(x="GDP per capita (constant 2010 USD)",
       y="Fertility rate, total (births per woman)")

Более того, с помощью функции ggplot мы можем создавать динамические графики. Для этого нужно установить и загрузить два дополнительных пакета функций с помощью приведенного ниже кода. (Не забудьте удалить знак "#" перед строчками кода ниже.)

#install.packages('gifski')
library("gifski")
#install.packages('tweenr')
library("tweenr")
# install.packages("devtools")
# devtools::install_github('thomasp85/transformr')
# devtools::install_github('thomasp85/gganimate')
library("gganimate")
#install.packages('png')
library("png")

Создадим динамический график и сохраним его в формате .gif:

animated_ggplot <- ggplot(df,
       aes(gdp_per_capita,fertility,size=pop,colour=iso3)) + # добавим аргумент "size=pop"
  geom_point(show.legend = FALSE) +
  scale_x_log10() +
  labs(x="GDP per capita (constant 2010 USD)",
       y="Fertility rate, total (births per woman)") +
  transition_time(year)

anim_save("animated_ggplot.gif", animated_ggplot)

Выполнив этот код в RStudio, обязательно откройте вкладку “Viewer” справа, чтобы просмотреть анимацию!

Функция transition_time(year) генерирует динамический график, который показывает значения интересующих нас переменных в разные годы.

Мы можем добавить меняющийся заголовок (Year) к графику с помощью аргумента labs - title="Year: {frame_time}":

animated_ggplot2 <- ggplot(df,
       aes(gdp_per_capita,fertility,size=pop,colour=iso3)) + # добавим аргумент "size=pop"
  geom_point(show.legend = FALSE) +
  scale_x_log10() +
  labs(x="GDP per capita (constant 2010 USD)",
       y="Fertility rate, total (births per woman)",
       title="Year: {frame_time}") +
  transition_time(year)

anim_save("animated_ggplot2.gif", animated_ggplot2)

Мы можем настроить высоту и ширину графика:

options(gganimate.dev_args = list(width = 780, height = 440))

Теперь мы можем взять графический объект gganimate и конвертировать его в анимацию с помощью функции animate. Чтобы просмотреть все возможные аргументы функции "animate", выполните код ниже:

?animate 

Также можно ввести название этой функции в строке поиска Google, добавив слово “stackoverflow”.

Код ниже позволит замедлить анимацию animated_ggplot2:

animate(animated_ggplot2, nframes=350, fps=20) 

# настраиваем скорость отображения анимации, увеличивая или уменьшая количество кадров в секунду

# Объедините объект "df" и объект "codes"
df_merged <- merge(df, codes)

# Создайте динамический график с помощью функции "ggplot". В аргументах функции укажите, что цвет точек обозначает регионы
animated_ggplot3 <- ggplot(df_merged,
       aes(gdp_per_capita,fertility,size=pop,colour=region)) +
  geom_point(show.legend = FALSE) +
  scale_x_log10() +
  labs(x="GDP per capita (constant 2010 USD)",
       y="Fertility rate, total (births per woman)",
       title="Year: {frame_time}") +
  transition_time(year)

# Настройте параметры анимации (такие же, как в предыдущей анимации)
animate(animated_ggplot3, nframes=350, fps=20) 

5.1 Линейная гравитационная модель

5.1.1 Линейная регрессия

Начнем с построения самой простой гравитационной модели - линейной регрессии - с использованием функции lm (lm расшифровывается как "линейная модель").

Сначала сгенерируем значения переменных x и y с помощью функции rnorm. (Вернитесь в Раздел 1.2, чтобы вспомнить назначение этой функции, или выполните код ?rnorm.)

y <- rnorm(100,1,5)
x <- y * 2 + rnorm(100,0,1)

В аргументе функции в круглых скобках слева от тильды (~) указывают зависимую переменную, а справа - одну или более независимых переменных. В нашем случае мы построили уравнение модели y = x + c, где c является константой, значение которой будет определено в результате выполнения кода. reg1 - имя объекта, в который мы сохраняем результаты регрессии.

reg1 <- lm(y~x)
reg1

##
## Call:
## lm(formula = y ~ x)
##
## Coefficients:
## (Intercept)            x
##    -0.03085      0.49544

Мы можем просмотреть результаты регрессии при помощи функции summary. Она отобразит все базовые параметры линейной регрессии: значения коэффициентов, стандартные погрешности, p-значения, и пр. (Функция summary может применяться в R как к объектам, содержащим наборы данных, так и к объектам, содержащим параметры моделей.)

Применим функцию summary к объекту reg1. В первой строчке полученного результата содержится описание модели. В нашем случае: Call: lm(formula = y ~ x). Далее представлена информация о распределении регрессионных остатков (Residuals): минимальном (Min) и максимальном значениях (Max), медиане (Median), первой квартили (1Q) и третьей квартили (3Q). Например, значение медианы (Median) в нашем случае составило 0,04081.

Далее идут расчетные коэффициенты (Coefficients). В столбце Estimate указаны значения коэффициентов, а в столбце Std. Error - их стандартные отклонения. t-значения приведены в столбце t value, а p-значения - в столбце Pr(>|t|). В объекте reg1 значение свободного члена (intercept) составило 0,017038, а значение углового коэффициента (x) - 0,486928. Значение коэффициента x является статистически значимым, согласно показателю в столбце Pr(>|t|).

Под таблицей коэффициентов регрессии (Coefficients) приведена легенда, описывающая возможные уровни статистической значимости (Signif. codes). Например, уровень 0.001 ‘**’, помеченный **, соответствует уровню статистической значимости в 0,1%.

Далее приведены значения стандартного отклонения остатков (Residual standard error), коэффициента детерминации (Multiple R-squared) и скорректированного коэффициента детерминации (Adjusted R-squared), F-критерия (F-statistic) и его p-значения (p-value). В случае регрессионной модели reg1 коэффициент детерминации равен 0,9892.

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

summary(reg1)

##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
##      Min       1Q   Median       3Q      Max
## -1.40649 -0.29387 -0.02649  0.33124  1.29846
##
## Coefficients:
##              Estimate Std. Error t value            Pr(>|t|)
## (Intercept) -0.030853   0.053370  -0.578               0.565
## x            0.495440   0.005375  92.171 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5049 on 98 degrees of freedom
## Multiple R-squared:  0.9886, Adjusted R-squared:  0.9885
## F-statistic:  8496 on 1 and 98 DF,  p-value: < 0.00000000000000022

С помощью функции attributes можно просмотреть полный лист атрибутов, доступных через функцию summary, которую выше мы применяли к объекту reg1:

attributes(summary(reg1))

## $names
##  [1] "call"          "terms"         "residuals"     "coefficients"
##  [5] "aliased"       "sigma"         "df"            "r.squared"
##  [9] "adj.r.squared" "fstatistic"    "cov.unscaled"
##
## $class
## [1] "summary.lm"

Вы также можете извлечь значение отдельного коэффициента, например, коэффициента детерминации \(R^2\):

summary(reg1)$r.squared

## [1] 0.988596

# Ответ верный. Супер!

#Ответ B неверный. В коде не заданы значения переменных x и y. Верный ответ A.

#Ответ C неверный. Код задает значения переменных x и y, но не задает уравнения регрессии, а также не запрашивает отобразить результат регрессии. Правильный ответ A.

gd <- read.csv("short_data.csv",
               stringsAsFactors = FALSE)

names(gd)

##  [1] "X"                 "reporter"          "partner"
##  [4] "year"              "export_ij"         "import_ij"
##  [7] "trade_ij"          "export_i"          "import_i"
## [10] "trade_i"           "export_j"          "import_j"
## [13] "trade_j"           "ctcij"             "ntcij"
## [16] "tariffrateij_sa"   "tariffrateij_wa"   "tariff_ij_sa"
## [19] "tariff_ij_wa"      "tariffrateji_sa"   "tariffrateji_wa"
## [22] "tariff_ji_sa"      "tariff_ji_wa"      "ga_tariff_ijji_sa"
## [25] "ga_tariff_ijji_wa" "dist"              "contig"
## [28] "comlang_off"       "comlang_ethno"     "colony"
## [31] "comcol"            "curcol"            "col45"
## [34] "smctry"            "landlocked_i"      "landlocked_j"
## [37] "landlocked_ij"     "rta"               "gdp_i"
## [40] "gdpgrowth_i"       "gdppc_i"           "gdppcgrowth_i"
## [43] "rgdp_i"            "rgdppc_i"          "gdp_j"
## [46] "gdpgrowth_j"       "gdppc_j"           "gdppcgrowth_j"
## [49] "rgdp_j"            "rgdppc_j"          "dtf_i"
## [52] "tab_i"             "startbiz_i"        "credit_i"
## [55] "invest_i"          "tax_i"             "contract_i"
## [58] "dtf_j"             "tab_j"             "startbiz_j"
## [61] "credit_j"          "invest_j"          "tax_j"
## [64] "contract_j"        "lsci_i"            "lsci_j"
## [67] "lpi_i"             "lpi_j"

lm(gd$export_ij ~ gd$dist)

##
## Call:
## lm(formula = gd$export_ij ~ gd$dist)
##
## Coefficients:
## (Intercept)      gd$dist
##  1094335800       -73723

reg1 <- lm(gd$export_ij ~ gd$dist)

summary(reg1)

##
## Call:
## lm(formula = gd$export_ij ~ gd$dist)
##
## Residuals:
##          Min           1Q       Median           3Q          Max
##  -1090653864   -758408424   -507452840   -200191922 409695392687
##
## Coefficients:
##               Estimate Std. Error t value            Pr(>|t|)
## (Intercept) 1094335800   16216752   67.48 <0.0000000000000002 ***
## gd$dist         -73723       1932  -38.16 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5307000000 on 374236 degrees of freedom
##   (137806 observations deleted due to missingness)
## Multiple R-squared:  0.003875,   Adjusted R-squared:  0.003872
## F-statistic:  1456 on 1 and 374236 DF,  p-value: < 0.00000000000000022

nrow(gd)

## [1] 512044

nrow(gd[gd$export_ij==0,])

## [1] 120815

length(which(gd$export_ij==0))

## [1] 6409

length(gd)

## [1] 68

length(gd$export_ij)

## [1] 512044

length(gd$export_ij[is.na(gd$export_ij)])

## [1] 114406

gd[is.na(gd$export_ij),"export_ij"] <- 0

length(gd$export_ij[is.na(gd$export_ij)])

## [1] 0

reg2 <- lm(gd$export_ij ~ gd$dist)
summary(reg2)

##
## Call:
## lm(formula = gd$export_ij ~ gd$dist)
##
## Residuals:
##          Min           1Q       Median           3Q          Max
##   -939688372   -624180945   -404594794   -132356225 409761225877
##
## Coefficients:
##              Estimate Std. Error t value            Pr(>|t|)
## (Intercept) 943289085   13155595   71.70 <0.0000000000000002 ***
## gd$dist        -65972       1497  -44.07 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4698000000 on 478830 degrees of freedom
##   (33212 observations deleted due to missingness)
## Multiple R-squared:  0.004039,   Adjusted R-squared:  0.004037
## F-statistic:  1942 on 1 and 478830 DF,  p-value: < 0.00000000000000022

reg3 <- lm(log(gd$export_ij) ~ log(gd$dist))

gd[gd$export_ij==0,"export_ij"] <- 0.0001
reg3 <- lm(log(gd$export_ij) ~ log(gd$dist))
summary(reg3)

##
## Call:
## lm(formula = log(gd$export_ij) ~ log(gd$dist))
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -32.756  -3.341   3.907   7.381  19.156
##
## Coefficients:
##              Estimate Std. Error t value            Pr(>|t|)
## (Intercept)  36.37332    0.15809   230.1 <0.0000000000000002 ***
## log(gd$dist) -3.09398    0.01816  -170.4 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.51 on 478830 degrees of freedom
##   (33212 observations deleted due to missingness)
## Multiple R-squared:  0.05717,    Adjusted R-squared:  0.05717
## F-statistic: 2.904e+04 on 1 and 478830 DF,  p-value: < 0.00000000000000022

reg4 <- lm(log(export_ij) ~ log(dist), data = gd)
summary(reg4)

##
## Call:
## lm(formula = log(export_ij) ~ log(dist), data = gd)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -32.756  -3.341   3.907   7.381  19.156
##
## Coefficients:
##             Estimate Std. Error t value            Pr(>|t|)
## (Intercept) 36.37332    0.15809   230.1 <0.0000000000000002 ***
## log(dist)   -3.09398    0.01816  -170.4 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.51 on 478830 degrees of freedom
##   (33212 observations deleted due to missingness)
## Multiple R-squared:  0.05717,    Adjusted R-squared:  0.05717
## F-statistic: 2.904e+04 on 1 and 478830 DF,  p-value: < 0.00000000000000022

#Переменная `trade_ij` содержит нули. Выполним те же действия, которые мы выполнили для переменной "export_ij":
gd[gd$trade_ij==0,"trade_ij"] <- 0.0001

#Теперь оценим следующее гравитационное уравнение:
lm(log(trade_ij) ~ log(dist) + log(export_ij), data = gd)

#Согласно полученному результату, значение коэффициента для log(export_ij) равно 0,2976 или 0,30, а значит правильным является ответ А.

reg5 <- lm(log(export_ij) ~ log(dist) +
             log(gdp_i) +
             log(gdp_j) +
             rta +
             log(ga_tariff_ijji_sa)
             , data = gd)
summary(reg5)

##
## Call:
## lm(formula = log(export_ij) ~ log(dist) + log(gdp_i) + log(gdp_j) +
##     rta + log(ga_tariff_ijji_sa), data = gd)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -33.334  -1.688   1.713   4.409  19.737
##
## Coefficients:
##                          Estimate Std. Error  t value            Pr(>|t|)
## (Intercept)            -62.352427   0.273433 -228.036 <0.0000000000000002 ***
## log(dist)               -2.516327   0.018142 -138.703 <0.0000000000000002 ***
## log(gdp_i)               2.489102   0.006209  400.867 <0.0000000000000002 ***
## log(gdp_j)               1.381105   0.006204  222.598 <0.0000000000000002 ***
## rta                      1.322291   0.042079   31.424 <0.0000000000000002 ***
## log(ga_tariff_ijji_sa)  -0.224746   0.242981   -0.925               0.355
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.107 on 301412 degrees of freedom
##   (210626 observations deleted due to missingness)
## Multiple R-squared:  0.4224, Adjusted R-squared:  0.4224
## F-statistic: 4.409e+04 on 5 and 301412 DF,  p-value: < 0.00000000000000022

write.table(
  summary(reg5)$coefficients #выберем только коэффициенты — это тоже сработает
  , "clipboard", sep="\t", row.names=FALSE)

reg6 <- lm(log(export_ij) ~ log(dist) +
             log(gdp_i) +
             log(gdp_j) +
             rta +
             I(rta*log(gdp_j)) + #включаем компоненту взаимодействия с помощью оператора "|"
             I(dist^2) + #полиномиальное преобразование первого порядка
             log(ga_tariff_ijji_sa)
           , data = gd)
summary(reg6)

##
## Call:
## lm(formula = log(export_ij) ~ log(dist) + log(gdp_i) + log(gdp_j) +
##     rta + I(rta * log(gdp_j)) + I(dist^2) + log(ga_tariff_ijji_sa),
##     data = gd)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -34.160  -1.726   1.697   4.412  20.456
##
## Coefficients:
##                                 Estimate        Std. Error  t value
## (Intercept)            -63.9832459652895   0.3271008450554 -195.607
## log(dist)               -2.9782272151947   0.0263709439899 -112.936
## log(gdp_i)               2.5488480880571   0.0062645188343  406.871
## log(gdp_j)               1.5342760391309   0.0067795539273  226.309
## rta                     23.4357915076500   0.4081369407485   57.421
## I(rta * log(gdp_j))     -0.8966090092956   0.0163094395377  -54.975
## I(dist^2)                0.0000000062601   0.0000000002558   24.469
## log(ga_tariff_ijji_sa)  -1.1094688930275   0.2419365488330   -4.586
##                                    Pr(>|t|)
## (Intercept)            < 0.0000000000000002 ***
## log(dist)              < 0.0000000000000002 ***
## log(gdp_i)             < 0.0000000000000002 ***
## log(gdp_j)             < 0.0000000000000002 ***
## rta                    < 0.0000000000000002 ***
## I(rta * log(gdp_j))    < 0.0000000000000002 ***
## I(dist^2)              < 0.0000000000000002 ***
## log(ga_tariff_ijji_sa)           0.00000452 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.063 on 301410 degrees of freedom
##   (210626 observations deleted due to missingness)
## Multiple R-squared:  0.4295, Adjusted R-squared:  0.4295
## F-statistic: 3.242e+04 on 7 and 301410 DF,  p-value: < 0.00000000000000022

#install.packages("estimatr")
library(estimatr)

reg_robust <- lm_robust(log(export_ij) ~ log(dist) + log(gdp_i) + log(gdp_j) + rta + log(ga_tariff_ijji_sa),
          data = gd, clusters = dist, se_type = "stata")
summary(reg_robust)

##
## Call:
## lm_robust(formula = log(export_ij) ~ log(dist) + log(gdp_i) +
##     log(gdp_j) + rta + log(ga_tariff_ijji_sa), data = gd, clusters = dist,
##     se_type = "stata")
##
## Standard error type:  stata
##
## Coefficients:
##                        Estimate Std. Error  t value
## (Intercept)            -62.3524    0.86939 -71.7201
## log(dist)               -2.5163    0.05655 -44.4980
## log(gdp_i)               2.4891    0.01948 127.7587
## log(gdp_j)               1.3811    0.01766  78.2119
## rta                      1.3223    0.10906  12.1240
## log(ga_tariff_ijji_sa)  -0.2247    0.64297  -0.3495
##                                                      Pr(>|t|) CI Lower CI Upper
## (Intercept)            0.000000000000000000000000000000000000  -64.057  -60.648
## log(dist)              0.000000000000000000000000000000000000   -2.627   -2.405
## log(gdp_i)             0.000000000000000000000000000000000000    2.451    2.527
## log(gdp_j)             0.000000000000000000000000000000000000    1.346    1.416
## rta                    0.000000000000000000000000000000001213    1.109    1.536
## log(ga_tariff_ijji_sa) 0.726687161813804882726230971456971020   -1.485    1.036
##                           DF
## (Intercept)            12607
## log(dist)              12607
## log(gdp_i)             12607
## log(gdp_j)             12607
## rta                    12607
## log(ga_tariff_ijji_sa) 12607
##
## Multiple R-squared:  0.4224 ,    Adjusted R-squared:  0.4224
## F-statistic:  4266 on 5 and 12607 DF,  p-value: < 0.00000000000000022

summary(reg5)

##
## Call:
## lm(formula = log(export_ij) ~ log(dist) + log(gdp_i) + log(gdp_j) +
##     rta + log(ga_tariff_ijji_sa), data = gd)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -33.334  -1.688   1.713   4.409  19.737
##
## Coefficients:
##                          Estimate Std. Error  t value            Pr(>|t|)
## (Intercept)            -62.352427   0.273433 -228.036 <0.0000000000000002 ***
## log(dist)               -2.516327   0.018142 -138.703 <0.0000000000000002 ***
## log(gdp_i)               2.489102   0.006209  400.867 <0.0000000000000002 ***
## log(gdp_j)               1.381105   0.006204  222.598 <0.0000000000000002 ***
## rta                      1.322291   0.042079   31.424 <0.0000000000000002 ***
## log(ga_tariff_ijji_sa)  -0.224746   0.242981   -0.925               0.355
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.107 on 301412 degrees of freedom
##   (210626 observations deleted due to missingness)
## Multiple R-squared:  0.4224, Adjusted R-squared:  0.4224
## F-statistic: 4.409e+04 on 5 and 301412 DF,  p-value: < 0.00000000000000022

# Выполните код:
lm_robust(log(trade_ij) ~ log(dist) + log(export_ij), data = gd, clusters = dist, se_type = "stata")

# Полученное значение стандартных ошибок для log(export_ij) равно 0,001470578 или 0,0015. Значит, ответ А верный.

5.2 Метод квази-максимального правдоподобия Пуассона (PPML)

Метод квази-максимального правдоподобия Пуассона (PPML) часто используется при анализе торговли, так как он работает и при наличии нулевых значений переменных, позволяет учесть фиксированные эффекты и облегчает интерпретацию результатов. Рассмотрим базовые компоненты модели PPML в R. Дополнительная информация о методе доступна в руководстве ЭСКАТО по выполнению гравитационного анализа в R.

Функция PPML в R входит в пакет gravity:

#install.packages("gravity")
library(gravity)

В качестве зависимой переменной выбираем export_ij. В качестве независимой переменной указываем расстояние между экономиками (переменная dist), а также включаем любые дополнительные независимые переменные (например, ВВП страны-партнера по экспорту gdp_j и ВВП страны-экспортера gdp_i). Отметим важный момент: когда мы выполняем в R регрессию для переменной dist по методу PPML, мы трансформируем её с помощью натурального логарифма, аналогично использованию функции log в R.

ppml("export_ij", dist="dist", c('gdp_j' ,'gdp_i' ), data=gd)

##
## Call:  y_ppml ~ dist_log + gdp_j + gdp_i
##
## Coefficients:
##         (Intercept)             dist_log                gdp_j
## 28.0022547089188549  -1.0674026797574185   0.0000000000003008
##               gdp_i
##  0.0000000000002836
##
## Degrees of Freedom: 435085 Total (i.e. Null);  435082 Residual
##   (43746 observations deleted due to missingness)
## Null Deviance:       1423000000000000
## Residual Deviance: 724400000000000   AIC: NA

Как было отмечено выше, одним из достоинств PPML является то, что он может работать с нулевыми значениями переменных. Поэтому вернем нулевые значения переменной export_ij:

gd[gd$export_ij==0.0001,"export_ij"] <- 0
ppml("export_ij", dist="dist", c('gdp_j' ,'gdp_i' ), data=gd)

##
## Call:  y_ppml ~ dist_log + gdp_j + gdp_i
##
## Coefficients:
##         (Intercept)             dist_log                gdp_j
## 28.0022547089189793  -1.0674026797574607   0.0000000000003008
##               gdp_i
##  0.0000000000002836
##
## Degrees of Freedom: 435085 Total (i.e. Null);  435082 Residual
##   (43746 observations deleted due to missingness)
## Null Deviance:       1423000000000000
## Residual Deviance: 724400000000000   AIC: NA

Во-первых, обратите внимание, что переменная dist выражена через натуральный логарифм.

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

# Выполните код ниже, если вы еще этого не сделали:
gd[gd$export_ij==0.0001,"export_ij"] <- 0

# Затем укажите PPML модель:
ppml("export_ij", dist="dist", c('gdp_j' ,'gdp_i' ), data=gd)

# Выполнив весь приведенный выше код, вы получите значение коэффициента `rta` равное 1,103045. Таким образом, ответ D верный. 

6.1 Коэффициент покрытия

Рассчитаем сначала коэффициент покрытия, который описывает, какая доля всего объема импорта конкретной страны регулируется нетарифными мерами. Ниже приведена формула для расчета коэффициента покрытия:

где \(NTM_{ik}\) - бинарная переменная, которая показывает, регулируется ли импортируемый товар \(k\) хотя бы одной НТМ; а \(X_{ik}\) - стоимостный объем импорта товара \(k\) в страну \(i\).

Для вычисления индикатора по этой формуле нам нужны значения только двух переменных: \(NTM_{ik}\) и \(X_{ik}\). У нас уже есть значения стоимостного объёма импорта каждого вида товаров \(X_{ik}\) (столбец value в объекте lao_imp). В отношении значений этой переменной не нужно проводить никаких дополнительных расчетов. Чтобы получить значения переменной \(NTM_{ik}\) (указывает, регулируется ли импорт отдельного товара нетарифной мерой), извлечем все коды товаров (hs6) из объекта lao_ntm. Мы затем посчитаем количество НТМ, применяемых к каждому отдельному коду товара. Мы присвоим значение 1 переменной \(NTM_{ik}\), если к одному коду товара применяется хотя бы одна НТМ： (\(NTM_{ik} = 1\)), и значение 0, если к одному коду товара не применяется ни одной НТМ： (\(NTM_{ik} = 0\)).

Сделаем это пошагово.

Сначала нам необходимо создать переменную \(NTM_{ik}\), которая будет показывать наличие/отсутствие НТМ, применяемой к коду товара. Сначала возьмем объект lao_ntm. Код ниже позволяет посчитать количество НТМ, установленных для каждого кода товара в столбце hs6:

lao_ntm %>% group_by(hs6) %>% count()

## # A tibble: 5,205 x 2
## # Groups:   hs6 [5,205]
##      hs6     n
##    <int> <int>
##  1 10121    10
##  2 10129    13
##  3 10130    13
##  4 10190    13
##  5 10221    10
##  6 10229    13
##  7 10231    10
##  8 10239    13
##  9 10290    13
## 10 10310    10
## # ... with 5,195 more rows

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

Таким образом, приведенный выше код lao_ntm %>% group_by(hs6) %>% count()выполняет следующие операции (при этом код не меняет ничего в объекте, так как в данном случае мы не используем функцию <-):
1. берет объект с таблицей данных lao_ntm
2. group_by(hs6): группирует данные на основе переменных в столбце hs6 (то есть несколько строк с одинаковым кодом товара(hs6) группируются в одну строку)
3. count(): считает количество строк в оригинальной таблице данных, в которых были указаны одинаковые коды товара в столбце (hs6)

Мы используем функцию count, чтобы суммировать количество упоминаний одного и того же кода товара. Обратите внимание, что в данном случае в аргументе функции count мы не указываем названия объекта-источника данных. Это связано с тем, что этот источник данных - столбец hs6 - уже указан в аргументе предыдущей функции.

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

ntm_hs_freq <- lao_ntm %>% group_by(hs6) %>% count()

head(ntm_hs_freq)

## # A tibble: 6 x 2
## # Groups:   hs6 [6]
##     hs6     n
##   <int> <int>
## 1 10121    10
## 2 10129    13
## 3 10130    13
## 4 10190    13
## 5 10221    10
## 6 10229    13

На этом этапе выборка данных в объекте ntm_hs_freq не верная, так как в ней содержатся меры под кодом “P” - меры, применяемые к экспорту товаров экспортирующими странами. (Как было упомянуты выше, нас интересует импорт в Лаосскую НДР, а не экспорт, а значит эти меры нам необходимо исключить из выборки.)

Мы можем извлечь из объекта lao_ntm выборку, которая не содержит меры с кодом “P”, а затем пересчитать количество НТМ, применяемых к индивидуальным кодам товаров с помощью функции count:

ntm_hs_freq <- lao_ntm[lao_ntm$ntm_1_digit!="P",] %>%  #конвейерный оператор в конце строки говорит о том, что следующая строка содержит продолжение кода
  group_by(hs6) %>%
  count()

head(ntm_hs_freq)

## # A tibble: 6 x 2
## # Groups:   hs6 [6]
##     hs6     n
##   <int> <int>
## 1 10121     6
## 2 10129     9
## 3 10130     9
## 4 10190     9
## 5 10221     6
## 6 10229     9

Если вы поставите оператор %>% в начале строки, то код не выполнится. Но если вы поставите оператор %>% в конце строки, он сделает переход на следующую строку кода.

Кратко обобщим материал данного подраздела. Мы создали объект ntm_hs_freq c таблицей данных, которая в столбце hs6содержит все коды товаров, к которым применяются НТМ. В столбце n содержится число, отражающее количество НТМ, применяемых к конкретным кодам товаров. Теперь добавим в этот объект данные по импорту товаров в Лаосскую НДР. Для этого объединим объекты lao_imp и ntm_hs_freq.

Объединение объектов lao_imp и ntm_hs_freq осуществим по тем столбцам каждого их них, которые содержат коды товаров. Поскольку названия этих столбцов в каждом объекте совпадают, используем аргумент by:

coverage <- merge(lao_imp,ntm_hs_freq,by="hs6",
                  all.x = TRUE)
head(coverage)

##     hs6 partner    value n
## 1 10121     WLD     1430 6
## 2 10190     WLD     5100 9
## 3 10221     WLD   243284 6
## 4 10229     WLD 32125452 9
## 5 10231     WLD    22480 6
## 6 10239     WLD  8630620 9

Использование аргумента all.x = TRUE позволяет сохранить все строки в объекте lao_imp, иначе мы потеряем строки с кодами тех товаров, к которым не применяются НТМ.

Для расчета значения коэффициента покрытия воспользуемся данными в столбцах объекта coverage. Как вы помните, для расчёта нам нужны значения переменной value (\(X_{ik}\)) и фиктивной переменной \(NTM_{ik}\), указывающей на то, регулируется ли отдельный код товара хотя бы одной НТМ. Чтобы задать значения фиктивной переменной, воспользуемся значениями в столбце n.

Создадим новый столбец NTMik (\(NTM_{ik}\)), в ячейки которого сохраним нули:

coverage$NTMik <- 0

Напомним, что \(NTM_{ik} = 0\), если к коду товара не применяется ни одной НТМ. \(NTM_{ik} = 1\), если к коду товара применяется хотя бы одна НТМ.

Поэтому, если \(n > 0\), это значит, что к коду товара применятся хотя бы одна НТМ. А значит в столбце NTMik в соответствующую ячейку нужно записать значение 1:

coverage[coverage$n>0,"NTMik"] <- 1

В первой части формулы для расчёта индикатора мы умножаем NTMik на стоимостный объем импорта (переменная \(X_{ik}\) в формуле):

head(coverage$NTMik*coverage$value)

## [1]     1430     5100   243284 32125452    22480  8630620

Теперь суммируем полученные значения:

sum(coverage$NTMik*coverage$value)

## [1] 4107068464

Разделим результат на сумму всех значений в столбце value (переменная \(X_{ik}\)), а затем умножим на 100, чтобы получить результат в процентах:

sum(coverage$NTMik*coverage$value)/sum(coverage$value)*100

## [1] 100

Согласно этому результату, 100% всех импортируемых в Лаосскую НДР товаров регулируются нетарифными мерами.

После того, как мы разобрались с кодом, повторим выполненные нами шаги:
1. Извлекли нужные данные и записали их в объект ntm_hs_freq.
2. Объединили объекты ntm_hs_freq и lao_imp.
3. Создали переменную NTMik и присвоили ей нужные значения.
4. Вычислили с помощью формулы значение индикатора.

ntm_hs_freq <- lao_ntm[lao_ntm$ntm_1_digit=="A",] %>%  #конвейерный оператор в конце строки говорит о том, что следующая строка содержит продолжение кода
  group_by(hs6) %>%
  count()

ntm_hs_freq

## # A tibble: 1,139 x 2
## # Groups:   hs6 [1,139]
##      hs6     n
##    <int> <int>
##  1 10121     3
##  2 10129     4
##  3 10130     4
##  4 10190     4
##  5 10221     3
##  6 10229     4
##  7 10231     3
##  8 10239     4
##  9 10290     4
## 10 10310     3
## # ... with 1,129 more rows

coverage <- merge(lao_imp,ntm_hs_freq,
                  by="hs6",
                  all.x = TRUE)
head(coverage)

##     hs6 partner    value n
## 1 10121     WLD     1430 3
## 2 10190     WLD     5100 4
## 3 10221     WLD   243284 3
## 4 10229     WLD 32125452 4
## 5 10231     WLD    22480 3
## 6 10239     WLD  8630620 4

coverage$NTMik <- 0

coverage[coverage$n>0,"NTMik"] <- 1

coverage[is.na(coverage$n),"n"] <- 0
head(coverage)

##     hs6 partner    value n NTMik
## 1 10121     WLD     1430 3     0
## 2 10190     WLD     5100 4     0
## 3 10221     WLD   243284 3     0
## 4 10229     WLD 32125452 4     0
## 5 10231     WLD    22480 3     0
## 6 10239     WLD  8630620 4     0

coverage[is.na(coverage$n), ]

## [1] hs6     partner value   n       NTMik
## <0 rows> (or 0-length row.names)

coverage[coverage$n>0,"NTMik"] <- 1

sum(coverage$NTMik*coverage$value)/
  sum(coverage$value)*100

## [1] 15.57878

#Чтобы извлечь только те строки объекта "lao_ntm", которые содержат меры из главы В, нам нужно установить фильтр, выбирающий только те строки, которые в ячейке столбца "ntm_1_digit" содержат код В. 22% - неверное значение индикатора. Выполнив правильный код полностью, получаем ответ 32%. Значит, ответ В верный!

#Правильный код:
ntm_hs_freq <- ntm_hs_freq <- lao_ntm[lao_ntm$ntm_1_digit=="B",] %>%
     group_by(hs6) %>%
     count()

coverage <- merge(lao_imp,ntm_hs_freq,
                  by="hs6",
                  all.x = TRUE)

coverage$NTMik <- 0
coverage[is.na(coverage$n),"n"] <- 0
coverage[coverage$n>0,"NTMik"] <- 1

sum(coverage$NTMik*coverage$value)/
  sum(coverage$value)*100

## [1] 32.20915

#Ответ B верный! Супер!   

6.2 Индекс частоты

Индекс частоты похож на коэффициент покрытия. Только в этом случае значение фиктивной переменной присутствия/отсутствия НТМ мы умножаем не на объем импорта, а на значения фиктивной переменной наличия/отсутствия импорта продукции, относящейся к отдельному коду ГС. Поэтому значение этого коэффициента описывает, какая доля импортируемых товарных позиций регулируется нетарифными мерами. Ниже приведена формула для расчёта данного показателя:

где \(NTM_{ik}\) - бинарная переменная, которая показывает, регулируется ли импортируемый товар \(k\) хотя бы одной НТМ; \(D_{ik}\) - бинарная переменная, которая показывает, импортируется ли товар \(k\), относящийся к конкретному коду ГС, в страну.

Создаем соответствующую выборку данных из "lao_ntm", группируем строки по коду ГС и считаем количество строк с отдельными кодами: В этом случае мы также исключаем строки с экспортными мерами (“P”).

ntm_hs_freq <- lao_ntm[lao_ntm$ntm_1_digit!="P",] %>%
  group_by(hs6) %>%
  count()

Чтобы получить индекс частоты, объединяем объекты lao_imp и ntm_hs_freq.

freq_index <- merge(lao_imp,
                    ntm_hs_freq,
                    by="hs6",
                    all.x = TRUE)

Создаем новую переменную Dik, которая равна 1, если страна \(i\) импортирует любое количество товара \(k\), и равна 0, если рассматриваемый товар не импортируется:

freq_index$Dik <- 0

Как сгенерировать фиктивную переменную, которая равна 1, если в Лаосскую НДР импортируется любое количество отдельного товара, и нулю, если импорт товара равен нулю?

Переменная Dik должна быть равна 1, если в Лаосскую НДР было импортировано товара хотя бы на 1 доллар.

Начнем с проверки противоположного критерия и найдем все строки, где значение переменной value является NA (отсутствующее значение). Это значит, что Лаосская НДР не импортировала эти товары:

head(freq_index[is.na(freq_index$value),])

## [1] hs6     partner value   n       Dik
## <0 rows> (or 0-length row.names)

Чтобы выбрать строки, где значение переменной "value" не является отсутствующим, добавим “!” перед is.na:

head(freq_index[!is.na(freq_index$value),])

##     hs6 partner    value n Dik
## 1 10121     WLD     1430 6   0
## 2 10190     WLD     5100 9   0
## 3 10221     WLD   243284 6   0
## 4 10229     WLD 32125452 9   0
## 5 10231     WLD    22480 6   0
## 6 10239     WLD  8630620 9   0

Зададим нашей переменной значение 1, если в соответствующей ячейке столбца "value" нет отсутствующих значений:

freq_index[!is.na(freq_index$value),"Dik"] <- 1

Строго говоря, мы не можем использовать критерий \(n > 0\), так как мы имеем дело с отсутствующими значениями.

Сгенерируем теперь переменную \(NTM_ik\). Напомним, что значение переменной NTMik (\(NTM_ik\)) показывает, регулируется ли отдельный код товара хотя бы одной НТМ. Снова создадим новую переменную NTMik. Сначала заполним этот столбец нулями, а затем сохраним значение 1 в те ячейки, для которых соблюдается условие \(n > 0\).

freq_index$NTMik <- 0

freq_index[freq_index$n>0,"NTMik"] <- 1

Затем разделим полученное значение на сумму всех значений переменной Dik, и умножим на 100.

sum(freq_index$NTMik*freq_index$Dik)/
  sum(freq_index$Dik)*100

## [1] 100

Это значит, что 100% импортируемых товаров регулируются хотя бы одной НТМ.

ntm_hs_freq <- lao_ntm[lao_ntm$ntm_1_digit=="A",] %>%
     group_by(hs6) %>%
     count()

freq_index <- merge(lao_imp,
                    ntm_hs_freq,
                    by="hs6",
                    all.x = TRUE)

freq_index$Dik <- 0

freq_index[!is.na(freq_index$value),"Dik"] <- 1

freq_index$NTMik <- 0

freq_index[is.na(freq_index$n),"n"] <- 0

freq_index[freq_index$n>0,"NTMik"] <- 1

sum(freq_index$NTMik*freq_index$Dik)/
  sum(freq_index$Dik)*100

## [1] 18.3908

#Измените строку кода, с помощью которой создается объект "ntm_hs_freq"  следующим образом:
ntm_hs_freq <- lao_ntm[lao_ntm$ntm_1_digit=="B"|lao_ntm$ntm_1_digit=="C",] %>%
     group_by(hs6) %>%
  count()

#Затем выполните весь код в RStudio и получите результат 31,41762 или 31%. Это значит, что ответ А верный.
freq_index <- merge(lao_imp,
                    ntm_hs_freq,
                    by="hs6",
                    all.x = TRUE)

freq_index$Dik <- 0
freq_index[!is.na(freq_index$value),"Dik"] <- 1
freq_index$NTMik <- 0
freq_index[is.na(freq_index$n),"n"] <- 0
freq_index[freq_index$n>0,"NTMik"] <- 1

sum(freq_index$NTMik*freq_index$Dik)/
  sum(freq_index$Dik)*100

6.3 Показатель распространенности для СФС мер

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

где \(\#NTM_{ik}\) - число НТМ, относящихся к товару \(k\); а \(D_{ik}\) - бинарная переменная, которая показывает, импортируется ли товар \(k\) в страну.

Воспользуемся объектом freq_index, созданным в предыдущем подразделе, где мы вычисляли индекс частоты, так как этот объект уже содержит переменную n (\(\#NTM_{ik}\)) и переменную Dik (\(D_{ik}\)). Напомним, что СФС меры относятся к главе "А" классификации НТМ.

# создадим объект, отражающий число СФС мер, применяемых к каждому коду товара
ntm_hs_freq <- lao_ntm[lao_ntm$ntm_1_digit=="A",] %>%
     group_by(hs6) %>%
     count()

# объединим его с данными по импорту товаров страну
freq_index <- merge(lao_imp,
                    ntm_hs_freq,
                    by="hs6",
                    all.x = TRUE)
#создадим переменную Dik
freq_index$Dik <- 0
freq_index[!is.na(freq_index$value),"Dik"] <- 1

#заменим отсутствующие значения в столбце "n" на 0
freq_index[is.na(freq_index$n),"n"] <- 0

Ниже приведена формула для расчёта данного показателя:

sum(freq_index$n*freq_index$Dik)/
  sum(freq_index$Dik)

## [1] 0.8125798

Полученный результат говорит о том, что в среднем 0,8 СФС меры применяется к одному коду импортируемого в Лаосскую НДР товара.

Здесь можно посмотреть значения всех показателей для отдельных стран: Ссылка на сайт ЮНКТАД.

#С помощью года ниже создайте выборку данных `ntm_hs_freq`, исключающую экспортные меры:
ntm_hs_freq <- lao_ntm[lao_ntm$ntm_1_digit!="P",] %>%
     group_by(hs6) %>%
  count()

#А теперь используйте код ниже, чтобы получить значение показателя распространенности:
freq_index <- merge(lao_imp,
                    ntm_hs_freq,
                    by="hs6",
                    all.x = TRUE)

freq_index$Dik <- 0
freq_index[!is.na(freq_index$value),"Dik"] <- 1
freq_index$NTMik <- 0
freq_index[is.na(freq_index$n),"n"] <- 0
freq_index[freq_index$n>0,"NTMik"] <- 1


sum(freq_index$n*freq_index$Dik)/sum(freq_index$Dik)

#Полученное значение 3,818008 или 3,8 говорит о том, что 3,8 нетарифных меры приходятся в среднем на 1 товарную позицию импортируемой в Лаос продукции.

7.1 Загрузка из WITS и подготовка данных о торговле Вьетнама

Загружаем данные в R. (Эти данные мы скачали с сайта WITS.)

trade_data_VNM <- read.csv("trade_data_VNM_WITS_2015-2017.csv",
                     stringsAsFactors = FALSE)
head(trade_data_VNM)

##   reporter partner   hs6 year flow   value
## 1      VNM     AUS 10129 2016    5 211.878
## 2      VNM     AUS 10129 2017    5 681.523
## 3      VNM     MYS 10129 2017    5  49.167
## 4      VNM     WLD 10129 2015    5 249.600
## 5      VNM     WLD 10129 2016    5 899.994
## 6      VNM     WLD 10129 2017    5 730.690

Затем мы выбираем только те товары, которые имеют код ГС не выше 980000, так как главы 98 и 99 ГС не гармонизированы. Они предназначены для использования странами по своему усмотрению. Поэтому они несопоставимы между странами. Для этого конвертируем код ГС в числовой формат.

trade_data_VNM$hs6 <- as.numeric(trade_data_VNM$hs6)
trade_data_VNM <- trade_data_VNM[!(trade_data_VNM$hs6 > 980000),]

В столбце "reporter" меняем код для Румынии: вместо “ROM” указываем “ROU”, чтобы соблюсти единообразие с таблицей данных по НТМ. Позже в Разделе 7.4 мы объединим обе таблицы данных и рассчитаем коэффициент покрытия, индекс частоты и показатель распространенности.

trade_data_VNM$reporter[trade_data_VNM$reporter=="ROM"] <- "ROU"

7.2 Загрузка и подготовка данных по нетарифным мерам (НТМ) Вьетнама

Загружаем данные по НТМ и выбираем только те строки, в которых в столбце "reporter" указан Вьетнам. (Оригинальные данные по НТМ можете найти здесь.) Выбираем столбцы с данными о стране, применяющей меру, стране, по отношению к которой применяется мера, классификационном коде меры (четырехзначный код и буква главы классификации), и коде товара (код ГС). Переименовываем столбец, содержащий код НТМ, из ntmcode в NTM и меняем название столбца StartDate на year для удобства.

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

ntms_VNM <- read.dta13("ntm_hs6_2016 v.12.dta")
ntms_VNM <- ntms_VNM[ntms_VNM$reporter=="VNM", ] ## выбираем строки таблицы с данными по Вьетнаму
ntms_VNM <-  ntms_VNM %>% select("partner", "ntmcode",  "hs6", "nbr", "StartDate", "ntm_1_digit") %>% # выбираем интересующие нас столбцы
  rename(NTM = ntmcode, year = StartDate) #задаем новые имена переменных
head(ntms_VNM)

##      partner  NTM    hs6 nbr year ntm_1_digit
## 2525     WLD A140 010121   1 2004           A
## 2754     WLD A150 010121   1 2004           A
## 3219     WLD A310 010121   1 2004           A
## 3514     WLD A810 010121   2 2004           A
## 6095     WLD A890 010121   1 2004           A
## 7193     WLD P130 010121   1 2004           P

Исключаем коды ГС выше 980000, как мы это сделали для выборки данным о торговле. Мы делаем это, так как главы 98 и 99 номенклатуры не гармонизированы. Они предназначены для использования странами по своему усмотрению, поэтому они несопоставимы между странами.

ntms_VNM$hs6 <- as.numeric(ntms_VNM$hs6)
ntms_VNM <- ntms_VNM[!(ntms_VNM$hs6 > 980000),]

EUN <- c('AUT', 'BEL', 'BGR', 'HRV', 'CYP', 'CZE', 'DNK', 'EST', 'FIN', 'FRA', 'DEU', 'GRC', 'HUN', 'IRL', 'ITA', 'LVA', 'LTU', 'LUX', 'MLT', 'NLD', 'POL', 'PRT', 'ROM', 'SVK', 'SVN', 'ESP', 'SWE', 'GBR')
EUN <- paste(EUN,collapse=",")

All_members <- unique(trade_data_VNM$partner)
All_members <- All_members[!(All_members %in% c("WLD","EUN"))]
All_members <- paste0(All_members, collapse=",")

ntms_VNM$partner  <- gsub("WLD", All_members, ntms_VNM$partner)
ntms_VNM$partner  <- gsub("EUN", EUN, ntms_VNM$partner)

ntms_VNM <- ntms_VNM %>% separate_rows(partner, sep = ",")
nrow(ntms_VNM)

## [1] 2726605

head(ntms_VNM)

##   partner  NTM   hs6 nbr year ntm_1_digit
## 1     AUS A140 10121   1 2004           A
## 2     MYS A140 10121   1 2004           A
## 3     ZAF A140 10121   1 2004           A
## 4     LAO A140 10121   1 2004           A
## 5     NZL A140 10121   1 2004           A
## 6     THA A140 10121   1 2004           A

ntms_VNM <- ntms_VNM[ntms_VNM$partner !="VNM",]

ntms_VNM <- ntms_VNM[ntms_VNM$ntm_1_digit !="P",]
head(ntms_VNM)

##   partner  NTM   hs6 nbr year ntm_1_digit
## 1     AUS A140 10121   1 2004           A
## 2     MYS A140 10121   1 2004           A
## 3     ZAF A140 10121   1 2004           A
## 4     LAO A140 10121   1 2004           A
## 5     NZL A140 10121   1 2004           A
## 6     THA A140 10121   1 2004           A

ntms_VNM <- ntms_VNM %>% group_by(NTM, hs6, partner) %>% count() %>%
  select(NTM, hs6, partner)
head(ntms_VNM)

## # A tibble: 6 x 3
## # Groups:   NTM, hs6, partner [6]
##   NTM     hs6 partner
##   <chr> <dbl> <chr>
## 1 A140  10121 AGO
## 2 A140  10121 AND
## 3 A140  10121 ARE
## 4 A140  10121 ARG
## 5 A140  10121 AUS
## 6 A140  10121 AUT

7.3 Подготовка объекта с данными по импорту товаров

Применяем функцию subset к объекту "trade_data_VNM" для создания новой выборки данных в объекте "VNM_import", в которой значение в столбце "flow" равно 5 - только импорт. Сохраняем в новом объекте только те строки, которые содержат уникальные коды стран из столбца "partners" объекта "ntms_VNM". Проверяем количество строк.

VNM_import <- subset(trade_data_VNM,
                     trade_data_VNM$flow==5&
                     trade_data_VNM$partner %in% unique(ntms_VNM$partner),
                     select = c(partner, hs6, year, value))
nrow(VNM_import)

## [1] 208941

head(VNM_import)

##   partner   hs6 year    value
## 1     AUS 10129 2016  211.878
## 2     AUS 10129 2017  681.523
## 3     MYS 10129 2017   49.167
## 7     ZAF 10129 2015  249.600
## 8     ZAF 10129 2016  688.116
## 9     AUS 10221 2015 8301.276

Этот объект мы используем для расчета среднего значения стоимостного объема импорта для каждой комбинации страны-партнера, кода товара и года. Используем функцию summarize(value = mean(value)), для которой указаны такие аргументы, которые позволяют рассчитать среднеарифметическое значение для значений переменной value, относящихся к одинаковым комбинациям страны-партнера, кода ГС и года. Полученные результаты сохраняются в столбец “value” переформатированной таблицы данных. Мы это делаем, чтобы избавиться от строк, содержащих совпадающие комбинации кодов стран, кодов ГС и лет, но отличающиеся по стоимостному объему импорта.

VNM_import <- VNM_import %>% group_by(partner, hs6, year) %>% summarize(value = mean(value))
nrow(VNM_import)

## [1] 208941

head(VNM_import)

## # A tibble: 6 x 4
## # Groups:   partner, hs6 [5]
##   partner    hs6  year   value
##   <chr>    <dbl> <int>   <dbl>
## 1 AGO      30749  2017  154.
## 2 AGO      50790  2017    7.84
## 3 AGO     120799  2017   12.6
## 4 AGO     150420  2016  207.
## 5 AGO     230120  2015  601
## 6 AGO     230120  2016 2308.

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

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

VNM_import_years <- unique(VNM_import$year)
head(VNM_import_years)

## [1] 2017 2016 2015

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

VNM_import <- spread(VNM_import, year, value)
VNM_import[is.na(VNM_import)] <- 0
head(VNM_import)

## # A tibble: 6 x 5
## # Groups:   partner, hs6 [6]
##   partner    hs6 `2015` `2016`   `2017`
##   <chr>    <dbl>  <dbl>  <dbl>    <dbl>
## 1 AGO      30749      0     0    154.
## 2 AGO      50790      0     0      7.84
## 3 AGO     120799      0     0     12.6
## 4 AGO     150420      0   207.     0
## 5 AGO     230120    601  2308. 13380.
## 6 AGO     230400      0     0    544.

Теперь у нас нет отсутствующих данных ни для каких пар "код ГС-год". Снова конвертируем в длинный формат переменные "year" и "value" с помощью функции gather. Функция gather принимает несколько столбцов и преобразует их в пары "ключ - значение", в результате чего "широкие" таблицы данных превращаются в "длинные"..

VNM_import <- gather(VNM_import, year, value, as.character(VNM_import_years), factor_key=FALSE)
head(VNM_import)

## # A tibble: 6 x 4
## # Groups:   partner, hs6 [6]
##   partner    hs6 year     value
##   <chr>    <dbl> <chr>    <dbl>
## 1 AGO      30749 2017    154.
## 2 AGO      50790 2017      7.84
## 3 AGO     120799 2017     12.6
## 4 AGO     150420 2017      0
## 5 AGO     230120 2017  13380.
## 6 AGO     230400 2017    544.

Теперь вычислим среднегодовое значение.

VNM_import <- VNM_import %>% group_by(partner, hs6) %>% summarize(value=mean(value))
head(VNM_import)

## # A tibble: 6 x 3
## # Groups:   partner [1]
##   partner    hs6   value
##   <chr>    <dbl>   <dbl>
## 1 AGO      30749   51.4
## 2 AGO      50790    2.61
## 3 AGO     120799    4.18
## 4 AGO     150420   69.1
## 5 AGO     230120 5429.
## 6 AGO     230400  181.

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

7.4 Рассчитаем коэффициент покрытия, индекс частоты и показатель распространенности

Рассчитаем теперь коэффициент покрытия, индекс частоты и показатель распространенности для Вьетнама.

Создадим объект с итоговой таблицей, в которую сохраним следующие данные: последний год, для которого в WITS имеются данные по импорту во Вьетнам, а также значения коэффициентов охвата, индексов частоты и показателей распространенности для всего импорта и по секторам. Импорт товаров мы разобьем на такие группы, как продукция сельского хозяйства, природные ресурсы и продукция обрабатывающей промышленности. Обозначение “agro” используем для продукции сельского хозяйства, “NR” - для природных ресурсов, и “manu” - для продукции обрабатывающей промышленности. Например, с помощью CR_agro мы обозначим коэффициент охвата нетарифными мерами импортируемой во Вьетнам продукции сельского хозяйства.

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

summary_data_VNM <- data.frame(matrix(vector(),1,13,
                                  dimnames = list(c(),
                                                  c("latest_year",
                                                    "CR",
                                                    "FI",
                                                    "PS",
                                                    "CR_agro",
                                                    "CR_NR",
                                                    "CR_manu",
                                                    "FI_agro",
                                                    "FI_NR",
                                                    "FI_manu",
                                                    "PS_agro",
                                                    "PS_NR",
                                                    "PS_manu"
                                                    )
                                  )),
                           stringsAsFactors=F)

summary_data_VNM[1,"latest_year"] <- max(VNM_import_years) 

VNM_ntm_import_CR <- data.frame(ntms_VNM) %>% select(hs6, partner)

VNM_ntm_import_CR <- VNM_ntm_import_CR %>% group_by(hs6, partner) %>% summarise(count = n()) %>% select(hs6, partner)

VNM_ntm_import_CR$NTM <- 1

head(VNM_ntm_import_CR)

## # A tibble: 6 x 3
## # Groups:   hs6 [1]
##     hs6 partner   NTM
##   <dbl> <chr>   <dbl>
## 1 10121 AGO         1
## 2 10121 AND         1
## 3 10121 ARE         1
## 4 10121 ARG         1
## 5 10121 AUS         1
## 6 10121 AUT         1

CR_ntm_data <- merge(VNM_import,VNM_ntm_import_CR, by=c("partner", "hs6"), all.x = TRUE )

CR_ntm_data$NTM[is.na(CR_ntm_data$NTM)] <- 0

head(CR_ntm_data)

##   partner    hs6       value NTM
## 1     AGO  30749   51.361000   1
## 2     AGO  50790    2.613333   0
## 3     AGO 120799    4.183333   1
## 4     AGO 150420   69.090000   1
## 5     AGO 230120 5429.471667   1
## 6     AGO 230400  181.450000   1

CR <- sum(CR_ntm_data$NTM*CR_ntm_data$value)/sum(CR_ntm_data$value)*100
CR

## [1] 62.06449

summary_data_VNM[1, "CR"] <- CR

FI_ntm_data <- CR_ntm_data
FI_ntm_data$valueDummy <- 0
FI_ntm_data[!is.na(FI_ntm_data$value),"valueDummy"] <- 1
head(FI_ntm_data)

##   partner    hs6       value NTM valueDummy
## 1     AGO  30749   51.361000   1          1
## 2     AGO  50790    2.613333   0          1
## 3     AGO 120799    4.183333   1          1
## 4     AGO 150420   69.090000   1          1
## 5     AGO 230120 5429.471667   1          1
## 6     AGO 230400  181.450000   1          1

FI <- sum(FI_ntm_data$NTM*FI_ntm_data$valueDummy)/sum(FI_ntm_data$valueDummy)*100
FI

## [1] 50.81422

summary_data_VNM[1, "FI"] <- FI

VNM_ntm_PS <- ntms_VNM %>% group_by(hs6, partner) %>% summarise(count = n())
head(VNM_ntm_PS)

## # A tibble: 6 x 3
## # Groups:   hs6 [1]
##     hs6 partner count
##   <dbl> <chr>   <int>
## 1 10121 AGO         5
## 2 10121 AND         5
## 3 10121 ARE         5
## 4 10121 ARG         5
## 5 10121 AUS         5
## 6 10121 AUT         5

PS_ntm_data <- merge(VNM_ntm_PS, VNM_import, by=c("hs6", "partner"), all.y = TRUE)

head(PS_ntm_data)

##     hs6 partner count      value
## 1 10129     AUS    NA  297.80033
## 2 10129     MYS    NA   16.38900
## 3 10129     ZAF    NA  312.57200
## 4 10221     AUS     5 5002.01867
## 5 10221     NZL     5   17.00333
## 6 10221     THA     5  393.93000

PS_ntm_data$count[is.na(PS_ntm_data$count)] <- 0

PS_ntm_data$value[!(is.na(PS_ntm_data$value))&(PS_ntm_data$value>0)] <- 1
PS_ntm_data$value[is.na(PS_ntm_data$value)|PS_ntm_data$value==0] <- 0
head(PS_ntm_data)

##     hs6 partner count value
## 1 10129     AUS     0     1
## 2 10129     MYS     0     1
## 3 10129     ZAF     0     1
## 4 10221     AUS     5     1
## 5 10221     NZL     5     1
## 6 10221     THA     5     1

PS <- sum(PS_ntm_data$value*PS_ntm_data$count)/sum(PS_ntm_data$value)
PS

## [1] 1.86705

summary_data_VNM[1, "PS"] <- PS

summary_data_VNM[1, "CR_agro"] <- sum(CR_ntm_data$NTM[CR_ntm_data$hs6<250000]*CR_ntm_data$value[CR_ntm_data$hs6<250000])/sum(CR_ntm_data$value[CR_ntm_data$hs6<250000])*100
summary_data_VNM[1, "CR_NR"] <- sum(CR_ntm_data$NTM[CR_ntm_data$hs6>=250000&CR_ntm_data$hs6<280000]*CR_ntm_data$value[CR_ntm_data$hs6>=250000&CR_ntm_data$hs6<280000])/sum(CR_ntm_data$value[CR_ntm_data$hs6>=250000&CR_ntm_data$hs6<280000])*100
summary_data_VNM[1, "CR_manu"] <- sum(CR_ntm_data$NTM[CR_ntm_data$hs6>=280000]*CR_ntm_data$value[CR_ntm_data$hs6>=280000])/sum(CR_ntm_data$value[CR_ntm_data$hs6>=280000])*100

summary_data_VNM[1, "FI_agro"] <- sum(FI_ntm_data$NTM[FI_ntm_data$hs6<250000]*FI_ntm_data$valueDummy[FI_ntm_data$hs6<250000])/sum(FI_ntm_data$valueDummy[FI_ntm_data$hs6<250000])*100
summary_data_VNM[1, "FI_NR"] <- sum(FI_ntm_data$NTM[FI_ntm_data$hs6>=250000&FI_ntm_data$hs6<280000]*FI_ntm_data$valueDummy[FI_ntm_data$hs6>=250000&FI_ntm_data$hs6<280000])/sum(FI_ntm_data$valueDummy[FI_ntm_data$hs6>=250000&FI_ntm_data$hs6<280000])*100
summary_data_VNM[1, "FI_manu"] <- sum(FI_ntm_data$NTM[FI_ntm_data$hs6>=280000]*FI_ntm_data$valueDummy[FI_ntm_data$hs6>=280000])/sum(FI_ntm_data$valueDummy[FI_ntm_data$hs6>=280000])*100

rm(CR, FI) # удаляем объекты "CR" и "FI"
write.csv(summary_data_VNM, "summary_data_VNM.csv", row.names = FALSE)

# В качестве основы возьмем код выше, но воспользуемся формулой для расчета показателя распространенности
summary_data_VNM[1, "PS_agro"] <- sum(PS_ntm_data$value[PS_ntm_data$hs6<250000]*PS_ntm_data$count[PS_ntm_data$hs6<250000])/sum(PS_ntm_data$value[PS_ntm_data$hs6<250000])
summary_data_VNM[1, "PS_NR"] <- sum(PS_ntm_data$value[PS_ntm_data$hs6>=250000&PS_ntm_data$hs6<280000]*PS_ntm_data$count[PS_ntm_data$hs6>=250000&PS_ntm_data$hs6<280000])/sum(PS_ntm_data$value[PS_ntm_data$hs6>=250000&PS_ntm_data$hs6<280000])
summary_data_VNM[1, "PS_manu"] <- sum(PS_ntm_data$value[PS_ntm_data$hs6>=280000]*PS_ntm_data$count[PS_ntm_data$hs6>=280000])/sum(PS_ntm_data$value[PS_ntm_data$hs6>=280000])

# Значит Ответ D верный.

8.1 Расчёт значений индекса расхождения

В данном разделе мы познакомимся с методикой расчёта значения индекса расхождения НТМ для разных пар стран. Индекс расхождения - это индикатор, который измеряет степень расхождения в применении НТМ между различными государствами. С помощью этого индикатора можно увидеть, какие страны участвуют в одних и тех же региональных торговых соглашения (РТС). Чтобы получить значение индекса расхождения для пары стран, необходимо суммировать все случаи, где для двух стран существуют одинаковые комбинации кодов ГС и кодов НТМ, и разделить полученное число на количество всех имеющихся уникальных комбинаций кодов ГС и кодов НТМ:

В формуле выше \(N\) - это количество всех уникальных сочетаний кодов ГС и кодов НТМ, \(n_{ilk}\) - фиктивная переменная, которая показывает, применяет ли страна \(i\) нетарифную меру \(l\) к товару \(k\). \(n_{jlk}\) также является фиктивной переменной, которая показывает, применяет ли страна \(j\) нетарифную меру \(l\) к товару \(k\). (Более подробную информацию об индексе расхождения НТМ можно найти в этой публикации ЮНКТАД.)

В данном разделе мы построим матрицу значений индикатора расхождения для разных пар государств, которые являются членами ЭСКАТО. Для этого выполним следующие шаги:

Шаг 1: Создадим список государств-членов ЭСКАТО, для которых имеются данные по НТМ.
Шаг 2: Создадим пустую матрицу, использовав в качестве названий столбцов и строк код ISO3 для государств-членов ЭСКАТО.
Шаг 3: Рассчитаем \(N\) (количество всех уникальных комбинаций кодов ГС и кодов НТМ в регионе ЭСКАТО).
Шаг 4: Создадим цикл cо вторым вложенным for-циклом, который позволит вычислить значения индекса расхождения для каждой пары стран.

Но перед этим загрузим необходимые пакеты функций. Если вы еще не установили некоторые пакеты функций, выполните соответствующий код, убрав знак решетки.

#install.packages("dplyr")
library("dplyr")
#install.packages("readstata13")
library("readstata13")

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

Загрузим в R список всех товаров, к которым применяются НТМ. Напомним, что это займет несколько минут.

ntm_data <- read.dta13("ntm_hs6_2016 v.12.dta")

После загрузки данных в объект R исключим из него строки, которые в столбце "ntmcode" не содержат отсутствующие значения. Нас интересуют только три переменные: "reporter", "ntmcode" и "hs6".

ntm_data <- ntm_data[!is.na(ntm_data$ntmcode)&ntm_data$ntmcode!="",]
ntm_data <- ntm_data[,c("reporter", "ntmcode", "hs6")]

Сгруппируем данные по переменным "reporter", "ntmcode" и "hs6" . С помощью функции "count" посчитаем количество уникальных комбинаций и сохраним результат в столбец count.

ntm_data <- ntm_data %>% group_by(reporter, ntmcode, hs6) %>%
  summarise(count = n())
head(ntm_data)

## # A tibble: 6 x 4
## # Groups:   reporter, ntmcode [1]
##   reporter ntmcode hs6    count
##   <chr>    <chr>   <chr>  <int>
## 1 AFG      A110    010110     1
## 2 AFG      A110    010190     1
## 3 AFG      A110    010210     1
## 4 AFG      A110    010290     2
## 5 AFG      A110    010410     1
## 6 AFG      A110    010420     1

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

Мы выполним анализ только для государств-членов ЭСКАТО, так как выполнение кода может занять много времени. Если вы хотите выполнить анализ для всех стран мира, для которых имеются данные по НТМ, мы рекомендуем использовать код avail_iso3s <- unique(ntm_data$reporter).

Загрузим файл с кодами стран по ISO3 и извлечем только коды государств-членов ЭСКАТО.

countries <- read.csv("Country Categories.csv", stringsAsFactors = FALSE)

ESCAP <- countries$iso3[countries$escap==1]

Затем извлечем все уникальные коды стран из столбца "reporter" таблицы с данными по НТМ. Таким образом мы получим список все стран, для которых имеются данные по НТМ.

avail_iso3s <- unique(ntm_data$reporter)
head(avail_iso3s)

## [1] "AFG" "ARE" "ARG" "ATG" "AUS" "BEN"

Используем функцию intersect, чтобы извлечь только те коды стран, которые содержатся в объектах avail_iso3s и ESCAP одновременно. Обратите внимание, что в последнем случае мы получили список кодов государств-членов ЭСКАТО.

avail_iso3s <- intersect(avail_iso3s,ESCAP)
head(avail_iso3s)

## [1] "AFG" "AUS" "BRN" "CHN" "HKG" "IDN"

Создадим пустую матрицу значений индекса расхождения. Зададим аргумент length(avail_iso3s)+1, чтобы задать число столбцов, равное количеству элементов в объекте "avail_iso3s", плюс один столбец для столбца "reporter". Теперь зададим названия столбцов с помощью аргумента dimnames, а также переменной reporter и кодов стран ISO3 .

regulatory_distance_matrix <- data.frame(matrix(vector(),0,length(avail_iso3s)+1,
                                   dimnames = list(c(), c("reporter", avail_iso3s )
                                   )),
                            stringsAsFactors=F)

Теперь рассчитаем значения индексов расхождения в соответствии с формулой:

Так как переменная \(N\) является константой, сначала рассчитаем её. Как вы помните, \(N\) - это количество всех уникальных комбинаций кодов НТМ и кодов ГС во всем наборе данных. Значит, мы можем использовать функцию group_by. Создадим объект N, в который сохраним это значение.

N <- ntm_data %>% group_by(ntmcode, hs6) %>% count()
N <- nrow(N)

Теперь мы можем заполнить матрицу значений индекса расхождения. Для этого создадим цикл, который выполнит итерации кода для каждого государства-члена ЭСКАТО и сохранит полученные значения индексов в итоговой матрице regulatory_distance_matrix, которую мы создали ранее. Внимательно изучим код, содержащий вложенный цикл.

Сначала извлечем все комбинации кодов НТМ (ntmcode) и кодов ГС (hs6) для страны из столбца "reporter", обрабатываемой в итерации g (далее "страна А"): country_a <- ntm_data[ntm_data$reporter==avail_iso3s[g],c("ntmcode", "hs6")]. Создадим дополнительный столбец, в котором сохраним значение 1 для всех этих комбинаций : country_a$country_a <- 1. Сохраним код ISO3 страны А в ячейку столбца "reporter" с помощью кода regulatory_distance_matrix[g,"reporter"] <- avail_iso3s[g].

Далее введем вложенный цикл, так как нам нужно пройти через все комбинации страны А со всеми другими странами. Используем вложенный цикл, чтобы извлечь все комбинации кодов НТМ (ntmcode) и кодов ГС (hs6) страны-партнёра (далее "страна В"): country_b <- ntm_data[ntm_data$reporter==avail_iso3s[k],c("ntmcode", "hs6")]. Создадим дополнительный столбец, в котором сохраним значение 1 для всех этих комбинаций : country_b$country_b <- 1. Объединим два новых объекта с данными по странам А и B, сохранив все строки каждого из объектов: merged <- merge(country_a, country_b, by=c("ntmcode", "hs6"), all = TRUE). В объединенном объекте появились отсутствующие значения для тех комбинаций кодов НТМ и ГС, которые присутствовали в одном из объединяемых объектов, но отсутствовали в другом. Заменим отсутствующие значение на 0: merged[is.na(merged)] <- 0. Создадим переменную (в новом столбце), которая отражает абсолютную разницу между значениями бинарных переменных обеих стран: merged$abs_diff <- abs(merged$country_a-merged$country_b). Суммируем все значения в новом столбце и разделим результат на число в объекте N (количество всех уникальных комбинаций кодов НТМ и ГС в регионе): rd <- sum(merged$abs_diff)/N. Таким образом мы получаем значение индекса расхождения rd, которое мы сохраняем в матрицу: regulatory_distance_matrix[g,avail_iso3s[k]] <- rd.

Чтобы повысить эффективность цикла, мы включаем в него выражение с условным оператором if, которое позволяет не вычислять значение показателя для пары стран А и В, если значение для стран В и А уже было получено в предыдущей итерации цикла: if (!is.na(regulatory_distance_matrix[k,avail_iso3s[g]])){next }. (В результате мы получаем треугольную матрицу.)

Когда будете готовы, выполните код ниже. Это может занять несколько минут.

for (g in 1:length(avail_iso3s)){
  country_a <- ntm_data[ntm_data$reporter==avail_iso3s[g],c("ntmcode", "hs6")]
  country_a$country_a <- 1
  regulatory_distance_matrix[g,"reporter"] <- avail_iso3s[g]

  for (k in 1:length(avail_iso3s)){

     if (!is.na(regulatory_distance_matrix[k,avail_iso3s[g]])){next }

  country_b <- ntm_data[ntm_data$reporter==avail_iso3s[k],c("ntmcode", "hs6")]
  country_b$country_b <- 1
  merged <- merge(country_a, country_b, by=c("ntmcode", "hs6"), all = TRUE)
  merged[is.na(merged)] <- 0
  merged$abs_diff <- abs(merged$country_a-merged$country_b)
  rd <- sum(merged$abs_diff)/N
  regulatory_distance_matrix[g,avail_iso3s[k]] <- rd

   }
}

Теперь заполним отсутствующие значения и создадим файл .csv. Например, если значение индекса для пары стран A и B отсутствует в матрице (условное выражение if (is.na(regulatory_distance_matrix[k,avail_iso3s[g]]))), мы используем значение индекса, рассчитанное для пары стран B и A (код regulatory_distance_matrix[k,avail_iso3s[g]] = regulatory_distance_matrix[g,avail_iso3s[k]]).

for (g in 1:length(avail_iso3s)){
  for (k in 1:length(avail_iso3s)){
    if (is.na(regulatory_distance_matrix[k,avail_iso3s[g]])){
      regulatory_distance_matrix[k,avail_iso3s[g]] <- regulatory_distance_matrix[g,avail_iso3s[k]]
      }
  }
}

write.csv(regulatory_distance_matrix, "regulatory_distance_matrix_ESCAP.csv")

# Вам нужно выполнить код выше и изменить путь в аргументе функции "write.csv". Откройте сохраненный файл и найдите значение индекса для KAZ и CHN. Значение индекса расхождения для этой пары стран 0,32280137848764 или 0,32, а значит, ответ В верный. Полностью заполненную матрицу можно скачать здесь: https://r.tiid.org/regulatory_distance_matrix_ESCAP.csv

8.2 Многомерное шкалирование

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

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

#install.packages("magrittr")
library(magrittr)
#install.packages("dplyr")
library(dplyr)
#install.packages("ggpubr")
library(ggpubr)

Сначала нужно соответствующим образом подготовить матрицу значений индекса расхождения. Создадим копию матрицы значений индекса расхождения НТМ в объекте regulatory_distance_matrix_MDS. Затем с помощью данных, содержащихся в первом столбце, зададим названия строк таблицы данных.

regulatory_distance_matrix_MDS <- regulatory_distance_matrix
rownames(regulatory_distance_matrix_MDS) <- regulatory_distance_matrix_MDS$reporter
regulatory_distance_matrix_MDS <- regulatory_distance_matrix_MDS %>% select (-c(reporter))

Теперь выполним анализ главных координат с помощью функции cmdscale и сохраним полученные значения в таблицу. Созданные в ней столбцы назовем Dim.1 и Dim.2.

regulatory_distance_matrix_MDS <- regulatory_distance_matrix_MDS %>% cmdscale() %>% as_tibble()

colnames(regulatory_distance_matrix_MDS) <- c("Dim.1", "Dim.2")

head(regulatory_distance_matrix_MDS)

## # A tibble: 6 x 2
##     Dim.1    Dim.2
##     <dbl>    <dbl>
## 1  0.0168  0.0158
## 2 -0.0163 -0.00568
## 3  0.0149  0.0111
## 4 -0.271  -0.0340
## 5  0.0159  0.0107
## 6  0.0123  0.0171

Пары значений в столбцах Dim.1 и Dim.2 отражают положение стран в двухмерной системе координат, которые размещены так, чтобы расстояние между ними примерно соответствовало значениям индекса расхождения между всеми парами стран. Значения в этих столбцах сами по себе не несут никакой информации. Для нас важно только расстояние между парами стран.

Теперь мы можем построить график и визуализировать относительные расстояния между странами.

Воспользуемся функцией ggscatter для построения графика рассеяния. Добавим горизонтальную geom_hline и вертикальную линии geom_vline на наш график, чтобы показать горизонтальную и вертикальную оси соответственно.

scatterplot <- ggscatter(regulatory_distance_matrix_MDS, x = "Dim.1", y = "Dim.2",
          label = regulatory_distance_matrix$reporter,
          size = 1,
          repel = TRUE)
scatterplot + geom_hline(yintercept=0, linetype="dashed", color = "red") +
  geom_vline(xintercept = 0, linetype="dotted",
                color = "blue")

Из графика видно, что большинство стран АСЕАН расположены ближе друг к другу (Мьянма, Борнео, Малайзия, Таиланд и Сингапур). Россия, Казахстан и Кыргызстан тоже расположены достаточно близко друг другу, будучи членами Евразийского экономического сообщества (ЕАЭС). Также видно, что Китай расположен далеко от всех стран. Это ожидаемо, так как применяемое Китаем количество НТМ одно из самых высоких в мире.