Październik 2011 – SmarterPoland.pl

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.

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)
          }
)

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.

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

# 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.

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)

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.

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

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.

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

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.

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")

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.

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)
          }
)

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.

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(...)
          }
)

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.

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)
          }
)

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().

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)
          }
)

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.

#
# 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)
          }
)

# 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.

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(...)
})

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.

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(...)
})

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.

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(...)
})

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.

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"),]

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}.

#
# 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")

# 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”).

#
# 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=".")

# 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ą.

#
# 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=".")

# 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.

#
# 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))

# 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().

#
# 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(...)
})

# 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.

#
# 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(...)
})

# 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.

Gartner chart

Widziałem kilka dni temu świetną wizualizację. Nie wykonaną w R, nie dotyczącą polskich danych ale wciąż wystarczająco świetną by ją tu wrzucić. Świetna oznacza że jest w niej na tyle interesującej treści by przez przynajmniej kilka minut analizować ją bez znudzenia.

Firma Gatner przygotowuje co roku wizualizacje jaka technologia jest w jakim stadium rozwoju, jakie są oczekiwania i jak daleko jej do wielkoskalowych rzowiązań.

Diagnoza Społeczna 2011

Diagnoza społeczna to badanie prowadzone przez radę monitoringu społecznego od roku 2000. Więcej informacji o tym badaniu można znaleźć na stronie http://diagnoza.com/. Jest to badanie panelowe, dane zbierane są co 2-3 lata. Niedawno pojawiły się dane z edycji 2011. Badane jest bardzo wiele parametrów, można naprawdę prześledzić co ciekawego działo się w Polsce przez ostatnie 11 lat. Te dane nadają się świetnie na ćwiczenia ze statystycznej analizy danych dla studentów i nie tylko. Tydzień temu Paweł Teisseyre z IPIPANu używał tego zbioru danych do demonstrowania regularyzowanej wersji regresji logistycznej w R na WZUR 4.0.

Dane są publicznie dostępne. Niestety na stronie projektu dane są w postaci plików programu SPSS. Na potrzeby tego bloga zostały przekonwertowane do formatu programu R.

Katalog z danymi znajduje się tutaj.

Dane podzielone są na dwa zbiory, z opisem gospodarstw domowych i opisem osób o wieku ponad 16 lat zamieszkujących w tych gospodarstwach.

Dane o gospodarstwach można ściągnąć w postaci pliku RData, pliku w formacie csv oraz pliku z opisami kolumn, w zbiorze danych jest 20655 wierszy i 1820 kolumn.

Dane o osobach można ściągnąć w postaci pliku RData, pliku w formacie csv oraz pliku z opisami kolumn, w zbiorze danych jest 65373 wierszy i 2427 kolumn.

Skrypt wczytujący dane dostępny jest tutaj.

Na stronach projektu znaleźć można obszerne raporty które na kilkuset stronach prezentują tysiące wniosków i dziesiątki rysunków. Postaram się w najbliższej przyszłości umieścić kilka celowanych wizualizacji tak by na jednym rysunku upakować całą historię. Jeżeli studenci coś ciekawego na tym zbiorze danych zrobią to też dodam do bloga.

Cytowanie: Rada Monitoringu Społecznego (2011). Diagnoza społeczna: zintegrowana baza danych. www.diagnoza.com 20-X-2011;

Ceny mieszkań spadły o 20%!!! albo wzrosły o 5%!!! albo się nie zmieniły!!!

Ostatnio w gazetach pojawiło się kilka artykułów w których aby przyciągnąć więcej czytelników zastosowano bardzo głupie tytułu w stylu ,,Realne ceny mieszkań spadły/ mieszkania potaniały o 25% także w Warszawie”. Artykuły np. [gazeta na podstawie open finance], [raport oferty net], [wyborcza biz finanse] mają tą cechę wspólną, że podają liczby bez żadnej refleksji na temat tego jak liczby uzyskano i jakie ma to znaczenie. W każdym z przypadków badano zmianę ceny w trakcie ostatnich czterech lat (czyli od tzw górki 2007) tak więc i poniższe rozważania oprzemy na tym okresie.

Celem tego wpisu jest nie polemika z tym co to znaczy realna cena mieszkania, ani z tym czy spadek wynosił 25% czy nie, ale z tym, że jeżeli nie porównujemy ceny tego samego mieszkania to mówiąc o zmianie wartości należy dokładnie opisać zmianę wartości czego liczymy, bo wyniki mogą być bardzo różne. Dla każdej dzielnicy Warszawy na 8 sposobów policzymy jak zmieniały się ceny przez ostatnie 4 lata. I w każdym przypadku otrzymamy 8 różnych liczb opisujących zmianę ceny.

