[R] Club convergence of regional wage across Indonesian provinces in R

Club convergence of regional wage across Indonesian provinces in R and conditioning factors

library(ConvergenceClubs) # For implementing log-t-tes and identify club convergence library(mFilter) # For filtering series using HP filter library(tidyverse) # Modern data manipulation library(tibble) # Modern data manipulation library(data.table) # Modern data manipulation library(ggpubr) # Provides some easy-to-use functions for creating and customizing ggplot2 library(plyr) # Modern data manipulation library(plotly) # For interactive visualization library(sf) # Encodes spatial vector data # Adjustments knitr::opts_chunk$set(fig.align = 'center') options(warn = -1) options(scipen = 10000) options(repr.plot.width = 6, repr.plot.height = 4)

── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

✔ ggplot2 3.3.4     ✔ purrr   0.3.4
✔ tibble  3.1.2     ✔ dplyr   1.0.7
✔ tidyr   1.1.3     ✔ stringr 1.4.0
✔ readr   1.4.0     ✔ forcats 0.5.1

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()


Attaching package: ‘data.table’


The following objects are masked from ‘package:dplyr’:

    between, first, last


The following object is masked from ‘package:purrr’:

    transpose


------------------------------------------------------------------------------

You have loaded plyr after dplyr - this is likely to cause problems.
If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
library(plyr); library(dplyr)

------------------------------------------------------------------------------


Attaching package: ‘plyr’


The following object is masked from ‘package:ggpubr’:

    mutate


The following objects are masked from ‘package:dplyr’:

    arrange, count, desc, failwith, id, mutate, rename, summarise,
    summarize


The following object is masked from ‘package:purrr’:

    compact



Attaching package: ‘plotly’


The following objects are masked from ‘package:plyr’:

    arrange, mutate, rename, summarise


The following object is masked from ‘package:ggplot2’:

    last_plot


The following object is masked from ‘package:stats’:

    filter


The following object is masked from ‘package:graphics’:

    layout


Linking to GEOS 3.9.0, GDAL 3.2.2, PROJ 7.2.1

Part 1. Identifying club convergence

Load the data

nominal.wage <- read_delim("https://raw.githubusercontent.com/haginta/Project_Club_convergence_wage_Indonesia/main/nominal_wage_A.csv", ";", escape_double = FALSE, trim_ws = TRUE) nominal.wage


── Column specification ────────────────────────────────────────────────────────
cols(
  Region = col_character(),
  Province = col_character(),
  `2005` = col_character(),
  `2006` = col_double(),
  `2007` = col_double(),
  `2008` = col_double(),
  `2009` = col_double(),
  `2010` = col_double(),
  `2011` = col_double(),
  `2012` = col_double(),
  `2013` = col_double(),
  `2014` = col_double(),
  `2015` = col_double(),
  `2016` = col_double(),
  `2017` = col_double(),
  `2018` = col_double(),
  `2019` = col_double(),
  `2020` = col_double()
)

glimpse(nominal.wage) head(nominal.wage)

