Tag: kod w R

Ceny używanych aut po nowym roku a prosty pająk sieciowy

Trzy miesiące temu opisywałem zbiór danych o ofertach sprzedaży samochodów z serwisu otomoto.pl (zobacz tutaj). Po kwartale ponownie zebrałem informacje o ofertach sprzedaży. Porównując te dwa zbiory będzie można sprawdzić jak zmieniają się ceny różnych modeli w odpowiedzi na podwyższoną akcyzę na olej napędowy. Kolejne interesujące pytanie to globalna zmiana cen aut używanych w ,,czasach kryzysu’’. Interesujący będize trend zmian w liczbie ofert samochodów sprowadzanych z zagranicy. Efekt starzenia się auta, z dnia na dzień auto staje sie o rok starsze, ciekawe czy dla aut używanych wiek w latach ma duzy wpływ na cenę auta.

Nowy zbiór danych jest dostępny w postaci csv i Rdata pod tym adresami [csv 80Mb], [Rdata 5,5 MB].

Skrypt wczytujący dane bezpośrednio z Internetu do R znajduje się tutaj.

 

Ponieważ ostatnio pojawiały się pytania o skrypty użyte do parsowania stron interntowych i zbierania automatycznego danych, więc do tych danych dołączam wszystkie skrypty potrzebne do ich zebrania i przetworzenia. Poniżej krótko je opiszę. Wykorzystuję tutaj głównie Perla i wyrażenia regularne pozwalające na łatwe i szybkie parsowanie tekstu. Z pewnością istnieją lepsze (łatwiejsze) sposoby parsowania stron HTML, chętnie usłyszę sugestie.

W moim rozwiązaniu proces zbierania danych jest dwuetapowy. Pierwszy etap to przejrzenie hierarchicznej struktury ofert sprzedaży samochodów oraz zapisanie z każdej oferty wybranych pól. Dodatkowo należy zatroszczyć się o to by to samo auto nie zostało dwukrotnie dodane do listy aut. Ponieważ parsowanie 200 tys ofert trochę trwa, a tym czasie do serwisu dodawane są nowe oferty, więc potrzebna jest dodatkowa tablica haszująca pamiętająca, które auta zostały już przetworzone. Skrypt perlowy, który wykorzystałem w tym celu znajduje się tutaj [skrypt perl]. Plik wynikowy, z zebranymi danymi znajduje się tutaj [surowe dane po pierwszym etapie 140MB].

Ponieważ dla różnych ofert sprzedaży  przedstawiane są różne charakterystyki, niektórych opisów brakuje, niektóre parametry są w różnych formatach, dlatego też w drugim kroku wykorzystuję skrypt R obrabiający surowe dane i przekształcający je do postaci tabelarycznej, łatwiejszej w obróbce. Przy okazji cena aut w walucie przeliczana jest na PLN, liczona jest pojemność skokowa, moc i inne parametry, które w surowych danych sa zapisane w formacie tekstowym. Skrypt użyty do transformacji znajduje sie tutaj [skrypt r], a plik wynikowy znajduje się tutaj [plik csv].

Smacznego!

 

O czym się mówi w Sejmie i Senacie

