“Эконометрика. Домашнее задание №1”

Илья Тищенко, 12/06/2021

Введение, вопрос

Поскольку данные RLMS временно надоели, то захотелось посмотреть в другую сторону. В данном исследовании я решил посмотреть на зависимость между уровнем рождаемости и другими показателями развития стран (ВВП, ожидаемая продолжительность жизни женщин, доступность питьевой воды, рента от природных ресурсов в доле ВВП, урбанизация). Данные были взяты из world bank open data.

Обзор данных

В качестве объясняющих переменных я использовал следующие:

• gdp ВВП в ценах 2015 выраженный в US$, конвертированный из местной валюты по курсу 2015
  

• trade сумма импорта и экспорта товаров деленная на ВВП, в процентах
  

• urban_ratio процент городского населения. Эту переменную я превратил в качественную, при значении больше или равно 70% - “high”, меньше 70% - “low”
  

• water процент людей имеющих доступ к питьевой воде
  

• resour общая рента от природных ресурсов, в процентах от ВВП
  

• life_exp_female ожидаемая продолжительность жизни для женщин
  

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

Объясняемая переменная:

birth_rate количество рождений на 1000 человек

Посмотрим на данные. Средний по миру уровень рождаемости 20, минимальный 7 максимальный 47.

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

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

Данные ВВП и доля импорта, экспорта имеют несколько очень далеко отстоящих значений (см. США, Китай, Япония, Сингапур, Гонконг, Джибути)

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

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

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

Также можно графически, а потом и точно оценить сколько точек имеется для каждого региона.

На этом графике видны следующие особенности:

• в Европе зависимость между уровнем рождаемости и ожидаемой продолжительностью жизни женщин практически остутствует

• в Америке есть одна точка выброс (Гаити), хотя скорее всего это не нарушит корреляцию

• точек для Океании довольно мало.  

Ну вообщем тут есть над чем работать. 

Можно это немного формализовать и подсчитать корреляцию между этими двумя параметрами (ожидаемая продолжительность жизни и уровнь рождаемости) по региону. Как и следовало ожидать в Европе эти два параметра практически не коррелируют. На графике ниже показана корреляция по регионам.

Мультиколлинеарность

Мультиколлинеарность между некоторыми параметрами можно увидеть на кроссплот-матрице ниже. Например параметр доступ к воде (water) коррелирует с параметром ожидаемая продолжительность жизни женщин(life_exp_female) - корреляция 0.832 (без разделения по урбанизации).  

Главные компоненты

Разложение регрессоров на главные компоненты показывает, что переменные (water, life_exp_female) - доступность к питьевой воде и ожидаемая продолжительность жизни женщин, вместе с общей рентой от природных ресурсов описывают 44.5% дисперсии в данных. ВВП и торговая доля в ВВП отвечают за оставшиеся 22.4%, таким образом первые две главные компоненты описывают 66.9% дисперсии исходных данных.

Также на графике ниже видны страны, которые сильно отличаются по показателю ВВП (США, Китай, Япония) и по торговой доле от ВВП (Гонконг, Сингапур, Джибути).  

Модели на всех данных

Сначала я построил модель на всех данных, ниже показано сравнение модели с обычными предпосылками (слева) и робастными(справа) стандартными ошибками.  

Сравнение простых моделей

	simple		simple robust
(Intercept)	79.732 ***	(4.247)	79.732 ***	(3.867)
water	-0.498 ***	(0.025)	-0.498 ***	(0.031)
log(gdp)	-0.640 ***	(0.182)	-0.640 ***	(0.164)
N. obs.	181		181
R squared	0.738		0.738
F statistic	250.083		201.074
P value	0.000		0.000
*** p < 0.001; ** p < 0.01; * p < 0.05.

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

Сравнение сложных моделей

	complex		complex robust
(Intercept)	83.09 ***	(71.13 : 95.05)	83.09 ***	(69.48 : 96.70)
water	-0.18 ***	(-0.25 : -0.12)	-0.18 ***	(-0.27 : -0.09)
resour	0.11 **	(0.04 : 0.19)	0.11 *	(0.03 : 0.20)
regionAmericas	-4.24 ***	(-6.39 : -2.10)	-4.24 **	(-6.95 : -1.53)
regionAsia	-3.45 ***	(-5.39 : -1.50)	-3.45 **	(-5.90 : -1.00)
regionEurope	-6.62 ***	(-9.01 : -4.23)	-6.62 ***	(-9.37 : -3.86)
regionOceania	-1.15	(-3.82 : 1.52)	-1.15	(-5.17 : 2.87)
log(gdp)	-0.12	(-0.43 : 0.19)	-0.12	(-0.45 : 0.20)
log(trade)	-1.47 *	(-2.60 : -0.35)	-1.47 **	(-2.55 : -0.40)
life_exp_female	-0.48 ***	(-0.63 : -0.33)	-0.48 ***	(-0.65 : -0.30)
urban_qlow	0.55	(-0.87 : 1.98)	0.55	(-0.91 : 2.02)
N. obs.	181		181
R squared	0.87		0.87
F statistic	118.30		126.64
P value	0.00		0.00
*** p < 0.001; ** p < 0.01; * p < 0.05.

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

Сравнение сложных моделей

	complex		complex robust
