Tag: kod w R

Bilion miliardów a Tufte

Europa gasi pożar olbrzymią ilością wirtualnych pieniędzy (niestety pożar nie jest wirtualny, ale to temat dla innego bloga).
Rzeczpospolita (nomen omen, której infografiki dosyć lubię) artykuł o ratowaniu sytuacji finansowej w Europie (tutaj link do artykułu) okrasiła takim oto wykresem

 

Jaki jest problem z tym wykresem? To, że zgodnie z opisem dane są w miliardach a na rysunku pojawia się słowo bilion, które sugeruje że potencjalne zasoby będą wynosiły miliard bilionów, to już pomijam.

Klasyczny problem z takimi wykresami polega na tym, że nie jest jasne czy podanym liczbom odpowiada wysokość, szerokość czy pole graficzki.

Pomiar średnicy bilionowej monetki daje 190px, pomiar średnicy drugiej co do wielkości monetki daje 125px. A więc to pole odpowiada podanej liczbie.

Szkoda, bo jak się okazuje ludzie bardzo niedokładnie potrafią porównywać pola figur, gorzej jest już tylko z objętościami. Najlepiej porównuje się długości odcinków w poziomie (przynajmniej zgodnie z badaniami Tufty’ego).

Więc zróbmy wykres w R tych samych liczb, ale bez fajerwerków.

I kod R który ten wykres wygenerował

 

1
2
3
4
5
6
7
8
9
10

x <- c(106.4, 376, 440, 1000, 74.1)
par(mar=c(5,20,3,3))
plot(1,type="n", las=1, yaxt="n",ylab="", xlab="nakłady w mld euro",xlim=c(0,1000), ylim=c(0.5,5.5), bty="n")
sapply(1:5, function(i) lines(c(0,x[i]),c(i,i)))
abline(v=seq(100,1000,100),col="grey95")
abline(v=0)
points(x, 1:5, pch=19)
mtext(side=2,line=-1.5,at=1:5, c("potrzeby kapitalowe bankow w Europie", 
     "pozyczki udzielone Portugali Irlandi Grecji", "dotychczasowe zasoby Europejskiego Funduszu Stabilnosci", 
     "potencjalne zasoby Europejskiego Funduszu Stabilnosci", "rezerwy walutowe Polski na koniec X 2011"), las=1, cex=0.9)

Rysujemy rozkład cen krok po kroku, część 4

Czas na ostatnią część wyjaśnień krok po kroku jak konstruowane były wykresy o cenach mieszkań.
Tym razem wykorzystamy wykres pudełkowy pokazany na wpisie tutaj do pokazania rozkładów cen w dzielnicach Warszawy.

Wczytujemy pierwsze 33 linie kodu z poprzedniego wpisu a następnie uruchamiamy linie 142-187. Wyjaśnijmy od razu po co była funkcja nazwyIprocenty(). Otóż w pakiecie lattice dosyć łatwo narysować wykres w podziale na poziomy pewnej zmiennej grupującej. Grupa obserwacji odpowiadająca poszczególnym poziomom rysowana jest na kolejnym panelu. Nazwy poziomów znajdują się w nagłówku panelu. W naszym przykładzie funkcja nazwyIprocenty() zmieniła nazwy wszystkich poziomów w ten sposób, że do nazw dzielnic dodano cztery liczby określające procentową zmianę ceny w określonej dzielnicy (zmianę liczoną na różne sposoby, zobacz komentarze wewnątrz tej funkcji).

Dzięki temu warunkując po zmiennej dzielnica2 powinniśmy uzyskać zbiór wykresów pudełkowych w rozbiciu na dzielnicę.
Poniższy kod od kodu z poprzedniego wpisu różni się praktycznie wyłącznie formułą cenam2~dataF|dzielnica2.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

bwplot(cenam2~dataF|dzielnica2,
         data=mieszkaniaKWW2011Warszawa3,
         scales=list(y=list(log=T, at=at), x=list(rot=90, cex=0.6)),
         ylim=c(5000, 20000),
         panel =  function(...) {
            tmp <- trellis.par.get("plot.symbol")
              tmp$pch=19
              tmp$col="grey20"
              tmp$cex=1/2
            trellis.par.set("plot.symbol",tmp)
            tmp <- trellis.par.get("box.rectangle")
              tmp$col="grey20"
            trellis.par.set("box.rectangle",tmp)
            tmp <- trellis.par.get("box.umbrella")
              tmp$col="grey20"
            trellis.par.set("box.umbrella",tmp)
 
            kolory <- brewer.pal(4, "Set1")
            panel.abline(h=log(at,10), col="grey85")
            panel.abline(v=4.5 + c(0,12,24,36), col="grey85")
 
            panel.bwplot(..., col="grey20")
            args <- list(...)
            mod  <- rlm(args$y~as.numeric(args$x))
            panel.abline(mod, col=kolory[1], lwd=4, alpha=0.85)
 
            mediany <- sapply(miesiace, function(x) median(args$y[args$x == x], na.rm=T))
            mod2    <- rlm(mediany~seq_along(miesiace))
            panel.abline(mod2, col=kolory[2], lwd=4, alpha=0.85)
            indtmp  <- c(1, length(mediany))
            llines(indtmp, mediany[indtmp], col=kolory[4], lwd=4, alpha=0.85)
 
            panel.loess(..., col=kolory[3], lwd=4, alpha=0.85)
          }
)

