#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>>>>>>>>>>>>>ROZDZIAŁ 6>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>>>>>>> MODELOWANIE PRZESTRZENNE>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>


# 6.3

######################
###### Przykład 13 ######
######################
# regresja liniowa
euro.lm <-lm(PZ2~PZ4+PZ5+PZ6, data=euro.dane)
summary(euro.lm)

# statystyka Morana dla reszt, test jednostronny
euro.reg.moran1 <-lm.morantest(euro.lm, euro.listw)
print(euro.reg.moran1)

# tworzenie obiektu reszt modelu
res<-residuals(euro.lm)
is.data.frame(res)
res.frame<-as.data.frame(res)
is.data.frame(res.frame)
print(res.frame)

# rysowanie reszt
brks<-c(round(min(res.frame$res),digits=2),0,round(max(res.frame$res),digits=2))
cols<-c("grey80","grey45")
plot(euro.poly, col=cols[findInterval(res.frame$res, brks)],forcefill=FALSE)
legend(x=c(-600000,-600000), y=c(-400000,-400000), legend=c("reszty ujemne", "reszty dodatnie"), fill=cols, bty="n")

# drukowanie reszt wg regionów
id<-order(euro.dane$REG)
printCoefmat(data.frame(res.frame[id,], row.names=euro.dane$REG[id]),check.names=FALSE)

# test join-count dla reszt
resj<-cut(res.frame$res, breaks=c(min(res.frame$res),0 ,max(res.frame$res)), labels=c("ujemne", "dodatnie"), include.lowest=TRUE)
resj.f<-factor(resj)
joincount.test(resj.f, euro.listw)
joincount.mc(resj.f, euro.listw, 99)

######################

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# 6.4

######################
###### Przykład 14 ######
######################
euro.testLM<-lm.LMtests(euro.lm, euro.listw, test="all")
print(euro.testLM)

######################

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# 6.7
# modele przestrzenne

euro.lag<-lagsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.listw, method="eigen", quiet=TRUE)
euro.error<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.listw, method="eigen", quiet=TRUE)
euro.lm<-lm(PZ1~PZ4+PZ5+PZ6, data=euro.dane)

# kryteria informacyjne dla pojedynczego modelu
AIC(euro.lag) #kryterium informacyjne Akaike
AIC(euro.lag, k=log(nrow(euro.dane))) # kryterium informacyjne bayesowskie
logLik(euro.lag) # kyterium logLik


# kryteria informacyjne dla kilku modeli
AIC(euro.lag, euro.error, euro.lm) # kryterium informacyjne AIC
AIC(euro.lag, euro.error, euro.lm, k=log(nrow(euro.dane))) # kryterium informacyjne BIC
anova.sarlm(euro.lag, euro.error, euro.lm) # porównanie kryt.info.AIC oraz logLik

# łączenie wyników w jeden
out1<-anova.sarlm(euro.lag, euro.error, euro.lm)
out2<-AIC(euro.lag, euro.error, euro.lm, k=log(nrow(euro.dane)))
colnames(out2)<-c("df","BIC") # nadawanie nazw wierszom
out3<-cbind(out1, out2) # połączenie obiektów
out3 # wyświetlenie wyniku

######################
###### Przykład 15 ######
######################
# model opóźnienia przestrzennego
euro.lag<-lagsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.listw, method="eigen", quiet=TRUE, tol.solve=1.0e-10)
summary(euro.lag) # w modelu może być użyta metoda "eigen" lub "sparse"

# model błędu przestrzennego
euro.error<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.listw, method="eigen", quiet=TRUE)
summary(euro.error)

######################

######################
###### Przykład 16 ######
######################
# tworzenie macierzy wag
euro.W.listw<-nb2listw(euro.nb, style="W")
euro.B.listw<-nb2listw(euro.nb, style="B")
euro.S.listw<-nb2listw(euro.nb, style="S")

coord<-as.matrix(euro.dane[,7:8])
euro.knn<- knearneigh(coord, k=10) # powstaje obiekt klasy knn
euro.k.nb<- knn2nb(euro.knn)
euro.k.sym.nb<-make.sym.nb(euro.k.nb)
euro.k.sym.listw<-nb2listw(euro.k.sym.nb)

maxidist<-900.90*200 #900.90*...km
conti200<- dnearneigh(coord, 0, maxidist, row.names = rownames(euro.dane), lonlat = FALSE)
euro.d.listw<-nb2listw(conti200)

dist <- nbdists(euro.nb, coord) # powstaje obiekt klasy nbdist, który łatwo przerobić na klasę nb
dist1 <- lapply(dist, function(x) 1/x) # obiekt klasy nb
euro.dist.listw<-nb2listw(euro.nb, glist=dist1)


