Libraries

library(tidyverse)
library(dplyr)
# library(stargazer)
library(writexl)

# geo
library(sp)
library(spdep)
library(rgdal)
library(maptools)
library(rgeos)
library(sf)
library(spatialreg)
library(modelsummary)
# mass retype()
library(hablar)

# multicollinearity test omcdiag imcdiag
library(mctest)

library(GGally)

#tables
library(xtable)
# library(texreg)
library(remotes)
library(modelsummary)

#spatial panels
library(splm)

library(plm)
Panel_Ru_shape = readOGR(
  'C:/Users/Kolya/Documents/New_studies/GeoAnalisys/DataGeo/Term_paper/Panel_regions.shp')
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\Kolya\Documents\New_studies\GeoAnalisys\DataGeo\Term_paper\Panel_regions.shp", layer: "Panel_regions"
## with 79 features
## It has 121 fields
## Integer64 fields read as strings:  educ_2006 educ_2007 educ_2008 Inv_2008 educ_2009 Inv_2009 educ_2010 Inv_2010 educ_2011 Inv_2011 educ_2012 Inv_2012 educ_2013 Inv_2013 educ_2014 Inv_2014 educ_2015 Inv_2015 educ_2016 Inv_2016 educ_2017 Inv_2017
Panel_Ru <- Panel_Ru_shape@data

# to numeric
Panel_Ru_shape@data[, c(-1)] <-  Panel_Ru_shape@data[, c(-1)] %>% mutate_if(is.character, as.numeric)
# Panel_Ru_shape@data %>% glimpse() #everything is character

#rename
names(Panel_Ru_shape@data)[c(1, 7, 17, 27, 37, 47, 57, 67, 77, 87, 97, 107, 117)] <- c('region',  paste("EmpWF", c(2006:2017), sep="_"))


### LONG vs.WIDE format
long_spatial_data <- Panel_Ru_shape@data %>% pivot_longer(cols = !region,
                                            names_to = c(".value", "year"),
                                          names_pattern = "([A-Za-z]+)_(\\d+)")


#where is NA
which(is.na(long_spatial_data), arr.ind=TRUE)
##      row col
## [1,] 937   3
## [2,] 938   3
## [3,] 939   3
## [4,] 940   3
## [5,] 941   3
pdim(long_spatial_data)
## Balanced Panel: n = 79, T = 12, N = 948
#impute missing data - only 5/948 observations - it's appropriate, because without it- the panel is understood as unbalanced by the model
# long_spatial_data[, 3] = apply(long_spatial_data[, 3], 2, function(x) ifelse(is.na(x), mean(x,na.rm=TRUE), x))

#we can also replace all NA with the mean value (was calculated in General_dataset.rmd)
long_spatial_data[is.na(long_spatial_data)] = 51

#check again
which(is.na(long_spatial_data), arr.ind=TRUE)
##      row col
pdim(long_spatial_data)
## Balanced Panel: n = 79, T = 12, N = 948
#projection
proj4string(Panel_Ru_shape) <- CRS("+init=epsg:4326")
## Warning in proj4string(obj): CRS object has comment, which is lost in output; in tests, see
## https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html
#listw
map_crd <- coordinates(Panel_Ru_shape)
      Points <- SpatialPoints(map_crd)
      Panel_Ru_shape.knn_8 <- knearneigh(Points, k=8)
      K8_nb <- knn2nb(Panel_Ru_shape.knn_8)
      Ru_weight_sq <- nb2listw(K8_nb, style="W")

#Specifications

###Panel models

#Specifications V
grp_V <- log(GRP) ~ V + log(Inv) + educ + Wage
labpr_V <- LabPr ~ V + log(Inv) + educ + Wage
emp_V <- EmpWF ~ V + log(Inv) + educ + Wage
unemp_V <- Unemp ~ V + log(Inv) + educ + Wage