Wadą tego wykresu są ponownie dzielnice w których mało jest oferowanych mieszkań. Usuńmy dzielnice w których jest mniej niż 1000 mieszkań średniej wielkości oferowanych do sprzedaży w ostatnich 4 latach. Poniżej prezentujemy tylko kod usuwający odpowiednie wiersze, następnie używamy tego samego kodu co powyżej aby wygenerować wykres dla dzielnic, tym razem już tylko 12.

37
38
39
40

# usun male dzielnice
usun <- names(which(table(mieszkaniaKWW2011Warszawa3$dzielnica)<1000))
mieszkaniaKWW2011Warszawa3 <- mieszkaniaKWW2011Warszawa3[!(mieszkaniaKWW2011Warszawa3$dzielnica %in% usun),]
mieszkaniaKWW2011Warszawa3$dzielnica <- factor(mieszkaniaKWW2011Warszawa3$dzielnica)

Rysujemy rozkład cen krok po kroku, część 3

Dzisiaj kontynuujemy rozpisywanie krok po kroku wizualizacji cen mieszkań. Ten wpis poświęcony jest wykresowi pudełkowemu. Kory w programie R użyte poniżej można znaleźć na stronie tutaj.

Podobnie jak poprzedni pierwsza część to przygotowanie danych, druga to ich wizualizacja.

Częśc 1.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

library(lattice)
library(RColorBrewer)
library(MASS)
 
# wczytujemy dane o mieszkaniach
mieszkaniaKWW2011 <-
     read.table("http://tofesi.mimuw.edu.pl/~cogito/smarterpoland/mieszkaniaKWW2011/mieszkaniaKWW2011.csv",
         row.names=NULL, sep=";", header=TRUE,
         colClasses=c("factor", "factor", "numeric", "numeric", "factor", "numeric", "numeric", "factor", "Date"))
 
# wybieramy tylko mieszkania z Warszawy o powierzchni do 300m2
mieszkaniaKWW2011Warszawa <- mieszkaniaKWW2011[mieszkaniaKWW2011$miasto =="Warszawa", ]
mieszkaniaKWW2011Warszawa$dzielnica <- factor(mieszkaniaKWW2011Warszawa$dzielnica)
mieszkaniaKWW2011Warszawa <- mieszkaniaKWW2011Warszawa[mieszkaniaKWW2011Warszawa$powierzchnia<300, ]
 
# do analiz trafią tylko oferty złożone po 13 września 2007
dataStart <- as.Date("13-09-2007","%d-%m-%Y")
mieszkaniaKWW2011Warszawa2 <- mieszkaniaKWW2011Warszawa
mieszkaniaKWW2011Warszawa2 <- mieszkaniaKWW2011Warszawa2[mieszkaniaKWW2011Warszawa2$data >= dataStart,]
#
# dodajemy zmienną określającą miesiąc w którym złożona została oferta
mieszkaniaKWW2011Warszawa2$dataF <- factor(as.character(mieszkaniaKWW2011Warszawa2$data, "%Y-%m"))
 
#
# z wszystkich meiszkań do obiektu mieszkaniaKWW2011Warszawa3 wybieramy tylko te o średnim rozmiarze, tz o powierzchni d 49 do 68m2
mieszkaniaKWW2011Warszawa2$rozmiar <- cut(mieszkaniaKWW2011Warszawa2$powierzchnia, c(0,49,68,200), labels=c("do 49 m2", "od 49 do 68 m2", "powyzej 68 m2"))
mieszkaniaKWW2011Warszawa3 <- mieszkaniaKWW2011Warszawa2[which(mieszkaniaKWW2011Warszawa2$rozmiar == "od 49 do 68 m2"),]
 
#
# na potrzeby tego wykresu wybieramy tylko mieszkania z Zoliborza
mieszkaniaKWW2011Warszawa3Zoliborz    <- mieszkaniaKWW2011Warszawa3[which(mieszkaniaKWW2011Warszawa3$dzielnica == "Zoliborz"),]
 