Rows: 35
Columns: 18
$ Region   <chr> "SUMATERA", "SUMATERA", "SUMATERA", "SUMATERA", "SUMATERA", "…
$ Province <chr> "Aceh", "North Sumatera", "West Sumatera", "Riau", "Riau Isla…
$ `2005`   <chr> "1102.687", "808.136", "1112.456", "1218.72", "0", "1012.558"…
$ `2006`   <dbl> 1196.115, 889.477, 1116.158, 1213.140, 1327.662, 1182.479, 87…
$ `2007`   <dbl> 1260.6, 1009.7, 1064.5, 1233.7, 1544.9, 928.0, 1125.6, 983.9,…
$ `2008`   <dbl> 1348.9, 1093.0, 1152.6, 1307.9, 1632.5, 1002.3, 1211.8, 1067.…
$ `2009`   <dbl> 1420.4, 1301.7, 1430.9, 1359.5, 1860.2, 1244.0, 1400.6, 1195.…
$ `2010`   <dbl> 1456.8, 1344.0, 1488.1, 1422.8, 1897.9, 1300.5, 1441.8, 1222.…
$ `2011`   <dbl> 1497.5, 1425.4, 1640.6, 1739.9, 2178.2, 1349.6, 1549.2, 1418.…
$ `2012`   <dbl> 1567.9, 1477.5, 1716.5, 1817.7, 2277.3, 1395.5, 1587.2, 1475.…
$ `2013`   <dbl> 1660.6, 1579.6, 1770.1, 1929.2, 2360.0, 1527.9, 1699.2, 1585.…
$ `2014`   <dbl> 1731.2, 1676.0, 1881.1, 2134.1, 2744.9, 1764.2, 1867.6, 1596.…
$ `2015`   <dbl> 1744.6, 1741.6, 1899.5, 2099.4, 3411.1, 1860.0, 1931.3, 1807.…
$ `2016`   <dbl> 1918.7, 1921.4, 2013.9, 2307.7, 3503.4, 1976.8, 2048.4, 1958.…
$ `2017`   <dbl> 2402.1, 2297.2, 2473.2, 2545.8, 3446.4, 2204.6, 2497.4, 2113.…
$ `2018`   <dbl> 2192.7, 2152.6, 2394.4, 2280.1, 3559.8, 2093.3, 2427.8, 2050.…
$ `2019`   <dbl> 2400.3, 2349.6, 2564.3, 2688.4, 3646.4, 2313.2, 2494.2, 2180.…
$ `2020`   <dbl> 2549.4, 2436.0, 2769.3, 2776.2, 3972.6, 2433.1, 2548.6, 2256.…

Select complete data only from 2008 onwards and remove "Region" column, and remove the observation of Indonesia

nominal.wage.complete <- nominal.wage %>% select(-c(1,3:5)) %>% subset(Province!="National") nominal.wage.complete

Transform the wage series to log

nominal.wage.log <- log(nominal.wage.complete[,-1]) nominal.wage.log

Remove short-run noise using HP filter

f.nominal.wage.log <- apply(nominal.wage.log, 1, function(x){hpfilter(x, freq=6.25, type="lambda")$trend} ) f.nominal.wage.log <- data.frame(Province = nominal.wage.complete[,1], t(f.nominal.wage.log), stringsAsFactors=FALSE ) colnames(f.nominal.wage.log) <- colnames(nominal.wage.complete)

Inspect the filtered data

head(f.nominal.wage.log)

Run the log-t-test

H.nominal <- computeH(f.nominal.wage.log[,-1], quantity = "H") round(estimateMod(H.nominal, time_trim=0.333, HACmethod = "FQSB"), 3)

Identify convergence clubs

clubs.nominal <- findClubs(f.nominal.wage.log, dataCols=2:14, unit_names = 1, refCol=14, time_trim=0.333, cstar=0, HACmethod = 'FQSB')

summary(clubs.nominal)

Number of convergence clubs: 2
Number of divergent units: 0

        | # of units  | beta     | std.err  | tvalue   | cstar  
-------- ------------- ---------- ---------- ---------- --------
 club1  |   12        |  -0.378  |   0.087  |  -4.316  |   0    
 club2  |   22        |  -0.058  |   0.104  |  -0.558  |   0

print(clubs.nominal)

======================================================================== 
club 1
------------------------------------------------------------------------ 
Banten, Papua, East Kalimantan, Riau Islands, North Kalimantan, 
Jakarta, West Java, North Sulawesi, Central Kalimantan, West Papua, 
Bali, South Sulawesi

beta:     -0.3775 
std.err:   0.0875 
tvalue:   -4.3162 
pvalue:    0 
cstar:     0 

======================================================================== 
club 2
------------------------------------------------------------------------ 
South Kalimantan, Bangka Belitung, North Maluku, Maluku, West Sumatera, 
Riau, Southeast Sulawesi, Bengkulu, Aceh, West Kalimantan, Central 
Sulawesi, Gorontalo, North Sumatera, West Nusa Tenggara, East Java, 
Jambi, Lampung, Yogyakarta, Central Java, East Nusa Tenggara, South 
Sumatera, West Sulawesi

