# Libraries

library(spdep)

library(spatialreg)

library(rgdal)

library(rgeos)

library(geoR)

library(gvlma)

# Data & Map

file.shape <- readOGR(dsn = "C:...", layer = "X")

plot(file.shape)

names(file.shape)

## K neighbors

coords <- coordinates(file.shape)

NEIGH <- row.names(as(file.shape, "data.frame"))

fileshape_kn10 <- knn2nb(knearneigh(coords, k = 10), row.names = NEIGH)

fileshape10.listw=nb2listw(fileshape_kn10)

listw1= fileshape10.listw

## Check distances

dist <- unlist(nbdists(fileshape_kn10, coords))

summary(dist)

## Weight matrix (row standard y binary)

cursomodelos_kn10_w <- nb2listw(fileshape_kn10, zero.policy = T)

# Descriptives

summary(Y)

summary(X)

plot(X, Y)

abline(lm(Y~X), col="red")

lines(lowess(Y,X), col="blue")

## Global Moran I

moran.test(Y, listw = fileshape10.listw, zero.policy = T)

# OLS: y=XB+e

reg.eq1=Y~X

ols=lm(reg.eq1,data=file.shape)

summary(ols)

## Diagnostics

dols <- gvlma(ols)

summary(dols)

# SAR: y=pWy+XB+e

sar=lagsarlm(reg.eq1,data=file.shape, listw1)

summary(sar)

(est <- cbind(Estimate = coef(sar), confint(sar)))

impacts(sar,listw=listw1)

summary(impacts(sar,listw=listw1,R=500),zstats=TRUE)

sar$deviance/sar$null.deviance

mean(fitted(sar))

hist(fitted(sar), main="SAR")

anova(sar)

logLik(sar)

LR1.sarlm(sar)

Wald1.sarlm(sar)

bptest.sarlm(sar,studentize=TRUE)

# SEM: y=XB+u, u=LWu+e

sem=errorsarlm(reg.eq1,data=file.shape, listw1)

summary(sem)

(est <- cbind(Estimate = coef(sem), confint(sem)))

sem$deviance/sem$null.deviance

mean(fitted(sem))

hist(fitted(sem), main="SEM")

bptest.sarlm(sem,studentize=TRUE)

anova(sem)

logLik(sem)

LR1.sarlm(sem)

Wald1.sarlm(sem)

bptest.sarlm(sem,studentize=TRUE)

Hausman.test(sem)