at = seq(1000,24000,1000)

Dane są wczytane, czas na wykres. Wykorzystamy funkcję bwplot() z pakietu lattice.
Narysujemy jak zmieniają się ceny m2 średniej wielkości mieszkań w kolejnych miesiącach.

34
35
36

bwplot(cenam2~dataF,
         data=mieszkaniaKWW2011Warszawa3Zoliborz,
         main="Zoliborz")

Nie wygląda to najlepiej, zajmijmy się na początek osiami. Ponieważ cena jest zmienną silnie prawo skośną, przedstawimy ją na osi logarytmicznej. Ponieważ etykiety na osi OX zachodzą na siebie to zmniejszymy je i pokażemy pionowo.

37
38
39
40

bwplot(cenam2~dataF,
         data=mieszkaniaKWW2011Warszawa3Zoliborz,
         scales=list(y=list(log=T, at=at), x=list(rot=90, cex=0.6)),
         main="Zoliborz")

Przy takiej rozpiętości na osi OY trudno analizować delikatne zmiany w cenach mieszkań, więc w dalszej części zawęzimy zainteresowania do przedziału osi OY od 7 do 14k.

41
42
43
44
45

bwplot(cenam2~dataF,
         data=mieszkaniaKWW2011Warszawa3Zoliborz,
         scales=list(y=list(log=T, at=at), x=list(rot=90, cex=0.6)),
         ylim=c(7000, 14000),
         main="Zoliborz")

Trudno zauważyć jakiś trend. Dorysujmy więc krzywą trendu liniowego przedefiniowując funkcję rysującą panel. Nowa funkcja rysująca panel rysuje linie pomocnicze siatki, rysuje wykresy pudełkowe używając panel.bwplot() oraz dorysowuje linię odpornej regresji liniowej.

46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

bwplot(cenam2~dataF,
         data=mieszkaniaKWW2011Warszawa3Zoliborz,
         scales=list(y=list(log=T, at=at), x=list(rot=90, cex=0.6)),
         ylim=c(7000, 14000),
         main="Zoliborz",
         panel =  function(...) {
            panel.abline(h=log(at,10), col="grey85")
            panel.abline(v=4.5 + c(0,12,24,36), col="grey85")
 
            panel.bwplot(..., col="grey20")
            args = list(...)
            mod = rlm(args$y~as.numeric(args$x))
            panel.abline(mod)
          }
)

Trend liniowy liczony na wszystkich punktach to tylko jedno z możliwych podejść do zagadnienia oceny trendu. Dodajemy trend lokalnie ważony wielomianami stopnia pierwszego oraz trend liniowy wyznaczony tylko na podstawie median cen w kolejnych miesiącach.

61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83

bwplot(cenam2~dataF,
         data=mieszkaniaKWW2011Warszawa3Zoliborz,
         scales=list(y=list(log=T, at=at), x=list(rot=90, cex=0.6)),
         ylim=c(7000, 14000),
         main="Zoliborz",
         panel =  function(...) {
            panel.abline(h=log(at,10), col="grey85")
            panel.abline(v=4.5 + c(0,12,24,36), col="grey85")
 
            panel.bwplot(..., col="grey20")
            args <- list(...)
            mod  <- rlm(args$y~as.numeric(args$x))
            panel.abline(mod)
 
            mediany <- sapply(miesiace, function(x) median(args$y[args$x == x], na.rm=T))
            mod2    <- rlm(mediany~seq_along(miesiace))
            panel.abline(mod2)
            indtmp  <- c(1, length(mediany))
            llines(indtmp, mediany[indtmp])
 
            panel.loess(...)
          }
)

Mało czytelne są te linie trendu. Narysujemy je grubszą kreską i dodatkowo użyjemy kolorów z pakietu RColorBrewer, które powinny być przyjemniejsze dla oka.

84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107

bwplot(cenam2~dataF,
         data=mieszkaniaKWW2011Warszawa3Zoliborz,
         scales=list(y=list(log=T, at=at), x=list(rot=90, cex=0.6)),
         ylim=c(7000, 14000),
         main="Zoliborz",
         panel =  function(...) {
            kolory  <- brewer.pal(4, "Set1")
            panel.abline(h=log(at,10), col="grey85")
            panel.abline(v=4.5 + c(0,12,24,36), col="grey85")
 
            panel.bwplot(..., col="grey20")
            args   <- list(...)
            mod    <- rlm(args$y~as.numeric(args$x))
            panel.abline(mod, col=kolory[1], lwd=4, alpha=0.85)
 
            mediany <- sapply(miesiace, function(x) median(args$y[args$x == x], na.rm=T))
            mod2    <- rlm(mediany~seq_along(miesiace))
            panel.abline(mod2, col=kolory[2], lwd=4, alpha=0.85)
            indtmp  <- c(1, length(mediany))
            llines(indtmp, mediany[indtmp], col=kolory[4], lwd=4, alpha=0.85)
 
            panel.loess(..., col=kolory[3], lwd=4, alpha=0.85)
          }
)