#Specifications RV, UV
grp_RV_UV <- log(GRP) ~ Rv + Uv + log(Inv) + educ + Wage
labpr_RV_UV <- LabPr ~ Rv + Uv + log(Inv) + educ + Wage
emp_RV_UV <- EmpWF ~ Rv + Uv + log(Inv) + educ + Wage
unemp_RV_UV <- Unemp ~ Rv + Uv + log(Inv) + educ + Wage


###Extra models with time-effects

#Specifications V, t

grp_V_t<- log(GRP) ~ V + log(Inv) + educ + Wage + factor(year)
labpr_V_t <- LabPr ~ V + log(Inv) + educ + Wage + factor(year)
emp_V_t <- EmpWF ~ V + log(Inv) + educ + Wage  + factor(year)
unemp_V_t <- Unemp ~ V + log(Inv) + educ + Wage  + factor(year)

#Specifications RV, UV, t
grp_RV_UV_t<- log(GRP) ~ Rv + Uv + log(Inv) + educ + Wage + factor(year)
labpr_RV_UV_t <- LabPr ~ Rv + Uv + log(Inv) + educ + Wage + factor(year)
emp_RV_UV_t <- EmpWF ~ Rv + Uv + log(Inv) + educ + Wage  + factor(year)
unemp_RV_UV_t <- Unemp ~ Rv + Uv + log(Inv) + educ + Wage  + factor(year)

#Models

fe_ind_grp_RV_UV_t_lag <- spml(formula = grp_RV_UV_t, data = long_spatial_data, index = NULL, listw = Ru_weight_sq,
lag = TRUE, spatial.error= "b", model = "within", effect = "individual", method = "eigen", na.action = na.fail )
sum_fe_ind_grp_RV_UV_t_lag <- summary(fe_ind_grp_RV_UV_t_lag)
fe_ind_labpr_RV_UV_t_lag <- spml(formula = labpr_RV_UV_t, data = long_spatial_data, index = NULL, listw = Ru_weight_sq,
lag = TRUE, spatial.error= "b", model = "within", effect = "individual", method = "eigen", na.action = na.fail )
sum_fe_ind_labpr_RV_UV_t_lag <- summary(fe_ind_labpr_RV_UV_t_lag)
fe_ind_emp_RV_UV_t_lag <- spml(formula = emp_RV_UV_t, data = long_spatial_data, index = NULL, listw = Ru_weight_sq,
lag = TRUE, spatial.error= "b", model = "within", effect = "individual", method = "eigen", na.action = na.fail )
sum_fe_ind_emp_RV_UV_t_lag <- summary(fe_ind_emp_RV_UV_t_lag)
fe_ind_unemp_RV_UV_t_lag <- spml(formula = unemp_RV_UV_t, data = long_spatial_data, index = NULL, listw = Ru_weight_sq,
lag = TRUE, spatial.error= "b", model = "within", effect = "individual", method = "eigen", na.action = na.fail )
sum_fe_ind_unemp_RV_UV_t_lag <- summary(fe_ind_unemp_RV_UV_t_lag)
fe_ind_grp_V_t_lag <- spml(formula = grp_V_t, data = long_spatial_data, index = NULL, listw = Ru_weight_sq,
lag = TRUE, spatial.error= "b", model = "within", effect = "individual", method = "eigen", na.action = na.fail )
sum_fe_ind_grp_V_t_lag <- summary(fe_ind_grp_V_t_lag)
fe_ind_labpr_V_t_lag <- spml(formula = labpr_V_t, data = long_spatial_data, index = NULL, listw = Ru_weight_sq,
lag = TRUE, spatial.error= "b", model = "within", effect = "individual", method = "eigen", na.action = na.fail )
sum_fe_ind_labpr_V_t_lag <- summary(fe_ind_labpr_V_t_lag)
fe_ind_emp_V_t_lag <- spml(formula = emp_V_t, data = long_spatial_data, index = NULL, listw = Ru_weight_sq,
lag = TRUE, spatial.error= "b", model = "within", effect = "individual", method = "eigen", na.action = na.fail )
sum_fe_ind_emp_V_t_lag <- summary(fe_ind_emp_V_t_lag)
fe_ind_unemp_V_t_lag <- spml(formula = unemp_V_t, data = long_spatial_data, index = NULL, listw = Ru_weight_sq,
lag = TRUE, spatial.error= "b", model = "within", effect = "individual", method = "eigen", na.action = na.fail )
sum_fe_ind_unemp_V_t_lag <- summary(fe_ind_unemp_V_t_lag)
#summary(fe_ind_unemp_V_t_lag)$CoefTable

