#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>>>>>>>>>>>>>>>>>>>>>>>>ROZDZIAŁ 2>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>>>>>>>>>PODSTAWOWE OPERACJE W PROGRAMIE R>>>>>>>>>>>>>>>>>>>>
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#2.2
# wyszukiwanie
?moran.test
apropos("mean")
library(help="spdep")
example("mean")
moran.test

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#2.3
# sprawdzanie pamięci
library(spdep)
search()

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#2.4
# Język programu R
X <- Y+Z

search()

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#2.5
# Wprowadzanie danych

PKB.dane<-read.table("c:/euro/PKB.txt", header=TRUE)
euro.dane<-read.csv("c:/euro/euro.csv",header=T, sep=";",dec="." )

names(euro.dane) #pokazuje nagłówki

names(euro.dane)
is.data.frame(euro.dane)
euro.dane.frame<-as.data.frame(euro.dane)
is.data.frame(euro.dane.frame)
names(euro.dane.frame)

######################
###### Przykład 1 #######
######################
names(euro.dane)

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

# eksport obiektu euro.dane na dysk c
utworzenie pliku nowedane.txt rozdzielanego przecinkami.
write.table(euro.dane, file="c:/nowedane.txt", sep=",")

is.na(euro.dane)<-euro.dane="."

wektor1<-c(1,2,3,4,5) # wektor liczb 1,2,3,4,5
1:20 # liczby od 1 do 20, co 1
wektor2<-seq(1,20, 0.1) # wektor liczb od 1 do 20, z krokiem 0.1
wektor3<-seq(1,20, length=100) # wektor liczb od 1 do 20 zawierający w sumie 100 liczb (razem z pierwszą i ostatnią).

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#2.6
#Operacje na obiektach

print(euro.dane) #wyrzuca wszystkie dane z obiektu na ekran
print(euro.dane$Z1_1995) # wyświetlenie jednej zmiennej
euro.dane # wyświetlenie całego zbioru, można pominąć komendę print()
euro.dane$Z1_1995 # tu też można pominąć komendę print()

ls()
class(euro.dane)
class(euro.nb)

euro.dane$Z1_1995[euro.dane$Z1_2000 & euro.dane$ID>15] # warunek koniunkcji (oraz)
euro.dane$Z1_1995[euro.dane$ID<10 | euro.dane$ID>30] #warunek alternatywy (lub)

Z1_1995<-euro.dane$Z1_1995 # wyodrębnienie zmiennej z macierzy
print(Z1_1995) # wyświetlenie nowego obiektu

######################
###### Przykład 2 #######
######################
euro.dane<-read.table("c:/euro/euro.csv",header=T, sep=";", dec=".") #wczytujemy nowe dane
names(euro.dane)
Z1<-euro.dane$Z1_1995
Z2<-log(Z1)
Z3<-cbind(Z1, Z2)

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

subset(euro.dane, KRAJ==2, select=PZ1:PZ3)

z1_1995<-euro.dane$Z1_1995 # wyodrębnianie zmiennej z obiektu bazy danych
z1_1995.cut<-cut(z1_1995,10, include.lowest=TRUE) # dzielenie zmiennej na 10 przedziałów
levels(z1_1995.cut) # wyświetlanie samych przedziałów
tz1<-table(z1_1995.cut)
is.matrix(tz1)
tz1m<-as.matrix(tz1) # formatowanie danych
tz1m # wyświetlanie szeregu rozdzielczego przedziałowego dla danej zmiennej

# przedziały długości 10000 od 0 do 120000
table(cut(euro.dane$Z1_1995, br=10000*(0:12)))

x<-5 #tworzymy obiekt x o wartości 5
x #wyświetlenie zawartości obiektu x
nowy_x<-x #zmiana nazwy obiektu z x na nowy_x
nowy_x #wyświetlenie zawartości obiektu nowy_x