Solą w oku są już tylko te niebieskie wykresy pudełkowe, lepiej wyglądałyby one w kolorze szarym, mają być tłem dla linii trendu. Aby zmienić kolor tych punktów używamy funkcji trellis.par.set() i trellis.par.get().

106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141

bwplot(cenam2~dataF,
         data=mieszkaniaKWW2011Warszawa3Zoliborz,
         scales=list(y=list(log=T, at=at), x=list(rot=90, cex=0.6)),
         ylim=c(7000, 14000),
         main="Zoliborz",
         panel =  function(...) {
            tmp <- trellis.par.get("plot.symbol")
              tmp$pch=19
              tmp$col="grey20"
              tmp$cex=1/2
            trellis.par.set("plot.symbol",tmp)
            tmp <- trellis.par.get("box.rectangle")
              tmp$col="grey20"
            trellis.par.set("box.rectangle",tmp)
            tmp <- trellis.par.get("box.umbrella")
              tmp$col="grey20"
            trellis.par.set("box.umbrella",tmp)
 
            kolory <- brewer.pal(4, "Set1")
            panel.abline(h=log(at,10), col="grey85")
            panel.abline(v=4.5 + c(0,12,24,36), col="grey85")
 
            panel.bwplot(..., col="grey20")
            args <- list(...)
            mod  <- rlm(args$y~as.numeric(args$x))
            panel.abline(mod, col=kolory[1], lwd=4, alpha=0.85)
 
            mediany <- sapply(miesiace, function(x) median(args$y[args$x == x], na.rm=T))
            mod2    <- rlm(mediany~seq_along(miesiace))
            panel.abline(mod2, col=kolory[2], lwd=4, alpha=0.85)
            indtmp  <- c(1, length(mediany))
            llines(indtmp, mediany[indtmp], col=kolory[4], lwd=4, alpha=0.85)
 
            panel.loess(..., col=kolory[3], lwd=4, alpha=0.85)
          }
)

Już jest nieźle. Ostatnia modyfikacja wykresu to dodanie opisu wykresu z liczbą procent o który zmieniła się cena mieszkania. Konstruujemy funkcję, która te procenty wyliczy i doklei do nazwy dzielnicy. Przy da się to do kolejnego przykładu, tutaj wyglądać może trochę sztucznie.

142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222

#
# ta funkcja doda do nazw dzielnic informacje o zmianie wartosci mieszkan w danej dzielnicy
# w badanym okresie czasu uzywajac czterech roznych metod
#
nazwyIProcenty <- function(mieszkaniaTmp) {
  # wektor z nazwami dzielnic ze zbioru danych i wektor nazw miesiecy
  ndzielnic  <- levels(mieszkaniaTmp$dzielnica2)
  miesiace   <- levels(mieszkaniaTmp$dataF)
  # anonimowa funkcja wyznacza cztery liczby
  # opisujace procentowa zmianę ceny
  ndzielnic <- paste(ndzielnic," ",sapply(ndzielnic, function(dzielnica) {
     tmp   <- mieszkaniaTmp[mieszkaniaTmp$dzielnica2==dzielnica,]
     skala <- diff(range(as.numeric(tmp$dataF)))
     # odporna (robust) regresja liniowa
     # szukamy trendu liniowego biorąc pod uwagę wszystkie ceny
     mod   <- rlm(tmp$cenam2 ~ I(as.numeric(tmp$dataF)/skala))$coeff
     proc1 <- round(1000*mod[2]/mod[1])
 
     mediany <- sapply(miesiace, function(miesiac) median(tmp[tmp$dataF == miesiac,"cenam2"]))
     # ponownie odporna regresja liniowa,
     # tym razem liczymy wspóczynniki trendu tylko na podstawie median dla kadego miesiąca
     mod2    <- rlm(mediany~I(seq_along(miesiace)/length(mediany)))$coeff
     proc2   <- round(1000*mod2[2]/mod2[1])
 
     # zmianę liczymy porównując medianę ceny z wrzenia 2007 z ceną z wrzesnia 2011
     # ceny w okresie przejciowym są ignorowane
     mediany <- na.omit(mediany)
     proc4   <- round(1000*(mediany[length(mediany)]-mediany[1])/mediany[1])
 
     # regresja lokalnie wygadzana wielomianami pierwszego stopnia
     # dodatkowe argumenty sa wymienione by zachowac zgodnosc z funkcj panel.loess
     tmp2  <- loess(tmp$cenam2 ~ I(as.numeric(tmp$dataF)/skala), span = 2/3, degree = 1, family = "symmetric" )$fitted
     # zmiane liczymy porówując oszacowany trend na wrzesien 2007
     # z oszacowanym trendem na wrzesien 2011
     tmp2  <- tmp2[order(I(as.numeric(tmp$dataF)/skala))]
     proc3 <- round(1000*(tmp2[length(tmp2)] - tmp2[1])/tmp2[1])
 
     paste(proc1/10, "% / ", proc2/10, "% / ", proc3/10, "% / ", proc4/10, "% ", sep="")
  }), sep="")
  ndzielnic
}
 