#Common table creating

#Rv, UV common table
common_table_spatial_panels_RV_UV <- cbind(sum_fe_ind_grp_RV_UV_t_lag$CoefTable[c(3:7), c(1, 2, 4)],
                                    sum_fe_ind_labpr_RV_UV_t_lag$CoefTable[c(3:7), c(1, 2, 4)],
                                    sum_fe_ind_emp_RV_UV_t_lag$CoefTable[c(3:7), c(1, 2, 4)],
                                    sum_fe_ind_unemp_RV_UV_t_lag$CoefTable[c(3:7), c(1, 2, 4)])

#V table                                   
common_table_spatial_panels_V <- cbind(sum_fe_ind_grp_V_t_lag$CoefTable[c(3:6), c(1, 2, 4)], 
                                    sum_fe_ind_labpr_V_t_lag$CoefTable[c(3:6), c(1, 2, 4)],
                                    sum_fe_ind_emp_V_t_lag$CoefTable[c(3:6), c(1, 2, 4)],
                                    sum_fe_ind_unemp_V_t_lag$CoefTable[c(3:6), c(1, 2, 4)])
#matrix into dataframe
common_table_spatial_panels_RV_UV <- as.data.frame(common_table_spatial_panels_RV_UV)
common_table_spatial_panels_V <- as.data.frame(common_table_spatial_panels_V)

#renaming
names(common_table_spatial_panels_RV_UV) <- c('Est_grp', 'SE_grp', 'Pvalue_grp',
                                              'Est_labpr', 'SE_labpr','Pvalue_labpr',
                                              'Est_emp', 'SE_emp', 'Pvalue_emp',
                                              'Est_unemp', 'SE_unemp', 'Pvalue_unemp')
names(common_table_spatial_panels_V) <- c('Est_grp', 'SE_grp', 'Pvalue_grp',
                                              'Est_labpr', 'SE_labpr','Pvalue_labpr',
                                              'Est_emp', 'SE_emp', 'Pvalue_emp',
                                              'Est_unemp', 'SE_unemp', 'Pvalue_unemp')

common_table_spatial_panels_RV_UV
common_table_spatial_panels_V
library(readxl)
#install.packages('modelsummary')
library(modelsummary)



ti_grp_V <- data.frame(
  term = c('V', 'Log(Inv)', 'Educ', 'Wage'),
  estimate = common_table_spatial_panels_V$Est_grp,
  std.error = common_table_spatial_panels_V$SE_grp,
  p.value = common_table_spatial_panels_V$Pvalue_grp)

ti_labpr_V <- data.frame(
  term = c('V', 'Log(Inv)', 'Educ', 'Wage'),
  estimate = common_table_spatial_panels_V$Est_labpr,
  std.error = common_table_spatial_panels_V$SE_labpr,
  p.value = common_table_spatial_panels_V$Pvalue_labpr)

ti_emp_V <- data.frame(
  term = c('V', 'Log(Inv)', 'Educ', 'Wage'),
  estimate = common_table_spatial_panels_V$Est_emp,
  std.error = common_table_spatial_panels_V$SE_emp,
  p.value = common_table_spatial_panels_V$Pvalue_emp)

