#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>>>>>>>>>>>>>>>ROZDZIAŁ 4>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>>>>>>>>>>>>>>>STATYSTYKI LOKALNE i GLOBALNE>>>>>>>>>>
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

$ 4.1
# globalna statystyka Morana
moran(euro.dane$Z1_1995, euro.listw, length(euro.nb), Szero(euro.listw))
moran.test(euro.dane$Z1_1995, euro.listw)
moran.test(euro.dane$Z1_1995, euro.listw, randomisation=FALSE)

# wartość p-value statystyki Morana
pval.norm<-1-pnorm(2.5895, mean=0, sd=1) # p-value dla I Morana (rozkład normalny)
pval.norm

pval.rand<-1-pnorm(2.6534, mean=0, sd=1) # p-value dla I Morana (podejście randomizacji)
pval.rand

# testowanie statystyki Morana metodą Monte Carlo
moran.mc(euro.dane$Z1_1995, euro.listw, 99)
moran.mc(euro.dane$Z1_1995, euro.listw, 999)

# wykres Morana dla zmiennej standaryzowanej
Z1_1995.std<-((euro.dane$Z1_1995-mean(euro.dane$Z1_1995))/(sd(euro.dane$Z1_1995)))
moran.plot(Z1_1995.std, euro.listw, labels=as.character(euro.dane$REG), pch=19)

######################
###### Przykład 8 #######
######################

x<-euro.dane$Z1_1995 # wyodrębnienie zmiennej ze zbioru
zx<-(x-mean(x))/sd(x) #standaryzacja zmiennej
mean(zx)
sd(zx)
wzx<-lag.listw(euro.listw, zx) # opóźnienie przestrzenne zmiennej x
morlm<-lm(wzx~zx) # regresja
slope<-morlm$coefficients[2]
intercept<-morlm$coefficients[1]
par(pty="s") # kwadratowe okno wykresu
plot(zx,wzx,xlab="zx",ylab="opóźnienie przestrzenne zx")
abline(intercept, slope) # linia regresji
abline(h=0, lty=2) # linia pozioma w y=0
abline(v=0, lty=2) #linia pionowa w x=0

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

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# 4.2
# statystykia globalna Geary'ego

geary(euro.dane$Z1_1995, euro.listw, length(euro.nb), length(euro.nb)-1, Szero(euro.listw))
geary.test(spNamedVec("Z1_1995",euro.dane),euro.listw)
geary.mc(euro.dane$Z1_1995, euro.listw,99)

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# 4.3
#statystyki i testy join-count

# podział graficzny obserwacji badanej zmiennej na trzy grupy wartości
plot(euro.dane$Z2_1995)
abline(h=13000, lty=2)
abline(h=25000, lty=2)

#przygotowanie danych - podział danych na trzy grupy
PKB1995<-cut(euro.dane$Z2_1995, breaks=c(0,14000,25000,39000), labels=c("malo", "srednio","duzo"))
PKB1995.f<-factor(PKB1995)

joincount.test(PKB1995.f, euro.listw)
joincount.mc(PKB1995.f, euro.listw,99)
joincount.multi(PKB1995.f, euro.listw)

######################
###### Przykład 9 #######
######################

summary(euro.dane$Z4_1995)
sd(euro.dane$Z4_1995) # odchylenie standardowe

plot(euro.dane$Z4_1995)

PKB1995<-cut(euro.dane$Z4_1995, breaks=c(0,500,1000,1700), labels=c("malo", "duzo", "outliersy"))
names(PKB1995)<-rownames(euro.dane)
PKB1995.f<-factor(PKB1995)

joincount.test(PKB1995.f, euro.listw)

joincount.mc(PKB1995.f, euro.listw, 99)

# przestrzenny rozkład wartości na trzy wyróżnione grupy.
brks1<- c(0,500,1000,1700)
cols<-c("grey80", "grey60", "grey30")
plot(euro.poly, col=cols[findInterval(euro.dane$Z4_1995,brks1)],forcefill=FALSE)
title(main="PKB per capita w 1995 w grupach")
legend(x=c(-700000,-700000), y=c(-300000,-300000), legend=c("malo", "duzo", "outliersy"), leglabs(brks1), fill=cols, bty="n")

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

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# 4.5
#statystyka lokalna Morana I

# kod generujący listę lokalnych statystyk Morana z przypisaniem nazwy regionu
locM<-localmoran(spNamedVec("Z1_1995", euro.dane),euro.listw)
oid1<-order(euro.dane$ID)
printCoefmat(data.frame(locM[oid1,], row.names=euro.dane$REG[oid1]), check.names=FALSE)

######################
###### Przykład 10 ######
######################

# wyznaczanie statystyk lokalnych Morana
locM<-localmoran(spNamedVec("Z1_1995", euro.dane), euro.listw)
locMdf<-as.data.frame(locM)

# wykres kolorystyczny statystyki lokalnej Morana na mapie
brks<-c(0.000000000000000000001, 0.05, 0.95, 0.99999999999999999999999)
cols<-c("grey30", "grey90", "grey60")