We wczorajszym wpisie dyskutowaliśmy różnice pomiędzy cenami m2 w dużych i małych mieszkaniach różnych dzielnic. Pisząc o strukturze koszyka pokazaliśmy, że badanie ,,tak po prostu” średniej mieszkań sprzedanych/oferowanych nie ma sensu jeżeli nie uwzględnimy jakie mieszkania (powierzchnia/lokalizacja/stan itp) wchodzą w skład koszyka. Poniżej będziemy podawać różne liczby ale nie po to by dowieść że cena zmieniła się o X%, ale by zaznaczyć jakie problemy wiążą się z analiza zmiany ceny. Wszystkie ceny to ceny ofertowe, nie pokrywają one całego rynku i pewnie nie są całkowicie reprezentacyjne. Nie mniej powinno być fajnie popatrzeć na obrazki.

Przyjrzyjmy się dzielnicy Żolibórz. Oś OY jest w skali logarytmicznej, pionowe szare linie oznaczają koniec roku. Na poniższym wykresie dla każdego miesiąca począwszy od września 2007 wyznaczamy wykres pudełkowy opisujący zmiany cen w danym miesiącu. W środku pudełka czarną kropką oznaczamy medianę cen w danym miesiącu. Następnie czterema liniami oznaczamy następujące trendy:

* linią czerwoną zaznaczamy trend liniowy wyznaczony na bazie wszystkich obserwacji używając regresji odpornej (funkcja rlm() w R), regresja odporna ma nas uchronić przed wpływem pojedynczych obserwacji odstających, nie uchroni nas jednak przed obciążeniem wynikającym z tego, że większość ofert sprzedaży dotyczy ostatniego roku.\ i to ten okres będzie głównie wpływał na ocenę trendu liniowego. W nagłówku wykresu podano cztery procenty, pierwszy z nich opisuje o ile ten trend liniowy się zmienił prze ostatnie cztery lata (-15.6%).

* linią niebieską zaznaczamy trend liniowy wyznaczony na bazie median cen używając regresji odpornej. Czyli najpierw wyznaczamy mediany cen w każdym miesiącu a później wyznaczamy trend. Tym razem trend nie jest już dominowany przez ceny z ostatniego roku. W nagłówku wykresu druga liczba to zmiana tego trendu przez ostatnie 4 lata (-7.6%, czyli około połowa poprzedniej liczby),

* linia zielona i trzecia liczba odpowiadają trendowi wyznaczonemu za pomocą lokalnego wygładzania. Ponieważ używanie trendu liniowego wydaje się być słabo uzasadnione, porównujemy trend liniowy z trendem wyznaczonym przez funkcję loess z parametrami rodzina = symmetric, span = 2/3, stopień wielomianu = 1. Oczywiście dla innych parametrów otrzymalibyśmy inne wygładzenie, te jednak wygląda w miarę dobrze. Trzecia liczba nagłówku opisuje o ile się zmieniła ocena trendu we wrześniu 2011 w stosunku do września 2007 (-11.1%).

* linia fioletowa łączy medianę ceny m2 mieszkań ferowanych we wrześniu 2007 z medianą ceny m2 mieszkań oferowanych we wrześniu 2011. Ponieważ danych o sprzedaży we wrześniu 2007 jest niewiele, ta linia wykazuje spora zmienność pomiędzy dzielnicami. Dla Żoliborza zmiana medianowej ceny m2 pomiędzy wrześniem 2007 a wrześniem 2011 to (-3.6%).

Mamy więc cztery liczby, począwszy od spadku niecałych 4% w cztery lata do spadku o ponad 15%. A to jeszcze nie wszystko. Te same analizy możemy zrobić uwzględniając strukturę wielkości mieszkań. Odrzućmy więc mieszkania najmniejsze i największe (te największe najwięcej straciły) i zajmujmy się tylko cenami m2 mieszkań o powierzchni od 49 do 68m2. Poniżej odpowiadający im wykres z oznaczeniami jak powyżej,

Jak widzimy ograniczając się tylko do segmentu mieszkań o średniej powierzchni otrzymujemy 4 inne trendy i cztery inne % opisujące zmianę tych trendów.