ti_unemp_V <- data.frame(
  term = c('V', 'Log(Inv)', 'Educ', 'Wage'),
  estimate = common_table_spatial_panels_V$Est_unemp,
  std.error = common_table_spatial_panels_V$SE_unemp,
  p.value = common_table_spatial_panels_V$Pvalue_unemp)

gl_V <- data.frame(
  Observations = "948",
  Model = "F + T + SP")

mod_grp_V <- list(tidy = ti_grp_V, glance = gl_V)
mod_labpr_V <- list(tidy = ti_labpr_V, glance = gl_V)
mod_emp_V <- list(tidy = ti_emp_V, glance = gl_V)
mod_unemp_V <- list(tidy = ti_unemp_V, glance = gl_V)

#Model class
class(mod_grp_V) <- "modelsummary_list"
class(mod_labpr_V) <- "modelsummary_list"
class(mod_emp_V) <- "modelsummary_list"
class(mod_unemp_V) <- "modelsummary_list"

#created named list
models_V <- list()
models_V[['GRP (1)']] <- mod_grp_V
models_V[['LabPr (2)']] <- mod_labpr_V
models_V[['Emp (3)']] <- mod_emp_V
models_V[['Unemp (4)']] <- mod_unemp_V

#Model
modelsummary(models_V, stars = T, title = 'Model with fixed and time effects, considering spatial weights', fmt = 3) #output = 'latex'
Model with fixed and time effects, considering spatial weights
GRP (1) LabPr (2) Emp (3) Unemp (4)
V 0.020** 0.097 0.001* 0.026
(0.007) (0.285) (0.001) (0.137)
Log(Inv) 0.137*** 2.975*** 0.002* 1.362***
(0.013) (0.569) (0.001) (0.272)
Educ 0.000 0.001 0.000 0.004**
(0.000) (0.003) (0.000) (0.001)
Wage 0.000 0.075** 0.000* 0.010
(0.001) (0.026) (0.000) (0.012)
Observations 948 948 948 948
Model F + T + SP F + T + SP F + T + SP F + T + SP
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 
# modelsummary(mod)
ti_grp_RV_UV <- data.frame(
  term = c('RV', 'UV', 'Log(Inv)', 'Educ', 'Wage'),
  estimate = common_table_spatial_panels_RV_UV$Est_grp,
  std.error = common_table_spatial_panels_RV_UV$SE_grp,
  p.value = common_table_spatial_panels_RV_UV$Pvalue_grp)

ti_labpr_RV_UV <- data.frame(
  term = c('RV', 'UV', 'Log(Inv)', 'Educ', 'Wage'),
  estimate = common_table_spatial_panels_RV_UV$Est_labpr,
  std.error = common_table_spatial_panels_RV_UV$SE_labpr,
  p.value = common_table_spatial_panels_RV_UV$Pvalue_labpr)

ti_emp_RV_UV <- data.frame(
  term = c('RV', 'UV', 'Log(Inv)', 'Educ', 'Wage'),
  estimate = common_table_spatial_panels_RV_UV$Est_emp,
  std.error = common_table_spatial_panels_RV_UV$SE_emp,
  p.value = common_table_spatial_panels_RV_UV$Pvalue_emp)

ti_unemp_RV_UV <- data.frame(
  term = c('RV', 'UV', 'Log(Inv)', 'Educ', 'Wage'),
  estimate = common_table_spatial_panels_RV_UV$Est_unemp,
  std.error = common_table_spatial_panels_RV_UV$SE_unemp,
  p.value = common_table_spatial_panels_RV_UV$Pvalue_unemp)

gl_RV_UV <- data.frame(
  Observations = "948",
  Model = "F + T + SP")

mod_grp_RV_UV <- list(tidy = ti_grp_RV_UV, glance = gl_RV_UV)
mod_labpr_RV_UV <- list(tidy = ti_labpr_RV_UV, glance = gl_RV_UV)
mod_emp_RV_UV <- list(tidy = ti_emp_RV_UV, glance = gl_RV_UV)
mod_unemp_RV_UV <- list(tidy = ti_unemp_RV_UV, glance = gl_RV_UV)