levels(mieszkaniaKWW2011Warszawa3$dzielnica2) = nazwyIProcenty(mieszkaniaKWW2011Warszawa3)
mieszkaniaKWW2011Warszawa3Zoliborz = mieszkaniaKWW2011Warszawa3[which(mieszkaniaKWW2011Warszawa3$dzielnica == "Zoliborz"),]
 
bwplot(cenam2~dataF,
         data=mieszkaniaKWW2011Warszawa3Zoliborz,
         scales=list(y=list(log=T, at=at), x=list(rot=90, cex=0.6)),
         ylim=c(7000, 14000),
         main=as.character(mieszkaniaKWW2011Warszawa3Zoliborz$dzielnica2)[1],
         panel =  function(...) {
            tmp <- trellis.par.get("plot.symbol")
              tmp$pch=19
              tmp$col="grey20"
              tmp$cex=1/2
            trellis.par.set("plot.symbol",tmp)
            tmp <- trellis.par.get("box.rectangle")
              tmp$col="grey20"
            trellis.par.set("box.rectangle",tmp)
            tmp <- trellis.par.get("box.umbrella")
              tmp$col="grey20"
            trellis.par.set("box.umbrella",tmp)
 
            kolory <- brewer.pal(4, "Set1")
            panel.abline(h=log(at,10), col="grey85")
            panel.abline(v=4.5 + c(0,12,24,36), col="grey85")
 
            panel.bwplot(..., col="grey20")
            args <- list(...)
            mod  <- rlm(args$y~as.numeric(args$x))
            panel.abline(mod, col=kolory[1], lwd=4, alpha=0.85)
 
            mediany <- sapply(miesiace, function(x) median(args$y[args$x == x], na.rm=T))
            mod2    <- rlm(mediany~seq_along(miesiace))
            panel.abline(mod2, col=kolory[2], lwd=4, alpha=0.85)
            indtmp  <- c(1, length(mediany))
            llines(indtmp, mediany[indtmp], col=kolory[4], lwd=4, alpha=0.85)
 
            panel.loess(..., col=kolory[3], lwd=4, alpha=0.85)
          }
)

Wykres gotowy, w kolejnym odcinku pokażemy jak wygenerować taki wykres dla wszystkich dzielnic.

Rysujemy rozkład cen krok po kroku, część 2

Kontynuując temat z wczoraj, narysujemy rozkład cen mieszkania dla każdej z  dzielnic Warszawy.

Punktem wyjścia jest przygotowanie danych, wykonajmy pierwsze 32 linie tak jak w poprzednim wpisie.

Aby wyświetlić na rożnych panelach dane dla kolejnych dzielnic, wystarczy zmodyfikowac formułę na cenam2~data|dzielnica, oraz za zabiór danych wskazać mieszkaniaKWW2011Warszawa2.

33
34
35
36
37
38
39
40
41
42
43
44
45
46

at = seq(1000,24000,1000)
 
xyplot(cenam2~data|dzielnica, group=rozmiar,
        data=mieszkaniaKWW2011Warszawa2,
        scales=list(y=list(log=T, at=at)), ylim=c(6000, 16000),
        type=c("smooth"), pch='.', lwd=3,
        auto.key=list(space="top", columns=3, pch=19),
        panel=function(...) {
  panel.abline(h=log(at,10), col="grey85")
  panel.abline(v=365.25 * (38:41), col="grey85")
  panel.xyplot(list(...)$x, list(...)$y, type="p",col="grey", cex=1/4)
  panel.xyplot(list(...)$x, list(...)$y, type="r",col="black", lwd=2,lty=2)
  panel.xyplot(...)
})

