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.

Dodaj komentarz

Twój adres email nie zostanie opublikowany. Pola, których wypełnienie jest wymagane, są oznaczone symbolem *