x<-seq(1, 20, 0.4) # generowanie zmiennej x
y<-seq(2,21,0.4) # generowanie zmiennej y
z<-seq(3,22,0.4) # generowanie zmiennej z
xyz<-cbind(x,y,z) # łączenie zmiennych w jeden obiekt xyz
xyz # wyświetlenie obiektu xyz
is.data.frame(xyz) # sprawdzenie czy xyz jest obiektem klasy data frame
is.matrix(xyz) # sprawdzenie czy xyz jest macierzą
names(xyz) # żądanie wyświetlenia nagłówków
colnames(xyz)<-c("zmienna x", "zmienna y", "zmienna z") # nadanie nowych nazw kolumnom
rownames(xyz)<-rownames(xyz, do.NULL = FALSE, prefix = "Obs.") # nadanie nazw wierszom
xyz # wyświetlenie obiektu xyz po zmianach
xyz.df<-as.data.frame(xyz) # konwersja xyz na data frame
names(xyz.df) # żądanie wyświetlenia nagłówków kolumn
colnames(xyz.df) # żądanie wyświetlenia nagłówków kolumn
rownames(xyz.df) # żądanie wyświetlenia nagłówków wierszy

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#2.7
# Podstawowe statystyki

mean(euro.dane$Z1_1995)

summary(euro.dane) # pokazuje podstawowe statystyki dla wszystkich zmiennych

# podsumowanie zmiennej Z1_1995 w grupach krajów - zmienna KRAJ
by(euro.dane$Z1_1995, euro.dane$KRAJ, summary)

var1<-rnorm(1000)
quantile(var1)

# własne określenie progów prawdopodobieństwa
quantile(var1, probs=c(1,2,5,95,99)/100)

nlevels(factor(euro.dane$ID))

levels(euro.dane$KOD)

summary(factor(euro.dane$KOD))

mean(euro.dane$Z1_1995[euro.dane$KRAJ==1]) #średnia z Z1_1995 gdy zmienna KRAJ wynosi 1

fivenum(euro.dane$Z1_1995) # liczy statystyki pozycyjne dla zmiennej Z1_1995

x<- euro.dane$Z1_1995
zx<-(x-mean(x))/sd(x) # standaryzacja
mean(zx) # powinna wynosić 0
var(zx) # powinna wynosić 1

x1<- euro.dane$Z1_1995
zx1<-scale(x1, mean(x1), sd(x1))
mean(zx1)
var(zx1)

Z1_1995<-euro.dane$Z1_1995 # tworzenie obiektu ze zmiennej Z1_1995
par(mfrow=c(1,2)) # podział okna na dwie kolumny i jeden wiersz
boxplot(Z1_1995, col="bisque") # wykres pudełkowy ze zmiennej ogółem
boxplot(Z1_1995~euro.dane$KRAJ, col="blue") # wykresy pudełkowy względem zmiennej KRAJ

par(mfrow=c(1,1))
qqnorm(log(euro.dane$Z1_1995), main="wykres kwantylowy dla Z1_1995")
qqline(log(euro.dane$Z1_1995))

par(mfrow=c(1,2))
hist(euro.dane$Z1_1995, col="salmon1", labels=TRUE, br=10, border="brown", xlab="PKB wg PPP w mln euro w 1995") #histogram dla zmiennej Z1_1995
hist(euro.dane[,5], col="salmon2", labels=TRUE, br=10, border="brown", xlab="PKB wg PPP w mln euro w 2000") #histogram dla zmiennej Z1_2000, która jest w piątej kolumnie

######################
###### Przykład 3 #######
######################

# Narysuj histogramy logarytmów zmiennych Z1_1995 i Z1_2000.
par(mfrow=c(1,2))
hist(log(euro.dane$Z1_1995), xlab="log PKB wg PPP w mln euro w 1995")
hist(log(euro.dane[,5]), xlab="log PKB wg PPP w mln euro w 2000")

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

