W poprzedni piątek, w tym wpisie, opisałem trzecią wersję zbioru danych o cenach aut używanych.  W trzech najbliższych wpisach przyjrzymy się bliżej temu zbiorowi danych. Wpisy będą znacznie bardziej technicznie niż poprzednie, ale mam nadzieję że ciekawie dla osób zaczynających przygodę z modelowaniem statystycznym.

Dzisiaj przedstawię kilka umiarkowanie złożonych podejść do modelowania ceny samochodu z użyciem modeli liniowych. Będzie kod w R i trochę narzekań na niezrównoważony zbiór danych. Kolejny wpis przedstawiać będzie przykładowe wizualizacje zbioru danych o cenach aut. A w trzecim wpisie przedstawię krótką prezentację wykonaną z użyciem narzędzia dostępnego online na stronie prezi.com. Narzędzie to jest bezpłatne w wersji średniej dla studentów i nauczycieli akademickich. Ciekawe jestem jak wypadnie prezentacja w nim wykonana, czy lepiej niż w Beamerze?

Pomimo iż większość (żmudnej) pracy związanej z wyborem modelu nie zostanie przedstawiona, i tak będzie dzisiaj bardzo technicznie.

Zaznaczę jeszcze, że modelowana jest cena używanego auta a nie zmiana ceny czy spadek wartości. Modelując cenę mogę dodać efekt wieku auta i mówić, że auta o rok starszy są o x PLN tańsze, ale to są inne auta! W innej wersji i innym wyposażeniem. Efekt wyposażenia można jeszcze usunąć ale efektu związanego z wersja nie. Można porównywać średnią cenę Passata z roku 2006 z ceną Passata z roku 2005. Ale z tego porównania nie można wyciągać wniosków co do utraty wartości tego auta. Są to różne wersje, z różnymi problemami estetyki wersji B5 czy problemów wieku młodzieńczego wersji B6. Modelowanie utraty wartości jest tematem ciekawym i wrócimy do tego gdy będziemy mieli ceny aut zebranych rok wcześniej i rok później.

Wybór podzbioru danych – filtrowanie

Zbiór danych o cenach aut z marca 2012 ma ponad 220 tys wierszy. Kusząca jest możliwość wykorzystania wszystkich tych wierszy do modelowania. Kuszące rzeczy często jednak są złe.

Zgodnie ze zÅ‚otÄ… zasadÄ… ,,garbage in, garbage out” jeżeli do modelu wÅ‚ożymy Å›mieciowe dane, otrzymamy Å›mieciowe wyniki. Ponieważ bÄ™dziemy wykorzystywać dosyć wrażliwe metody estymacji, w pierwszym kroku wyczyÅ›cimy zbiór danych tak by zawieraÅ‚ jak najmniej kÅ‚opotliwych zmiennych.

Nawet jeżeli czyszczenie danych jest dosyć żmudne, mało widowiskowe i potrafi zabrać więcej czasu niż cała reszta wizualizacji i modelowania to i tak jest to jeden z ważniejszych kroków. Co z tego, że użyjemy super wyrafinowanej metody estymacji skoro w danych mamy śmieci?

W tym przypadku proces czyszczenia danych polegał na wybraniu aut, które:

• należą do segmentu C, jest to najpopularniejszy segment aut. WybieraÅ‚em z tego segmentu modele dla których znalazÅ‚em przynajmniej 100 ofert sprzedaży, co mam nadzieje wystarczy do rozsÄ…dnej estymacji. Wybrane modele to: Golf, Astra, A3, C4, Focus, Civic, i30, Cee’d, 308, Octavia.
• majÄ… mniej niż 12 lat, czyli rok produkcji to przynajmniej rok 2000,
• ma nadwozie Kombi lub Hatchback, pozostaÅ‚e wersje nadwozia sÄ… znacznie mniej popularne,
• jest zarejestrowane w Polsce,
• pochodzi z jednego z krajów: Polska, Wlochy, Czechy, Francja, Holandia, Belgia, Niemcy,
• rodzaj paliwa to diesel, benzyna lub benzyna+LPG, inne bardziej egzotyczne źródÅ‚a energii pomijam,
• pojemność silnika jest z przedziaÅ‚u 1150 – 2200 cm3,
• skrzynia biegów jest manualna,
• auto nie jest uszkodzone a cena jest cenÄ… brutto.

Po zastosowaniu tych filtrów pozostaje jedynie 13 419 aut, czyli mniej niż 10%. Wybrana grupa jest jednak bardziej jednorodna, mniej jesteÅ›my wiÄ™c narażeni na ,,niespodziewane” problemy zwiÄ…zane z brakiem zrównoważenia danych. Np. nie musimy siÄ™ przejmować artefaktem zwiÄ…zanym z tym, że najstarsze golfy majÄ… po kilkadziesiÄ…t lat, a modele i30 sÄ… co najwyżej kilkuletnie.

