library(spdep)
library(ggplot2)
library(tmap)
library(readxl)
library(RColorBrewer)
library(spatialreg)
library(lmtest)
library(rgdal)
library(car)

setwd('.')

# nacitanie dat
data <- read_excel('data.xlsx')
map <- readOGR(dsn = '.', 'SPH_OBEC')
names(map) <- c('uzkod', 'KOD_LAU1', 'KOD_ORP', 'KOD_OPU', 'KOD_LAU2', 'NAZEV_LAU2')
spatial_tab <- merge(map, data, by = 'uzkod')


# LINEARNY MODEL

OLS <- lm(podielNezam_T ~ podielVS_T + podielCudzincov_T + podielRomov_T 
        + podielVeriacich_T + podielRozvedenych_T  + 
          podielVdovcov_T + + podielZien_T + podiel20az29_T + podielZSbezVZD_T,
        data = spatial_tab)

summary(OLS)

vif(OLS) # variance inflation factor

shapiro.test(sample(residuals(OLS), 5000))
bptest(OLS) # Breusch-Paganov test

# vytvorenie matice vah
coords <- coordinates(map)
IDs <- row.names(as(map, "data.frame"))
nb <- poly2nb(map, queen = TRUE)
listw <- nb2listw(nb, style = 'W')

# Moranov test autokorelacie rezidui linearneho modelu
lm.morantest(OLS, listw = listw)

# skorove a robustne skorove testy pre SEM a SLM model
lm.LMtests(OLS, listw, test = c("LMerr","LMlag","RLMerr","RLMlag")) 

# mapa rezidui linearneho modelu bez premennych podielZien_T a podiel20az29_T
OLS2 <- lm(podielNezam_T ~ podielVS_T + podielCudzincov_T + podielRomov_T 
          + podielVeriacich_T + podielRozvedenych_T  + 
            podielVdovcov_T + podielZSbezVZD_T,
          data = spatial_tab)

spatial_tab$OLS2_rezidua <- residuals(OLS2)
tm_shape(spatial_tab) + tm_fill("OLS2_rezidua", palette = "Reds") +
  tm_borders(alpha=0) + tm_layout(legend.outside = T, frame = F) 



# SLM model
SLM <- lagsarlm(podielNezam_T ~ podielVS_T + podielCudzincov_T + podielRomov_T 
                + podielVeriacich_T + podielRozvedenych_T  + podielVdovcov_T + 
                  podielZSbezVZD_T,
                data = spatial_tab, listw = listw, method = 'LU') 
summary(SLM)

moran.test(residuals(SLM), listw = listw) # Moranov test autokorelacie rezidui

## SEM model
SEM <- errorsarlm(podielNezam_T ~ podielVS_T + podielCudzincov_T + podielRomov_T 
                    + podielVeriacich_T + podielRozvedenych_T  + podielVdovcov_T + podielZSbezVZD_T,
                    data = spatial_tab, listw = listw, method = 'LU') 

summary(SEM)

moran.test(residuals(SEM), listw = listw) # Moran test pre rezidua SEM modelu



# SAC MODEL
SAC <- sacsarlm(podielNezam_T ~ podielVS_T + podielCudzincov_T + podielRomov_T 
                + podielVeriacich_T + podielRozvedenych_T  + podielVdovcov_T + podielZSbezVZD_T,
                data = spatial_tab, listw = listw, method = 'LU') 

summary(SAC)

moran.test(residuals(SAC), listw = listw)


# SPATIAL DURBIN MODEL
SDM <- lagsarlm(podielNezam_T ~ podielVS_T + podielCudzincov_T + podielRomov_T 
                + podielVeriacich_T + podielRozvedenych_T  + podielVdovcov_T + podielZSbezVZD_T,
                data = spatial_tab, listw = listw, type="mixed", method = 'LU', tol.solve=1.0e-20)
summary(SDM)

moran.test(residuals(SDM), listw = listw) 



# VAZENE MODELY

# LINEARNY VAZENY MODEL
OLS_weights <- lm(podielNezam_T ~ podielVS_T + podielCudzincov_T + podielRomov_T 
        + podielVeriacich_T + podielRozvedenych_T  + podielVdovcov_T + podielZSbezVZD_T,
        data = spatial_tab, weights = ekonomicky_aktivni)

summary(OLS_weights)
sigma_hat <- sqrt(sum(residuals(OLS_weights)^2) / 6250) # 6250 = pocet stupnov volnosti
lm.morantest(OLS_weights, listw = listw) # moranov test autokorelacie rezidui


# VAZENY SEM MODEL
SEM_weights <- errorsarlm(podielNezam_T ~ podielVS_T + podielCudzincov_T + podielRomov_T
                          + podielVeriacich_T + podielRozvedenych_T  + podielVdovcov_T + podielZSbezVZD_T,
                          data = spatial_tab, listw = listw, weights = ekonomicky_aktivni, method = 'LU')

summary(SEM_weights)
sigma_hat <- sqrt(sum(residuals(SEM_weights)^2) / 6248)
moran.test(residuals(SEM_weights), listw = listw)

# vysledny vazeny SEM model s maticou vah definovanou na zaklade k = 25 najblizsich susedov
nb2 <- knn2nb(knearneigh(coords, k = 25), row.names = IDs)
listw2 <- nb2listw(nb2, style = 'W')
SEM_weights_vysledny <- errorsarlm(podielNezam_T ~ podielVS_T + podielCudzincov_T + podielRomov_T 
                                    + podielVeriacich_T + podielRozvedenych_T  + podielVdovcov_T + podielZSbezVZD_T,
                                    data = spatial_tab, listw = listw2, weights = ekonomicky_aktivni, method = 'LU')

summary(SEM_weights_vysledny)
sigma_hat <- sqrt(sum(residuals(SEM_weights_vysledny)^2) / 6248)


# mapa rezidui vysledneho vazeneho SEM modelu
spatial_tab$SEM_weights_vysledny_rezidua <- residuals(SEM_weights_vysledny)
tm_shape(spatial_tab) + tm_fill("SEM_weights_vysledny_rezidua", palette = "Reds") +
  tm_borders(alpha=0) + tm_layout(legend.outside = T, frame = F) 

# QQPLOT rezidua vysledneho vazeneho SEM modelu
p <- data.frame(residuals(SEM_weights_vysledny))
p$rezidua_SEM_weights_vysledny <- residuals(SEM_weights_vysledny)
  
ggplot(p, aes(sample = rezidua_SEM_weights_vysledny)) +
  stat_qq()  + theme_minimal() + 
  stat_qq_line()