# R jako kalkulator
3+5

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#2.8
# Grafika
par(mfrow=c(1,2))
split.screen(c(1,2))

cols<-grey(4:1/5) # wektor czterech szarych kolorów
od 1 do 5
cols<-rev(topo.colors(4)) # wektor czterech kolorów z palety topo

# wyświetlanie nazw wbudowanych kolorów
cl<-colors()
length(cl);cl[1:657]

######################
###### Przykład 4 #######
######################

##------ Some palettes ------------
demo.pal <-
function(n, border = if (n<32) "light gray" else NA,
main = paste("color palettes; n=",n),
ch.col = c("rainbow(n, start=.7, end=.1)", "heat.colors(n)",
"terrain.colors(n)", "topo.colors(n)", "cm.colors(n)"))
{
nt <- length(ch.col)
i <- 1:n; j <- n / nt; d <- j/6; dy <- 2*d
plot(i,i+d, type="n", yaxt="n", ylab="", main=main)
for (k in 1:nt) {
rect(i-.5, (k-1)*j+ dy, i+.4, k*j,
col = eval(parse(text=ch.col[k])), border = border)
text(2*j, k * j +dy/4, ch.col[k])
}
}
n <- if(.Device == "postscript") 64 else 16
# Since for screen, larger n may give color allocation problem
demo.pal(n)

####################
par(mfrow=c(1,2))
# generowanie zmiennych
x<-seq(1,25,1)
y1<-0.5*x^2
# wykresy ramek
plot(y1, type="l", bty="7", lty=2, main="Rys.a")
text(10,200, "bty 7")
text(11,250, "ramka górna i prawa")
plot(y1, type="l", bty="l", lty=3, main="Rys.b")
text(10,200, "bty l")
text(11, 250, "ramka dolna i lewa")

####### kody rysunku 2.6
par(mfrow=c(3,2))
# generowanie zmiennych
x<-seq(1,25,1)
y1<-0.5*x^2
y2<-0.5*x^2+50
y3<-0.5*x^2+100
y4<-0.5*x^2+150
y5<-0.5*x^2+200

# rysunek na grubość linii
plot(y1, type="l", main="Rys.a")
lines(y2, lwd=2)
lines(y3, lwd=3)
lines(y4, lwd=4)
lines(y5, lwd=5)
text(4,-0.5,"lwd=1", cex=0.8)
text(3,37,"lwd=2", cex=0.8)
text(3,92,"lwd=3", cex=0.8)
text(3,142,"lwd=4", cex=0.8)
text(3,192,"lwd=5", cex=0.8)
text(15, 35, "grubość linii")


# rysunek na różne znaczniki
plot(y1, , main="Rys.b")
points(y2, pch="*")
points(y3, pch="&")
text(18,35, "wykresy punktowe")

# rysunek różnych linii
plot(y1, type="l", main="Rys.c")
lines(y2, type="o")
lines(y3, type="b")
lines(y4, type="s")
lines(y5, type="S")
text(4,0, "typ l", cex=0.8)
text(4,45, "typ o", cex=0.8)
text(4,95, "typ b", cex=0.8)
text(4,145, "typ s", cex=0.8)
text(4,195, "typ S", cex=0.8)
text(17,35, "wykresy liniowe")

# wykres linii przerywanych lty
plot(y1, type="l", lty=1, main="Rys.d")
lines(y2, lty=2)
lines(y3, lty=3)
lines(y4, lty=4)
text(5,0, "lty 1", cex=0.8)
text(4,43, "lty 2", cex=0.8)
text(4,93, "lty 3", cex=0.8)
text(4,143, "lty 4", cex=0.8)
text(18,35, "wykresy liniowe przerywane")

plot(y1, type="h", main="Rys.e")
text(5, 50, "typ h")
text(6,250, "wykres słupkowy")