beta:     -0.0579 
std.err:   0.1037 
tvalue:   -0.5585 
pvalue:    0.2883 
cstar:     0

plot(clubs.nominal)

plot(clubs.nominal, clubs=NULL, avgTP = TRUE, legend=TRUE)

Merge clubs

mclubs.nominal <- mergeClubs(clubs.nominal, mergeMethod='PS') summary(mclubs.nominal)

Number of convergence clubs: 2
Number of divergent units: 0

        | merged clubs  | # of units  | beta     | std.err  | tvalue   
-------- --------------- ------------- ---------- ---------- ----------
 club1  | clubs: 1      |   12        |  -0.378  |   0.087  |  -4.316  
 club2  | clubs: 2      |   22        |  -0.058  |   0.104  |  -0.558

print(mclubs.nominal)

======================================================================== 
club 1
------------------------------------------------------------------------ 
Banten, Papua, East Kalimantan, Riau Islands, North Kalimantan, 
Jakarta, West Java, North Sulawesi, Central Kalimantan, West Papua, 
Bali, South Sulawesi

beta:     -0.3775 
std.err:   0.0875 
tvalue:   -4.3162 
pvalue:    0 

======================================================================== 
club 2
------------------------------------------------------------------------ 
South Kalimantan, Bangka Belitung, North Maluku, Maluku, West Sumatera, 
Riau, Southeast Sulawesi, Bengkulu, Aceh, West Kalimantan, Central 
Sulawesi, Gorontalo, North Sumatera, West Nusa Tenggara, East Java, 
Jambi, Lampung, Yogyakarta, Central Java, East Nusa Tenggara, South 
Sumatera, West Sulawesi

beta:     -0.0579 
std.err:   0.1037 
tvalue:   -0.5585 
pvalue:    0.2883

plot(mclubs.nominal)

plot(mclubs.nominal, clubs=NULL, avgTP = TRUE, legend=TRUE)

Part 2. Generate relative values for transition path plot

Convert "m.clubs.nom" from list to table with name "df.mclubs.nom"

table.mclubs.nom <- map(mclubs.nominal, as.data.table) df.mclubs.nom <- rbindlist(table.mclubs.nom, fill = T, idcol = T) df.mclubs.nom

colnames(df.mclubs.nom)[c(1,5)] <- c("Club","Province") df.mclubs.nom

Merge the "df.mclubs.nom" with "f.nominal.wage.log", which is the log of wage after HP filtered and save as dataframe with name "nominal.wage.club"

nominal.wage.club <- as.data.frame( inner_join(f.nominal.wage.log, df.mclubs.nom, by="Province") %>% select(-contains(c("clubs","id","model"))) ) nominal.wage.club

Generate relative values of wage for each year by dividing the wage of province i at time t to the mean wage of all provinces at each time.

relative <- list() for(a in 2:14) { relative[[a]] <- data.frame(nominal.wage.club$Province, nominal.wage.club[,a]/mean(nominal.wage.club[,a])) colnames(relative[[a]])[1]<- "Province" colnames(relative[[a]])[2]<- paste("rel",colnames(nominal.wage.club)[a],sep="_") } for (x in 2:13) { relative[[x+1]]<-left_join(relative[[x]],relative[[x+1]], by="Province") }

table.relative <- map(relative, as.data.table) df.relative <- rbindlist(table.relative, fill = T, idcol = T) %>% drop_na() df.relative <- join(df.relative,nominal.wage.club,by=c('Province')) %>% select(Province,contains("rel"),Club) df.relative <- as.data.frame(df.relative)

colnames(df.relative) <- gsub("rel_","",colnames(df.relative)) df.relative

Transform "df.relative" to long format dataframe so the ggplot would be happy to plot.