(Intercept)	82.81 ***	(74.44 : 91.18)	82.81 ***	(73.45 : 92.17)
water	-0.18 ***	(-0.25 : -0.12)	-0.18 ***	(-0.27 : -0.09)
resour	0.10 **	(0.03 : 0.17)	0.10 *	(0.02 : 0.19)
regionAmericas	-4.06 ***	(-6.19 : -1.94)	-4.06 **	(-6.67 : -1.45)
regionAsia	-3.42 ***	(-5.33 : -1.51)	-3.42 **	(-5.81 : -1.03)
regionEurope	-6.63 ***	(-8.99 : -4.27)	-6.63 ***	(-9.32 : -3.95)
regionOceania	-0.76	(-3.33 : 1.81)	-0.76	(-4.66 : 3.14)
log(trade)	-1.30 *	(-2.38 : -0.21)	-1.30 *	(-2.33 : -0.26)
life_exp_female	-0.52 ***	(-0.66 : -0.38)	-0.52 ***	(-0.68 : -0.36)
N. obs.	181		181
R squared	0.87		0.87
F statistic	148.00		155.35
P value	0.00		0.00
*** p < 0.001; ** p < 0.01; * p < 0.05.

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

	GVIF	Df	GVIF^(1/(2*Df))
water	3.59	1	1.895
resour	1.31	1	1.144
region	3.257	4	1.159
log(trade)	1.22	1	1.105
life_exp_female	4.301	1	2.074

Выбор финальной модели

Тесты  

С помощью F-теста можно посмотреть, действительно ли более сложная модель лучше. Гипотеза о том что более простая модель верна отвергается на 5% уровне.

Analysis of Variance Table

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
178	4636	NA	NA	NA	NA
172	2241	6	2396	30.65	7.232e-25

Тест Рамсея показывает, что в сложной модели есть пропущенные переменные, нулевая гипотеза что коэффициенты \(\gamma\) при вспомогательных регрессорах равны нулю не отвергается. (\(p-value < \alpha=5\%\))

RESET test: fit.complex_hc

Test statistic	df1	df2	P value
5.964	2	170	0.003139 * *

Интерпретация коэффициентов

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

Для параметра:  

• ожидаемая продолжительность жизни женщин, при увеличении на один год ожидаемая рождаемость (\(Y\)) уменьшается на 0.52  
• при увеличении на один процент (переменная в процентах) доступа к воде, \(Y\) уменьшается на 0.18
• при увеличении на один процент (переменная в процентах) ренты от природных ресурсов, \(Y\) увеличивается на 0.1
• при изменении доли торговли от ВВП на один процент (т.к. переменная взята с логарифмом), \(Y\) уменьшается на 1.3  
• базовый регион в регрессии Aфрика  
• в Европе рождаемость ниже в среднем на 6.63  
• в Америках рождаемость ниже в среднем на 4.06  
• в Азии рождаемость ниже в среднем на 3.42  
• в Океании недостаточно данных (коэффициент не значим) чтобы регрессия отличалась от Африки  

Тестовая и обучающая выборка

Случайным образом я поделил данные на 70:30%. Первую часть (70% = 126 точек) использовал как обучающую, вторую (55 точек) как тестовую.

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

Сравнение сложных моделей

	simple robust		complex robust
(Intercept)	82.78 ***	(73.49 : 92.06)	81.24 ***	(69.20 : 93.27)
water	-0.50 ***	(-0.56 : -0.43)	-0.20 ***	(-0.29 : -0.11)
log(gdp)	-0.76 ***	(-1.13 : -0.39)
resour			0.13 **	(0.04 : 0.23)
regionAmericas			-2.99	(-6.03 : 0.05)
regionAsia			-2.72	(-5.58 : 0.14)
regionEurope			-5.93 ***	(-9.13 : -2.74)
regionOceania			1.95	(-1.62 : 5.52)
log(trade)			-0.62	(-1.85 : 0.61)
life_exp_female			-0.53 ***	(-0.73 : -0.33)
N. obs.	126		126
R squared	0.74		0.88
F statistic	149.93		108.10
P value	0.00		0.00
*** p < 0.001; ** p < 0.01; * p < 0.05.

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

region	Количество наблюдений
Africa	35
Americas	24
Asia	35
Europe	25
Oceania	7

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

LASSO

Для поиска оптимального значения лямбды для LASSO регрессии я использовал кросс-валидацию и построил стандартный график по оси Х - логарифм лямбда, по оси Y - средне квадратичная ошибка.

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

lambda.1se	lambda.min
0.843	0.0904

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

term	estimate
(Intercept)	79.7
water	-0.186
resour	0.0613
Europe_q	-2.45
Africa_q	0.913
life_exp_female	-0.578

Прогноз моделей на тестовую часть

Теперь можно подсчитать прогноз на тестовую часть по всем трем моделям. В качестве оценки моделей, я подсчитал средне квадратичную ошибку и среднюю абсолютную ошибку в процентах (MAPE) по формулам: \[ MSE=\frac{1}{n} \displaystyle\sum_{i=1}^{n}(Y_i-\hat{Y_i})^2 \]

\[ MAPE=\frac{100\%}{n} \displaystyle\sum_{n-1}^{n}\left|\frac{Y_i-\hat{Y_i}}{Y_i}\right| \]

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

RMSE lasso	RMSE simple	RMSE complex	MAPE lasso	MAPE simple	MAPE complex
4.02	5.46	4.04	16.3	22.8	16.2

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

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