#
# Filtrujemy dane. Pozostawiamy tylko auta pasujące do powyższego opisu.
# W wyniku otrzymujemy trochÄ™ ponad 13 tys aut
#
load("cenyAutIII2012.Rdata")
modele <- c("Golf", "Astra", "A3", "C4", "Focus", "Civic", "i30", "Cee'd", "308", "Octavia") 
segmentC <- cenyAutIII2012[(cenyAutIII2012$Model %in% modele) & 
                 (cenyAutIII2012$Nadwozie %in% c("Kombi", "Hatchback")) & 
                 (cenyAutIII2012$Kraj.aktualnej.rejestracji == "Polska" | cenyAutIII2012$Kraj.aktualnej.rejestracji == "") & 
                 (cenyAutIII2012$Kraj.pochodzenia %in% c("Wlochy", "Czechy", "Francja", "Holandia", "Belgia", "Polska", "Niemcy", "")) & 
                 (cenyAutIII2012$Rodzaj.paliwa %in% c("benzyna+LPG", "benzyna", "olej napedowy (diesel)")) &
                 (cenyAutIII2012$Pojemnosc.skokowa > 1150 & cenyAutIII2012$Pojemnosc.skokowa < 2201) &
                 cenyAutIII2012$Rok.produkcji > 2000 & 
                 cenyAutIII2012$Skrzynia.biegow == "manualna" & 
                 cenyAutIII2012$Brutto.netto == "brutto" & (cenyAutIII2012$Liczba.drzwi %in% c("2/3", "4/5")) & cenyAutIII2012$Pojazd.uszkodzony == "",]
segmentC <- segmentC[!is.na(segmentC$Cena),]

Kodowanie zmiennych i wstępne transformacje

Zanim zaczniemy cokolwiek modelować, trzeba zastanowić się jak zakodować zmienne w modelu. Po serii prób wybrałem następujące transformacje przygotowujące dane

• zamiast roku produkcji analizowany bÄ™dzie wiek auta, czyli 2012 – rok produkcji. UÅ‚atwi to interpretacjÄ™ ceny efektu wieku auta (porównywaÅ‚em też wyniki z wiekiem traktowanym jako zmienna jakoÅ›ciowa, ale nie byÅ‚y lepsze, wiÄ™c wiek pozostaÅ‚ zmiennÄ… iloÅ›ciowÄ…)
• pojemność silnika zostaÅ‚a zamieniona na zmiennÄ… jakoÅ›ciowÄ… z poziomami co 100cm3.
• z adresu auta wyciÄ…gnÄ…Å‚em informacjÄ™ o kodzie pocztowym, odpowiednio pierwszej cyfrze kodu pocztowego, dwóch pierwszych cyfrach i wszystkich piÄ™ciu cyfrach.
• z kolumn wyposażenie dodatkowe i informacje dodatkowe wyciÄ…gnÄ…Å‚em elementy wyposażenia pojawiajÄ…ce siÄ™ w przynajmniej 100 ofertach. DodaÅ‚em kolumny kodujÄ…ce binarnie czy auto posiada: niezalezne ogrzewanie, instalacja gazowa, szyberdach, bagażnik na dach, blokada skrzyni biegow, skorzana tapicerka, ksenony, EDS, system nawigacji, hak, podgrzewane fotele, pod. przednia szyba, przyciemniane szyby, czujnik deszczu, czujnik parkowania, tempomat, kierownica wielofunkcyjna, welurowa tapicerka, ASR, garażowany, ESP, alufelgi, autoalarm, lwiatla przeciwmglowe, pierwszy wlasciciel, serwisowany w ASO, bezwypadkowy, komputer, el. lusterka, radio / CD, immobiliser, el. szyby, poduszka powietrzna, klimatyzacja, centralny zamek, ABS, wspomaganie kierownicy.
• cenÄ™ wyrażonÄ… w różnych walutach zamieniÅ‚em na cenÄ™ w PLN.
• zmiennÄ… do analizy bÄ™dzie logarytm dwójkowy ceny. Analiza ceny z użyciem modeli gaussowskich nie ma sensu, majÄ…c do wyboru transformacjÄ™ ceny albo rozważanie bardziej zÅ‚ożonych klas modeli wybraÅ‚em transformacjÄ™. Z rodziny transformacji Boxa Coxa najlepiej wypadaÅ‚a transformacja y^0.2, ale nie różniÅ‚a siÄ™ ona znaczÄ…co od logarytmu wiÄ™c wybraÅ‚em logarytm bo Å‚atwiej go interpretować.

#
# Dodajemy do zbioru danych zmienne po transformacji.
#
# Zamieniamy rok produkcji na wiek
segmentC$wiek <- 2012 - segmentC$Rok.produkcji
 
