#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>>>>>>>>>>>>>>>ROZDZIAŁ 5>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>>>>>>>WARIOGRAM I KORELOGRAM>>>>>>>>>>>>>>>>>>>>>>>

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#5.2
#Przygotowanie bibliotek

# wczytywanie bibliotek
library(shapefiles)
library(maptools)
library(spdep)
library(geoR)

# liczenie koordynat
euro.shape1<-read.shape("c:/euro/euro.shp", dbf.data=TRUE, verbose=TRUE)
centr<-get.Pcent(euro.shape1)
centr1<-as.data.frame(centr)
names(centr1)<-c(�xcoords","ycoords")
print(centr1)
names(euro.dane) # zostały dodane dwie ostatnie zmienne

#tworzenie obiektu klasy geodata
euro.geoR<-as.geodata(euro.dane, coords.col = 7:8, data.col = 9)
is.geodata(euro.geoR) # sprawdzenie, czy obiekt jest klasy geodata

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#5.3
# liczenie i rysowanie wariogramu
vario1<-variog(euro.geoR)
plot(vario1, type="l", main="Wariogram wielokierunkowy dla zmiennej Z1_1995
PKB wg PPP w 1995r.")
text(2.3e+05, 5.2e+08, "nugget variance", adj=c(1,1))
text(1.085e+06,1.8e+09,"sill")
text(1.085e+06,0,"range")

# Podsumowanie wariogramu
print(vario1)

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#5.4
#wariogram skalowany

# porównanie wariogramów bez skalowania i ze skalowaniem.
euro.geoR.Z1_1995<-as.geodata(euro.dane, coords.col = 7:8, data.col = 9)
euro.geoR.PZ1<-as.geodata(euro.dane, coords.col = 7:8, data.col = 27)
par(mfrow=c(1,2))
vario1<-variog(euro.geoR.Z1_1995)
vario2<-variog(euro.geoR.PZ1)
plot(vario1, type="l", main="nieskalowany dla Z1_1995")
plot(vario2, type="l", main="nieskalowany dla PZ1")
plot(vario1, type="l", scaled=TRUE, main="skalowany dla Z1_1995", var.lines=TRUE)
plot(vario2, type="l", scaled=TRUE, main="skalowany dla PZ1", var.lines=TRUE)

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# 5.5
# przetwarzanie wyników wariogramu
# modyfikacja wydruku
vario.val1<-vario1$v # tworzenie obiektu z wartości wariogramu
vario.val<-as.data.frame(vario.val1) # zmiana formatu na liczbowy obiektu wartości wariogramu
names(vario.val)<-c(�variog.values") # zmiana nazwy zmiennej
is.data.frame(vario.val) # sprawdzenie, czy obiekt jest klasy data frame
vario.val # wyświetlenie zmiennej wartości wariogramu jako danych liczbowych
var(euro.dane$Z1_1995) # policzenie wariancji próbkowej
variog.scaled<-vario.val$variog.values/var(euro.dane$Z1_1995) # wartości wariogramu skalowanego
vario.dist1<-vario1$u # obiekt wartości odległości dla jakich liczony jest wariogram
vario.dist<-as.data.frame(vario.dist1)
names(vario.dist)<-c(�variog.distance") # zmiana nazwy zmiennej
vario.dist.km<-vario.dist$variog.distance/900.90 # przeliczenie punktów na mapie na km
vario.all<-cbind(vario.dist, vario.dist.km, vario.val, variog.scaled) # połączenie wszystkich zmiennych w jeden wydruk
vario.all


#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# 5.6
# porównania graficzne wariogramów

# dwa wariogramy na jednym wykresie
plot(vario1, type="l", scaled=TRUE, var.lines=TRUE)
lines(vario2, type="l", scaled=TRUE)


#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# 5.7
# obwiednie wariogramu

par(mfrow=c(1,2))
# rysowanie wariogramu i testowanie autokorelacji przestrzennej
euro.geoR<-as.geodata(euro.dane, coords.col = 7:8, data.col = 10)
vario1<-variog(euro.geoR)
vario1.env<-variog.mc.env(euro.geoR, obj.var=vario1)
plot(vario1, envelope=vario1.env, type="l")
title(main="wariogram dla zmiennej Z2_1995")

euro.geoR<-as.geodata(euro.dane, coords.col = 7:8, data.col = 28)
vario1<-variog(euro.geoR)
vario1.env<-variog.mc.env(euro.geoR, obj.var=vario1)
plot(vario1, envelope=vario1.env, type="l")
title(main="wariogram dla zmiennej PZ2")

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# 5.8
# ręczne dopasowanie wariogramu
# dopasowanie wzrokowe funkcji rozkładu do empirycznego wariogramu
eyefit(vario1, silent = FALSE)

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# 5.9
# estymacja wariogramu

euro.geoR.Z1_1995<-as.geodata(euro.dane, coords.col = 7:8, data.col = 9)

# dopasowanie przez MNW do danych empirycznych funkcja
likfit()
vario.likfit<-likfit(euro.geoR.Z1_1995, ini=c(0.5, 0.5), fix.nugget=TRUE)
print(vario.likfit)
summary(vario.likfit)


######################
###### Przykład 12 ######
######################
# dopasowanie przez MNW do danych empirycznych - funkcja likfit()
vario.likfit<-likfit(euro.geoR.Z1_1995, ini=c(0.5, 0.5), fix.nugget=TRUE)
print(vario.likfit)
summary(vario.likfit)

#dopasowanie przez MNK do danych empirycznych - funkcja variofit()
vario.variofit<-variofit(vario1, ini.cov.pars=c(100000000,1000000))
print(vario.variofit)
summary(vario.variofit)


# rysunek wariogramów
plot(vario1)
lines(vario1.fit)
lines(vario1.fitvario)

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

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
# 5.11
# wyliczenie i rysowanie korelogramu

# korelogram oparty na klasycznym współczynniku korelacji
corr.classic<- sp.correlogram(euro.nb, euro.dane$Z1_1995, order=6, method="corr")

# korelogram oparty na współczynniku autokorelacji I Morana
corr.moran<-sp.correlogram(euro.nb, euro.dane$Z1_1995, order=6, method="I")

print(corr.classic)
print(corr.moran)

par(mfrow=c(1,2))
plot(corr.moran, main="Rys.a")
plot(corr.classic, main="Rys.b")

# testy istotności dla I Morana
corr.moran<-sp.correlogram(euro.nb, euro.dane$Z1_1995, order=6, method="I")
# tworzenie statystyki testowej
stat <- apply(corr.moran$res, 1, function(x) (x[1]-x[2])/sqrt(x[3]))
# tworzenie tabeli z p-value
corr.moran<-cbind(corr.moran$res, stat, p.adjust(2*pnorm(abs(stat), lower.tail=FALSE)))
colnames(corr.moran)<-c("I Morana", "E(I)", "Var(I)", "ZI", "p-value")
corr.moran

# testy istotności dla współczynnika korelacji liniowej
corr.classic<- sp.correlogram(euro.nb, euro.dane$Z1_1995, order=6, method="corr")
# pętla licząca liczbę obserwacji n w każdym opóźnieniu
for(i in 1:6)
{obs[i]<-sum(unlist(corr.classic$cardnos[i])*as.numeric(names(unlist(corr.classic$cardnos[i]))))}
corr<-cbind(corr.classic$res, obs)
corr<-as.data.frame(corr)
colnames(corr)<-c("r", "n")
# statystyka testująca i skorygowane p-value
stat<-(corr$r/sqrt(1-(corr$r)^2))*sqrt(corr$n-2)
pval<-p.adjust(1-pnorm(abs(stat)), method="bonferroni", n=length(stat))
corr<-cbind(corr, stat, pval)
corr