#Model class
class(mod_grp_RV_UV) <- "modelsummary_list"
class(mod_labpr_RV_UV) <- "modelsummary_list"
class(mod_emp_RV_UV) <- "modelsummary_list"
class(mod_unemp_RV_UV) <- "modelsummary_list"

#created named list
models_RV_UV <- list()
models_RV_UV[['GRP (5)']] <- mod_grp_RV_UV
models_RV_UV[['LabPr (6)']] <- mod_labpr_RV_UV
models_RV_UV[['Emp (7)']] <- mod_emp_RV_UV
models_RV_UV[['Unemp (8)']] <- mod_unemp_RV_UV

#Model
modelsummary(models_RV_UV, stars = T, title = 'Model with fixed and time effects, considering spatial weights', fmt = 3)
Model with fixed and time effects, considering spatial weights
GRP (5) LabPr (6) Emp (7) Unemp (8)
RV 0.016 0.421 0.000 0.555+
(0.014) (0.611) (0.001) (0.294)
UV 0.023* 0.103 0.002** 0.385+
(0.010) (0.436) (0.001) (0.209)
Log(Inv) 0.137*** 2.974*** 0.002* 1.357***
(0.013) (0.569) (0.001) (0.270)
Educ 0.000 0.001 0.000 0.004***
(0.000) (0.003) (0.000) (0.001)
Wage 0.000 0.074** 0.000* 0.008
(0.001) (0.026) (0.000) (0.012)
Observations 948 948 948 948
Model F + T + SP F + T + SP F + T + SP F + T + SP
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 
models_total <- list()
models_total[['GRP (1)']] <- mod_grp_V
models_total[['LabPr (2)']] <- mod_labpr_V
models_total[['Emp (3)']] <- mod_emp_V
models_total[['Unemp (4)']] <- mod_unemp_V
models_total[['GRP (5)']] <- mod_grp_RV_UV
models_total[['LabPr (6)']] <- mod_labpr_RV_UV
models_total[['Emp (7)']] <- mod_emp_RV_UV
models_total[['Unemp (8)']] <- mod_unemp_RV_UV

modelsummary(models_total, stars = T, title = 'Model with fixed and time effects, considering spatial weights', fmt = 3, coef_map = c('V', 'RV', 'UV', 'Log(Inv)', 'Educ', 'Wage')) # output = 'latex'
Model with fixed and time effects, considering spatial weights
GRP (1) LabPr (2) Emp (3) Unemp (4) GRP (5) LabPr (6) Emp (7) Unemp (8)
V 0.020** 0.097 0.001* 0.026
(0.007) (0.285) (0.001) (0.137)
RV 0.016 0.421 0.000 0.555+
(0.014) (0.611) (0.001) (0.294)
UV 0.023* 0.103 0.002** 0.385+
(0.010) (0.436) (0.001) (0.209)
Log(Inv) 0.137*** 2.975*** 0.002* 1.362*** 0.137*** 2.974*** 0.002* 1.357***
(0.013) (0.569) (0.001) (0.272) (0.013) (0.569) (0.001) (0.270)
Educ 0.000 0.001 0.000 0.004** 0.000 0.001 0.000 0.004***
(0.000) (0.003) (0.000) (0.001) (0.000) (0.003) (0.000) (0.001)
Wage 0.000 0.075** 0.000* 0.010 0.000 0.074** 0.000* 0.008
(0.001) (0.026) (0.000) (0.012) (0.001) (0.026) (0.000) (0.012)
Observations 948 948 948 948 948 948 948 948
Model F + T + SP F + T + SP F + T + SP F + T + SP F + T + SP F + T + SP F + T + SP F + T + SP
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001