df.relative.long <- df.relative %>% pivot_longer(-c(Province,Club), names_to = "Time", values_to="Rel_Wage") df.relative.long

Calculate the mean by clubs

df.relative.path <- aggregate(Rel_Wage ~ Club + Time, df.relative.long, mean) %>% arrange(Club) df.relative.path

Plot all clubs

path_all <- df.relative.path %>% ggplot(aes(x=Time,y=Rel_Wage, group=Club, col=Club)) + geom_line() + labs(title = "Transition path of all clubs", caption = "Rel_Wage = Relative Wage") + theme_bw() path_all

Plot transition path within club

Generate relative values Club 1

club1 <- nominal.wage.club %>% filter(Club == "club1")

relative.club1 <- list() for(a in 2:14) { relative.club1[[a]] <- data.frame(club1$Province, club1[,a]/mean(club1[,a])) colnames(relative.club1[[a]])[1]<- "Province" colnames(relative.club1[[a]])[2]<- paste("rel",colnames(club1)[a],sep="_") } for (x in 2:13) { relative.club1[[x+1]]<-left_join(relative.club1[[x]],relative.club1[[x+1]], by="Province") }

table.relative.club1 <- map(relative.club1, as.data.table) df.relative.club1 <- rbindlist(table.relative.club1, fill = T, idcol = T) %>% drop_na() library(plyr) df.relative.club1 <- join(df.relative.club1,club1,by=c('Province')) %>% select(Province,contains("rel"),Club) df.relative.club1 <- as.data.frame(df.relative.club1)

colnames(df.relative.club1) <- gsub("rel_","",colnames(df.relative.club1)) df.relative.club1

df.relative.club1.long <- df.relative.club1 %>% pivot_longer(-c(Province,Club), names_to = "Time", values_to="Rel_Wage") df.relative.club1.long

df.relative.club1.path <- aggregate(Rel_Wage ~ Province + Time, df.relative.club1.long, mean) %>% arrange(Province)

Plot transition path of regions in Club 1

path_club1 <- df.relative.club1.path %>% ggplot(aes(x=Time,y=Rel_Wage, group=Province, col=Province)) + geom_line() + labs(title = "Transition path of regions in Club 1", caption = "Rel_Wage = Relative Wage") + theme_bw() + theme(legend.title = element_text(size = 8), legend.text = element_text(size = 8), axis.text.x = element_text(angle = 90)) path_club1

ggplotly(path_club1)

Generate relative values Club 2

club2 <- nominal.wage.club %>% filter(Club == "club2")

relative.club2 <- list() for(a in 2:14) { relative.club2[[a]] <- data.frame(club2$Province, club2[,a]/mean(club2[,a])) colnames(relative.club2[[a]])[1]<- "Province" colnames(relative.club2[[a]])[2]<- paste("rel",colnames(club2)[a],sep="_") } for (x in 2:13) { relative.club2[[x+1]]<-left_join(relative.club2[[x]],relative.club2[[x+1]], by="Province") }

table.relative.club2 <- map(relative.club2, as.data.table) df.relative.club2 <- rbindlist(table.relative.club2, fill = T, idcol = T) %>% drop_na() library(plyr) df.relative.club2 <- join(df.relative.club2,club2,by=c('Province')) %>% select(Province,contains("rel"),Club) df.relative.club2 <- as.data.frame(df.relative.club2)

colnames(df.relative.club2) <- gsub("rel_","",colnames(df.relative.club2)) df.relative.club2

df.relative.club2.long <- df.relative.club2 %>% pivot_longer(-c(Province,Club), names_to = "Time", values_to="Rel_Wage") df.relative.club2.long

df.relative.club2.path <- aggregate(Rel_Wage ~ Province + Time, df.relative.club2.long, mean) %>% arrange(Province)

Plot transition path of regions in Club 2