Прогноз для финальной модели

Поскольку в задании нам нужно построить прогноз с предиктивным интервалом то пришлось взять модель МНК. Для прогноза я взял данные по 2018 году, все предсказания до этого строились по 2017. То есть это данные которые модель “не видела”. Прогноз подсчитан с предиктивным интервалом и показан на карте и кроссплоте ниже.

Гетероскедастичность

В качестве визуального теста на гетероскедастичность я отобразил остатки для этого прогноза. По оси Х прогноз, по оси Y остатки (прогноз-исходное значение). Серые линии нарисованы мной - они призваны показать, что некоторое неравномерное распределение дисперсии наблюдается, т.е. гетероскедастичность есть.  

В качестве формального теста на гетероскедастичность я подсчитал тест Бройш-Пагана. В данном случае нулевая гипотеза (о том что квадраты остатков не предсказываются регрессорами) отвергается на 5% уровне значимости (p-value=2.8 < 5%), т.е. в модели есть гетероскедастичность.  

## 
##  studentized Breusch-Pagan test
## 
## data:  fit.complex_hc
## BP = 17.182, df = 8, p-value = 0.02827

Прогноз на карте

Ну а поскольку я работал со странами, то мне стало интересно посмотреть на результаты на карте. Если навести курсор на карту, то показывается исходное значение (Birth Rate:), прогноз (Fit:), отклонение прогноза от истины в процентах (Resid:) и предиктивный интервал (lwr:up).  

Но если честно, карта красивая, время я на неё потратил, но каких-то прям очевидных закономерностей не увидел :) Ну разве что при такой разбивке на бины, в Африке и Западной Европе (за исключением Италии и Финляндии) моя модель в основном ошибается в меньшую сторону. У меня, всё-таки довольно широкий предиктивный интервал.  

Но если нужен более традиционный прогноз, то вот кроссплот ниже. Серыми линиями показан предиктивный интервал, синими - доверительный, красная линия y=x.  

Скрипт

## ----setup, include=FALSE-------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = FALSE,message=FALSE,warning = FALSE, fig.align = 'left')
library(readr)
library(tidyverse)
library(skimr)
library(corrplot)
library(sandwich)
library(lmtest)
library(estimatr)
library(glmnet) 
library(tidymodels)
library(car) 
library(factoextra)
library(huxtable)
library(cowplot)
library(plotly)
library(sjPlot) 
library(GGally)
library(vcd)
library(DT)
library(sf)
library(leaflet)
library(ggrepel)
library(pander)
library(ggridges)



## ----data input, cache=TRUE-----------------------------------------------------------------------------------
world <- read.csv("World_dev_indicators.csv",check.names = FALSE)
code_region <- read.csv("Data/country_region.txt")
code_region_sel <- code_region %>% 
  select(name,alpha.3,region,sub.region) 
#filter data
world_flt <-  world %>% 
  filter(`Series Name` %in% c("People using at least basic drinking water services (% of population)",
                              "Access to electricity (% of population)",
                              "Urban population (% of total population)",
                              "Birth rate, crude (per 1,000 people)",
                              "GDP (constant 2015 US$)",
                              "Merchandise trade (% of GDP)",
                              "Total natural resources rents (% of GDP)",
                              "Life expectancy at birth, female (years)")) %>% 
  # View()
select(-c(`Series Code`,`2018 [YR2018]`)) %>%
    # View()
pivot_wider(id_cols = c(1:2),names_from = `Series Name`,values_from = `2017 [YR2017]`) %>% 
  # View()
      rename(birth_rate="Birth rate, crude (per 1,000 people)",
                urban_ratio="Urban population (% of total population)",
                electr_acc="Access to electricity (% of population)",
                gdp="GDP (constant 2015 US$)",
                water="People using at least basic drinking water services (% of population)",
             trade="Merchandise trade (% of GDP)",
             resour="Total natural resources rents (% of GDP)",
             life_exp_female="Life expectancy at birth, female (years)") %>%
  slice_head(n=217) %>% 
#create quantiative variable from urban
  mutate(urban_q=case_when(
    urban_ratio < 70 ~ "low",
    urban_ratio >= 70 ~ "high")) %>%
  mutate(urban_q = factor(urban_q)) 

#data for regression
df <- world_flt%>% 
  left_join(code_region_sel,by=c("Country Code"="alpha.3")) %>% 
  select(`Country Name`,birth_rate,gdp, trade, urban_q, electr_acc, water,resour, region, sub.region,life_exp_female) %>% 
  mutate(region=factor(region),sub.region=factor(sub.region)) %>% 
  na.omit() %>% 
  data.frame()



## ----import 2018 data-----------------------------------------------------------------------------------------
world_flt_2018 <-  world %>% 
  filter(`Series Name` %in% c("People using at least basic drinking water services (% of population)",
                              "Access to electricity (% of population)",
                              "Urban population (% of total population)",
                              "Birth rate, crude (per 1,000 people)",
                              "GDP (constant 2015 US$)",
                              "Merchandise trade (% of GDP)",
                              "Total natural resources rents (% of GDP)",
                              "Life expectancy at birth, female (years)")) %>% 
  # View()
select(-c(`Series Code`,`2017 [YR2017]`)) %>%
    # View()