Kontynuując wpis sprzed dwóch dni, dziś będzie o Sejmie i Senacie. Interesującym portalem nawiązującym do tematu smart voting jest Sejmometr (http://sejmometr.pl/). Umożliwia on obserwowanie prac Sejmu i Senatu. Niebawem dostępne będzie API pozwalające na automatyczny dostęp do zbieranych w ramach tego projektu danych. W bardzo estetyczny i łatwy w nawigacji sposób przedstawione są informacje o posłach i senatorach, ich wypowiedzi, wyniki głosowań itp. Bardzo ciekawy portal dla osób chcących być naprawdę na bieżąco.

Rzecz której mi brakuje to spojrzenie na prace sejmu ,,z lotu ptaka’’. Agregaty pozwalające na orientację co się działo przez ostatnie pół roku/dwa lata. Rozwiązania w stylu chmura tagów, z informacjami jakie tematy są najczęściej poruszane albo analiza częstości słów byłoby mile widziane. Tagi mogłyby być generowane automatycznie lub jeszcze lepiej, użytkownicy portalu mogliby otagowywać wystąpienia posłów. Poczekajmy więc na API i zobaczymy co w tym temacie można zrobić.

A w międzyczasie pokażę przykład analizy danych pochodzących ze stenogramów, z lotu bardzo wysoko latającego ptaka. Punktem wyjścia są dane z Korpusu Języka Polskiego (http://korpus.pl/), projektu rozwijanego przez IPIPAN. Wiele tekstów zostało poddanych analizie w ramach tego korpusu, między innymi stenogramy z posiedzeń Sejmu i Senatu. Dla każdego z posiedzeń, dla każdego (ok., dla większości) wypowiedzianego słowa przypisano odpowiedni fleksem, czyli jedną z ponad dwudziestu klas gramatycznych. Więcej informacji o fleksemach znaleźć można np. tutaj. Mamy więc zbiór danych dla 65 milionów słów wypowiedzianych przez 4 kadencje Sejmu i 4 kadencje Senatu z informacją do której klasy fleksyjnej (których klas) należy to słowo. Możemy teraz z takiego zbioru danych policzyć tablice rozdzielczą (krzyżową, kontyngencji, zwał jak zwał) z informacją w której kadencji Sejmu i Senatu padło ile słów z odpowiedniej klasy fleksyjnej (opisującej formę i znaczenie słowa). Tę tablicę liczb (27 x 8) umieściłem w pakiecie PBImisc programu R w zbiorze danych SejmSenat. Jak znaleźć wzorce w takiej tablicy 217 liczb? Użyjemy do tego celu analizy odpowiedniości / analizy korespondencji. Analiza ta pozwala na określenie, które wiersz (fleksemy) i kolumny (kadencje Sejmu i Senatu) mają podobne profile używalności słów a również które wiersze / kolumny współwystępują częściej niż wskazywałaby na to częstość występowania w języku. Trochę więcej o stronie technicznej później, a na razie zobaczmy wyniki analizy korespondencji na zbiorze danych o używalności klas fleksyjnych w stenogramach  Sejmu i Senatu.

[Rys. 1. Czerwone strzałki odpowiadają profilom stenogramów z posiedzeń Sejmu Senatu, niebieskie punkty odpowiadają profilom używalności fleksemów. Im bliższe zwroty tym większa zależność pomiędzy profilami. W analizie usunięto fleksem interp, ponieważ bardzo odstawał od pozostałych. Wersja png znajduje się tutaj. Warto ten wykres powiększyć by zobaczyć gdzie są jakie fleksemy, na szczęście jest to grafika wektorowa.]

Osie pozioma i pionowa odpowiadają dwóm automatycznie znalezionym komponentom. Tak się jednak składa, że te komponenty mają naturalną interpretację, którą łatwo odczytać z wykresu. Pierwszy komponent (oś pozioma) odpowiada za zmiany w używalności fleksemów pomiędzy Sejmem a Senatem. Im wyższa wartość pierwszej składowej, tym profil używalności bardziej charakterystyczny dla stenogramów z Senatu. Drugi komponent odpowiada za zmianę w używalności fleksemów z czasem, im wyższa wartość drugiej składowej tym profil bardziej charakterystyczny dla starszych  posiedzeń.

Pięknie. Skoro osie mają taką naturalną interpretację, to zobaczmy jakie fleksemy są częściej używane w Sejmie, a które w Senacie, które były częściej używane kiedyś a które obecnie.

Na osi poziomej dwa interesujące fleksemy to np. num i depr. Num to skrót od ,liczebnik główny’ a depr to skrót od rzeczownik deprecjatywny (najczęściej używany do oceny negatywnej).  To co można więc z wykresu łatwo odczytać (i sprawdzić ręcznie w tablicy kontyngencji) to, że w stenogramach z Senatu częściej występują liczebniki niż w stenogramach z Sejmu, widać więcej mówi się o liczbach, konkretach. W stenogramach z Sejmu częściej występują rzeczowniki w znaczeniu deprecjatywnym, widać atmosfera jest gorętsza.

Na osi pionowej interesujące fleksemy to np. winien i będzie. Skrót ‘winien’ odpowiada słowom typu ,winien’, ‘powinien’ itp., skrót ‘bedzie’ odpowiada przyszłej formie czasownika być. Kiedyś jak widać częściej w Sejmie i Senacie mówiło się o tym jak być powinno, tendencja jest tak aby częściej mówić o tym jak będzie.

Podsumowując. Można z danych o stenogramach szukać trendów widocznych w większej skali czasowej. Potrzebujemy tylko dostępu do przetworzonych stenogramów i pomysłu na to czego w tych stenogramach szukać.

W programie R jest kilka pakietów do analizy korespondencji, np., pakiet ca i anacor, można o nich przeczytać np. tutaj (pakiet anacor, pakiet ca). W tych artykułach przedstawione są zarówno  podstawy matematyczne jak i przykłady zastosowań. Technicznie, podobnie jak dla PCA, bazuje się na dekompozycji SVD, ale oczywiście ważne jest co chcemy dekomponować.

Kod R wraz z dokładniejszymi wynikami analizy korespondencji znajduje się tutaj.

 

 

Czym jeżdżą nasi wybrańcy, czyli rzut okiem na oświadczenia majątkowe połów i senatorów

Dzisiejszy wpis ma zawierać będzie trzy elementy, obywatelski, humorystyczny i techniczny.

Element obywatelski:

Jakiś czas temu przedstawiłem kilka wizualizacji danych dotyczących głosowań posłów z poprzedniej kadencji (dot glosowań,dot głosowań 2). Podsumowania przebiegu obrad sejmu to ważna sprawa, podsumowanie pracy posłów również. Zamiast pokazywać kto ile pieniędzy wydał na diety i paliwo może lepiej byłoby zbadać czy dany poseł na to paliwo i diety zapracował. Mało kto ma tyle czasu i zaparcia by czytać stenogramy z obrad, nawet podsumowania każdego z głosowań ilu za ilu przeciw to również nudnawa lektura. Z drugiej strony warto wiedzieć nad czym pracują nasi wybrańcy i kto jakie poglądy reprezentuje. Nie tylko od święta gdy podgrzewa się opinię publiczną debatą na temat definicji życia. Są na szczęście portale i stowarzyszenia działające w temacie tzw. smart voting. Jak by tego nie tłumaczyć tego terminu, ich celem jest przybliżenie działalności sejmu i senatu zwykłym ludziom po to by Ci lepiej wybierali w kolejnych wyborach lub po prostu tworzyli presję. Jedna z takich organizacji to stowarzyszenie art.61 (http://art61.pl/) prowadzące obywatelski portal informacyjny ,,Mam prawo wiedzieć’’ (http://mamprawowiedziec.pl). Jakiś czas temu otrzymałem od nich dwa zbiory danych, jeden dotyczył poglądów kandydatów a drugi oświadczeń majątkowych złożonych przez posłów i senatorów z poprzedniej kadencji. Ok, zastanówmy się co ciekawego można z takich danych pokazać.

Element humorystyczny:

Oświadczenia majątkowe nie są zbyt pasjonującą lekturą same w sobie. Można oczywiście skupić się na pokazaniu najbogatszych lub najbiedniejszych posłów ale nie ma w takim zestawieniu zabawy. Więc dziś zrobimy coś mniej poważnego ale bardziej zabawnego.

Zobaczymy jakie samochody posłowie i senatorowie zadeklarowali w ostatnim złożonym oświadczeniu majątkowym. Ponieważ oświadczenia wypisuje się ręcznie, portal ,,Mam prawo wiedzieć’’ poświęcił sporo czasu by te dane zdigitalizować do postaci nadającej się do dalszej obróbki. Nie sposób jednak uniknąć błędów i niespójności, poniższe wizualizacje dotyczą więc tylko dych oświadczeń w których udało się automatycznie zidentyfikować markę i model auta.

Na poniższej wizualizacji zaznaczono ile sumarycznie samochodów określonego modelu i marki zadeklarowali w oświadczeniach posłowie i senatorowie.

[Rys 1. Liczba aut określonych marek i modeli wskazanych w ostatnim oświadczeniu majątkowym przez posłów z 6 kadencji. Każdy prostokąt oznacza jedno auto. Szerokość i wysokość prostokąta jest proporcjonalna do rozmiarów auta. Jeżeli nie wyświetla się wersje svg zobacz tutaj wersję png.]

[Rys 2. Liczba aut określonych marek i modeli wskazanych w ostatnim oświadczeniu majątkowym przez senatorów z 7 kadencji. Każdy prostokąt oznacza jedno auto. Szerokość i wysokość prostokąta jest proporcjonalna do rozmiarów auta. Jeżeli nie wyświetla się wersje svg zobacz tutaj wersję png.]

 

Element techniczny:

Kod źródłowy użyty do wygenerowania tych rysunków znajduje się tutaj.

Dane do analizy zostały przekazane jako pliki Excelowe, po jednym na jednego posła/senatora (łącznie ponad 570 plików Excelowych). Pliki takie są czytelne dla człowieka ale do automatycznego przetwarzania nadają się tak sobie (różna liczba wierszy i kolumn, posklejane komórki, kolory, różne standardy oznaczania samochodów, różne informacje podawane nt. samochodu). Zazwyczaj do odczytania danych z plików w formacie xls używa się funkcji read.xls(). Ta funkcja jest dostępna w pakiecie xlsReadWrite, który wykorzystuje zewnętrzną bibliotekę dll napisaną w Delphi do odczytywania plików xls, oraz w pakiecie gdata, który wykorzystuje moduły Perla do konwersji formatu xls do formatu csv a później odczytuje dane przez funkcję read.csv. Niestety w przypadku otrzymanych plików żadna z tych funkcji nie dawała zadowalających wyników. W większości przypadków funkcje te odczytywały wyłącznie wartości NA. Najlepsze rozwiązanie udało się otrzymać z użyciem funkcji sqlFetch() i odbcConnectExcel() z pakietu RODBC. Do tego kilka różnych warunków i wyszukiwanie wzorców z użyciem funkcji grep() wystarczyło by zidentyfikować pole w którym w oświadczeniach majątkowych wpisane są posiadane samochody.

Kolejnym problemem jest niesamowita liczba literówek dotyczących wpisywanych nazw modeli aut. Np. w oświadczeniach posłów model Oktavia występował w 6 różnych formach łącząc różne możliwe wielkości liter oraz nazwę modelu zamiennie stosując litery c/k i w/v. Podobną liczbę wariantów uzyskała Corolla. Poza tym nazwy modeli i marek aut otoczone były różnymi informacjami, czasem o przebiegu, czasem o roku produkcji, czasem o wartości. Aby z tego morza tekstu wyłowić nazwy modeli wykorzystałem dane z allegro, które zebrałem jakiś czas temu (link). Można teraz szukać w tekście oświadczeń majątkowych słów, które przypominają nazwy modeli i marek znane z danych z allegro. Dzięki temu można odsiać większość śmieciowych słów a aby poprawnie odczytać większość modeli wystarczy juz tylko około 100 dodatkowych ręcznie wprowadzonych reguł.

Ostatnia rzecz dotyczy wielkości aut. Otóż stwierdziłem, że fajnie byłoby zaznaczyć jak duże są pokazywane auto. Nie każdy musi wiedzieć czy VW Polo jest większy czy mniejszy od Passata. Poza tym o ile dotąd czytelnik mógł zadawać sobie pytanie czy ta automatyzacja ma sens, czy nie lepiej byłoby przeczytać dane z oświadczeń i z użyciem ,,ludzkiego przetwornika’’ zbudować nową tabelkę z informacją o tym kto jakie auto ma, tak w sposób ręczny szukanie informacji o parametrach technicznych aut to prawdziwie syzyfowa praca. Więc napisałem mały programik parsujący strony www.motofakty.pl i odczytujący parametry techniczne różnych wersji i modeli aut. Na przedstawionej wizualizacji szary prostokąt oznacza wymiary 2m szerokości i 5m długości, a kolorowy prostokąt oznacza średnie wymiary danego modelu (parametry te różnią się dla rocznika i typu nadwozia).

Udało się uzyskać pełną automatyzację, dzięki której było łatwo wygenerować podobny wykres dla senatorów i będzie łatwo powtórzyć razem z nowymi oświadczeniami majątkowymi.

I żyli długo i szczęśliwie…

Kilka dni temu popełniłem wpis opisujący prawdopodobieństwo dożycia wieku emerytalnego. Przy okazji pojawiła się dyskusja nt. tego czy oczekiwana długość życia będzie się w Polsce wydłużała czy nie (zobacz ten wpis).
Łatwo być adwokatem optymistycznych jak i pesymistycznych scenariuszy, gdy opiera się wyłącznie na przypuszczeniach. Temat nie dawał mi spokoju, więc znalazłem dane na podstawie których zobaczmy jak wygląda  oczekiwana długość życia w Polsce i w innych krajach. Może zobaczymy czy istnieje i czy dotarliśmy do maksymalnej średniej życia a jeżeli tak to gdzie ona jest.

Dane z których korzystam pochodzą z bazy danych http://www.mortality.org/. W tej bazie danych dostępne są tablice życia i inne pochodne miary zebrane dla 37 krajów. Tablice życia są dostępne dla pewnej liczby ostatnich lat, dla różnych krajów długość tej historii jest różna. Najdłuższa historia jest dla Szwecji i sięga ponad 200 lat, dla Polski mniej więcej 50 lat.

Mając te dane zobaczmy co się dzieje z oczekiwaną długością życia w różnych krajach w ostatnich 50 latach (zbiór danych life expectancy). Zobaczmy wykres poniżej. W tym zestawieniu Polska charakteryzuje się najniższą oczekiwaną długością życia. Oczywiście są kraje w których żyje się krócej, ale nie znalazły się w tym zestawieniu. Ma to tę zaletę, że przed nami prawdopodobnie wydłużająca się średnia długość życia, nie widać na razie w tych prognozach sufitu. W większości krajów współczynnik wzrostu wynosi około 3 lat średniej życia na dekadę.

[Wersja png rysunku. Rys 1. Oczekiwane średnie życia mężczyzn dla wybranych 11 krajów w ostatnich 50 latach. ]
Dla niektórych krajów mamy dane ze znacznie większej liczby lat. Zobaczmy dla mniejszej grupy krajów jak wygląda zmiana oczekiwanej długości życia w szerszym przedziale czasu. Dynamika zmian długości życia jest różna w różnych krajach, w Szwajcarii czy Szwecji widać mniej więcej stały wzrost oczekiwanej długości życia.


[Wersja png rysunku. Rys 2. Oczekiwane średnie życia mężczyzn dla wybranych 5 krajów w ostatnich 140 latach. Kolory jak na Rys 1.]

Nie jestem zwolennikiem średnich i pracowania na wartościach oczekiwanych. Najchętniej zobaczyłbym do jakiego wieku dożywa 50% mężczyzn. To znacznie ciekawszy współczynnik, który można łatwiej zinterpretować. Problem jest tylko taki, że ponieważ taki ,,połowiczny rozpad’’ dla mężczyzn będzie wynosił około 70 lat więc by go policzyć dokładnie potrzebujemy danych z tablic życia wstecz o ponad 70 lat. Dla Polski takich danych nie mam, ale mam dla Szwajcarii. Na poniższym rysunku porównuję trzy współczynniki, mogące opisywać długość życia:

  1. Wiek jakiego dożywa 50% mężczyzn urodzonych w roku X
  2. Średni czas życia chłopców urodzonych w roku X

Co ciekawe jednak, współczynnik 1 jest wyższy niż pozostałe dwa. Dlaczego? Argument, który przychodzi do głowy jest taki, że długość życia to zmienna lewostronnie skośna. Tz. jest okres w którym najwięcej osób umiera (okolice 70 roku), ale znacznie więcej jest osób które umarło 65 lat wcześniej, niż 65 lat później. Dla skośnych rozkładów średnia nie pokrywa się z medianą, a w tym przypadku co jest ciekawe wiek którego dożyje 50% mężczyzn jest wyższy niż oczekiwany czas życia w chwili narodzin. Oczywiście w międzyczasie wydarzyły się dwie wojny światowe, które jakoś wpłynęły na skośność rozkładu czasu życia. Dla Szwajcarii jednak w mniejszym stopniu niż dla sąsiednich krajów. Do tego dochodzi wysoka śmiertelność najmłodszych. W prezentowanym okresie czasu w pierwszych dwóch latach życia umierało około 15% chłopców.

Pointa.: Nie dość że średnia życia rośnie, to ponad połowa mężczyzn będzie żyła dłużej niż ich oczekiwana długość życia.

Wreszcie coś optymistycznego na święta.


[Wersja png rysunku. Rys 3. Oczekiwany lub połówkowy czas życia mężczyzn w Szwajcarii. Na osi OX przedstawiłem dane dla których połówkowy czas życia czy średni czas życia mogłem policzyć na podstawie danych historycznych a nie szacować. Dlatego oś OX kończy się w okolicach roku 1920.]

 

Wykresy radarowe a kobieta manager

Kilka dni temu przy okazji wpisu http://smarterpoland.pl/index.php/2011/12/kobieta-menedzer-a-szansa-na-sukces/ w komentarzu od @Analfabeta pojawiła się sugestia by przedstawić te dane na wykresie radarowym.

Zobaczmy więc jak wyglądać będą te dane na wykresie radarowym, poniżej informacja o tym które elementy są ważne do osiągnięcia sukcesu zaprezentowana w kolorach na wykresie radarowym.

Kod w programie R za pomocą którego można taki wykres wygenerować znaleźć można tutaj.

Który element został uznany za ważny a procent wskazań wykonanych przez kobiety.

Elementy uznane za ważne do osiągnięcia sukcesu w opinii mężczyzn i kobiet.

Jak wyżej, ale nie łączymy punktów linią, poprawniej ale trudniej się czyta.

Tutaj zamiast funkcji lines używamy polygon dzięki czemu mamy wypełnione ,,radary”.

Kobieta menedżer a szansa na sukces

Andrzej P. podesłał mi artykuł zatytułowany ,,Kobieta menedżer ma mniejsze szanse na awans” (artykuł tutaj). Artykuł ten jest wyjątkowo ciekawym przykładem jak nie pokazywać danych. W artykule autorka stara się nas przekonać, że kobiety menedżerki (to słowo jest już nawet w SJP) mają mniejsze szanse na awans. Przekonać ma nas o tym niezbicie pierwszy wykres.

Już nawet nie czepiam się wykresu kołowego, ani tego że jest on 3D, ani że odpowiedź która ma się najbardziej rzucać w oczy jest na czerwono. Najbardziej zdziwiony jestem, że pytanie które zostało zadane to ,,czy szanse na awans są TAKIE SAME?”. To już autorka zadecydowała że nierówność musi oznaczać faworyzowanie mężczyzn.

 

Ciekawy jest też drugi wykres prezentowany w tym artykule.

Teoretycznie z takich danych można by się dowiedzieć, które elementy są częściej wskazywane przez mężczyzn a które przez kobiety. Teoretycznie, ponieważ sposób prezentacji to uniemożliwia, trudno porównywać iloczyny długości słupków pomiędzy sobą.

Również teoretycznie można by odczytać z takich danych które elementy są uznawane za najważniejsze w sumie. Ale ponownie tylko teoretycznie, ponieważ pochyłość słupków utrudnia określenie który słupek jest dłuższy. A liczby odpowiedzi nie są podane w sumie, więc by dowiedzieć się ile osób wybrało daną odpowiedź trzeba szybko dodawać trzycyfrowe liczby.

 

Postarajmy się jednak być konstruktywni w tej krytyce. Czy można inaczej przedstawić te dane? Kod w programie R użyty do wygenerowania poniższego wykresu znajduje się tutaj.

I ten sam obrazek obrócony o 45 stopni.

Używając wykresu punktowego/rozrzutu przedstawiliśmy te same liczby, ale tym razem odczytując położenie punktów możemy porównać elementy decydujące o awansie pomiędzy sobą. Im wyżej jest kropka (dotyczy drugiego wykresu) tym częściej ten element jest wskazywany przez mężczyzn, im niżej tym częściej przez kobiety. Im bardziej na prawo jest kropka tym więcej osób w sumie uznało dany element za istotny.

 

Sugerując się komentarzami dodałem kolory. Wrzosowy i piaskowy kolor oznaczają obszary na którym jedna płeć wybiera określone elementy o ponad 20% częściej niż druga płeć. Mam nadzieję, że dzięki temu widać że niektóre elementy są preferowane przez jedną z płci.

Polskie ogonki a iconv()

Napisał do mnie maila Krzysztof T. z informacją, że strona kodowa windows-1250, którą zakodowałem polskie znaczki w zbiorze danych Diagnoza Społeczna źle wygląda pod Linuxami.

Zmieniłem więc pliki z danymi usuwając znaki diaktrytyczne. Można też było zmienić kodowanie na UTF-8, ale usunięcie ogonków gwarantuje zgodność z każdym systemem operacyjnym.

Do zmiany kodowania w programie R można użyć funkcji iconv(), która wykorzystuje specyficzne dla systemu narzędzia do konwersji. Listę obsługiwanych stron kodowych wyświetla funkcja iconvlist().

Poniżej przykładowy kod R który usuwa znaki diaktrytyczne ze zbioru danych diagnozaOsoby2011. Podanie argumentu to=”UTF-8″ spowodowałoby konwersje do formatu UTF-8. W systemie Windows od wersji R 2.11 aby usunąć ogonki należy podać argument to=”ASCII//TRANSLIT”, pod innymi systemami wystarczy to=”ASCII”.

 

# konwertujemy nazwy kolumn
colnames(diagnozaOsoby2011) <- iconv(colnames(diagnozaOsoby2011), from="windows-1250", to="ASCII//TRANSLIT")
 
# konwertujemy nazwy poziomów w zmiennych jakościowych
for (i in 1:ncol(diagnozaOsoby2011)) 
   if ("factor" %in% class(diagnozaOsoby2011[,i])) 
      levels(diagnozaOsoby2011[,i]) <- iconv(levels(diagnozaOsoby2011[,i]), from="windows-1250", to="ASCII//TRANSLIT")

Co jest ważne w pracy?

Dwa  dni temu pokazywaliśmy przykład analizy gradacyjnej w badaniu co jest ważne w życiu. Dziś zobaczymy co dla ankietowanych jest ważne w pracy. W latach 2007 i 2011 zadano respondentom pytanie o to co jest ważne w pracy. Podobnie jak w przypadku wartości ważnej w życiu, można było wybrać maksymalnie trzy cechy dorej pracy (z listy: Brak napięć i stresów, Duza samodzielnosc, Możliwość rozwoju osobistego, Praca zgodna z umiejetnosciami, Możliwość szybkiego awansowania, Stabilnosc zatrudnienia, Dogodne godziny pracy, Możliwość wykonywania pracy w domu, Dlugi urlop, Zajecie powazane przez ludzi, Odpowiednia płaca, Inne czynniki).

Używając tych samych technik co ostatnio, sprawdzimy czy oczekiwania w stosunku do pracy sie zmienily.

 

Po prawej stronie przedstawiono dla każdej cechy dotyczącej pracy informacje jaka frakcja osób uznała tę cechę za ważną. Po lewej stronie mamy wynik jednowymiarowej analizy gradacyjnej.

Zauważmy na początek że odległość tej krzywej od przekątnej, jest dużo większa niz w przypadku pytan o to co ważne w życiu. Wydaje sie to zgodne z intuicja ze pogląd dotyczący wartości waznych w zyciu zmienia sie wolniej niz dotyczacy wartosci waznych w pracy.

Największe zmiany dotyczyły wzrostu liczby osob uwazajacych ze wazna jest stabilnosc zatrudnienia (z 11.8% do 19% a więc zmiana o ponad 60%), duża samodzielnośc w pracy, brak napiec i stresow. Mniej osób za najważniejsze wymienia odpowiednia place czy prace zgodna z umiejętnościami. Mam nadzieje ze jest to zwiazane z tym ze podstawowe potrzeby zwiazane z wystarczająca placa i zatrudnieniem w odpowiednim miejscu zostaly zaspokojone i teraz osoby mogą sie skupic na wyzszych potrzebach. Moze to tez byc związane z rosnacym wiekiem respondentów, sa o 4 lata starsi moga juz cenic inne rzeczy.

Warto zrobic taka analize w podziale na grupy wiekowe, moze wiec wrocimy do tego tematu nastepnym razem.

 

1
2
3
4
5
6
7
8
9
10
11
12

# zbieramy dane kto co uwazal za wazne w pracy w roku 2007 i 2011
zb1 <- colSums(diagnozaOsoby2011[,paste("dp106_",1:12,sep="")]=="TAK zaznaczone",na.rm=T)
zb2 <- colSums(diagnozaOsoby2011[,paste("dp106_",1:12,sep="")]=="NIE zaznaczone",na.rm=T)
zb3 <- colSums(diagnozaOsoby2011[,paste("fp113",c("_1","_2","_3","_4","_5","_6","_7","_8","_9","10","11","12"),sep="")]=="TAK",na.rm=T)
zb4 <- colSums(diagnozaOsoby2011[,paste("fp113",c("_1","_2","_3","_4","_5","_6","_7","_8","_9","10","11","12"),sep="")]=="NIE",na.rm=T)
etykiety <- c("Brak napiec i stresow", "Duza samodzielnosc", "Mozliwosc rozwoju osobistego", "Praca zgodna z umiejetnosciami", "Mozliwosc szybkiego awansowania", "Stabilnosc zatrudnienia", "Dogodne godziny pracy", "Mozliwosc wykonywania pracy w domu", "Dlugi urlop", "Zajecie powazane przez ludzi", "Odpowiednia placa", "Inne czynniki")
dane <- data.frame(TAK2007 = zb1, NIE2007 = zb2, TAK2011=zb3, NIE2011=zb4)
rownames(dane) <- etykiety
zm1 <- dane[,1,drop=F]/dane[,2]
zm2 <- dane[,3,drop=F]/dane[,4]
 
plotGradeStat(zm1, zm2, osX="rok 2007", osY="rok 2011", skala=c(0.002,0.6), cutoff=0.011)

Co jest w życiu ważne?

Ostatnio moi magistranci na mini-seminarium prezentowali jednowymiarową analizę gradacyjną. Służyć może ona między innymi do porównania czy pomiędzy dwoma wektorami obserwacji zmieniła się struktura odpowiedzi. Wygląda to na ciekawą metodę, więc warto ją zaimplementować w R i zobaczyć jak dziala.

Kilka dni temu pisaliśmy o zbiorze Diagnoza Społeczna (http://smarterpoland.pl/index.php/2011/10/diagnoza-spoleczna-2011/), już dołączony do repozytorium. Wykorzystamy go na potrzeby badania analizy gradacyjnej.

W latach 2005 i 2009 w Diagnozie Społecznej ankieterzy pytali respondentów o wskazanie wartości ważnych w ich życiu (zmienne cp2.1-cp2.14 i ep2.1-ep2.14) . Badany mógł wybrać maksymalnie trzy odpowiedzi ze zbioru 14 możliwych (PIENIADZE, DZIECI, UDANE MALZENSTWO, PRACA, PRZYJACIELE, OPATRZNOSC, BOG, POGODA DUCHA, OPTYMIZM, UCZCIWOŚĆ, ŻYCZLIWOŚĆ I SZACUNEK OTOCZENIA, WOLNOSC, SWOBODA, ZDROWIE, WYKSZTALCENIE, SILNY CHARAKTER, INNE). Wykorzystamy analizę gradacyjną by sprawdzić czy zmieniła się struktura wartości w badanej grupie respondentów w przeciągu czterech lat.

Zaczniemy od analizy dwóch czternastoelementowych wektorów. Każdy wektor określi jaka frakcja osób uznała daną wartość za ważną w ich życiu. Porównamy oba wektory, by sprawdzić które wartości zyskały, a które straciły na znaczeniu pomiędzy rokiem 2009 a 2005.

 

 

Kod generujący powyższy rysunek znajduje się poniżej. Po lewej prezentowane są wyniki analizy gradacyjnej, po prawej zwykły wykres rozrzutu. Oba wykresy prezentują te same dane.

Zacznijmy od prawego wykresu. Frakcje osób uznających daną wartośc za ważną unormowano tak, by po zsumowaniu wszystkich wartości otrzymać 1. Osobno dla roku 2005 osobno dla 2009. Każdy punkt opisuje jedną wartość. Współrzędne punktu odpowiadają unormowanej frakcji osób uznających tą wartość za ważną w roku 2005 i 2009. Dorysowano przekątną, dzięki temu punkty pod przekątną odpowiadają wartościom których znaczenie spadło do roku 2009, punkty nad odpowiadają wartosciom których znaczenie wzrosło.

Po lewej stronie przedstawiono te frakcje w sposób skumulowany. Kolejność odpowiada procentowej zmianie ważności w stosunku do roku 2009. Na początku wykresu, przy punkcie 0,0 znajdują się wartości, które zyskały na znaczeniu. Pod koniec wartości, ktore stracily na znaczeniu. Długość kroku odpowiada frakcji osob uznających daną wartość za ważną. Odległość wyrysowanej łamanej od przekątnej obrazuje jak bardzo zmieniła się struktura wartości. W tym przypadku łamana jest blisko przekątnej, więc ludzie nie zmienili istotnie swojego systemu wartości. Dzieci i zdrowie zyskały na ważności. Pieniądze i praca straciły, choć w obu przypadkach nie są to duże zmiany.

 

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
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

plotGradeStat <- function(zmienna1, zmienna2, uporzadkujMalejaco = TRUE, osX = "", osY = "", skala=c(0.005,0.5), cex.text=0.8, cutoff = 0.01) {
# normalizacja obu cech
  zm1r   <- zmienna1/sum(zmienna1)
  zm2r   <- zmienna2/sum(zmienna2)
  iloraz <- zm1r/zm2r
# jezeli zachodzi taka potrzeba to zmienna sa porzadkowane w kolejn  
  if (uporzadkujMalejaco) {
    zm1r   <- zm1r[order(iloraz, decreasing=FALSE), 1, drop=FALSE]
    zm2r   <- zm2r[order(iloraz, decreasing=FALSE), 1, drop=FALSE]
    iloraz <- zm1r/zm2r
  }
# dwa wykresy w poziomie  
  par(mfrow=c(1,2))
  par(xpd=F)
# pierwszy wykres to analiza gradacyjna
# dla jednowymiarowych danych  
  plot(c(0,cumsum(zm1r[,1])),c(0,cumsum(zm2r[,1])),type="b",pch=19,xlab=osX,ylab=osY)
  abline(0,1,col="grey")
  par(xpd=NA)
# korekta na zachodzace etykiety
  odleglosci <- sqrt(diff(c(0,cumsum(zm1r[,1])))^2+diff(c(0,cumsum(zm2r[,1])))^2)
  korekta    <- numeric(length(odleglosci))
  for (i in seq_along(korekta)) {
      if (odleglosci[i] < cutoff) 
          korekta[i] <- cutoff + korekta[i-1]
  }
  text(cumsum(zm1r[,1])+korekta+2*cutoff,cumsum(zm2r[,1])+korekta-2*cutoff,rownames(zm1r), srt=-45, adj=c(0,0),cex=cex.text)
  text(cumsum(zm1r[,1])+korekta-2*cutoff,cumsum(zm2r[,1])+korekta+2*cutoff,paste(round((1/iloraz[,1]-1)*1000)/10," %",sep=""), srt=-45, adj=c(1,1),cex=cex.text)
  par(xpd=F)
# drugi wykres to klasyczny wykres rozrzutu
  plot(1,type="n",log="xy",xlim=skala,ylim=skala, las=1, cex.axis=0.8, xlab=osX, ylab=osY)
  abline(0,1,col="grey")
  abline(h=c(0.0001*c(1,2,5),0.001*c(1,2,5),0.01*c(1,2,5),0.1*c(1,2,5)),col="grey95")
  abline(v=c(0.0001*c(1,2,5),0.001*c(1,2,5),0.01*c(1,2,5),0.1*c(1,2,5)),col="grey95")
  points(zm1r[,1],zm2r[,1],pch=19)
  par(xpd=NA)
  text(zm1r[,1],zm2r[,1],rownames(zm1r), srt=-45, adj=c(-0.1,-0.1),cex=cex.text)
  par(xpd=F)
}
 
# zbieramy dane kto co uwazal za istotne w zyciu w roku 2005 i 2009
zb1 <- colSums(diagnozaOsoby2011[,paste("cp2_",1:14,sep="")]=="TAK",na.rm=T)
zb2 <- colSums(diagnozaOsoby2011[,paste("cp2_",1:14,sep="")]=="NIE",na.rm=T)
zb3 <- colSums(diagnozaOsoby2011[,paste("ep2_",1:14,sep="")]=="TAK",na.rm=T)
zb4 <- colSums(diagnozaOsoby2011[,paste("ep2_",1:14,sep="")]=="NIE",na.rm=T)
# etykiety, co jest wazne w zyciu
etykiety <- c("PIENIADZE", "DZIECI", "UDANE MALZENSTWO", "PRACA", "PRZYJACIELE", "OPATRZNOSC, BOG", "POGODA DUCHA, OPTYMIZM", "UCZCIWOSC", "ZYCZLIWOSC I SZACUNEK OTOCZENIA", "WOLNOSC, SWOBODA", "ZDROWIE", "WYKSZTALCENIE", "SILNY CHARAKTER", "INNE")
# tabela opisujaca ktora wartosc ile osob zaznaczylo lub nie w wymienionych powyzej latach
dane <- data.frame(TAK2005 = zb1, NIE2005 = zb2, TAK2009=zb3, NIE2009=zb4)
rownames(dane) <- etykiety
zm1 <- dane[,1,drop=F]/dane[,2]
zm2 <- dane[,3,drop=F]/dane[,4]
 
plotGradeStat(zm1, zm2, osX="rok 2005", osY="rok 2009", skala=c(0.001,0.5),cutoff=0.01)

Parsowanie stron HTML, meta-analiza, rak jelita i oczywiście obrazek

Dostałem dzisiaj pytanie od Macieja B. o kod użyty do wyciągania danych z portalu otomoto.pl.
Jak będę miał chwilę to ten kod wygładzę i opiszę na blogu, ale przy okazji dziś wpadłem na ciekawą funkcję służącą do parsowania danych, więc się nią podzielę.

Chodzi o funkcję readHTMLTable() z pakietu XML. Pozwala ona na wyciągnięcie danych z pliku HTML i wczytanie ich automatycznie do R.
Cool!
Jako przykład wykorzystamy zbiór danych o zachorowalności na nowotwór jelit w Wielkiej Brytanii, więcej o tym zbiorze danych i jego analizie przeczytać można na stronie http://blog.ouseful.info/2011/10/31/power-tools-for-aspiring-data-journalists-r/.

Poniższy fragment kodu wczytuje dane bezpośrednio ze strony HTML, dodaje nazwy kolumn i zmienia typy na liczbowe.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

library(XML)
#
# funkcja readHTMLTable wyciaga tabele z danymi z pliku HTML
cancerdata <- readHTMLTable('http://www.guardian.co.uk/commentisfree/2011/oct/28/bad-science-diy-data-analysis', which=1)
colnames(cancerdata) <- c('Area','Rate','Population','Number')
for (i in 2:4) cancerdata[,i] <- as.numeric(as.character(cancerdata[,i]))
 
> # zobaczmy podsumwoanie wczytanych danych
> summary(cancerdata)
            Area          Rate         Population          Number      
 Aberdeen     :  1   Min.   : 9.16   Min.   :  31332   Min.   :  6.00  
 Aberdeenshire:  1   1st Qu.:15.22   1st Qu.: 140110   1st Qu.: 24.00  
 Adur         :  1   Median :17.44   Median : 189202   Median : 32.00  
 Allerdale    :  1   Mean   :17.64   Mean   : 224741   Mean   : 39.61  
 Amber Valley :  1   3rd Qu.:19.89   3rd Qu.: 267794   3rd Qu.: 45.00  
 Angus        :  1   Max.   :31.09   Max.   :1268959   Max.   :251.00  
 (Other)      :373

Skoro już ten zbiór danych wczytaliśmy to może jeszcze słowo komentarza skąd meta-analiza w nazwie tego wpisu. Zacznijmy od przedstawienia częstość zachorowań na nowotwór jelit na 100 tys mieszkańców a liczbę osób zamieszkałych na obszarze w którym ta częstość jest liczona.

1
2
3
4
5
6
7
8

library(lattice)
#
# wykres czestosci zachorowan od rozmiaru populacji
ktoreZNazwy <- with(cancerdata, Population^0.25 * Rate > 560 | Rate >= 27)
xyplot(Rate~Population, data = cancerdata, pch=19, col="red", panel=function(...) {
	panel.xyplot(...)
	with(cancerdata, ltext(Population[ktoreZNazwy], Rate[ktoreZNazwy], Area[ktoreZNazwy], cex=0.8, adj=c(0,-0.5)))
})

Dla małych miejscowości ocena częstości zachorowań obarczona jest większą przypadkowością, jeżeli mamy miasto o 100 mieszkańcach i jeden zachoruje to unormowana częstość skacze do 1000 na 100 tys. nawet jeżeli nie jest to obszar szczególnie narażony na podwyższone ryzyko. Dla zaludnionych obszarów takie losowe fluktuacje mają mały wpływ. Zmierza to w kierunku meta-analizy w której na podstawie pomiarów z wielu obszarów ocenilibyśmy oczekiwaną zmienność dla obszaru o zadanym zaludnieniu i porównywali ją z obserwowaną wartością, tutaj zachorowalności, na danym obszarze.

Na powyższym wykresie widać, że patrząc na częstość zachorowań Glasgow ma podobną zachorowalność jak Orkney Islands, ale jeżeli dodatkowo uwzględni się liczbę osób zamieszkałych na obu obszarach to Orkney Islands ma zachorowalność mieszczącą się w granicach losowych fluktuacji, a dla Glasgow zachorowalność ta jest znacząco powyżej oczekiwanej na bazie pomiarów z całego kraju. Ciekawe prawda. Kiedyś o meta-analizie napiszę więcej, bo warto. Co ciekawe o wykresie tunelowym (funnel-plot) bez skrępowania piszą w Wielkiej Brytanii takie gazety jak Guardian (http://www.guardian.co.uk/commentisfree/2011/oct/28/bad-science-diy-data-analysis). Jak widać nawet duże dzienniki mogą serwować rzetelne informacje a nie tylko plotki o tym co nowego u celebrytów.