Która z tych 8 liczb lepiej opisuje zmianę cen mieszkań? Oczywiście to zależy jaki aspekt zmiany chcemy analizować. Ceny mieszkań zmieniają się w sposób złożony. Mnie najbardziej podoba się trend lokalnie wygładzany, nie zakłada on liniowości tak jak linie czerwone i niebieskie a jednocześnie czerpie informacje z prawie wszystkich danych (w przeciwieństwie do linii fioletowej która patrzy tylko na skrajne miesiące).

W przypadku trendu nieliniowego nadużyciem jest stwierdzenie, że ceny spadały o X%, ponieważ w różnych okresach zachowywały się różnie.

Poniżej wyniki dla wszystkich dzielnic, pierwszy wykres bazuje na wszystkich obserwacjach, drugi dotyczy tylko mieszkań o średniej wielkości.

Czym się różnią ceny mieszkań na Żoliborzu od cen mieszkań na Bemowie

Dzisiejszy wpis to przygotowanie gruntu do jutrzejszego, w ktorym rozliczymy sie ze zmianami cen mieszkan. Od jakiegos czasu na różnych stronach można znaleźć artykuły o nazwach ”realne ceny mieszkań spadły o X %’’ (dobrze że realne, co to by byly gdyby spadaly ceny nierealne). Celem tego i kolejnego wpisu jest pokazanie że zmianę ceny można liczyć na wiele różnych sposobów otrzymując wiele różnych wyników. Więc podanie zmiany nie mówiąc dokładnie jak była ona liczona to zwykłą propaganda.

Będziemy korzystać ze zbioru danych mieszkaniaKWW2011, przeanalizujemy tylko ceny z miasta Warszawa dla mieszkan o powierzchni do 300m2. Na początek przyjrzymy się dwóm dzielnicom. Jednej w której mieszkam (Bemowo) i jednej w której ceny się ciekawie zachowują (Żoliborz).

Przedstawimy zmianę trendu ceny mieszkań w czasie w rozbiciu na trzy grupy wielkości mieszkań. Kwantyle rzędu 1/3 i 2/3 z rozkładu wielkości mieszkań wynoszą 49m2 i 68m2, wiec będziemy analizować cenę m2 w grupie mieszkań w trzech grupach: do 49m2, w grupie mieszkań dużych o powierzchni powyżej 68m2 i w grupie pozostałych – średnich co do wielkości mieszkań.

Na poniższym wykresie każdy punkt to jedna oferta sprzedaży. Linia zielona, różowa i niebieska oznaczają lokalne wygładzenie wielomianowe dla ceny m2 dla różnych grup wielkości mieszkań, a czarna linia przerywana to globalny trend w cenach mieszkań bez uwzględnienia struktury mieszkań. Struktura jest istotna ponieważ hipotetycznie, jeżeli mieszkania duże maja niższa cene za m2, i przed czterema laty mieszkań dużych bylo sprzedawanych mniej niż teraz to zaobserwujemy zmiane w średniej cenie nie ze względu na rzeczywistą zmianę ceny ale ze wzgledu na zmianę struktury sprzedawanych mieszkań. Podobnie z analizą cen w Warszawie. Jeżeli w tym roku sprzedaje sie więcej mieszkań na obrzeżach Warszawy niż przed czterema laty (z roku na rok buduje się średnio coraz dalej bo tam jeszcze jest miejsce) to zmiana średniej ceny będzie zwiazana z tym ze w koszyku jest coraz więcej mieszkań dalekich od centrum a wiec tańszych.

No dobrze, to tyle tytułem teorii a teraz obrazki dla dwoch obiecanych dzielnic.

Dla Żoliborza jak widzimy ceny m2 dla mieszkań dużych sa wyzsze niz dla malych, moze byc to zwiazane ze duze mieszkania to juz apartamenty dlatego cenę ich m2 ciezko porównywać z cena m2 mieszkania o standardowym wykonczeniu. Pdobnie beda zachowywaly sie mieszkania w Śródmieściu. Te apartamenty tez najbardziej straciły na wartości. Na Bemowie jest inaczej. Nie buduje sie raczej apartamentow w dzielnicy sypialni wiec tutaj to male mieszkania maja wyższą cenę za m2. W obu przypadkach ceny maja tendencje do spadania ale w kazdej kategorii tempo spadku jest inne.

Poniżej prezentujemy wykresy dla wszystkich dzielnic Warszawy.

Zmiany zamożności Polaków na kolorowo