# Zamieniamy ilościową zmienną pojemność na zmienną jakościową
segmentC$pojemnosc <- factor(((segmentC$Pojemnosc.skokowa - 50) %/% 100) * 100 + 50)            
 
# konstruujemy zmiennÄ… opisujÄ…ca lokalizacjÄ™ geograficznÄ…, 
# odpowiadającą pierwszej cyfrze kodu (mniej więcej województwo), 
# dwóm pierwszym cyfrom kodu pocztowego lub wszystkim pięciu cyfrom kodu pocztowego.
segmentC$kodPocztowy5 <- factor(substr(segmentC$Adres, 1,6))
segmentC$kodPocztowy1 <- factor(substr(segmentC$Adres, 1,1))
segmentC$kodPocztowy2 <- factor(substr(segmentC$Adres, 1,2))
 
# usuwamy zbędne puste poziomy
segmentC$Model = factor(segmentC$Model)
segmentC$Nadwozie = factor(segmentC$Nadwozie)
segmentC$Rodzaj.paliwa = factor(segmentC$Rodzaj.paliwa)
 
# Konstruujemy macierz cech dodatkowych, najpierw identyfikujemy listę nazw wyposażenia dodatkowego 
# a następnie budujemy odpowiednią ramkę danych
tmp <- paste(as.character(segmentC$Wyposazenie.dodatkowe),as.character(segmentC$Informacje.dodatkowe), sep=", ")
tmps <- strsplit(tmp, split=", *")
cechy <- names(sort(table(factor(unlist(tmps))))[10:46])
ncechy <- sapply(cechy, function(x) grepl(x, tmp))
 
# Składamy nowe zmienne w ramkę danych segmentC2
segmentC2 <- data.frame(segmentC[,c("Cena.w.PLN", "Model", "Nadwozie", "Pojemnosc.skokowa", "Rodzaj.paliwa", 
  "wiek", "kodPocztowy2", "kodPocztowy1")], ncechy)

#
# Pomocnicza funkcja, przedstawiająca tabelę efektów w bardziej zwartej i czytelnej postaci. Przekształcamy efekty liczone na logarytmach tak by przedstawiały ,,efekt multiplikatywny'' czyli porównanie procentowej ceny aut.
#
opisz <- function(x, filtr="") {
 tt <- data.frame(round(2^x[grepl(filtr, rownames(x)),1,drop=F]*1000)/10, 
       signif = cut(x[grepl(filtr, rownames(x)),4,drop=F], c(-1,0.001,0.01,1), c("**","*","")))
 tt[order(tt[,1]),]
}
 
#
# Wyniki analizy wariancji
# Jak widzimy największy wpływ na cenę auta ma wiek, w drugiej kolejności model auta i typ nadwozia
#
> anova(model <- lm(log(Cena.w.PLN,2)~Model:Nadwozie+wiek, segmentC2))
Analysis of Variance Table
 
Response: log(Cena.w.PLN, 2)
                  Df Sum Sq Mean Sq  F value    Pr(>F)    
wiek               1 5873.4  5873.4 67373.01 < 2.2e-16 ***
Model:Nadwozie    18  706.9    39.3   450.47 < 2.2e-16 ***
Residuals      13399 1168.1     0.1                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
#
# Referencyjnym poziomem jest Skoda Octavia Kombi
# Wartość w kolumnie ,Estimate' odpowiada mniej więcej temu ile procent średniej ceny 
# Skody Octavii Kombi  stanowi średnia cena danego auta w określonej wersji nadwozia. 
# Czyli np. cee'd kombi jest o 15% tańszy od octavii kombi (przyjmując tę samą utratę 
# wartości z uwagi na wiek), za to Audi A3 kombi jest średnio o 35% droższe niż Octavia Kombi.
#
# Wartość przy zmiennej wiek oznacza że auta średnio tracą około 14% wartości na rok.
# 
> opisz(summary(model)$coefficients)
                                Estimate signif
ModelCivic:NadwozieKombi            69.7       
ModelFocus:NadwozieKombi            70.2     **
ModelC4:NadwozieHatchback           72.9     **
Modeli30:NadwozieHatchback          73.9     **
ModelCee'd:NadwozieHatchback        74.5     **
ModelAstra:NadwozieKombi            75.3     **
ModelFocus:NadwozieHatchback        75.4     **
Modeli30:NadwozieKombi              78.0     **
ModelAstra:NadwozieHatchback        78.5     **
Model308:NadwozieKombi              79.6     **
Model308:NadwozieHatchback          80.6     **
ModelCee'd:NadwozieKombi            86.3     **
ModelOctavia:NadwozieHatchback      93.4     **
ModelGolf:NadwozieKombi             93.6     **
ModelCivic:NadwozieHatchback        96.0     **
ModelGolf:NadwozieHatchback         98.8       
ModelA3:NadwozieHatchback          124.5     **
ModelA3:NadwozieKombi              135.4     **
wiek                                86.0     **

