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.