pivot_wider(id_cols = c(1:2),names_from = `Series Name`,values_from = `2018 [YR2018]`) %>% 
  # View()
      rename(birth_rate="Birth rate, crude (per 1,000 people)",
                urban_ratio="Urban population (% of total population)",
                electr_acc="Access to electricity (% of population)",
                gdp="GDP (constant 2015 US$)",
                water="People using at least basic drinking water services (% of population)",
             trade="Merchandise trade (% of GDP)",
             resour="Total natural resources rents (% of GDP)",
             life_exp_female="Life expectancy at birth, female (years)") %>%
  slice_head(n=217) %>% 
#create quantiative variable from urban
  mutate(urban_q=case_when(
    urban_ratio < 70 ~ "low",
    urban_ratio >= 70 ~ "high")) %>%
  mutate(urban_q = factor(urban_q)) 

#data for regression
df_2018 <- world_flt_2018%>% 
  left_join(code_region_sel,by=c("Country Code"="alpha.3")) %>% 
  select(`Country Name`,`Country Code`,birth_rate,gdp, trade, urban_q, electr_acc, water,resour, region, sub.region,life_exp_female) %>% 
  mutate(region=factor(region),sub.region=factor(sub.region)) %>% 
  na.omit() %>% 
  data.frame()

df_2018_dummy <- df_2018 %>% 
  mutate(Europe_q=ifelse(region == "Europe",1,0),
         Asia_q=ifelse(region == "Asia",1,0),
         Africa_q=ifelse(region == "Africa",1,0),
         Americas_q=ifelse(region == "Americas",1,0),
         urban_q_low=ifelse(urban_q==1,1,0)) 



## -------------------------------------------------------------------------------------------------------------
# GDP at purchaser's prices is the sum of gross value added by all resident producers in the economy plus any product taxes and minus any subsidies not included in the value of the products. It is calculated without making deductions for depreciation of fabricated assets or for depletion and degradation of natural resources. Data are in constant 2015 prices, expressed in U.S. dollars. Dollar figures for GDP are converted from domestic currencies using 2015 official exchange rates. For a few countries where the official exchange rate does not reflect the rate effectively applied to actual foreign exchange transactions, an alternative conversion factor is used.

# Merchandise trade as a share of GDP is the sum of merchandise exports and imports divided by the value of GDP, all in current U.S. dollars.

# Urban population refers to people living in urban areas as defined by national statistical offices. The data are collected and smoothed by United Nations Population Division.


# Access to electricity is the percentage of population with access to electricity. Electrification data are collected from industry, national surveys and international sources.

# The percentage of people using at least basic water services.  This indicator encompasses both people using basic water services as well as those using safely managed water services.  Basic drinking water services is defined as drinking water from an improved source, provided collection time is not more than 30 minutes for a round trip.  Improved water sources include piped water, boreholes or tubewells, protected dug wells, protected springs, and packaged or delivered water.

# Total natural resources rents are the sum of oil rents, natural gas rents, coal rents (hard and soft), mineral rents, and forest rents.

# Crude birth rate indicates the number of live births occurring during the year, per 1,000 population estimated at midyear. Subtracting the crude death rate from the crude birth rate provides the rate of natural increase, which is equal to the rate of population change in the absence of migration.




## ----skim table numeric---------------------------------------------------------------------------------------
df %>% 
  skim_without_charts() %>% 
  select(-c(n_missing,complete_rate)) %>% 
  yank("numeric") %>% 
  # mutate_if(is.numeric,round,2) %>% 
  datatable(rownames = FALSE,options = list(pageLength=10),autoHideNavigation = TRUE) %>% 
  formatSignif(digits=2,interval=3,columns=c("mean","sd","p50","p0","p25","p75","p100"))


## ----density plot---------------------------------------------------------------------------------------------
df %>% 
  ggplot(aes(x=birth_rate,fill=urban_q))+
  geom_density(alpha=0.4)+
  scale_y_continuous(expand=expansion(mult=c(0,0.05)))+
  scale_x_continuous(expand=expansion(mult=c(0,0.05)))+
  theme_minimal_hgrid()+
  labs(x="Уровень рождаемости", y="Плотность",fill="Урбанизация", subtitle="Распределение уровня рождаемости по урбанизации")


## ----density plot water---------------------------------------------------------------------------------------
df %>% 
  ggplot(aes(x=water,fill=urban_q))+
  geom_density(alpha=0.7,position = position_identity())+
  scale_y_continuous(expand=expansion(mult=c(0,0.05)))+
  scale_x_continuous(expand=expansion(mult=c(0,0.05)))+
  theme_minimal_hgrid()+
  labs(x="Доступ к питьевой воде, %", y="Плотность",fill="Урбанизация", subtitle="Доступ к питьевой воде в процентах")


## ----gdp, trade-----------------------------------------------------------------------------------------------
df %>% 
  select(Country.Name,urban_q,gdp,trade) %>% 
  pivot_longer(cols=3:4) %>% 
  ggplot(aes(x=urban_q,y=value,color=urban_q,fill=urban_q))+
  geom_jitter(width = 0.2, shape=21)+
  geom_text_repel(data=df %>%
                    select(Country.Name,urban_q,gdp,trade) %>% 
                    pivot_longer(cols=3:4) %>% 
                    filter(Country.Name %in% c("United States","China",
                                               "Hong Kong SAR, China",
                                               "Djibouti",
                                               "Japan")),
                  aes(label=Country.Name),
                  color="grey30",
                  nudge_x=0.5,
                  nudge_y=-2,
                  box.padding = 0.9,
                  segment.curvature=-0.1,
                  direction="both",
                  min.segment.length = 0,
                  max.overlaps = Inf)+
  geom_boxplot(alpha=0.4)+
  facet_wrap(.~name,scales = "free")+
  theme_cowplot()



