Jak wyglądają ceny aut używanych w segmencie C i problemy modelowania

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.