# estymacja modeli przy użyciu wszystkich macierz wag
meuro.W.listw<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.W.listw, method="eigen")
meuro.B.listw<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.B.listw, method="eigen")
meuro.S.listw<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.S.listw, method="eigen")
meuro.k.listw<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.k.sym.listw, method="eigen")
meuro.d.listw<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.d.listw, method="eigen")
meuro.dist.listw<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.dist.listw, method="eigen")

# połączenie wszystkich kryteriów informacyjnych w jednen wydruk
out11<-anova.sarlm(meuro.W.listw, meuro.B.listw, meuro.S.listw, meuro.k.listw, meuro.d.listw, meuro.dist.listw)
out12<-AIC(meuro.W.listw, meuro.B.listw, meuro.S.listw, meuro.k.listw, meuro.d.listw, meuro.dist.listw, k=log(nrow(euro.dane)))
colnames(out12)<-c("df","BIC") # nadawanie nazw wierszom
out13<-cbind(out11, out12) # połączenie obiektów
out13 # wyświetlenie wyniku

# model z macierzą wag według k najbliższych sąsiadów
meuro.k.listw<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.k.sym.listw, method="eigen")
summary(meuro.k.listw)

# model liniowy
meuro.lm<-lm(PZ1~PZ4+PZ5+PZ6, data=euro.dane)
summary(meuro.lm)

######################

######################
###### Przykład 17 ######
######################

# Przykład - siła dyfuzji

coords<-as.matrix(euro.dane[,7:8])

# macierze wag wg sąsiedztwa w promieniu d km

maxidist<-900.90*50 #900.90*...km
conti50<- dnearneigh(coords, 0, maxidist, row.names = rownames(euro.dane), lonlat = FALSE)
euro.d50.listw<-nb2listw(conti50, zero.policy=TRUE)

maxidist<-900.90*100 #900.90*...km
conti100<- dnearneigh(coords, 0, maxidist, row.names = rownames(euro.dane), lonlat = FALSE)
euro.d100.listw<-nb2listw(conti100, zero.policy=TRUE)

maxidist<-900.90*150 #900.90*...km
conti150<- dnearneigh(coords, 0, maxidist, row.names = rownames(euro.dane), lonlat = FALSE)
euro.d150.listw<-nb2listw(conti150, zero.policy=TRUE)

maxidist<-900.90*200 #900.90*...km
conti200<- dnearneigh(coords, 0, maxidist, row.names = rownames(euro.dane), lonlat = FALSE)
euro.d200.listw<-nb2listw(conti200)

maxidist<-900.90*250 #900.90*...km
conti250<- dnearneigh(coords, 0, maxidist, row.names = rownames(euro.dane), lonlat = FALSE)
euro.d250.listw<-nb2listw(conti250)

maxidist<-900.90*300 #900.90*...km
conti300<- dnearneigh(coords, 0, maxidist, row.names = rownames(euro.dane), lonlat = FALSE)
euro.d300.listw<-nb2listw(conti300)

maxidist<-900.90*350 #900.90*...km
conti350<- dnearneigh(coords, 0, maxidist, row.names = rownames(euro.dane), lonlat = FALSE)
euro.d350.listw<-nb2listw(conti350)

maxidist<-900.90*400 #900.90*...km
conti400<- dnearneigh(coords, 0, maxidist, row.names = rownames(euro.dane), lonlat = FALSE)
euro.d400.listw<-nb2listw(conti400)

maxidist<-900.90*450 #900.90*...km
conti450<- dnearneigh(coords, 0, maxidist, row.names = rownames(euro.dane), lonlat = FALSE)
euro.d450.listw<-nb2listw(conti450)

maxidist<-900.90*500 #900.90*...km
conti500<- dnearneigh(coords, 0, maxidist, row.names = rownames(euro.dane), lonlat = FALSE)
euro.d500.listw<-nb2listw(conti500)

maxidist<-900.90*550 #900.90*...km
conti550<- dnearneigh(coords, 0, maxidist, row.names = rownames(euro.dane), lonlat = FALSE)
euro.d550.listw<-nb2listw(conti550)