## ----gdp, trade log-------------------------------------------------------------------------------------------
df %>% 
  select(urban_q,gdp,trade) %>% 
  mutate(gdp=log(gdp),trade=log(trade)) %>% 
  pivot_longer(cols=2:3) %>% 
  ggplot(aes(x=urban_q,y=value,fill=urban_q,color=urban_q))+
  geom_boxplot(alpha=0.4)+
  geom_jitter(width = 0.2)+
  facet_wrap(.~name,scales = "free")+
  theme_cowplot()+
  theme(legend.position = 'none')+
  labs(x="Урбанизация",y="Значение, логарифмическая шкала", color="",fill="")



## ----boxplot female life exp----------------------------------------------------------------------------------
df %>% 
  select(urban_q,life_exp_female) %>% 
  ggplot(aes(x=urban_q,y=life_exp_female,color=urban_q,fill=urban_q))+
  geom_boxplot(alpha=0.4)+
  geom_jitter(width = 0.2)+
  theme_cowplot()+
  labs(x="Урбанизация",y="Ожидаемая продолжительность жизни женщин", color="")+
  theme(legend.position = 'none')



## ----exploratory data analysis, fig.width= 9------------------------------------------------------------------

df %>% 
  select(`Country.Name`,birth_rate,urban_q,gdp,trade,water,resour,life_exp_female) %>%
  mutate(gdp=log(gdp),trade=log(trade)) %>%
  filter(!is.na(urban_q)) %>%
  pivot_longer(cols=c(4:8)) %>% 
  ggplot(aes(y=birth_rate,x=value,color=urban_q))+
  geom_point()+
  facet_wrap(.~name,scales = "free")+
  theme_cowplot()+
  labs(color="Урбанизация", y="Уровень рождаемости", x="Значение параметра")


## ----exploratory data analysis be region, fig.width= 9--------------------------------------------------------

df %>% 
  select(`Country.Name`,birth_rate,region,life_exp_female, urban_q) %>%
  ggplot(aes(x=birth_rate,y=life_exp_female,color=urban_q, text= `Country.Name`))+
  geom_point()+
  facet_wrap(.~region,scales = "free")+
  theme_cowplot()+
  labs(color="Урбанизация", x="Уровень рождаемости", y="Ожидаемая продолжительность жизни у женщин")+
  geom_text_repel(data=df %>%
                      select(`Country.Name`,birth_rate,region,life_exp_female, urban_q) %>%
                    filter(!is.na(urban_q),region=="Americas") %>%
                    filter(Country.Name %in% c("Guatemala","Haiti")),
                  aes(label=Country.Name),
                  color="grey30",
                  nudge_x=0.1,
                  nudge_y=-0.1,
                  box.padding = 0.9,
                  segment.curvature=-0.1,
                  direction="both",
                  min.segment.length = 0,
                  max.overlaps = Inf)



## ----corr matrix region,fig.width=9---------------------------------------------------------------------------
cor_region <- df %>% 
  select(`Country.Name`,urban_q,birth_rate,region,life_exp_female) %>%
  filter(!is.na(urban_q)) %>%
  pivot_wider(names_from = region,values_from = birth_rate) %>% 
  select(-c(Country.Name,urban_q)) %>% 
  cor(use = "pairwise.complete.obs") 
  # View()
# ggpairs(columns = 3:8, ggplot2::aes(colour=urban_q,alpha=0.5))

cor_reg <-  cor_region[,1] %>% data.frame(cor=.) %>% rownames_to_column(var="region")
cor_reg %>% 
  # View()
  filter(cor != 1) %>% 
  mutate_if(is.numeric,round,2) %>% 
  ggplot(aes(y=region,x=cor))+
  geom_col(alpha=0.7)+
  theme_cowplot()+
  geom_text(aes(label=cor,hjust=-0.2),color='white')+
  labs(x="Корреляция",y="Регион", subtitle="Корреляция между уровнем рождаемости и\nожидаемой продолжительностью жизни у женщин по регионам")


## ----corr matrix,fig.width=9----------------------------------------------------------------------------------
df %>% 
  select(`Country.Name`,urban_q,birth_rate,gdp,trade,water,resour,life_exp_female) %>%
  mutate(gdp=log(gdp),trade=log(trade)) %>%
  rename(log_gdp=gdp,log_trade=trade) %>% 
  filter(!is.na(urban_q)) %>%
ggpairs(columns = 2:8, ggplot2::aes(colour=urban_q,alpha=0.5))



## ----pca------------------------------------------------------------------------------------------------------
df_pc <- select(df,-c(birth_rate,urban_q,region,sub.region,electr_acc)) %>% remove_rownames() 

df_pc <- column_to_rownames(df_pc,var="Country.Name")
# View(df_pc)
df_pca <-  prcomp(df_pc, scale = TRUE)
# fviz_eig(df_pca)
fviz_pca_biplot(df_pca, repel = TRUE)