path_club2 <- df.relative.club2.path %>% ggplot(aes(x=Time,y=Rel_Wage, group=Province, col=Province)) + geom_line() + labs(title = "Transition path of regions in Club 2", caption = "Rel_Wage = Relative Wage") + theme_bw() + theme(legend.title = element_text(size = 8), legend.text = element_text(size = 8), axis.text.x = element_text(angle = 90)) path_club2

ggplotly(path_club2)

Part 3. Analyzing conditioning factors of convergence clubs

ologit <- read_delim("https://raw.githubusercontent.com/haginta/Project_Club_convergence_wage_Indonesia/main/wage_convergence_ologit.csv", ";", escape_double = FALSE, trim_ws = TRUE) ologit


── Column specification ────────────────────────────────────────────────────────
cols(
  Region = col_character(),
  Province = col_character(),
  club = col_double(),
  clubs = col_double(),
  club_real = col_double(),
  clubs_real = col_double(),
  real_wage_2008 = col_double(),
  pmtb_gdrp = col_double(),
  tpak = col_double(),
  emp_manu = col_double(),
  grdp = col_double(),
  grdp_g = col_double()
)

Interactive data exploration

explore <- ggplot(ologit, aes(Region, grdp_g, label=Province, color=Region)) + geom_point() + xlab("Province") + ylab("GDP growth 2008-2020 (%, yoy)") + theme_bw() ggplotly(explore)

ggbarplot(ologit, x = "Region", y = "grdp_g", add = c("mean_se", "point"), color = "Region", alpha = 0.5) + xlab("Region") + ylab("GDP growth 2020 (%, yoy)") + ggrepel::geom_text_repel(aes(label = Province), size = 2) + theme(axis.text.y = element_text(size = 8), axis.text.x = element_text(size = 8, angle = 90, vjust = 0.5, hjust=1),legend.position = "none")

explore_box <- ggbarplot(ologit, x = "Region", y = "grdp_g", color="Region", add = c("mean_se","point")) + geom_point(aes(label=Province, color=Region)) + xlab("Region") + ylab("GDP growth 2008-2020 (%, yoy)") + theme(axis.text.y = element_text(size = 8), axis.text.x = element_text(size = 8, angle = 90, vjust = 0.5, hjust=1),legend.position = "none") ggplotly(explore_box)

scatter <- ggplot(ologit, aes(x=emp_manu, y=grdp_g, size=grdp)) + geom_point(aes(label=Province)) + geom_smooth(method = "lm", se = F) + labs(x = "Manufacturing employment share", y="Average GDP growth 2008-2020") ggplotly(scatter)

`geom_smooth()` using formula 'y ~ x'

Load shape file

IDNmap <- read_sf("https://gist.githubusercontent.com/haginta/22a1d27cf8ae440530758a5715f129c4/raw/1506205198c4337a6dbb6430b070fc24ba80e940/IDN_poly34_simple.JSON") colnames(IDNmap)

colnames(IDNmap)[4] <- 'Province' colnames(IDNmap)

ggplot(IDNmap) + geom_sf() + theme_bw() + ggtitle("Provinces level boundaries map of Indonesia", subtitle = paste0("(", length(unique(IDNmap$POLY_ID)), " provinces)")) + labs(x = "Longitude", y="Latitude")

Merge non-spatial and spatial data

df_map <- merge(IDNmap, ologit, by=c("Province")) df_map <- st_make_valid(df_map) glimpse(df_map)