Trzy dni temu pisaliśmy o tym jak zmienia się zamożność gospodarstw badanych w ramach Diagnozy Społecznej (ponad 20 tys gospodarstw).
Generalny wniosek jest taki, że coraz więcej gospodarstw domowych poprawiło swój standard przez ostatnie 6 lat (porównywaliśmy wyniki z lat 2005 i 2011). Zobaczmy jak ta sytuacja wygląda w rożnych województwach.
Na poniższych wykresach wykonamy analizę gradacyjną, kolorami zaznaczając wyniki różnych województw.

Aby było czytelniej wybraliśmy pięć województw, dla których wyniki były ciekawe. Są to województwa Zachodni-pomorskie, Mazowieckie, Dolnośląskie, Lubuskie i Świętokrzyskie. Etykiety zamiast przy punktach zostały umieszczone w prawej dolnej legendzie. Jak czytać te wykresy? Ponieważ odpowiedzi są w skali uporządkowanej możemy interpretować bezpośrednio położenie k-tego punktu, anie tylko ścieżkę do niego prowadzącą (jak na poprzednich przykładach).
Etykieta 5 oznacza, że wystarcza tylko na najtańsze jedzenie, ubranie, opłaty, kredyt. Ponieważ na wykresie współrzędne punktów to skumulowane częstości dla danego i niższych poziomów, więc współrzędne punktów z etykietą 5 oznaczają frakcję osób, którym starcza tylko na najtańsze jedzenie, kredyt, ubrania lub i na to nie. W województwie zachodnio-pomorskim w roku 2005 takich osób było około 40% (współrzędna OX brązowej 5), ale w roku 2011 było już takich osób niewiele ponad 20%. Praktycznie w każdym z narysowanych województw współrzędna OY dla cyfry 7 to 80%, co oznacza, że w roku 2011 80% gospodarstw w tych województwach deklarowało, że muszą żyć oszczędnie, bardzo oszczędnie a czasem i to nie wystarcza. W roku 2005 procent takich deklaracji był wyższy w każdym z województw, najwięcej spadł w Świętokrzyskim z około 90%.

Im dalej punkt od przekątnej tym większa zmiana do tego poziomu. Przykładowo różowa 6 ilustruje, że w województwie Lubuskim osób które żyją bardzo oszczędnie a i to czasem nie starcza było ponad 60% w roku 2005 a w 6 lat później było ich już tylko około 40%.

Mniejsze zmiany dotknęły województw Mazowieckiego i Dolnośląskiego.

Na zakończenie wykres dla wszystkich województw, dosyć gęsty, 16 krzywych zachodzi na sobie co utrudnia odcyfrowywanie wyników. W każdym województwie krzywa jest w większości pod przekątną co znaczy, że jest raczej lepiej niż przed 6 laty.

Jak pokazać 27 by wyglądało jak 20

Wczoraj studenci z proseminarium prezentowali różne znalezione grafiki dotyczące wyborów i partii politycznych.
Poniższa przypadła mi najbardziej do gustu, pochodzi ze strony szczecinek24.

Co jest w tej grafice ciekawego? Przyjrzyjmy się wysokościom słupków. Po starannym zmierzeniu w programie graficznym mamy:

słupek PO 99 pixeli

słupek PiS 65 pixeli

słupek SLD 25 pixeli

słupek PSL 11 pixeli

Gdybyśmy przeskalowali te pixele tak by PO miało 31 jednostek, to kolejne słupki miałyby następującą liczbę jednostek

PO 31, PiS 20, SLD 8, PSL 3

Ale nijak to nie pasuje do procentów z tego obrazka! Powinno być 31, 27, 16, 7.5

Zobaczmy jak wyglądałby ten wykres, gdyby wysokość słupków odpowiadała wpisanym procentom

Różnica pomiędzy Po i PiS nie wygląda tak samo. Poparcie dla SLD po prawej stronie też wygląda lepiej niż po lewej

Btw: dosyć często zdarza mi się krytykować wykresy kołowe za to, że oszukują percepcje. Na co najczęściej rozmówcy odpowiadają, co z tego że oszukują percepcję, skoro i tak podawane są procenty. Więc procenty można sobie porównać. Powyższy przykład powinien pokazać, że nawet jeżeli podamy procenty to i tak przy szybkim czytaniu opieramy się na grafice i a nie na liczbach. Patrząc na lewy wykres mamy wrażenie że poparcie dla PSL jest znikome, dla SLD małe a dla PO znacznie większe niż dla pozostałych partii. Gdybyśmy porównywali liczby albo słupki z prawego wykresu to już byśmy takiego wrażenia nie mieli.