Stavíme si dataset
Cílem tohoto cvičení je projít praktickou ukázku, ve které si ukážeme jak se staví dataset pro výzkum v mezinárodní migraci. Budeme používat pouze volně dostupná data.
- Ukážeme si několik nových balíků pro stahování dat z otevřených databází
- Provičíme praktikou aplikaci nástrojů z balíků tidyr a dplyr
- Ukážeme si, jaké problémy můžeme potkat v divočině
Trocha teorie…
V neoklasických modelech migrace se vychází z toho, že migrant (nebo jeho rodina) maximalizují užitek z migrace. Kontinuálně zvažují, do které země \(j \in J\) se přestěhují. Součástí množiny \(J\) je i země \(i\), ve které se momentálně nalézají. Užitek z migrace do zěmě \(j\) je přitom dán jako rozdíl mezi očekávanými přínosy z migrace a jejími náklady.
Na hraní můžeme použít následující empirickou specifikaci migračního modelu:
\[\log(\frac{M_{ijt}}{P_{it}})=\alpha + \beta_1 \log(GDP_{i,t-1}) + \beta_2 \log(GDP_{j,t-1}) + \beta_3 \log(U_{i,t-1}) + \beta_4 \log(U_{j,t-1})\] \[+ \beta_5\log(\frac{M_{i,t-1}}{P_{i,t-1}}) + \mathbf{\gamma} X_{ij} + \Theta_i + \Theta_j + \Theta_t + \varepsilon\] Kde \(M\) je migrační tok ze země \(i\) do země \(j\) v čase \(t\), \(P\) je velikost populace, \(GDP\) je GDP per capita, \(U\) je míra nezaměstnanosti a \(X\) je vektor bilaterálních gravotačních proměnných (vzdálenost, atp.).
Získávání dat z connections
Stahování dat ze vzdálených zdrojů
V divočině na síti najdete mnoho zdrojů dat. Pro jednoduchost se dají rozdělit do tří skupin:
- zdroj má API (Application Programming Interface)
- zdroj nemá API, ale umožňuje interakci například prostřednictvím simulace chování lidského uživatele
- data se válí v podobě připraveného datasetu (
dta– Stata,sav– SPSS,csv,xls) na webu - data jsou připravená v balíčku pro R
Zdroje s API
Pro mnoho poskytovatelů dat (statistických úřadů,…) existují možnosti jak stahovat aktuální data “přímo do R”. Zajímavá je tato možnost zejména v případě, že existuje balík specializovaný pro komunikaci s API konkrétního poskytovatele.
World bank
Balík WDI obsahuje World Development Indicators – tedy většinu dat veřejně poskytovaných World Bank
library(WDI)Balík obsahuje dvě základní funkce: WDIsearch a WDI. WDIsearch umožňuje vyhledávat indikátory z databáze WDI.
WDIsearch(string = "gdp", field = "name", short = TRUE, cache = NULL)Vrací matrix (help tvrdí, že data.frame) s celkem logickou strukturou:
WDIsearch() %>% head## indicator
## [1,] "BG.GSR.NFSV.GD.ZS"
## [2,] "BM.KLT.DINV.GD.ZS"
## [3,] "BN.CAB.XOKA.GD.ZS"
## [4,] "BN.CUR.GDPM.ZS"
## [5,] "BN.GSR.FCTY.CD.ZS"
## [6,] "BN.KLT.DINV.CD.ZS"
## name
## [1,] "Trade in services (% of GDP)"
## [2,] "Foreign direct investment, net outflows (% of GDP)"
## [3,] "Current account balance (% of GDP)"
## [4,] "Current account balance excluding net official capital grants (% of GDP)"
## [5,] "Net income (% of GDP)"
## [6,] "Foreign direct investment (% of GDP)"Vlastní stahování dat zařizuje funkce WDI():
WDI(country = "all", indicator = "NY.GNS.ICTR.GN.ZS",
start = 2005, end = 2011, extra = FALSE, cache = NULL)country– ISO-2 kódindicator– kód indikátoru – lze najít pomocíWDIsearch()nebo na webu World Bank (http://data.worldbank.org/). Tamtéž lze najít bližší popis indikátorů.startaend– udává první a poslední rok. WDI obsahuje pouze roční data.extra– pokud je nastaveno naTRUEvrací více proměnné – jméno regionu, iso3c kód,… (defaultněFALSE, doporučujiTRUE)
WDI(indicator = c("NY.GNS.ICTR.GN.ZS","BN.KLT.DINV.CD.ZS"), extra = TRUE) -> wdi_data
wdi_data %>% class## [1] "data.frame"head(wdi_data)## iso2c country year NY.GNS.ICTR.GN.ZS BN.KLT.DINV.CD.ZS iso3c region
## 1 Africa 2005 NA 2.816084 <NA> <NA>
## 2 Africa 2006 NA 2.504447 <NA> <NA>
## 3 Africa 2007 NA 3.028385 <NA> <NA>
## 4 Africa 2008 NA 3.183462 <NA> <NA>
## 5 Africa 2009 NA 2.896574 <NA> <NA>
## 6 Africa 2010 NA 2.382063 <NA> <NA>
## capital longitude latitude income lending
## 1 <NA> <NA> <NA> <NA> <NA>
## 2 <NA> <NA> <NA> <NA> <NA>
## 3 <NA> <NA> <NA> <NA> <NA>
## 4 <NA> <NA> <NA> <NA> <NA>
## 5 <NA> <NA> <NA> <NA> <NA>
## 6 <NA> <NA> <NA> <NA> <NA>Všiměte si, že balík WDI Vám dodá tidy data. Webové rozhraní nikoliv.
Nepříjemné je, že kromě dat pro jednotlivé země se vám vrátí i data pro aregáty (viz výše). Ty můžeme jednoduše odfiltrovat:
wdi_data %>%
as_tibble() %>%
filter(!is.na(iso3c)) %>%
print(n=5)## # A tibble: 1,638 × 12
## iso2c country year NY.GNS.ICTR.GN.ZS BN.KLT.DINV.CD.ZS iso3c
## <chr> <chr> <dbl> <dbl> <dbl> <fctr>
## 1 AD Andorra 2011 NA NA AND
## 2 AD Andorra 2007 NA NA AND
## 3 AD Andorra 2008 NA NA AND
## 4 AD Andorra 2010 NA NA AND
## 5 AD Andorra 2006 NA NA AND
## # ... with 1,633 more rows, and 6 more variables: region <fctr>,
## # capital <fctr>, longitude <fctr>, latitude <fctr>, income <fctr>,
## # lending <fctr>OECD
Balík OECD pro stahování dat z databáze OECD (http://stats.oecd.org/) funguje velmi podobně jako balík WDI. Poskytuje funkce pro hledání indikátorů a pro samotné stažení dat.
library(OECD)Funkce pro hledání indikátorů:
search_dataset(string, data = get_datasets(), ignore.case = TRUE)string– regulární výraz použitý pro vyhledání indikátoru
search_dataset("FDI") %>% head## # A tibble: 6 × 2
## id title
## <fctr> <fctr>
## 1 FDI_FLOW_INDUSTRY FDI flows by industry
## 2 FDI_POSITION_INDUSTRY FDI positions by industry
## 3 FDI_FLOW_PARTNER FDI flows by partner country
## 4 FDI_BOP_IIP FDI series of BOP and IIP aggregates
## 5 FDI_POSITION_PARTNER FDI positions by partner country
## 6 AEO11_INDEPTH_CHAPTER6_FIG5_EN Figure 6.5: African FDI inflows 1995-2008Samotná funkce pro stažení dat funguje podobně jako její ekvivalet z balíku WDI:
get_dataset(dataset, filter = NULL, start_time = NULL, end_time = NULL,
pre_formatted = FALSE, ...)dataset– string specifikující požadovaný dataset (idzsearch_dataset())filter– podmínka pro filtrování dat. Pokud není nastavena, stáhne se celý dataset.
get_dataset("EPL_OV",
filter = list(c("DEU", "FRA"),
c("EPRC_V1", "EPRC_V2")),
start_time = 2008, end_time = 2010) -> oecd_dataParametr filter může být velmi užitečný. Některé databáze OECD jsou velmi obsáhlé – klidně desítky MB – a jejich stahování pomocí balíků OECD bývá problémové. Nastavení filter může tento problém vyřešit. Pokud se tento problém vyskytne, tak spíše doporučuji stahovat data přes webové rozhraní. V takovém případě doporučuji vždy stahovat komprimované CSV. (Web OECD má nepříjemnou vlastnost. Při exportu do xls vrací pouze vyfiltrovanou část tabulky. Při exportu do CSV vrací celou. Tato vlastnost není, pokud vím, nikde dokumentovaná.)
Stahování dat přes balík OECD je obecně velmi pomalé a dost nestabilní.
oecd_data %>% print(n=5)## # A tibble: 12 × 5
## COUNTRY SERIES TIME_FORMAT obsTime obsValue
## * <chr> <chr> <chr> <chr> <dbl>
## 1 FRA EPRC_V1 P1Y 2008 2.468254
## 2 FRA EPRC_V1 P1Y 2009 2.384921
## 3 FRA EPRC_V1 P1Y 2010 2.384921
## 4 DEU EPRC_V1 P1Y 2008 2.678571
## 5 DEU EPRC_V1 P1Y 2009 2.678571
## # ... with 7 more rowsJe příjemné, že get_dataset() vrací tibble, který sice není tidy, ale obyvkle (jako v tomto případě) není konverze do tidy dat problém:
oecd_data %>%
spread(SERIES,obsValue) %>%
print(n=5)## # A tibble: 6 × 5
## COUNTRY TIME_FORMAT obsTime EPRC_V1 EPRC_V2
## * <chr> <chr> <chr> <dbl> <dbl>
## 1 DEU P1Y 2008 2.678571 2.948980
## 2 DEU P1Y 2009 2.678571 2.948980
## 3 DEU P1Y 2010 2.678571 2.948980
## 4 FRA P1Y 2008 2.468254 2.727324
## 5 FRA P1Y 2009 2.384921 2.667800
## # ... with 1 more rowsEurostat
Data z Eurostatu jsou dostupná přes balík eurostat.
library(eurostat)Ten poskytuje celou řadu funkcí, které mohou být velmi užitečné (umožňuje například stahovat geolokovaná – spatial – data atp.). Základní funkcionalitou je ovšem hledání a stažení dat.
Vhledat data lze pomocí funkce search_eurostat():
search_eurostat(pattern, type = "dataset")type– může nabývat hodnotdataset,folder,tableneboall
search_eurostat("migration") %>% names## [1] "title" "code"
## [3] "type" "last update of data"
## [5] "last table structure change" "data start"
## [7] "data end" "values"search_eurostat("migration") %>% head()## # A tibble: 6 × 8
## title
## <chr>
## 1 Main scenario - Net migration by age, sex and NUTS 2 regions
## 2 Main scenario - Net migration by age, sex and NUTS 3 regions
## 3 Immigration by age and sex
## 4 Immigration by five year age group, sex and citizenship
## 5 Immigration by five year age group, sex and country of birth
## 6 Immigration by age, sex and broad group of citizenship
## # ... with 7 more variables: code <chr>, type <chr>, `last update of
## # data` <chr>, `last table structure change` <chr>, `data start` <chr>,
## # `data end` <chr>, values <chr>search_eurostat("migration", type = "folder") %>% head()## # A tibble: 6 × 8
## title code type
## <chr> <chr> <chr>
## 1 Demography and migration pop folder
## 2 Immigration migr_immi folder
## 3 Emigration migr_emi folder
## 4 Asylum and managed migration migr folder
## 5 Enforcement of Immigration Legislation migr_eil folder
## 6 Population by migration characteristics cens_11rmc folder
## # ... with 5 more variables: `last update of data` <chr>, `last table
## # structure change` <chr>, `data start` <chr>, `data end` <chr>,
## # values <chr>search_eurostat("migration", type = "table") %>% head()## # A tibble: 4 × 8
## title code type
## <chr> <chr> <chr>
## 1 Crude rate of net migration plus adjustment tsdde230 table
## 2 Immigration tps00176 table
## 3 Emigration tps00177 table
## 4 Crude rate of net migration plus adjustment tsdde230 table
## # ... with 5 more variables: `last update of data` <chr>, `last table
## # structure change` <chr>, `data start` <chr>, `data end` <chr>,
## # values <chr>Pro vlastní satžení dat je klíčová funkce get_eurostat():
get_eurostat(id, time_format = "date", filters = "none", type = "code",
select_time = NULL, cache = TRUE, update_cache = FALSE,
cache_dir = NULL, compress_file = TRUE,
stringsAsFactors = default.stringsAsFactors(), keepFlags = FALSE, ...)id– kód datasetuselect_time– umožňuje zvolit si frekvenci pozorování
# Immigration by age and sex
get_eurostat("migr_imm8") -> eu_data## Table migr_imm8 cached at /tmp/RtmpxQquXy/eurostat/migr_imm8_date_code_TF.rdseu_data %>% print(n=5)## # A tibble: 158,914 × 7
## agedef age unit sex geo time values
## <fctr> <fctr> <fctr> <fctr> <fctr> <date> <dbl>
## 1 COMPLET TOTAL NR F AT 2014-01-01 52494
## 2 COMPLET TOTAL NR F BE 2014-01-01 58318
## 3 COMPLET TOTAL NR F BG 2014-01-01 11903
## 4 COMPLET TOTAL NR F CH 2014-01-01 75173
## 5 COMPLET TOTAL NR F CY 2014-01-01 6422
## # ... with 1.589e+05 more rowsget_eurostat() vrací tibble, který se dá relativně lehce překlopit do tidy formátu.
Ostatní
Balíků, které umožňují přístup k datům je celá řada, viz například:
- https://www.quandl.com/tools/r
- https://cran.r-project.org/web/packages/quantmod/index.html
- https://cran.r-project.org/web/views/Finance.html
- https://cran.r-project.org/web/views/TimeSeries.html
Existují i balíky, které umožňují načítat data z webových služeb a sociálních sítí (Wikipedia, Google Trends,…):
Web Scraping
Lokální soubory
Data obsažená v R packages
Některá užitečná data mají své vlastní balíky, které zpravidla neobsahují nic jiného než právě data. Mezi zajímavé kousky patří:
- Penn World Tables v balících
pwtapwt8 - World Population Prospects v balících
wpp2008,wpp2010,wpp2012awpp2015 - …
Tabulka postavená od základů
Pro přípravu tabulky je vhodné si rozmyslet, jaký bude její záběr. Chceme zahrnout všechna myslitelná pozorování, nebo nám stačí soustředit se na ty, u kterých můžeme získat data?
V tomto konkrétním případě je rozumné zvolit spíše druhou variantu. Ta v praxi znamená, že získáme maximální rozsah dat pro vysvětlovanou proměnnou a následně budeme doplňovat data pro pravou stranu rovnice.
Migrační data
Data o migraci jsou dostupná v několika databázích:
- World Bank (ne WDI!) – počty imigrantů dle země narození pro celá desetiletí 1960, 1970, …, 2000 (Özden et al. 200X)
- UN – počty imigrantů pro celá destiletí a nově jejich interpolace pro poloviny desetiletí
- Odhadované globální bilaterální migrační toky (Abel, Sander 2014, Abel 2013, 2015, 2016) dostupné v agregaci pro pětiletky – nástroje pro odhady jsou dostupné v balíku
migest - Databáze OECD zahrnující roční data o migraci do OECD – velmi reprezentativní zdroj pro migraci po roce 1990
- Databáze z Adsera & Pytliková (2015), která doplňuje chybějící hodnoty v datech OECD
- Databáze Eurostatu – překvapivě děravá. Skoro nikdo ji nepoužívá.
My budeme používat databázi OECD International Migration Database. Tu lze stáhnout pomocí balíku OECD nebo z webové aplikace. Výsledný soubor CSV má 89 MB a proto je celkem jasnou volbou stahování z webové aplikace. Výsledný soubor je MIG_OECD.CSV.
U práce s takto velkými soubory už se vyplatí přemýšlet o tom, jakou funkci použít pro načtení dat:
system.time(
read.csv("data/MIG_OECD.csv") -> data
)## user system elapsed
## 7.560 0.208 7.862system.time(
read_csv("data/MIG_OECD.csv",
col_types = "cccccccccnncc") -> data
)## user system elapsed
## 1.372 0.028 1.401read_csv() z balíku readr je nepřekvapivě jasný vítěz.
data %>% print(n=5)## # A tibble: 718,142 × 13
## `"CO2"` `Country of birth/nationality` VAR
## <chr> <chr> <chr>
## 1 AFG Afghanistan B11
## 2 AFG Afghanistan B11
## 3 AFG Afghanistan B11
## 4 AFG Afghanistan B11
## 5 AFG Afghanistan B11
## # ... with 7.181e+05 more rows, and 10 more variables: Variable <chr>,
## # GEN <chr>, Gender <chr>, COU <chr>, Country <chr>, YEA <chr>,
## # Year <dbl>, Value <dbl>, `Flag Codes` <chr>, Flags <chr>Data z OECD obsahují hodně sloupců, jejichž obsah je dost kriptický. Rozhodně stojí za bližší podívání:
data %>% distinct(VAR, Variable)## # A tibble: 10 × 2
## VAR Variable
## <chr> <chr>
## 1 B11 Inflows of foreign population by nationality
## 2 B12 Outflows of foreign population by nationality
## 3 B14 Stock of foreign-born population by country of birth
## 4 B13 Inflows of asylum seekers by nationality
## 5 B15 Stock of foreign population by nationality
## 6 B16 Acquisition of nationality by country of former nationality
## 7 B21 Inflows of foreign workers by nationality
## 8 B22 Inflows of seasonal foreign workers by nationality
## 9 B23 Stock of foreign-born labour by country of birth
## 10 B24 Stock of foreign labour by nationalityNa první pohled je vidět, že vstupní soubor není tidy a že dokonce obsahuje data, která nás nezajímají. Saoustředíme se na proměnné B11, B14 a B15. B11 slibuje nejlepší pokrytí daty. Proměnné B14 a B15 použijeme pro aproximaci diaspor. Ostatní proměnné (řádky) zahodíme:
data %<>%
filter(VAR %in% c("B11","B14","B15"))
print(data)## # A tibble: 321,684 × 13
## `"CO2"` `Country of birth/nationality` VAR
## <chr> <chr> <chr>
## 1 AFG Afghanistan B11
## 2 AFG Afghanistan B11
## 3 AFG Afghanistan B11
## 4 AFG Afghanistan B11
## 5 AFG Afghanistan B11
## 6 AFG Afghanistan B11
## 7 AFG Afghanistan B11
## 8 AFG Afghanistan B11
## 9 AFG Afghanistan B11
## 10 AFG Afghanistan B11
## # ... with 321,674 more rows, and 10 more variables: Variable <chr>,
## # GEN <chr>, Gender <chr>, COU <chr>, Country <chr>, YEA <chr>,
## # Year <dbl>, Value <dbl>, `Flag Codes` <chr>, Flags <chr>Další mysteriózní proměnná je “Flags”:
data %>%
rename(FlagCodes = `Flag Codes`) %>%
distinct(FlagCodes,Flags)## # A tibble: 3 × 2
## FlagCodes Flags
## <chr> <chr>
## 1 <NA> <NA>
## 2 b Break
## 3 e Estimated valueNyní je potřeba zjistit, co vlastně “flags” znamenají. Čisté řešení je nastudovat si metodiku zjišťování a publikování dat. My se podíváme na následující:
data %>%
filter(`Flag Codes`=="e", !is.na(Value))## # A tibble: 0 × 13
## # ... with 13 variables: "CO2" <chr>, Country of birth/nationality <chr>,
## # VAR <chr>, Variable <chr>, GEN <chr>, Gender <chr>, COU <chr>,
## # Country <chr>, YEA <chr>, Year <dbl>, Value <dbl>, Flag Codes <chr>,
## # Flags <chr>data %>%
filter(`Flag Codes`=="b", !is.na(Value))## # A tibble: 0 × 13
## # ... with 13 variables: "CO2" <chr>, Country of birth/nationality <chr>,
## # VAR <chr>, Variable <chr>, GEN <chr>, Gender <chr>, COU <chr>,
## # Country <chr>, YEA <chr>, Year <dbl>, Value <dbl>, Flag Codes <chr>,
## # Flags <chr>Zdá se, že z flags nás hlava bolet nemusí – u všechny hodnot jde o NA. Můžeme tedy celkem bez starostí vyčistit tabulku. Ponecháme jen relevantní sloupce:
data %<>%
select(
origin_code = 1,
origin_name = `Country of birth/nationality`,
gender = Gender,
destination_code = COU,
destination_name = Country,
year = Year,
VAR,
Value
)Pojmenovávání sloupců je věda sama o sobě. Každý postupuje jinak. V případě migračních toků se občas používá i značení pomocí dolních indexů: x_i pro charakteristiky země původu a x_ij pro charakteristiku bilaterálního vztahu (např. dataset k Adsera & Pytliková, 2015). Osobně preferuji ukecanější pojmenování, které je jen obtížně zaměnitelné. Nicméně je to věc vkusu.
Teď přišel čas data ještě trochu pročistit:
- provedeme konverzi roků na double/integer – v praxi je to celkem jedno
- převedeme data do tidy formátu
data %<>%
unite(variable,VAR,gender,sep="_") %>%
spread(variable,Value) %>%
rename(
flow = B11_Total,
flow_women = B11_Women,
stock_cob = B14_Total,
stock_cob_women = B14_Women,
stock_nat = B15_Total,
stock_nat_women = B15_Women
) %>%
mutate(
year = year %>% as.character() %>% as.numeric()
)
print(data, n=5)## # A tibble: 95,319 × 11
## origin_code origin_name destination_code destination_name year flow
## <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 AFG Afghanistan AUS Australia 1996 NA
## 2 AFG Afghanistan AUS Australia 1997 368
## 3 AFG Afghanistan AUS Australia 1998 620
## 4 AFG Afghanistan AUS Australia 1999 852
## 5 AFG Afghanistan AUS Australia 2000 887
## # ... with 9.531e+04 more rows, and 5 more variables: flow_women <dbl>,
## # stock_cob <dbl>, stock_cob_women <dbl>, stock_nat <dbl>,
## # stock_nat_women <dbl>Nyní máme v ruce základní verzi tabulky. Protože nemáme kontrolu nad zdrojem dat, vyplatí se trochu si je osahat. Zaměříme se zejména na hodnoty pro celkový tok flow, který nás zajímá jednoznačně nejvíce.
data %>%
filter(!is.na(flow)) %>%
summarise(
year_min = range(year)[1],
year_max = range(year)[2],
origins = unique(origin_code) %>% length(),
destinations = unique(destination_code) %>% length,
nobs = n()
)## # A tibble: 1 × 5
## year_min year_max origins destinations nobs
## <dbl> <dbl> <int> <int> <int>
## 1 1980 2014 206 35 81184Data by měla obsahovat $3520635=$2.523510^{5} pozorování. Sumarizační tabulka však ukazuje, že jich máme pouze něco přes 80.000. Je tedy zřejmé, že tabulka obsahuje implicitní chybějící hodnoty. To je něco, s čím se musí počítat.
Také je rozumné podívat se na hodnoty proměnných – dávají smysl?
Je užitečné dívat se na variabilitu v proměnných, no odhlehlé hodnoty, atp… Momentálně na to zatím nemáme všechny nástroje, ale zdá se, že hodnoty jsou rozděleny cca tak, jak je možné očekávat:
data %>%
gather(variable,value,flow:stock_nat_women) %>%
ggplot(
aes(x=value)
) +
geom_histogram(bins = 100) +
facet_wrap(~variable, scales = "free", ncol = 2) +
theme_bw()## Warning: Removed 250695 rows containing non-finite values (stat_bin).
Mrkneme se i na extrémní hodnoty pro flows:
data %>%
arrange(desc(flow))## # A tibble: 95,319 × 11
## origin_code origin_name destination_code destination_name year
## <chr> <chr> <chr> <chr> <dbl>
## 1 MEX Mexico USA United States 1991
## 2 MEX Mexico USA United States 1990
## 3 MEX Mexico USA United States 1989
## 4 YUCS Former Yugoslavia DEU Germany 1992
## 5 ROU Romania ITA Italy 2007
## 6 POL Poland DEU Germany 1989
## 7 YUCS Former Yugoslavia DEU Germany 1991
## 8 MEX Mexico USA United States 2002
## 9 MEX Mexico USA United States 1992
## 10 POL Poland DEU Germany 1988
## # ... with 95,309 more rows, and 6 more variables: flow <dbl>,
## # flow_women <dbl>, stock_cob <dbl>, stock_cob_women <dbl>,
## # stock_nat <dbl>, stock_nat_women <dbl>Zdá se, že vše dává elementární smysl. Samozřejmě tato část by si zasloužila větší pozornost. Zatím ale úplně nevíme jak na to.
Lagování
V rovnici, jejíž parametry chceme nakonec odhadovat jsou vysvětlující proměnné lagované a logaritmované. Na rozdíl od mnoha jiných software nás vůbec nemusí trápit logaritmování. R si umí data transformovat při samotném odhadu modelu – není potřeba data dopředu (pozdravujeme Gretl, Statu,…). V některých případech to platí i pro lagování – např. při odhadu na panelových datech pomocí balíku plm. Obecně však doporučuji provést lagování ručně. (Ne všechno co je potřeba je implementováno v plm/pglm – typicky zde není Poisson pseudo-ML estimátor. Navíc si můžeme být jisti, že bude lagování uděláno správně.)
Pro lagování můžeme použít jednu z tzv. window functions z balíku dplyr: lag(). (Nebylo na přednášce.)
Její fugnování je možné demostrovat na populačních datech pro Norsko a Švédsko:
load("data/POP.Rdata")
POP_implicit %>%
mutate(
poplation_lag1 = lag(population)
)## # A tibble: 7 × 4
## country year population poplation_lag1
## <chr> <int> <int> <int>
## 1 Norway 1995 4359788 NA
## 2 Norway 1998 4440109 4359788
## 3 Norway 2000 4491572 4440109
## 4 Sweden 1996 8849420 4491572
## 5 Sweden 1997 8859106 8849420
## 6 Sweden 1998 8861204 8859106
## 7 Sweden 2000 8872284 8861204Funkce funguje tak, že penáší hodnoty o n (defaultně n=1L) řádků níže. Co to znamená:
- výsledek je citlivý na řazení
- výsledek je je citlivý na implicitní chybějící hodnoty
Obojí plyne z toho, že R vlastně neví, že data jsou fakticky panel. (Balík plm má vlastní class pro data viz plm::plm.data(), která mu lagování umožňuje.) Doplnění lagovaných hodnot můžeme provést dvojí způsobem:
- přičtením jedničky k roku a joinem s původní tabulkou
system.time(
data %>%
select(
-origin_name,-destination_name,-flow,-flow_women
) %>%
mutate(
year = year + 1
) %>%
left_join(data,.,
by=c("origin_code","destination_code","year"),
suffix = c("_lag0","_lag1")
) -> xdata
)## user system elapsed
## 0.176 0.000 0.176- pomocí rozšíření tabulky tak, aby obsahovala explicitní chybějící proměnné a grupovaných operací
system.time(
data %>%
complete(origin_code,destination_code,year) %>%
group_by(origin_code,destination_code) %>%
arrange(year) %>%
mutate(
stock_cob_lag1 = lag(stock_cob),
stock_cob_women_lag1 = lag(stock_cob_women),
stock_nat_lag1 = lag(stock_nat),
stock_nat_women_lag1 = lag(stock_nat_women)
) -> data
)## user system elapsed
## 0.404 0.000 0.411data %<>% ungroup Oba způsoby vedou z praktického hlediska k identickým výsledkům bez zásadního rozdílu v náročnosti. V praxi bude záležet, který z postupů bude mít nákladnější režii – to se může lišit zejména (můj odhad) s náročností joinování.
Druhý způsob však vrací tabulku s explicitními chybějícími pozorováními. První způsob žádné řádky nepřidává. Na druhou stranu ovšem vyžaduje přejmenování původních proměnných.
Každému co jeho jest. Obecně velmi doporučuji mít v datsetu explicitní chybějící hodnoty. Proto spíše inklinuji k druhému řešení. Na druhou stranu nic nikomu nebrání v zavolání complete() po provedení prvního způsobu přidání lagovaných hodnot.
Kódování
Pro jednoznačné určení pozorování potřebujeme znát místo, ze kterého pochází. To je v našem případě země nebo dvojice zemí. Rozhodně není moudré pro joinování nebo cokoliv jiného používat jméno země. To se může dramaticky lišit mezi databázemi – podle množství členů, velikosti písmen (R je case-sensitive!), množství a umístění slov typu “demokratický”, “lidový”, “republika” a podobně. Proto se používají standardizované (kéžby) kódy zemí.
Různé databáze používají různé standardy kódování zemí, které mají – bohuželva – dosti odlišné vlastnosti. V zásadě exstují (a používají se!) následující standardy:
- ISO třípísmenné (často něco jako iso3c) a dvoupísmenné (často něco jako iso2c) kódy.
- ISO numerické kódy
- pochopitelně všechny mezinárodní instituce mají za účelem zdůraznění své vyjímečnosti vlastní systémy kódování: UN, WB, IMF, FAO,…
Míra shody v rámci těchot systémů je různá a často jen zdánlivá a rozhodně na ni nelze spoléhat. Například kódování WB je skoro stejné jako iso3c. Skoro. Nejvýznamnější rozdíl je, že Rumunsko je v jednom “ROU” a ve druhém “ROM”. Tento rozdíl už zachránil tisíce životů.
Slabinou všech výše uvedených systémů je častá absence kódů pro historické (zaniklé) země, popřípadě i ignorování některých států (typicky Taiwan v databázi WB). Tímto problémem netrpí více historicky orientované standardy:
- Correlates of War (CoC) – numerická nebo písmenná varianta
- Gledish and Ward (GW)
Práci s kódováním může ulehčit balík countrycode, který obsahuje dva objekty: tabulku countrycode_data a funkci countrycode().
countrycode_data obsahuje vlastně pravidla pro konverzi jednotlivých kódování:
## country.name cowc cown fao fips104 imf ioc iso2c iso3c iso3n un wb
## 1 Afghanistan AFG 700 2 AF 512 AFG AF AFG 4 4 AFG
## 2 Aland Islands <NA> NA NA <NA> NA <NA> AX ALA 248 248 ALA
## 3 Albania ALB 339 3 AL 914 ALB AL ALB 8 8 ALB
## 4 Algeria ALG 615 4 AG 612 ALG DZ DZA 12 12 DZA
## 5 American Samoa <NA> NA NA AQ 859 ASA AS ASM 16 16 ASM
## 6 Andorra AND 232 6 AN NA AND AD AND 20 20 ADO
## regex continent region
## 1 afghan Asia Southern Asia
## 2 \\b(a|å)land Europe Northern Europe
## 3 albania Europe Southern Europe
## 4 algeria Africa Northern Africa
## 5 ^(?=.*americ).*samoa Oceania Polynesia
## 6 andorra Europe Southern EuropeFunkce countrycode() je využívá pro konverzi vstupního vektoru:
countrycode(sourcevar, origin, destination, warn = FALSE)Funkce bere několik parametrů – vstupní vektor a kód vstupního a výstupního kódování (viz help). Zrádné je nastavení parametru warn = FALSE. Tento paramtr potlačí vypsání všech zpráv o problemtické konverzi. Vždy nastavujte warn = TRUE. Balíček je zrádný v tom, že navozuje falešný pocit jistoty. Doporučuji vždy detailně kontrolovat vstupní i výstupní hodnoty.
Kódování v migrační databázi OECD
Databáze OECD přímo neodkazuje na použité kódování, takže se musíme podívat přímo do dat:
data %>%
filter(!is.na(origin_name)) %>%
distinct(origin_code, origin_name) -> codes
print(codes)## # A tibble: 206 × 2
## origin_code origin_name
## <chr> <chr>
## 1 DZA Algeria
## 2 ESP Spain
## 3 ITA Italy
## 4 KHM Cambodia
## 5 LAO Laos
## 6 MAR Morocco
## 7 POL Poland
## 8 PRT Portugal
## 9 SEN Senegal
## 10 TUN Tunisia
## # ... with 196 more rowsPoznámka: volání doplnění implicitních chybějících pozorování pomocí complete() vytvořilo chybějicí pozorování v origin_name a destination_name.
Poznámka 2: správně bychom měli vzít i origin_code a origin_name. V této konkrétní databázi to však není nutné.
Takže, co se nám to schovává pod pokličkou:
codes$origin_code %>% sort## [1] "AFG" "AGO" "ALB" "AND" "ARE" "ARG" "ARM" "ATG" "AUS" "AUT"
## [11] "AZE" "BDI" "BEL" "BEN" "BFA" "BGD" "BGR" "BHR" "BHS" "BIH"
## [21] "BLR" "BLZ" "BMU" "BOL" "BRA" "BRB" "BRN" "BTN" "BWA" "CAF"
## [31] "CAN" "CIV" "CMR" "COD" "COG" "COK" "COL" "COM" "CPV" "CRI"
## [41] "CSK" "CUB" "CYP" "CZE" "DEU" "DJI" "DMA" "DNK" "DOM" "DZA"
## [51] "ECU" "EGY" "ERI" "ESP" "EST" "ETH" "FIN" "FJI" "FRA" "FSM"
## [61] "GAB" "GBR" "GEO" "GHA" "GIN" "GMB" "GNB" "GNQ" "GRC" "GRD"
## [71] "GTM" "GUM" "GUY" "HKG" "HND" "HRV" "HTI" "HUN" "CHE" "CHL"
## [81] "CHN" "IDN" "IND" "IRL" "IRN" "IRQ" "ISL" "ISR" "ITA" "JAM"
## [91] "JOR" "JPN" "KAZ" "KEN" "KGZ" "KHM" "KIR" "KNA" "KOR" "KWT"
## [101] "LAO" "LBN" "LBR" "LBY" "LCA" "LIE" "LKA" "LSO" "LTU" "LUX"
## [111] "LVA" "MAC" "MAR" "MCO" "MDA" "MDG" "MDV" "MEX" "MHL" "MKD"
## [121] "MLI" "MLT" "MMR" "MNE" "MNG" "MOZ" "MRT" "MUS" "MWI" "MYS"
## [131] "NAM" "NER" "NGA" "NIC" "NIU" "NLD" "NOR" "NPL" "NRU" "NZL"
## [141] "OMN" "PAK" "PAN" "PER" "PHL" "PLW" "PNG" "POL" "PRI" "PRK"
## [151] "PRT" "PRY" "PSE" "QAT" "ROU" "RUS" "RWA" "SAU" "SCG" "SDN"
## [161] "SEN" "SGP" "SLB" "SLE" "SLV" "SMR" "SOM" "SRB" "STP" "SUN"
## [171] "SUR" "SVK" "SVN" "SWE" "SWZ" "SYC" "SYR" "TCD" "TGO" "THA"
## [181] "TJK" "TKL" "TKM" "TLS" "TON" "TTO" "TUN" "TUR" "TUV" "TWN"
## [191] "TZA" "UGA" "UKR" "URY" "USA" "UZB" "VCT" "VEN" "VNM" "VUT"
## [201] "WSM" "YEM" "YUCS" "ZAF" "ZMB" "ZWE"Od pohledu se zdá, že by mohlo jít o iso3c kódy. Můžeme to ověřit proti referenční tabulce z balíku countrycode. Je naše doměnka správná?
all(codes$origin_code %in% countrycode::countrycode_data$iso3c)## [1] FALSETak ne. Vypíšeme hříšníky:
codes %>%
filter(!(origin_code %in% countrycode::countrycode_data$iso3c))## # A tibble: 3 × 2
## origin_code origin_name
## <chr> <chr>
## 1 YUCS Former Yugoslavia
## 2 SCG Serbia and Montenegro
## 3 SUN Former USSRKlasika. Problematické jsou historické státy. Občas můžeme mít štěstí. Vypíšeme seznam pozorování, u kterých je problém:
data %>% filter(origin_code=="SUN") %>% filter(!is.na(flow))## # A tibble: 278 × 15
## origin_code destination_code year origin_name destination_name flow
## <chr> <chr> <dbl> <chr> <chr> <dbl>
## 1 SUN USA 1980 Former USSR United States 10543
## 2 SUN USA 1981 Former USSR United States 9223
## 3 SUN USA 1982 Former USSR United States 15462
## 4 SUN USA 1983 Former USSR United States 5214
## 5 SUN USA 1984 Former USSR United States 6088
## 6 SUN GBR 1985 Former USSR United Kingdom 1000
## 7 SUN USA 1985 Former USSR United States 3521
## 8 SUN USA 1986 Former USSR United States 2588
## 9 SUN USA 1987 Former USSR United States 2384
## 10 SUN USA 1988 Former USSR United States 2949
## # ... with 268 more rows, and 9 more variables: flow_women <dbl>,
## # stock_cob <dbl>, stock_cob_women <dbl>, stock_nat <dbl>,
## # stock_nat_women <dbl>, stock_cob_lag1 <dbl>,
## # stock_cob_women_lag1 <dbl>, stock_nat_lag1 <dbl>,
## # stock_nat_women_lag1 <dbl>data %>% filter(origin_code=="YUCS") %>% filter(!is.na(flow))## # A tibble: 318 × 15
## origin_code destination_code year origin_name destination_name
## <chr> <chr> <dbl> <chr> <chr>
## 1 YUCS GBR 1985 Former Yugoslavia United Kingdom
## 2 YUCS NOR 1986 Former Yugoslavia Norway
## 3 YUCS NOR 1987 Former Yugoslavia Norway
## 4 YUCS BEL 1988 Former Yugoslavia Belgium
## 5 YUCS CHE 1988 Former Yugoslavia Switzerland
## 6 YUCS NOR 1988 Former Yugoslavia Norway
## 7 YUCS SWE 1988 Former Yugoslavia Sweden
## 8 YUCS BEL 1989 Former Yugoslavia Belgium
## 9 YUCS DEU 1989 Former Yugoslavia Germany
## 10 YUCS CHE 1989 Former Yugoslavia Switzerland
## # ... with 308 more rows, and 10 more variables: flow <dbl>,
## # flow_women <dbl>, stock_cob <dbl>, stock_cob_women <dbl>,
## # stock_nat <dbl>, stock_nat_women <dbl>, stock_cob_lag1 <dbl>,
## # stock_cob_women_lag1 <dbl>, stock_nat_lag1 <dbl>,
## # stock_nat_women_lag1 <dbl>Zdá se, že tyto dva státy “SUN” a “YUCS” slouží jako sběrné kategorie. Taková pozorování jsou nám k ničemu, protože je logicky nedokážeme spojit s žádnou zemí a jejími charakteristikami. Proto je ze vzorku vyloučíme.
“Serbia and Montenegro” je tvrdší oříšek. Kód “SCG” byl validním kódem, ale po zániku státu se stal rezervovanou zkratkou. Moje zkušenost říká, že pod zkratkami typu “YUG”, “SRB” nebo “SCG” se často skrývají stejné státy. Dost nenápadná noticka v http://www.oecd.org/els/mig/international-migration-outlook-2016-statistical-annex.pdf ukazuje, že to tak +/- bude i tentokrát. Mělo by tedy postačit nechat pozorování pro “SRB”.
data %<>% filter(!(origin_code %in% c("SUN","YUCS","SCG")))Nyní máme mezi codes jen validní iso3c kódy. Ostatní databáze však mohou používat jiná kódování. Obůlíbený trik je doplinit si ostatní kódy do tabulky tak, aby bylo následné spojování tabulek jednodušší. K tabulce data tedy připojíme další kódy pro origin i destination:
data %<>%
select(-origin_name,-destination_name) %>%
rename(
origin_iso3c = origin_code,
destination_iso3c = destination_code
)
countrycode::countrycode_data %>%
select(
name = country.name,
iso3c, iso3n, iso2c, un, wb
) -> codes
codes %>%
set_names(
str_c("origin",names(.),sep="_")
) %>% left_join(data,.) -> data## Joining, by = "origin_iso3c"codes %>%
set_names(
str_c("destination",names(.),sep="_")
) %>% left_join(data,.) -> data## Joining, by = "destination_iso3c"A nakonec přidáme trochu štábní kultury (seřadíme sloupce):
data %<>%
select(
starts_with("origin"),
starts_with("destination"),
everything()
)
print(data)## # A tibble: 255,780 × 23
## origin_iso3c origin_name origin_iso3n origin_iso2c origin_un origin_wb
## <chr> <chr> <int> <chr> <int> <chr>
## 1 AFG Afghanistan 4 AF 4 AFG
## 2 AFG Afghanistan 4 AF 4 AFG
## 3 AFG Afghanistan 4 AF 4 AFG
## 4 AFG Afghanistan 4 AF 4 AFG
## 5 AFG Afghanistan 4 AF 4 AFG
## 6 AFG Afghanistan 4 AF 4 AFG
## 7 AFG Afghanistan 4 AF 4 AFG
## 8 AFG Afghanistan 4 AF 4 AFG
## 9 AFG Afghanistan 4 AF 4 AFG
## 10 AFG Afghanistan 4 AF 4 AFG
## # ... with 255,770 more rows, and 17 more variables:
## # destination_iso3c <chr>, destination_name <chr>,
## # destination_iso3n <int>, destination_iso2c <chr>,
## # destination_un <int>, destination_wb <chr>, year <dbl>, flow <dbl>,
## # flow_women <dbl>, stock_cob <dbl>, stock_cob_women <dbl>,
## # stock_nat <dbl>, stock_nat_women <dbl>, stock_cob_lag1 <dbl>,
## # stock_cob_women_lag1 <dbl>, stock_nat_lag1 <dbl>,
## # stock_nat_women_lag1 <dbl>Přídávání dat z WDI
Modelová specifikace obsahuje tři proměnné, které můžeme získat z databáze WDI:
- velikost populace
- velikost GDP per capita
- míru nezaměstanosti
WDIsearch("unemployment")## indicator
## [1,] "SL.UEM.1524.FE.ZS"
## [2,] "SL.UEM.1524.MA.ZS"
## [3,] "SL.UEM.1524.ZS"
## [4,] "SL.UEM.LTRM.FE.ZS"
## [5,] "SL.UEM.LTRM.MA.ZS"
## [6,] "SL.UEM.LTRM.ZS"
## [7,] "SL.UEM.PRIM.FE.ZS"
## [8,] "SL.UEM.PRIM.MA.ZS"
## [9,] "SL.UEM.PRIM.ZS"
## [10,] "SL.UEM.SECO.FE.ZS"
## [11,] "SL.UEM.SECO.MA.ZS"
## [12,] "SL.UEM.SECO.ZS"
## [13,] "SL.UEM.TERT.FE.ZS"
## [14,] "SL.UEM.TERT.MA.ZS"
## [15,] "SL.UEM.TERT.ZS"
## [16,] "SL.UEM.TOTL.FE.ZS"
## [17,] "SL.UEM.TOTL.MA.ZS"
## [18,] "SL.UEM.TOTL.ZS"
## name
## [1,] "Unemployment, youth female (% of female labor force ages 15-24)"
## [2,] "Unemployment, youth male (% of male labor force ages 15-24)"
## [3,] "Unemployment, youth total (% of total labor force ages 15-24)"
## [4,] "Long-term unemployment, female (% of female unemployment)"
## [5,] "Long-term unemployment, male (% of male unemployment)"
## [6,] "Long-term unemployment (% of total unemployment)"
## [7,] "Unemployment with primary education, female (% of female unemployment)"
## [8,] "Unemployment with primary education, male (% of male unemployment)"
## [9,] "Unemployment with primary education (% of total unemployment)"
## [10,] "Unemployment with secondary education, female (% of female unemployment)"
## [11,] "Unemployment with secondary education, male (% of male unemployment)"
## [12,] "Unemployment with secondary education (% of total unemployment)"
## [13,] "Unemployment with tertiary education, female (% of female unemployment)"
## [14,] "Unemployment with tertiary education, male (% of male unemployment)"
## [15,] "Unemployment with tertiary education (% of total unemployment)"
## [16,] "Unemployment, female (% of female labor force)"
## [17,] "Unemployment, male (% of male labor force)"
## [18,] "Unemployment, total (% of total labor force)"Zajímat nás bude celková nezaměstnanost pro jednotlivá pohlaví – tedy “SL.UEM.TOTL.FE.ZS”, “SL.UEM.TOTL.MA.ZS”, “SL.UEM.TOTL.ZS”. Googlením kódu rychle zjistíme, že se ve skutečnosti jedná o odhady ILO. Sice to není skutečně měření, ale alespoň budeme mít nejkompletnější dostupné pokrytí.
WDIsearch("population, total")## indicator name
## "SP.POP.TOTL" "Population, total"Opět se vyplatí vyggoglit si indikátor “SP.POP.TOTL” a podívávat se na jeho popis. Údaje jsou v tisících.
WDIsearch("gdp per capita")## indicator
## [1,] "GDPPCKD"
## [2,] "GDPPCKN"
## [3,] "NV.AGR.PCAP.KD.ZG"
## [4,] "NY.GDP.PCAP.CD"
## [5,] "NY.GDP.PCAP.KD"
## [6,] "NY.GDP.PCAP.KD.ZG"
## [7,] "NY.GDP.PCAP.KN"
## [8,] "NY.GDP.PCAP.PP.CD"
## [9,] "NY.GDP.PCAP.PP.KD"
## [10,] "NY.GDP.PCAP.PP.KD.ZG"
## [11,] "SE.XPD.PRIM.PC.ZS"
## [12,] "SE.XPD.SECO.PC.ZS"
## [13,] "SE.XPD.TERT.PC.ZS"
## name
## [1,] "GDP per Capita, constant US$, millions"
## [2,] "Real GDP per Capita (real local currency units, various base years)"
## [3,] "Real agricultural GDP per capita growth rate (%)"
## [4,] "GDP per capita (current US$)"
## [5,] "GDP per capita (constant 2000 US$)"
## [6,] "GDP per capita growth (annual %)"
## [7,] "GDP per capita (constant LCU)"
## [8,] "GDP per capita, PPP (current international $)"
## [9,] "GDP per capita, PPP (constant 2005 international $)"
## [10,] "GDP per capita, PPP annual growth (%)"
## [11,] "Expenditure per student, primary (% of GDP per capita)"
## [12,] "Expenditure per student, secondary (% of GDP per capita)"
## [13,] "Expenditure per student, tertiary (% of GDP per capita)"Chceme srovnávat v čase a napříč státy – ergo “NY.GDP.PCAP.PP.KD” je naše volba.
Nyní si můžeme data stáhnout:
WDI(indicator = c("SL.UEM.TOTL.FE.ZS",
"SL.UEM.TOTL.MA.ZS",
"SL.UEM.TOTL.ZS",
"SP.POP.TOTL",
"NY.GDP.PCAP.PP.KD"),
extra = TRUE,
start = 1975,
end = 2015
) %>% as_tibble() -> wdi_data
save(wdi_data, file="data/wdi_data.Rdata")
print(wdi_data)## # A tibble: 10,824 × 15
## iso2c country year SL.UEM.TOTL.FE.ZS SL.UEM.TOTL.MA.ZS SL.UEM.TOTL.ZS
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 AD Andorra 2004 NA NA NA
## 2 AD Andorra 1992 NA NA NA
## 3 AD Andorra 1981 NA NA NA
## 4 AD Andorra 2007 NA NA NA
## 5 AD Andorra 1982 NA NA NA
## 6 AD Andorra 1980 NA NA NA
## 7 AD Andorra 1994 NA NA NA
## 8 AD Andorra 1979 NA NA NA
## 9 AD Andorra 1985 NA NA NA
## 10 AD Andorra 2005 NA NA NA
## # ... with 10,814 more rows, and 9 more variables: SP.POP.TOTL <dbl>,
## # NY.GDP.PCAP.PP.KD <dbl>, iso3c <fctr>, region <fctr>, capital <fctr>,
## # longitude <fctr>, latitude <fctr>, income <fctr>, lending <fctr>API WDI mívá docela dlouhé výpadky. Proto je moudré si data uložit. Stahování aktuálních dat je pěkná věc, ale lehko se může stát, že ojde k jejich revizi a my pak už nikdy neodhalíme, proč se – dejme tomu – změnil regresní koeficient.
Data si transformujeme:
- odstraníme kryptické názvy proměnných
- vybereme jen relevantní sloupce
- převedem velikost populace na jednotky
wdi_data %<>%
select(
iso3c,
year,
country,
unemployment_women = SL.UEM.TOTL.ZS,
unemployment_men = SL.UEM.TOTL.MA.ZS,
unemployment = SL.UEM.TOTL.FE.ZS,
population = SP.POP.TOTL,
GDPpc = NY.GDP.PCAP.PP.KD
) %>%
mutate(
iso3c = as.character(iso3c),
population = population*1000
)
print(wdi_data)## # A tibble: 10,824 × 8
## iso3c year country unemployment_women unemployment_men unemployment
## <chr> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 AND 2004 Andorra NA NA NA
## 2 AND 1992 Andorra NA NA NA
## 3 AND 1981 Andorra NA NA NA
## 4 AND 2007 Andorra NA NA NA
## 5 AND 1982 Andorra NA NA NA
## 6 AND 1980 Andorra NA NA NA
## 7 AND 1994 Andorra NA NA NA
## 8 AND 1979 Andorra NA NA NA
## 9 AND 1985 Andorra NA NA NA
## 10 AND 2005 Andorra NA NA NA
## # ... with 10,814 more rows, and 2 more variables: population <dbl>,
## # GDPpc <dbl>V dalším kroku vytvoříme přidáme lagované proměnné:
wdi_data %<>%
complete(iso3c,year) %>%
group_by(iso3c) %>%
arrange(year) %>%
mutate(
unemployment_lag1 = lag(unemployment),
unemployment_women_lag1 = lag(unemployment_women),
unemployment_men_lag1 = lag(unemployment_men),
GDPpc_lag1 = lag(GDPpc),
population_lag1 = lag(population)
) %>% ungroup()
print(wdi_data)## # A tibble: 9,471 × 13
## iso3c year country unemployment_women unemployment_men
## <chr> <dbl> <chr> <dbl> <dbl>
## 1 ABW 1975 Aruba NA NA
## 2 AFG 1975 Afghanistan NA NA
## 3 AGO 1975 Angola NA NA
## 4 ALB 1975 Albania NA NA
## 5 AND 1975 Andorra NA NA
## 6 ARB 1975 Arab World NA NA
## 7 ARE 1975 United Arab Emirates NA NA
## 8 ARG 1975 Argentina NA NA
## 9 ARM 1975 Armenia NA NA
## 10 ASM 1975 American Samoa NA NA
## # ... with 9,461 more rows, and 8 more variables: unemployment <dbl>,
## # population <dbl>, GDPpc <dbl>, unemployment_lag1 <dbl>,
## # unemployment_women_lag1 <dbl>, unemployment_men_lag1 <dbl>,
## # GDPpc_lag1 <dbl>, population_lag1 <dbl>Teď chceme přidat proměnné z wdi do tabulky s migračními daty. Musíme si ale ověřit, jestli je všechno v pořádku:
wdi_data %>% distinct(iso3c,year) %>% nrow## [1] 9471nrow(wdi_data)## [1] 9471Nyní se podíváme, zde je možné tabulky bezpečně spojit. Balík WDI se tváří, že vrací iso3c kódy, ale jak říkali mudrci od Východu: důvěřuj, ale prověřuj.
codes <- wdi_data$iso3c %>% unique
codes %in% data$origin_iso3c -> loc
all(loc)## [1] FALSEPřes naše nebe přešel mráček. Podíváme se, který kód je problematický:
codes[!loc]## [1] "ABW" "ARB" "ASM" "CSS" "CUW" "CYM" "EMU" "EUU" "GIB" "GRL" "HIC"
## [12] "HPC" "CHI" "IMN" "INX" "LDC" "LIC" "LMC" "LMY" "MAF" "MIC" "MNP"
## [23] "NAC" "NCL" "OED" "OSS" "PSS" "PYF" "SAS" "SSD" "SST" "SXM" "TCA"
## [34] "UMC" "VIR" "WLD"Tak naštěstí se nic neděje. Kódy, které nejsou v migrační databázi jsou některé agregáty (se smyšleným kódem) a mikrozemě. Můžeme tedy bez obav slučovat.
Joinování budeme muset udělat dvakrát – jednou musíme sloučit podle origin a podruhé podle destination. Nejprve ovšem vyhodíme promněnnou country, kterou už dále nepotřebujeme:
wdi_data %<>% select(-country) %>% select(year,everything())
nrow(data)## [1] 255780wdi_data %>%
set_names(
c("year",
str_c("origin_",names(.)[-1])
)
) %>%
left_join(data,.,
by=c("origin_iso3c","year")
) -> data
nrow(data)## [1] 255780wdi_data %>%
set_names(
c("year",
str_c("destination_",names(.)[-1])
)
) %>%
left_join(data,.,
by=c("destination_iso3c","year")
) -> data
nrow(data)## [1] 255780Kontrolní počty řádků po slučování sedí. Ještě zkontrolujeme, zda jsou všechny řádky unikátní.
nrow(distinct(data,origin_iso3c,destination_iso3c,year)) == nrow(data)## [1] TRUEUfff… :-)
Přidávání dat z databáze GeoDist
GeoDist je široce používaná databáze gravitačních proměnných. Je dostupná jako excelovský soubor (nebo dta). Viz: http://www.cepii.fr/CEPII/en/publications/wp/abstract.asp?NoDoc=3877
read_excel("data/dist_cepii.xls") -> geodist
print(geodist)## # A tibble: 50,176 × 14
## iso_o iso_d contig comlang_off comlang_ethno colony comcol curcol col45
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ABW ABW 0 0 0 0 0 0 0
## 2 ABW AFG 0 0 0 0 0 0 0
## 3 ABW AGO 0 0 0 0 0 0 0
## 4 ABW AIA 0 0 1 0 0 0 0
## 5 ABW ALB 0 0 0 0 0 0 0
## 6 ABW AND 0 1 0 0 0 0 0
## 7 ABW ANT 0 1 1 0 1 0 0
## 8 ABW ARE 0 0 0 0 0 0 0
## 9 ABW ARG 0 1 0 0 0 0 0
## 10 ABW ARM 0 0 0 0 0 0 0
## # ... with 50,166 more rows, and 5 more variables: smctry <dbl>,
## # dist <dbl>, distcap <dbl>, distw <chr>, distwces <chr>Hodnoty vGeodist se nemění s časem. Vybereme si jen relevantní sloupce:
- identifikátory zemí
- dummy pro společnou hranici
- dummy pro společný jazyk
- dummy pro společnou koloniální minulost
- dummy pro společnou hranici
- vzdálenost hlavních měst
geodist %<>%
select(
origin_iso3c = iso_o,
destination_iso3c = iso_d,
border = contig,
language = comlang_ethno,
colony,
distance = distcap
)
print(geodist)## # A tibble: 50,176 × 6
## origin_iso3c destination_iso3c border language colony distance
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 ABW ABW 0 0 0 5.225315
## 2 ABW AFG 0 0 0 13257.810000
## 3 ABW AGO 0 0 0 9516.913000
## 4 ABW AIA 0 1 0 983.268200
## 5 ABW ALB 0 0 0 9091.742000
## 6 ABW AND 0 0 0 7572.788000
## 7 ABW ANT 0 1 0 136.384800
## 8 ABW ARE 0 0 0 12735.010000
## 9 ABW ARG 0 0 0 5396.220000
## 10 ABW ARM 0 0 0 11107.780000
## # ... with 50,166 more rowsPodívejme se, zda jsou skutečně všechny řádky unikátní.
nrow(distinct(geodist,origin_iso3c,destination_iso3c)) == nrow(geodist)## [1] TRUEJsou kódy kompatibilní?
codes <- geodist$origin_iso3c %>% unique
codes %in% data$origin_iso3c -> loc
all(loc)## [1] FALSEcodes[!loc]## [1] "ABW" "AIA" "ANT" "CCK" "CXR" "CYM" "ESH" "FLK" "FRO" "GIB" "GLP"
## [12] "GRL" "GUF" "MNP" "MSR" "MTQ" "NCL" "NFK" "PAL" "PCN" "PYF" "REU"
## [23] "ROM" "SHN" "SPM" "TCA" "TMP" "VGB" "WLF" "YUG" "ZAR"Tak jistě. Je patrných několik záměn:
- “ROU” -> “ROM”
- “COD” -> “ZAR”
- “YUG” -> “SRB”
Ostatní země jsou pidizemě nebo závislá území.
Záměny musíme ručně opravit:
geodist$origin_iso3c[geodist$origin_iso3c=="ROM"] <- "ROU"
geodist$origin_iso3c[geodist$origin_iso3c=="ZAR"] <- "COD"
geodist$origin_iso3c[geodist$origin_iso3c=="YUG"] <- "SRB"
geodist$destination_iso3c[geodist$destination_iso3c=="ROM"] <- "ROU"
geodist$destination_iso3c[geodist$destination_iso3c=="ZAR"] <- "COD"
geodist$destination_iso3c[geodist$destination_iso3c=="YUG"] <- "SRB"Nyní může provést spojení:
nrow(data)## [1] 255780data %<>% left_join(.,geodist)## Joining, by = c("origin_iso3c", "destination_iso3c")nrow(data)## [1] 255780Ufff… :-)
Databáze GeoDist podruhé
Databáze GeoDist obsahuje i soubor, ve kterém jsou charakteristiky jednotlivých zemí (ne bilaterálních vztahů). I ty můžeme chtít zohlednit v regersní analýze.
read_excel("data/geo_cepii.xls") -> geodist2
print(geodist2)## # A tibble: 238 × 34
## iso2 iso3 cnum country pays area
## <chr> <chr> <dbl> <chr> <chr> <dbl>
## 1 AW ABW 533 Aruba Aruba 193
## 2 AF AFG 4 Afghanistan Afghanistan 652225
## 3 AO AGO 24 Angola Angola 1246700
## 4 AI AIA 660 Anguilla Anguilla 102
## 5 AL ALB 8 Albania Albanie 28748
## 6 AD AND 20 Andorra Andorre 453
## 7 AN ANT 530 Netherland Antilles Antilles néerlandaises 960
## 8 AE ARE 784 United Arab Emirates Emirats arabes unis 83657
## 9 AR ARG 32 Argentina Argentine 2766889
## 10 AM ARM 51 Armenia Arménie 29800
## # ... with 228 more rows, and 28 more variables: dis_int <dbl>,
## # landlocked <dbl>, continent <chr>, city_en <chr>, city_fr <chr>,
## # lat <dbl>, lon <dbl>, cap <dbl>, maincity <dbl>, citynum <chr>,
## # langoff_1 <chr>, langoff_2 <chr>, langoff_3 <chr>, lang20_1 <chr>,
## # lang20_2 <chr>, lang20_3 <chr>, lang20_4 <chr>, lang9_1 <chr>,
## # lang9_2 <chr>, lang9_3 <chr>, lang9_4 <chr>, colonizer1 <chr>,
## # colonizer2 <chr>, colonizer3 <chr>, colonizer4 <chr>,
## # short_colonizer1 <chr>, short_colonizer2 <chr>, short_colonizer3 <chr>Do tabulky budeme chtít přidat infromaci o tom, zda je země “landlocked”.
geodist2 %<>%
select(iso3,landlocked)Odstraníme problémy kódování detekované u souboru s bilaterálními charakteristikami:
geodist2$iso3[geodist2$iso3=="ROM"] <- "ROU"
geodist2$iso3[geodist2$iso3=="ZAR"] <- "COD"
geodist2$iso3[geodist2$iso3=="YUG"] <- "SRB"Provedeme kontrolu před joinováním:
nrow(distinct(geodist2,iso3)) == nrow(geodist2)## [1] FALSENěkteré pozorování jsou v databázi dvakrát. To by se nemělo stát! Vypíšeme duplicitní řádky:
geodist2 %>%
group_by(iso3) %>%
filter(n()>1)## Source: local data frame [26 x 2]
## Groups: iso3 [13]
##
## iso3 landlocked
## <chr> <dbl>
## 1 AUS 0
## 2 AUS 0
## 3 BEN 0
## 4 BEN 0
## 5 BOL 1
## 6 BOL 1
## 7 BRA 0
## 8 BRA 0
## 9 CAN 0
## 10 CAN 0
## # ... with 16 more rowsŘádky jsou zjevně identické. Problém ve skutečnosti vznikl při změnách hlavních měst. Nicméně tato vlastnost databáze není nikde popsána a při jejím přehlédnutí by vedla k vytvoření duplicitních řádků v tabulce data.
geodist2 %>%
distinct() %>%
left_join(data,.,by=c("origin_iso3c"="iso3")) -> data
print(data)## # A tibble: 255,780 × 48
## origin_iso3c origin_name origin_iso3n origin_iso2c origin_un origin_wb
## <chr> <chr> <int> <chr> <int> <chr>
## 1 AFG Afghanistan 4 AF 4 AFG
## 2 AFG Afghanistan 4 AF 4 AFG
## 3 AFG Afghanistan 4 AF 4 AFG
## 4 AFG Afghanistan 4 AF 4 AFG
## 5 AFG Afghanistan 4 AF 4 AFG
## 6 AFG Afghanistan 4 AF 4 AFG
## 7 AFG Afghanistan 4 AF 4 AFG
## 8 AFG Afghanistan 4 AF 4 AFG
## 9 AFG Afghanistan 4 AF 4 AFG
## 10 AFG Afghanistan 4 AF 4 AFG
## # ... with 255,770 more rows, and 42 more variables:
## # destination_iso3c <chr>, destination_name <chr>,
## # destination_iso3n <int>, destination_iso2c <chr>,
## # destination_un <int>, destination_wb <chr>, year <dbl>, flow <dbl>,
## # flow_women <dbl>, stock_cob <dbl>, stock_cob_women <dbl>,
## # stock_nat <dbl>, stock_nat_women <dbl>, stock_cob_lag1 <dbl>,
## # stock_cob_women_lag1 <dbl>, stock_nat_lag1 <dbl>,
## # stock_nat_women_lag1 <dbl>, origin_unemployment_women <dbl>,
## # origin_unemployment_men <dbl>, origin_unemployment <dbl>,
## # origin_population <dbl>, origin_GDPpc <dbl>,
## # origin_unemployment_lag1 <dbl>, origin_unemployment_women_lag1 <dbl>,
## # origin_unemployment_men_lag1 <dbl>, origin_GDPpc_lag1 <dbl>,
## # origin_population_lag1 <dbl>, destination_unemployment_women <dbl>,
## # destination_unemployment_men <dbl>, destination_unemployment <dbl>,
## # destination_population <dbl>, destination_GDPpc <dbl>,
## # destination_unemployment_lag1 <dbl>,
## # destination_unemployment_women_lag1 <dbl>,
## # destination_unemployment_men_lag1 <dbl>, destination_GDPpc_lag1 <dbl>,
## # destination_population_lag1 <dbl>, border <dbl>, language <dbl>,
## # colony <dbl>, distance <dbl>, landlocked <dbl>