## ----simple model,fig.align='left'----------------------------------------------------------------------------
fit.smpl <- lm(birth_rate ~ water+log(gdp), df)
fit.smpl_hc <- lm_robust(birth_rate ~ water+log(gdp), df,se_type = "HC3")
ht <- huxreg("simple"=fit.smpl,"simple robust"=fit.smpl_hc, 
       number_format = 3,
       ci_level=0.95,
       error_pos = "right",
       # error_format = "({conf.low} -- {conf.high})",
       statistics = c("N. obs." = "nobs", "R squared" = "r.squared", "F statistic" = "statistic","P value" = "p.value")) %>% set_caption("Сравнение простых моделей")
position(ht) <- "left"
ht


## ----complex model--------------------------------------------------------------------------------------------
fit.complex <- lm(birth_rate ~ water+resour+region+log(gdp)+log(trade)+life_exp_female+urban_q,df)
fit.complex_hc <- lm_robust(birth_rate ~ water+resour+region+log(gdp)+log(trade)+life_exp_female+urban_q,df,se_type = "HC3")
ht <- huxreg("complex"=fit.complex,"complex robust"=fit.complex_hc, 
       number_format = 2,
       ci_level=0.95,
       error_pos = "right",
       error_format = "({conf.low} : {conf.high})",
       statistics = c("N. obs." = "nobs", "R squared" = "r.squared", "F statistic" = "statistic","P value" = "p.value")) %>% 
  map_background_color(by_values("regionOceania" = "#ff9c91","log(gdp)" = "#ff9c91","urban_qlow" = "#ff9c91")) %>% 
  set_caption("Сравнение сложных моделей")
position(ht) <- "left"
ht


## ----complex model without gdp--------------------------------------------------------------------------------
fit.complex <- lm(birth_rate ~ water+resour+region+log(trade)+life_exp_female,df)
fit.complex_hc <- lm_robust(birth_rate ~ water+resour+region+log(trade)+life_exp_female,df,se_type = "HC3")
ht <- huxreg("complex"=fit.complex,"complex robust"=fit.complex_hc, 
       number_format = 2,
       ci_level=0.95,
       error_pos = "right",
       error_format = "({conf.low} : {conf.high})",
       statistics = c("N. obs." = "nobs", "R squared" = "r.squared", "F statistic" = "statistic","P value" = "p.value")) %>% 
  map_background_color(by_values("regionOceania" = "#ff9c91","log(gdp)" = "#ff9c91")) %>% 
  set_caption("Сравнение сложных моделей")
position(ht) <- "left"
ht


## ----vif------------------------------------------------------------------------------------------------------
pander(vif(fit.complex))


## ----anova----------------------------------------------------------------------------------------------------
pander(anova(fit.smpl,fit.complex))


## -------------------------------------------------------------------------------------------------------------
pander(resettest(fit.complex_hc))


## -------------------------------------------------------------------------------------------------------------
plot_model(fit.complex_hc,show.values = TRUE,vline.color = "darkred")+
    # scale_y_continuous()+
  theme_sjplot2()+
  labs(title="Коэффициенты для построенной модели",  y="Значения коэффициентов с доверительным интервалом (95%)")


## ----sample test----------------------------------------------------------------------------------------------
set.seed(123)
# dim(df)
# 181*0.7=126 - сколько нужно выбрать рядов для 70%
randrow <- sample(nrow(df),126) 
df.train <- df[randrow,]
df.test <- df[-randrow,]


## ----fit train------------------------------------------------------------------------------------------------
fit.train.smpl_hc <- lm_robust(birth_rate ~ water+log(gdp), df.train,se_type = "HC3")
fit.train.complex_hc <- lm_robust(birth_rate ~ water+resour+region+log(trade)+life_exp_female, df.train,se_type = "HC3")
# huxreg(fit.train.smpl_hc)
ht <- huxreg("simple robust"=fit.train.smpl_hc,"complex robust"=fit.train.complex_hc, 
       number_format = 2,
       ci_level=0.95,
       error_pos = "right",
       error_format = "({conf.low} : {conf.high})",
       statistics = c("N. obs." = "nobs", "R squared" = "r.squared", "F statistic" = "statistic","P value" = "p.value")) %>% 
  map_background_color(by_values("regionOceania" = "#ff9c91",
                                 "regionAmericas" = "#ff9c91",
                                 "regionAsia"="#ff9c91",
                                 "log(trade)"="#ff9c91")) %>%
  set_caption("Сравнение сложных моделей")
position(ht) <- "left"
ht


## ----train data scatter---------------------------------------------------------------------------------------
df.train %>% 
  select(`Country.Name`,birth_rate,region,life_exp_female, urban_q) %>%
  # mutate(gdp=log(gdp)) %>%
  # filter(!is.na(urban_q)) %>%
  # pivot_longer(cols=c(4:8)) %>% 
  ggplot(aes(x=birth_rate,y=life_exp_female,color=urban_q, text= `Country.Name`))+
  geom_point()+
  facet_wrap(.~region,scales = "free")+
  theme_cowplot()+
  labs(color="Урбанизация", x="Уровень рождаемости", y="Ожидаемая продолжительность жизни женщин")


## ----summary train,fig.cap="Количество наблюдений по регионам в обучающей выборке"----------------------------
df.train %>% 
  group_by(region) %>% 
  summarize(`Количество наблюдений`=length(birth_rate))