# wykres czcionki
plot(0, main="Rys.f")
text(1,0.7, "czcionka cex=1", cex=1)
text(1, 0.0, "czcionka cex=2", cex=2)
text(1,-0.5, "czcionka cex=0.5", cex=0.5)


#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#2.9
# Regresja

# rysunek dwóch zmiennych względem siebie
par(mfrow=c(1,1))
plot(log(euro.dane$Z1_1995),log(euro.dane$Z1_2000), xlab=" PKB wg PPP 1995",ylab="PKB wg PPP 2000") # rysunek logarytmów zmiennych x i y.

# Komenda regresji
mod1<-lm(log(Z1_2000)~log(Z1_1995), data=euro.dane)
mod1 # wynik regresji
wersja uboższa
summary(mod1) # wynik regresji
wersja bogatsza

# wykres dopasowania modelu do danych empirycznych
plot(euro.dane$Z1_1995,euro.dane$Z1_2000)
abline(mod1) # mod1 jako obiekt wyniku regresji

# wyodrębnienie wyników regresji
plot(euro.dane$Z1_1995, euro.dane$Z1_2000)
intercept<-mod1$coefficients[1] # wyodrębnienie wartości wyrazu wolnego
slope<-mod1$coefficients[2] # wyodrębnienie wartości współczynnika kierunkowego zmiennej Z1_1995
intercept
slope
abline(intercept, slope) # dorysowanie linii regresji do wykresu

# analiza reszt
par(mfrow=c(1,2))
res<-residuals(mod1)
plot(res)
abline(h=0, lty=2)
title(main="wykres punktowy reszt")
hist(res, breaks=30) # rysowanie histogramu
# tworzenie funkcji gęstości rozkładu normalnego
mean.res<-mean(res)
sd.res<-sd(res)
x<-seq(-1, 1,length=100)
y<-dnorm(x,mean.res, sd.res)
lines(x,y)

# rysowanie reszt regresji
mod1.fit<-predict(mod1) # wartości teoretyczne - dopasowanie modelu do danych
blad<-log(euro.dane$Z1_2000)-mod1.fit # błąd jako różnica wartości empirycznych i teoretycznych
par(mfrow=c(1,2))
plot(blad)
abline(h=0, lty=2)
hist(blad, breaks=30)

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#2.10
# Przetwarzanie wyników
names(mod1)

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#2.11
# Obserwacje nietypowe (outliersy)

# posiadamy obiekt danych euro.dane
# wyodrębnienie danych dotyczących PKB, obiektu PKB.dane
PKB.dane<-cbind(euro.dane$ID, euro.dane$KOD, euro.dane$KRAJ, euro.dane$Z1_1995, euro.dane$Z1_2000)
PKB.dane<-as.data.frame(PKB.dane)
colnames(PKB.dane)<-c("ID", "KOD", "KRAJ", "Z1_1995", "Z1_2000")
names(PKB.dane)

plot(PKB.dane[,4:5])

par(mfrow=c(1,2))
hist(PKB.dane$Z1_1995)
hist(PKB.dane$Z1_2000)
par(mfrow=c(1,1))

PKB.dane[PKB.dane$Z1_1995>80000, ] # wyświetlenie obserwacji, których wartości sa wyższe niż 8 000

library(scatterplot3d)
scatterplot3d(PKB.dane[,4:5]) # wykres trójwymiarowy wszystkich danych
outy<-c(-3,-7,-12,-64,-67) # zdefiniowanie outliersów
scatterplot3d(PKB.dane[outy,4:5]) # wykres trójwymiarowy bez outliersów

# model na danych kompletnych
regresja i jej podsumowanie
mod1<-lm(Z1_2000~Z1_1995, data=euro.dane)
summary(mod1)

# regresja bez outliersów
outy<-c(-3,-7,-12,-64,-67) # wektor obserwacji nietypowych
mod2<-lm(euro.dane[outy,10]~euro.dane[outy, 9]) # regresja bez outliersów
summary(mod2)

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#2.12
# Wprowadzanie danych przestrzennych