Kolejne panele przedstawiają kolejne dzielnice, ale ich kolejność jest alfabetyczna. Taka sama jak kolejność poziomów zmiennej czynnikowej dzielnica. Nie zawsze kolejność alfabetyczna będzie najlepsza. Użyjemy funkcji reorder by zmienić kolejność poziomów tak by odpowiadała medianie ceny metra kwadratowego w danej dzielnicy. Kod generujący obrazek będzie taki sam, zmieni się tylko kolejność dzielnic.

47
48
49
50
51
52
53
54
55
56
57
58
59
60

mieszkaniaKWW2011Warszawa2$dzielnica2 <- reorder(mieszkaniaKWW2011Warszawa2$dzielnica, mieszkaniaKWW2011Warszawa2$cenam2, median)
 
xyplot(cenam2~data|dzielnica2, group=rozmiar,
        data=mieszkaniaKWW2011Warszawa2,
        scales=list(y=list(log=T, at=at)), ylim=c(6000, 16000),
        type=c("smooth"), pch='.', lwd=3,
        auto.key=list(space="top", columns=3, pch=19),
        panel=function(...) {
  panel.abline(h=log(at,10), col="grey85")
  panel.abline(v=365.25 * (38:41), col="grey85")
  panel.xyplot(list(...)$x, list(...)$y, type="p",col="grey", cex=1/4)
  panel.xyplot(list(...)$x, list(...)$y, type="r",col="black", lwd=2,lty=2)
  panel.xyplot(...)
})

Dla niektórych dzielnic jest mało punktów, co powoduje, że trudno mieć zaufanie do wyznaczonego trendu. Tym razem usuniemy te dzielnice, dla których nie ma przynajmniej 2000 wierszy. Kod generujący wykres jest bez zmian, usuwamy tylko obserwacje z dzielnic w których obserwacji było mniej niż 2k.

61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

usun <- names(which(table(mieszkaniaKWW2011Warszawa2$dzielnica)<2000))
mieszkaniaKWW2011Warszawa2 <- mieszkaniaKWW2011Warszawa2[!(mieszkaniaKWW2011Warszawa2$dzielnica %in% usun),]
mieszkaniaKWW2011Warszawa2$dzielnica2 <- factor(mieszkaniaKWW2011Warszawa2$dzielnica2)
 
 
xyplot(cenam2~data|dzielnica2, group=rozmiar,
        data=mieszkaniaKWW2011Warszawa2,
        scales=list(y=list(log=T, at=at)), ylim=c(6000, 16000),
        main="",
        type=c("smooth"), pch='.', lwd=3,
        auto.key=list(space="top", columns=3, pch=19),
        panel=function(...) {
  panel.abline(h=log(at,10), col="grey85")
  panel.abline(v=365.25 * (38:41), col="grey85")
  panel.xyplot(list(...)$x, list(...)$y, type="p",col="grey", cex=1/4)
  panel.xyplot(list(...)$x, list(...)$y, type="r",col="black", lwd=2,lty=2)
  panel.xyplot(...)
})

Rysujemy rozkład cen krok po kroku

Dostałem maila od Grzegorza L. z pytaniem jak zrobić krok po kroku wykresy z wpisów http://smarterpoland.pl/?p=169 i http://smarterpoland.pl/?p=172. W sumie cztery wykresy, każdy z nich postaram się opisać jak powstawał krok po kroku. Dziś zajmiemy się pierwszym z nich czyli zmiana cen mieszkań w dzielnicy Zoliborz w rożnych grupach wielkości mieszkań.

Zaczniemy od opisu jak dane były przygotowane a następnie pokażemy jak zrobić wykres krok po kroku. Skrypt R z wszystkimi komendami poniżej opisanymi znajduje się tutaj.

Część 1: Przygotowanie danych
Poniższy skrypt wczytuje do R dane odczytując je bezpośrednio z Internetu, następnie wybiera tylko dane dla Warszawy, dla dzielnicy Zoliborz. Wybiera tylko mieszkania o powierzchni do 300m2 i zgłoszone w ostatnich 3 latach. Następnie dodajemy do danych nowa zmienna opisującą wielkość mieszkania w jednej z trzech grup, małe do 49 m2, średnie od 49 do 68m2 i duże, o powierzchni ponad 68m2. Następnie dodatkowo wybieramy tylko podzbiór mieszkań z Zoliborza.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

library(lattice)
#
# wczytujemy dane o mieszkaniach bezposrednio z internetu
# moze chwile potrwac
mieszkaniaKWW2011 <-
     read.table("http://tofesi.mimuw.edu.pl/~cogito/smarterpoland/mieszkaniaKWW2011/mieszkaniaKWW2011.csv",
         row.names=NULL, sep=";", header=TRUE,
         colClasses=c("factor", "factor", "numeric", "numeric", "factor", "numeric", "numeric", "factor", "Date"))
 