## ----train data  boxplot--------------------------------------------------------------------------------------
df.train %>% 
  select(`Country.Name`,birth_rate,region,life_exp_female, urban_q) %>%
  # mutate(gdp=log(gdp)) %>%
  # filter(!is.na(urban_q)) %>%
  # pivot_longer(cols=c(4:8)) %>% 
  ggplot(aes(x=birth_rate,y=region,fill=urban_q))+
  geom_boxplot(alpha=0.9,width=0.5)+
  # geom_jitter(alpha=0.3,height = 0.01,width = 0.1)+
  # facet_wrap(.~region,scales = "free")+
  theme_cowplot()+
  labs(fill="Урбанизация", x="Уровень рождаемости", y="Регионы")


## ----lasso model----------------------------------------------------------------------------------------------

y <- df.train$birth_rate %>% unlist()
# create dummy for region

df.train_regionq <- df.train %>% 
  mutate(Europe_q=ifelse(region == "Europe",1,0),
         Asia_q=ifelse(region == "Asia",1,0),
         Africa_q=ifelse(region == "Africa",1,0),
         Americas_q=ifelse(region == "Americas",1,0),
         urban_q_low=ifelse(urban_q==1,1,0)) 
X0 <- model.matrix(data=df.train_regionq,birth_rate ~ 0+log(trade)+urban_q_low+
                                              water+resour+Europe_q+Asia_q+
                                              Africa_q+Americas_q+life_exp_female)

cv_lambda = cv.glmnet(X0, y, alpha = 1)
plot(cv_lambda)

fit.train.lasso <- glmnet(X0,y,alpha=1,lambda = cv_lambda$lambda.1se,standardize = TRUE)


## ----lambdas comparison---------------------------------------------------------------------------------------
tibble(lambda.1se=cv_lambda$lambda.1se,lambda.min=cv_lambda$lambda.min)


## ----lambda coeff---------------------------------------------------------------------------------------------

tidy(fit.train.lasso) %>%
  select(term,estimate)


## ----test prediction and mse----------------------------------------------------------------------------------
predict.smpl <- predict(fit.train.smpl_hc,newdata = df.test)
predict.cmplx <- predict(fit.train.complex_hc,newdata = df.test)
df.test_regionq <- df.test %>% 
  mutate(Europe_q=ifelse(region == "Europe",1,0),
         Asia_q=ifelse(region == "Asia",1,0),
         Africa_q=ifelse(region == "Africa",1,0),
         Americas_q=ifelse(region == "Americas",1,0),
         urban_q_low=ifelse(urban_q==1,1,0)) 

X0_test <- model.matrix(data=df.test_regionq,birth_rate ~ 0+log(trade)+urban_q_low+
                                              water+resour+Europe_q+Asia_q+
                                              Africa_q+Americas_q+life_exp_female)
predict.lasso <- predict(fit.train.lasso,newx = X0_test)
predict <- tibble(y=df.test$birth_rate,simple=predict.smpl,complex=predict.cmplx,
                 lasso=predict.lasso,Country=df.test$Country.Name,country=df.test$Country.Name,
                 region=df.test$region,urban_q=df.test$urban_q)
# View(predict)
tibble(`MSE lasso`=mean((predict$y-predict$lasso)^2),
`MSE simple`=mean((predict$y-predict$simple)^2),
`MSE complex`=mean((predict$y-predict$complex)^2),
`MAPE lasso`=mean( abs((predict$y-predict$lasso)/predict$y)*100),
`MAPE simple`=mean( abs((predict$y-predict$simple)/predict$y)*100),
`MAPE complex`=mean( abs((predict$y-predict$complex)/predict$y)*100)) %>% round(2)


## ----scatter test vs y----------------------------------------------------------------------------------------
predict %>% 
  ggplot(aes(x=lasso,y=y))+
  geom_point(size=2,aes(color=urban_q))+
  theme_cowplot()+
  background_grid()+
  geom_abline(slope = 1,intercept = 0,linetype=2,color='red')+
  # geom_smooth(method = "lm")+
  labs(x="Прогноз (LASSO)",y="Исходная переменная (рождаемость)",color="Урбанизация", subtitle="Кроссплот прогноз-исходная переменная, \nкрасная линия y=x")



## ----residuals lasso------------------------------------------------------------------------------------------
predict %>% 
  mutate(resid_lasso=y-lasso) %>% 
  ggplot(aes(x=lasso,y=resid_lasso,color=region))+
  geom_point(size=2)+
  theme_cowplot()+
  background_grid()+
  geom_abline(slope = 0,intercept = 0,linetype=2,color='red')+
  labs(x="Прогноз", y="Остатки",color="Регион", subtitle="Кроссплот остатки-прогноз для модели LASSO")
  



## ----predict complex_hc---------------------------------------------------------------------------------------
predict_2018 <- predict(fit.complex_hc,newdata=df_2018,interval="predict") %>% data.frame()
predict_2018_conf <- predict(fit.complex_hc,newdata=df_2018,interval="confidence") %>% data.frame()


df_2018_pred <- df_2018 %>% 
  bind_cols(fit=predict_2018$fit.fit,lwr=predict_2018$fit.lwr,upr=predict_2018$fit.upr,
            lwr_conf=predict_2018_conf$fit.lwr,upr_conf=predict_2018_conf$fit.upr)