# ładowanie pakietów do pamięci
library(shapefiles)
library(maptools)
library(spdep)

# wczytywanie danych - metoda 1
euro.shape<-read.shapefile("c:/euro/euro")
euro.poly<-shape2poly(euro.shape)
euro.nb<-poly2nb(euro.poly)
euro.listw<- nb2listw(euro.nb)
# wczytywanie danych liczbowych z pliku *.dbf
euro.dane<-euro.shape$dbf # powstaje obiekt danych klasy data.frame

# wczytywanie danych - metoda 2
euro.shape<-read.shape("c:/euro/euro.shp") # z pakietu maptools
euro.dane.<-euro.shape$att.data
euro.poly<-Map2poly(euro.shape) # z pakietu maptools
euro.nb<-poly2nb(euro.poly)
euro.listw<- nb2listw(euro.nb)

# wczytywanie danych liczbowych z pliku *.csv
euro.dane<-read.csv("c:/euro/euro.csv", header=T, sep=";", dec="." )

names(euro.dane)

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#2.13
# rysunek mapy konturowej z wykorzystaniem komendy plot.Map()
euro.shape1<-read.shape("c:/euro/euro.shp", dbf.data=TRUE, verbose=TRUE)
plot.Map(euro.shape1)

# rysunek mapy konturowej z wykorzystaniem komendy plot()
plot(euro.poly, forcefill=FALSE) #rysunek

brks<- c(1,2,3,4,5,6) # umowna nazwa obiektu brks, liczby w nawiasie są górnymi granicami przedziałów
# cols<- rev(cm.colors(6)) # umowna nazwa obiektu cols
cols<- c("grey25", "grey35", "grey85", "grey50", "grey60", "grey70")
plot(euro.poly, col=cols[findInterval(euro.dane$KRAJ, brks)],forcefill=FALSE)
legend(x=c(-700000,-700000), y=c(-310000,-310000), fill=cols,legend=c("Austria","Czechy","Niemcy","Węgry","Polska","Słowacja"),bty="n", cex=0.9)
title("Regiony wg przynależnosci do krajów")


brks1<-round(fivenum(euro.dane$Z1_1995),digits=1)
cols<-c("grey80", "grey60", "grey40", "grey20")
# cols<-rev(heat.colors(4)) # inna możliwość zdefiniowania kolorów
plot(euro.poly, col=cols[findInterval(euro.dane$Z1_1995,brks1)],forcefill=FALSE)
title(main="odsetek PKB wg PPP w mln euro w 1995 w ukladzie regionalnym")
legend(x=c(-700000,-700000), y=c(-400000,-400000), legend=c("mniej niz Mediana-Q","od Mediana-Q do Mediana","od Mediana do Mediana+Q","wiecej niz Mediana+Q"), leglabs(brks1), fill=cols, bty="n", cex=0.9)


# przedziały punktowe dla liczb całkowitych 1,2,3,4,5,6
brks1<- c(min(euro.dane$KRAJ),2,3,4,5,max(euro.dane$KRAJ))

# przdziały wg miar klasycznych: średniej i odchylenia standardowego
brks2<-round(c(min(euro.dane$Z1_1995), mean(euro.dane$Z1_1995)-sd(euro.dane$Z1_1995), mean(euro.dane$Z1_1995), mean(euro.dane$Z1_1995)+sd(euro.dane$Z1_1995),digits=1)

# przedziały wg miar pozycyjnych: mediany i odchylenia ćwiartkowego
brks3<-round(fivenum(euro.dane$Z1_1995),digits=1)

# przedziały punktowe dla liczb całkowitych 1,2,3,4,5,6 ( to samo co brks1)
brks4<- c(1,2,3,4,5,6)

# pięć przedziałów o równej długości
brks5<-cut(euro.dane$Z1_1995, 5)