# brks<-c(round(min(locMdf$"Pr(z > 0)",digits=5), 0.05000, 0.95000, round(max(locMdf$"Prob),digits=5))
# cols<-c(�orchid1", �snow", �maroon")

plot(euro.poly, col=cols[findInterval(locMdf$"Pr(z > 0)", brks)],forcefill=FALSE)
legend(x=c(-700000,-700000), y=c(-453000,-453000), legend=c("otoczony podobnymi wartościami, locM>0", "nieistotne", "otoczony odmiennymi wartościami, locM<0"), fill=cols, bty="n", cex=0.9)
title(main="statystyka Local Moran dla zmiennej Z1_1995,
PKB wg PPP w roku 1995", cex=0.7)

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

######################
###### Przykład 11 ######
######################

# mapa przynależności do ćwiartek wykresu punktowego Morana dla PKB wg PPP z 1995r.

# tworzenie zmiennej do analizy
x<-euro.dane$Z1_1995 # wyodrębnienie zmiennej ze zbioru
zx<-(x-mean(x))/sd(x) #standaryzacja zmiennej
mean(zx)
sd(zx)
wzx<-lag.listw(euro.listw, zx) # opóźnienie przestrzenne zmiennej x
cond1<-ifelse(zx>=0 & wzx>=0, 1,0) # I cwiartka
cond2<-ifelse(zx>=0 & wzx<0, 2,0) # II cwiartka
cond3<-ifelse(zx<0 & wzx<0, 3,0) # III cwiartka
cond4<-ifelse(zx<0 & wzx>=0, 4,0) # IV cwiartka
cond.all<-cond1+cond2+cond3+cond4 # zmienna wskazująca przynależność do ćwiartek
cond.all
cond<-as.data.frame(cond.all)
is.data.frame(cond)

# wykres - mapa kolorystyczna
brks<-c(1,2,3,4)
cols<-c("grey25", "grey60", "grey85", "grey45")
plot(euro.poly, col=cols[findInterval(cond$cond.all, brks)],forcefill=FALSE)
legend(x=c(-700000,-700000), y=c(-412000,-412000), legend=c("I ćw - HH - wysokie otoczone wysokimi", "II ćw - LH - niskie otoczone wysokimi", "III ćw - LL - niskie otoczone niskimi", "IV ćw - HL - wysokie otoczone niskimi"), fill=cols, bty="n", cex=0.85)
title(main="Przynalezność regionów do ćwiartek
z wykresu punktowego Morana")

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


#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# 4.6
# statystyka local Gi

Z1locG<-localG(spNamedVec("Z1_1995", euro.dane), nb2listw(euro.nb))
Z1G<-as.vector(Z1locG)
is.vector(Z1G)
summary(Z1G)
Z1Gdf<-as.data.frame(Z1G)
is.data.frame(Z1Gdf)
# test istotności t-studenta dla n=100
sig<-ifelse(Z1Gdf$Z1G<=-3.289 | Z1Gdf$Z1G>=3.289, "*", " ")
sigdf<-as.data.frame(sig)
all<-cbind(Z1Gdf$Z1G, sigdf)
all
# wykres statystyk istotnych
brks<-c(-10,-3.289, 3.289, 10)
cols<-c("grey70", "grey90", "grey55")
plot(euro.poly, col=cols[findInterval(Z1Gdf$Z1G, brks)],forcefill=FALSE)
legend(x=c(-700000,-700000), y=c(-430000,-430000), legend=c("G<0 - klaster niskich wartości", "G=0", "G>0 - klaster wysokich wartości") , fill=cols, bty="n", cex=0.9)
title(main="Istotne statystyki lokalne Getisa Gi")


# statystyka local Gi*
Z1locGstar<-localG(spNamedVec("Z1_1995", euro.dane), nb2listw(include.self(euro.nb)))
Z1G<-as.vector(Z1locGstar)
is.vector(Z1G)
summary(Z1G)
Z1Gdf<-as.data.frame(Z1G)
is.data.frame(Z1Gdf)
# test istotności t-studenta dla n=100
sig<-ifelse(Z1Gdf$Z1G<=-3.289 | Z1Gdf$Z1G>=3.289, "*", " ")
sigdf<-as.data.frame(sig)
all<-cbind(Z1Gdf$Z1G, sigdf)
all
# wykres statystyk istotnych
brks<-c(-10,-3.289, 3.289, 10)
cols<-c("grey40", "grey85", "grey55")
plot(euro.poly, col=cols[findInterval(Z1Gdf$Z1G, brks)],forcefill=FALSE)
legend(x=c(-700000,-700000), y=c(-430000,-430000), legend=c("G*<0 - klaster niskich wartości", "G*=0", "G*>0 - klaster wysokich wartości") , fill=cols, bty="n", cex=0.9)
title(main="Istotne statystyki lokalne Getisa Gi*")