## ----map for 2018 data,fig.width=9----------------------------------------------------------------------------
df_2018_pred %>% 
  mutate(resid=fit-birth_rate) %>% 
  ggplot(aes(y=resid,x=fit))+
  geom_point(alpha=0.5,size=3,aes(color=urban_q))+
  theme_cowplot()+
  geom_abline(slope=0.2,intercept = 2,linetype=2,color='grey60')+
  geom_abline(slope=-0.2,intercept = -2,linetype=2,color='grey60')+
  geom_abline(slope=0,intercept = 0,linetype=2,color='red')+
  labs(x="Предсказанный уровень рождаемости", y="Остатки",color="Урбанизация")
# View(df_2018_pred)


## ----Breusch-Pagan test---------------------------------------------------------------------------------------
bptest(fit.complex_hc)


## ----import shape, cache=TRUE---------------------------------------------------------------------------------
#import shape files for world
world_shapefiles <- read_sf(dsn = "~/Documents/DataQ/RU_Ref_2020_07_01/data_polygons/World/ne_50m_admin_0_countries.shp") 


## -------------------------------------------------------------------------------------------------------------
#join with summary by country 
df_2018_results <- world_shapefiles %>%
    left_join(df_2018_pred, by= c( "GU_A3"="Country.Code")) %>% 
  select(NAME,birth_rate,region,fit,lwr,upr,lwr_conf,upr_conf,ISO_A3) %>% 
  mutate(resid=(fit-birth_rate)/birth_rate*100) %>% 
  mutate_if(is.numeric,round,2)


## # .leaflet-container { background: #F1FDFF; }


## ----map with residuals, fig.width=9--------------------------------------------------------------------------

mybins <- c(-90,-50,-17,-5,5,17,50,90)
mypalette_resid <- colorBin( palette="RdBu", domain=df_2018_results$resid, na.color="#a3a2a2", bins=mybins,reverse = TRUE)
# summary(df_2018_results$resid)
mybins_birth <- c(5,10,15,25,35,50)
mypalette_birth <- colorBin( palette="Greens", domain=df_2018_results$birth_rate, na.color="#a3a2a2", bins=mybins_birth)

mytext <- paste("Country:",df_2018_results$NAME, "<br/>","Birth Rate:", df_2018_results$birth_rate,
                           "<br/>","Fit:",df_2018_results$fit,
                           "<br/>","Resid:",round(df_2018_results$resid,0),"%",
                           "<br/>","[lwr : up]","[",df_2018_results$lwr,":", df_2018_results$upr,"]") %>% 
  lapply(htmltools::HTML)

leaflet(width = "100%") %>%  
  setView(1,35,zoom = 2) %>%
addPolygons( group="Residual",
  data= df_2018_results,
            weight = 1,
            stroke = TRUE,color="grey",
            label=mytext,
            labelOptions = labelOptions( 
      style = list("font-weight" = "normal", padding = "3px 8px"), 
      textsize = "13px", 
      direction = "auto"
    ),
            fillColor  = ~mypalette_resid(resid),
            fillOpacity = 1) %>% 
  
  addPolygons( group="Birth Rate",
  data= df_2018_results,
            weight = 1,
            stroke = TRUE,color="grey",
            label=mytext,
            labelOptions = labelOptions( 
      style = list("font-weight" = "normal", padding = "3px 8px"), 
      textsize = "13px", 
      direction = "auto"
    ),
            fillColor  = ~mypalette_birth(birth_rate),
            fillOpacity = 1) %>% 
  addLegend(pal = mypalette_resid, values = df_2018_results$resid, title= "Ошибка в %",opacity=0.9,
            group="Residual",position="bottomleft") %>%
  addLegend(pal = mypalette_birth, values = df_2018_results$resid, title= "Birth rate",opacity=0.9,
            group="Birth Rate",position="bottomleft") %>%
  addLayersControl(
    overlayGroups =c("Residual","Birth Rate"),
    options = layersControlOptions(collapsed = FALSE),position = "topright") %>% 
  hideGroup("Birth Rate")


## ----predict scattter-----------------------------------------------------------------------------------------
gg1 <- df_2018_results %>% 
  mutate_if(is.numeric,round,2) %>% 
  ggplot(aes(y=birth_rate,x=fit))+
  geom_point(aes(color=region,text=paste(
    "Страна:",df_2018_results$NAME,
    "<br>Уровень\nрождаемости: ", round(df_2018_results$birth_rate,0),
    "<br>Прогноз: ", round(df_2018_results$fit,0),
    "<br>Предиктивный\nинтервал: ", round(df_2018_results$lwr,0), ":", round(df_2018_results$upr,0))),size=2)+
  geom_smooth(aes(x=fit,y=upr,text=df_2018_results$upr),method="lm",se=FALSE,linetype=2,size=0.7,color="grey40")+
  geom_smooth(aes(x=fit,y=lwr),method="lm",se=FALSE,linetype=2,size=0.7,color="grey40")+
  geom_smooth(aes(x=fit,y=upr_conf),method="lm",se=FALSE,linetype=4,size=0.5,color="blue")+
  geom_smooth(aes(x=fit,y=lwr_conf),method="lm",se=FALSE,linetype=4,size=0.5,color="blue")+
  geom_abline(intercept = 0,slope = 1,color="red",linetype=4)+
  theme_cowplot()+
  background_grid()+
  labs(x="Прогноз", y="Уровень рождаемости", color="Регион")

ggplotly(gg1,tooltip="text")