Rows: 34
Columns: 16
$ Province       <chr> "Aceh", "Bali", "Bangka Belitung", "Banten", "Bengkulu"…
$ POLY_ID        <int> 1, 2, 3, 4, 5, 11, 15, 28, 12, 16, 23, 6, 8, 9, 19, 21,…
$ ENGTYPE_1      <chr> "Autonomous Province", "Province", "Province", "Provinc…
$ prov           <chr> "ACH", "BLI", "BBL", "BTN", "BKL", "CJA", "CKL", "CSI",…
$ Region         <chr> "SUMATERA", "BALNUSTRA", "SUMATERA", "JAWA", "SUMATERA"…
$ club           <dbl> 2, 1, 2, 1, 2, 2, 1, 2, 2, 1, 2, 2, 1, 2, 2, 2, 1, 2, 1…
$ clubs          <dbl> 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1…
$ club_real      <dbl> 3, 2, 3, 1, 3, 3, 2, 3, 3, 2, 3, 3, 1, 3, 3, 3, 2, 3, 2…
$ clubs_real     <dbl> 1, 2, 1, 3, 1, 1, 2, 1, 1, 2, 1, 1, 3, 1, 1, 1, 2, 1, 2…
$ real_wage_2008 <dbl> 1006.5969, 1001.0670, 784.6343, 987.9893, 955.0992, 668…
$ pmtb_gdrp      <dbl> 31.96295, 32.41899, 21.76467, 30.04630, 42.77610, 28.91…
$ tpak           <dbl> 64.75956, 77.47804, 68.08067, 66.40072, 73.76226, 70.46…
$ emp_manu       <dbl> 5.058010, 13.765413, 5.960861, 23.902952, 3.748330, 19.…
$ grdp           <dbl> 28744.308, 31623.104, 11249.875, 90272.712, 9307.114, 1…
$ grdp_g         <dbl> 2.975360, 6.328427, 4.749206, 5.986227, 5.621255, 5.307…
$ geometry       <MULTIPOLYGON [°]> MULTIPOLYGON (((97.11918 2...., MULTIPOLYG…

Show interactive map

club_map <- ggplot(df_map) + geom_sf(aes(fill=factor(club_real), label=Province), lwd = 0.05) + scale_fill_manual(values = c("red", "yellow", "blue"), name = "Clubs") + theme_bw() + labs(x = "Longitude", y="Latitude", title = "Club convergence of regional wage in Indonesia, 2008-2020") club_inter <- ggplotly(club_map) club_inter

gdp_map <- ggplot(df_map) + geom_sf(aes(fill=grdp_g, label=Province), lwd = 0.05) + scale_fill_viridis_c(name = "GDP growth") + theme_bw() + labs(x = "Longitude", y="Latitude", title = "GDP growth of 34 provinces in Indonesia, average 2008-2020") gdp_inter <- ggplotly(gdp_map) gdp_inter

Perform ordered logit model

ologit$clubs <- factor(ologit$clubs)

ologit <- ologit %>% mutate(lngrdp = log(grdp))

glimpse(ologit)

Rows: 34
Columns: 13
$ Region         <chr> "SUMATERA", "SUMATERA", "SUMATERA", "SUMATERA", "SUMATE…
$ Province       <chr> "Aceh", "North Sumatra", "West Sumatra", "Riau", "Jambi…
$ club           <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 2, 2…
$ clubs          <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0…
$ club_real      <dbl> 3, 3, 3, 3, 3, 3, 3, 3, 3, 1, 1, 2, 3, 3, 3, 1, 2, 3, 3…
$ clubs_real     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 2, 1, 1, 1, 3, 2, 1, 1…
$ real_wage_2008 <dbl> 1006.5969, 872.7189, 885.5141, 1040.4209, 755.7823, 818…
$ pmtb_gdrp      <dbl> 31.96295, 29.88666, 29.62215, 28.46564, 22.44053, 37.71…
$ tpak           <dbl> 64.75956, 71.45083, 69.34770, 66.96679, 68.34050, 71.18…
$ emp_manu       <dbl> 5.058010, 8.072072, 7.247634, 6.319503, 3.840253, 5.372…
$ grdp           <dbl> 28744.308, 108075.305, 34395.038, 111631.432, 30224.747…
$ grdp_g         <dbl> 2.975360, 5.576577, 5.655802, 2.771346, 5.712246, 5.562…
$ lngrdp         <dbl> 10.266195, 11.590584, 10.445668, 11.622958, 10.316416, …

mylogit <- glm(clubs ~ real_wage_2008 + emp_manu + pmtb_gdrp + + tpak + lngrdp, data = ologit, family = "binomial") summary(mylogit)