Opisując transformacje zmiennych napisałem że do zbioru danych dodałem 37 binarnych zmiennych opisujących wyposażenie auta. Zobaczmy, które elementy wyposażenia mają największy wpływ na cenę auta.

Aby zmniejszyć objętość poniższego przykładu przedstawiam wyniki tylko dla wybranych elementów wyposażenia. Cztery elementy wyposażenia najdroższych aut to kseony, ESP, klimatyzacja i elektryczne lusterka. Auta, które w opisie wyposażenia miały wymienione wspomaganie kierownicy lub poduszkę powietrzną są o kilka procent tańsze od aut bez tych elementów w opisie. Nie oznacza to oczywiście, że z dwóch aut podobnych to bez wspomagania kierownicy jest droższe. Oznacza to raczej, że w opisie auta jeżeli pojawia się wspomaganie kierownicy to nie pojawiają się inne elementy jeszcze bardziej zwiększające ceny auta. Należy więc być bardzo ostrożnym interpretując wyniki z modelowania cen.

Poniższy model wyjaśnia już ponad 90% wariancji.

#
# Auta z kseonowymi lampami, klimatyzacjÄ…, ESP, elektrycznymi lusterkami sÄ… 
# średnio o kilka procent droższe od aut bez tych elementów. 
# Większość elementów wyposażenia nie różnicuje istotnie ceny auta
#
> model <- lm(log(Cena.w.PLN,2)~Model*Nadwozie+wiek+.-kodPocztowy1-kodPocztowy2, segmentC2)
> opisz(summary(model)$coefficients, "TRUE")
                              Estimate signif
wspomaganie.kierownicy        91.4     **
poduszka.powietrzna           94.6     **
welurowa.tapicerka            96.2     **
...
pierwszy.wlasciciel          100.2       
serwisowany.w.ASO            100.5       
bezwypadkowy                 100.7       
przyciemniane.szyby          101.0       
czujnik.deszczu              101.3       
system.nawigacji             101.9      *
tempomat                     101.9     **
bagaznik.na.dach             102.1       
skorzana.tapicerka           102.3      *
podgrzewane.fotele           102.5     **
el..szyby                    102.5      *
alufelgi                     103.0     **
kierownica.wielofunkcyjna    103.1     **
komputer                     103.2     **
klimatyzacja                 103.3     **
el..lusterka                 103.8     **
ESP                          104.3     **
ksenony                      108.2     **

Kuszące jest dodać do modelu informację o kodzie pocztowym aby sprawdzić czy miejsce w którym oferowane jest auto istotnie wpływa na cenę ofertową.

#
# Różnice pomiędzy cenami aut są ,,istotne statystycznie''
# Nie jest to jednak duże osiągnięcie, biorąc pod uwagę 13 tysięcy ofert 
# na których przeprowadzamy testowanie.
# Wartość różnic nie jest duża w porównaniu do efektu wieku, modelu czy typu nadwozia.
#
> anova(model <- lm(log(Cena.w.PLN,2)~Model:Nadwozie+wiek+kodPocztowy1+kodPocztowy2+kodPocztowy5, segmentC2))
Analysis of Variance Table
 
Response: log(Cena.w.PLN, 2)
                  Df Sum Sq Mean Sq    F value    Pr(>F)    
wiek               1 5873.4  5873.4 72953.9503 < 2.2e-16 ***
kodPocztowy1      10    7.4     0.7     9.1992 2.489e-15 ***
kodPocztowy2      88   39.0     0.4     5.5087 < 2.2e-16 ***
kodPocztowy5    1091  240.8     0.2     2.7419 < 2.2e-16 ***
Model:Nadwozie    18  604.7    33.6   417.2709 < 2.2e-16 ***
Residuals      12210  983.0     0.1                         
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ciekawym eksperymentem może być samodzielne znalezienie innych czynników istotnie różnicujące cenę auta. Szybko jednak się okazuje, że wyniki rozjeżdżają się z intuicją. I jest to za każdym razem efekt niezrównoważonych zmiennych (porównując diesle i silniki benzynowe trudno uciec od tego że diesle mają większy przebieg). Drugi ciekawy problem to szybko rosnąca macierz modelu. Jeżeli uwzględnimy wiek jako zmienna jakościowo i będziemy analizować jego interakcje z kodem pocztowym to jedna interakcja generuje ponad 1000 kolumn do macierzy modeli.

Analiza takich danych jest znacznie ciekawsza niż analiza zależności pomiędzy szerokościami płatków irysów jest też znacznie trudniejsza. Ale też bliższa rzeczywistym problemom.

Więc życzę miłej zabawy!

Plan na kolejny wpis: przedstawić powyższe liczbowe wyniki graficznie.