#
# wybieramy mieszkania tylko z Warszawy
# usuwamy informacje o dzielnicach z innych miast (funkcja factor to zalatwi)
# wybrane mieszkania zapisujemy do zmiennej mieszkaniaKWW2011Warszawa
#
mieszkaniaKWW2011Warszawa <- mieszkaniaKWW2011[mieszkaniaKWW2011$miasto =="Warszawa", ]
mieszkaniaKWW2011Warszawa$dzielnica <- factor(mieszkaniaKWW2011Warszawa$dzielnica)
 
#
# zawezimy horyzont do ofert sprzedazy z ostatnich 4 lat
# pominiemy oferty sprzed 13 wrzesnia 2007
# usuniemy tez mieszkania o powierzchni powyzej 300m2 (jest ich niewiele a psuja niektore wizualizacje)
#
dataStart <- as.Date("13-09-2007","%d-%m-%Y")
mieszkaniaKWW2011Warszawa2 <- mieszkaniaKWW2011Warszawa
mieszkaniaKWW2011Warszawa2 <- mieszkaniaKWW2011Warszawa2[mieszkaniaKWW2011Warszawa2$data >= dataStart,]
mieszkaniaKWW2011Warszawa2 <- mieszkaniaKWW2011Warszawa2[mieszkaniaKWW2011Warszawa2$powierzchnia<300, ]
 
#
# powierzchnie mieszkania zamienimy na zmienna jakosciowa
# wystepujaca na trzech poziomach, do 49m2, od 49 do 68 m2, powyzej 68m2
#
mieszkaniaKWW2011Warszawa2$rozmiar <- cut(mieszkaniaKWW2011Warszawa2$powierzchnia, c(0,49,68,200), labels=c("do 49 m2", "od 49 do 68 m2", "powyzej 68 m2"))
 
#
# wybierzemy tylko mieszkania z jednej dzielnicy, Zoliborz
#
mieszkaniaKWW2011Warszawa2Zoliborz <- mieszkaniaKWW2011Warszawa2[which(mieszkaniaKWW2011Warszawa2$dzielnica == "Zoliborz"),]

Część 2: Generowanie obrazków.

Ok, mamy przygotowane dane, teraz czas na zrobienie rysunku. Zaczniemy od prostego użycia funkcji xyplot{lattice}.

38
39
40
41
42
43
44
45

#
# pierwsze podejscie
# podajemy zmienne ktore maja byc narysowane na osi X i Y
# podajemy zmienna grupujaca i tytul wykresu
#
xyplot(cenam2~data, group=rozmiar,
        data=mieszkaniaKWW2011Warszawa2Zoliborz,
        main="Zoliborz")

Miło ze strony funkcji xyplot(), że jest w stanie narysować na osi OX zmienna typu Date, ale sam wykres nie powala na kolana.
Za bardzo w oczy rzucają się punkty, których jest tak dużo że nic ciekawego nie widać. Zmniejszymy więc wielkość punktów zaznaczając je jednopikselową kropka (pch=”.”), dodamy linie pomocnicze siatki i wygładzoną linię trendu (type=”g” i „smooth”).

46
47
48
49
50
51
52

#
# dodajemy pomocnicze linie siatki i wygladzona linie trendu
# zmieniamy tez sposob rysowania punktow na jednopunktowe kropki
#
xyplot(cenam2~data, group=rozmiar,
        data=mieszkaniaKWW2011Warszawa2Zoliborz,
        main="Zoliborz", type=c("p","g","smooth"), pch=".")

Już trochę lepiej, ale wciąż niewiele widać. Skala OY jest liniowa a dla prawoskośnych zmiennych dobrym pomysłem jest skala logarytmiczna. Dodajemy argument scales, który pozwala na zmianę osi OY na logarytmiczną.

53
54
55
56
57
58
59

#
# zmieniamy os OY na logarytmiczna
#
xyplot(cenam2~data, group=rozmiar,
        data=mieszkaniaKWW2011Warszawa2Zoliborz,
        scales=list(y=list(log=T)),
        main="Zoliborz", type=c("p","g","smooth"), pch=".")

Co nam się nie podoba tym razem? Wiele rzeczy. Po pierwsze przydała by się legenda. Dodamy ją argumentem auto.key, przy okazji określimy gdzie i jak ma być legenda rysowana, choć wcale nie musimy tego określać, wystarczyłoby auto.key=T. Druga rzecz to widoczność linii trendu w każdej z grup. Są mało widoczne z uwagi na grubość linii i z uwagi na duży rozstrzał na osi OY. Zmienimy i to i to, podamy grubość linii lwd=3 i ustalimy zakres zmienności na osi OY na 6-16k/m2. ostatnia poważna zmienna będzie dotyczyła znaczników na osi OY. Poprzednio wyglądały one słabo, jakieś potęgi 10 (ile to 10^4.2?), tym razem jawnie wskażemy w którym punktach mają pojawiać się znaczniki osi podając argument at.