# estysmacja modeli błędu przestrzennego
meuro.d50<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.d50.listw, method="eigen", zero.policy=TRUE)
meuro.d100<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.d100.listw, method="eigen", zero.policy=TRUE)
meuro.d150<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.d150.listw, method="eigen", zero.policy=TRUE)
meuro.d200<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.d200.listw, method="eigen")
meuro.d250<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.d250.listw, method="eigen")
meuro.d300<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.d300.listw, method="eigen")
meuro.d350<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.d350.listw, method="eigen")
meuro.d400<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.d400.listw, method="eigen")
meuro.d450<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.d450.listw, method="eigen")
meuro.d500<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.d500.listw, method="eigen")
meuro.d550<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.d550.listw, method="eigen")

# podsumowanie lambda i kryteriów informacyjnych
out.lambda<-rbind(meuro.d50$lambda, meuro.d100$lambda, meuro.d150$lambda, meuro.d200$lambda, meuro.d250$lambda, meuro.d300$lambda, meuro.d350$lambda, meuro.d400$lambda, meuro.d450$lambda, meuro.d500$lambda,meuro.d550$lambda)
out.lambda
out.AIC<-anova.sarlm(meuro.d50, meuro.d100, meuro.d150, meuro.d200, meuro.d250, meuro.d300, meuro.d350, meuro.d400, meuro.d450, meuro.d500, meuro.d550)
out.AIC
out.BIC<-AIC(meuro.d50, meuro.d100, meuro.d150, meuro.d200, meuro.d250, meuro.d300, meuro.d350, meuro.d400, meuro.d450, meuro.d500, meuro.d550, k=log(nrow(euro.dane)))
colnames(out.BIC)<-c("df","BIC") # nadawanie nazw wierszom
out.BIC
out.all<-cbind(out.AIC,out.BIC,out.lambda)
out.all

# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
euro.knn.1<- knearneigh(coords, k=1) # powstaje obiekt klasy knn
euro.k1.nb<- knn2nb(euro.knn.1) # macierz klasy nb dla k najbliższych sąsiadów
euro.k1.sym.nb<-make.sym.nb(euro.k1.nb) # usymetrycznianie macierzy
euro.k1.sym.listw<-nb2listw(euro.k1.sym.nb) # tworzenie obiektu klasy listw

euro.knn.5<- knearneigh(coords, k=5) # powstaje obiekt klasy knn
euro.k5.nb<- knn2nb(euro.knn.5) # macierz klasy nb dla k najbliższych sąsiadów
euro.k5.sym.nb<-make.sym.nb(euro.k5.nb) # usymetrycznianie macierzy
euro.k5.sym.listw<-nb2listw(euro.k5.sym.nb) # tworzenie obiektu klasy listw

euro.knn.10<- knearneigh(coords, k=10) # powstaje obiekt klasy knn
euro.k10.nb<- knn2nb(euro.knn.10) # macierz klasy nb dla k najbliższych sąsiadów
euro.k10.sym.nb<-make.sym.nb(euro.k10.nb) # usymetrycznianie macierzy
euro.k10.sym.listw<-nb2listw(euro.k10.sym.nb) # tworzenie obiektu klasy listw

euro.knn.15<- knearneigh(coords, k=15) # powstaje obiekt klasy knn
euro.k15.nb<- knn2nb(euro.knn.15) # macierz klasy nb dla k najbliższych sąsiadów
euro.k15.sym.nb<-make.sym.nb(euro.k15.nb) # usymetrycznianie macierzy
euro.k15.sym.listw<-nb2listw(euro.k15.sym.nb) # tworzenie obiektu klasy listw

euro.knn.20<- knearneigh(coords, k=20) # powstaje obiekt klasy knn
euro.k20.nb<- knn2nb(euro.knn.20) # macierz klasy nb dla k najbliższych sąsiadów
euro.k20.sym.nb<-make.sym.nb(euro.k20.nb) # usymetrycznianie macierzy
euro.k20.sym.listw<-nb2listw(euro.k20.sym.nb) # tworzenie obiektu klasy listw

euro.knn.25<- knearneigh(coords, k=25) # powstaje obiekt klasy knn
euro.k25.nb<- knn2nb(euro.knn.25) # macierz klasy nb dla k najbliższych sąsiadów
euro.k25.sym.nb<-make.sym.nb(euro.k25.nb) # usymetrycznianie macierzy
euro.k25.sym.listw<-nb2listw(euro.k25.sym.nb) # tworzenie obiektu klasy listw

euro.knn.30<- knearneigh(coords, k=30) # powstaje obiekt klasy knn
euro.k30.nb<- knn2nb(euro.knn.30) # macierz klasy nb dla k najbliższych sąsiadów
euro.k30.sym.nb<-make.sym.nb(euro.k30.nb) # usymetrycznianie macierzy
euro.k30.sym.listw<-nb2listw(euro.k30.sym.nb) # tworzenie obiektu klasy listw