60
61
62
63
64
65
66
67
68
69
70
71

#
# zmieniamy zakres osi OY
# jawnie wskazujemy gdzie maja byc rysowane znaczniki na osi OY
# dodajemy tez legende umieszczona na gorze wykresu rysowana poziomo (w trzech kolumnach, ale dla trzech wartosci wychodzi na jedno)
#
at = seq(1000,24000,1000)
xyplot(cenam2~data, group=rozmiar,
        data=mieszkaniaKWW2011Warszawa2Zoliborz,
        scales=list(y=list(log=T, at=at)), ylim=c(6000, 16000),
        main="Zoliborz",
        type=c("p", "g", "smooth"), pch='.', lwd=3,
        auto.key=list(space="top", columns=3, pch=19))

Wykres nabiera kształtów. Nadpiszemy teraz funkcję rysującą panel, w taki sposób by dodatkowo rysować trend wyznaczony metodą regresji liniowej. Rysujemy go na czarno przerywaną linią wykorzystując funkcję panel.xyplot().

72
73
74
75
76
77
78
79
80
81
82
83
84
85

#
# zmieniamy funkcje rysujaca panel, teraz ta funkcja dorysowuje krzywa regresji liniowej 
# przerywana linia
#
xyplot(cenam2~data, group=rozmiar,
        data=mieszkaniaKWW2011Warszawa2Zoliborz,
        scales=list(y=list(log=T, at=at)), ylim=c(6000, 16000),
        main="Zoliborz",
        type=c("p", "g", "smooth"), pch='.', lwd=3,
        auto.key=list(space="top", columns=3, pch=19),
        panel=function(...) {
  panel.xyplot(list(...)$x, list(...)$y, type="r",col="black", lwd=2,lty=2)
  panel.xyplot(...)
})

Czas na ostatnie modyfikacje. Co zmienimy? Po pierwsze punkty w legendzie są rysowane pustymi okręgami, a lepiej wyglądać będą wypełnione. Poniżej używam funkcji trellis.par.set() by o zmienić.
Po drugie, porzednie linie pomocnicze siatki nie były rysowane w tych samych miejscach co znaczniki na osi OY, usuwam więc argument type=”g” a zamiast tego w funkcji rysującej panel ręcznie dorysowuje linie pomocnicze w odpowiednich miejscach i pionowo i poziomo funkcją panel.abline(). Po trzecie zmieniam sposób rysowania punktów. Używam do tego funkcji panel.xyplot() z type=”p” i cex=1/4 dzięki czemu każda oferta sprzedaży jest zaznaczona jednym niewielkim szarym punktem.

86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109

#
# funkcja trellis.par.set() zmieniamy ksztalt io kolor kropek w legendzie
# z argumentu type usuwamy rysowanie punktow i lini pomocnicznych
# z tego tez powodu dodajemy do funkcji rysujacej panel kilka dodatkowych instrukcji
# rysujacych punkty, linie regresji liniowej i linie pomocnicze siatki
#
tt = trellis.par.get("superpose.symbol")
tt$pch=19
tt$col = c( "#0080ff", "#ff00ff", "darkgreen")
trellis.par.set("superpose.symbol",tt)
 
xyplot(cenam2~data, group=rozmiar,
        data=mieszkaniaKWW2011Warszawa2Zoliborz,
        scales=list(y=list(log=T, at=at)), ylim=c(6000, 16000),
        main="Zoliborz",
        type=c("smooth"), pch='.', lwd=3,
        auto.key=list(space="top", columns=3, pch=19),
        panel=function(...) {
  panel.abline(h=log(at,10), col="grey85")
  panel.abline(v=365.25 * (38:41), col="grey85")
  panel.xyplot(list(...)$x, list(...)$y, type="p",col="grey", cex=1/4)
  panel.xyplot(list(...)$x, list(...)$y, type="r",col="black", lwd=2,lty=2)
  panel.xyplot(...)
})

I koniec.

Jeszcze mała autoreklama. Więcej o pakiecie lattice można przeczytać w drugim wydaniu ,,Przewodnika po pakiecie R” wydanego przez wydawnictwo GiS w tym roku (jeden podrozdział dostępny w Internecie) lub w książce Lattice Multivariate Data Visualization with R, w całości poświęconej temu pakietowi (zobacz tutaj). Oczywiście można też znaleźć sporo informacji korzystając po prostu z googla.