# estysmacja modeli błędu przestrzennego
meuro.k1<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.k1.sym.listw, method="eigen")
meuro.k5<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.k5.sym.listw, method="eigen")
meuro.k10<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.k10.sym.listw, method="eigen")
meuro.k15<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.k15.sym.listw, method="eigen")
meuro.k20<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.k20.sym.listw, method="eigen")
meuro.k25<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.k25.sym.listw, method="eigen")
meuro.k30<-errorsarlm(PZ1~PZ4+PZ5+PZ6, data=euro.dane, euro.k30.sym.listw, method="eigen")

# podsumowanie lambda i kryteriów informacyjnych
out.lambda1<-rbind(meuro.k1$lambda, meuro.k5$lambda, meuro.k10$lambda, meuro.k15$lambda, meuro.k20$lambda, meuro.k25$lambda, meuro.k30$lambda)
out.lambda1
out.AIC1<-anova.sarlm(meuro.k1, meuro.k5, meuro.k10, meuro.k15, meuro.k20, meuro.k25, meuro.k30)
out.AIC1
out.BIC1<-AIC(meuro.k1, meuro.k5, meuro.k10, meuro.k15, meuro.k20, meuro.k25, meuro.k30, k=log(nrow(euro.dane)))
colnames(out.BIC1)<-c("df","BIC") # nadawanie nazw wierszom
out.BIC1
out.all1<-cbind(out.AIC1,out.BIC1,out.lambda1)
out.all1

######################

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# 6.8

######################
###### Przykład 18 ######
######################

attach(euro.dane)
# obiekt danych do oszacowania modelu
euro1995<-cbind(ID, REG, KOD, Z1_1995, Z4_1995, Z5_1995, Z6_1995)
is.data.frame(euro1995)
euro1995.df<-as.data.frame(euro1995)
is.data.frame(euro1995.df)
names(euro1995.df)<-c("ID", "REG", "KOD", "PKB_PPP", "VA_AG", "VA_IND", "VA_SER")
print(euro1995.df)

# obiekt danych do prognozy
euro2000<-cbind(ID, REG, KOD, Z1_2000, Z4_2000, Z5_2000, Z6_2000)
euro2000.df<-as.data.frame(euro2000)
names(euro2000.df)<-c("ID", "REG", "KOD", "PKB_PPP", "VA_AG", "VA_IND", "VA_SER")
detach(euro.dane)

# estymacja modelu opóźnienia przestrzennego
model1.lag<-lagsarlm(PKB_PPP~VA_AG+VA_IND+VA_SER, data=euro1995.df, euro.listw, method="eigen", tol.solve=1.0e-20)

# prognoza modelu na nowych danych
pred1.lag<-predict.sarlm(model1.lag, newdata=euro2000.df, listw=euro.listw)
print(pred1.lag)

# prognoza modelu na starych danych
pred2.lag<- predict.sarlm(model1.lag)

#RMSE dla modelu na nowych danych
sqrt(sum((euro2000.df$PKB_PPP - as.vector(pred1.lag))^2)/length(euro.nb))

# RMSE dla modelu na "starych" danych
sqrt(sum((euro1995.df$PKB_PPP - as.vector(pred2.lag))^2)/length(euro.nb))

>>>>>>>>>>>>>>>>>>

# estymacja modelu opóźnienia przestrzennego
model1.lm<-lm(PKB_PPP~VA_AG+VA_IND+VA_SER, data=euro1995.df)

# prognoza modelu na nowych danych
pred1.lm<-predict.lm(model1.lm, newdata=euro2000.df)
print(pred1.lm)

# prognoza modelu na starych danych
pred2.lm<- predict.lm(model1.lm)

#RMSE dla modelu na nowych danych
sqrt(sum((euro2000.df$PKB_PPP - as.vector(pred1.lm))^2)/length(euro.nb))

# RMSE dla modelu na "starych" danych
sqrt(sum((euro1995.df$PKB_PPP - as.vector(pred2.lm))^2)/length(euro.nb))

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# stworzenie obiektu klasy data.frame z wyniku prognozy
fore1<-as.vector(pred1.lag)
fore.df<-as.data.frame(fore1)
is.data.frame(fore.df)

# połączenie w jeden obiekt y z yhat
theo<-cbind(fore.df, euro2000.df$PKB_PPP)
names(theo)<-c("y_theo", "y_empir")
theo

# procentowe różnice pomiędzy prognozą a danymi rzeczywistymi
attach(theo)
error1<-((y_empir-y_theo)/y_empir)
detach(theo)
error1

summary(error1) # podsumowanie statystyczne różnic