shapper is on CRAN, it’s an R wrapper over SHAP explainer for black-box models

Written by: Alicja Gosiewska

In applied machine learning, there are opinions that we need to choose between interpretability and accuracy. However in field of the Interpretable Machine Learning, there are more and more new ideas for explaining black-box models. One of the best known method for local explanations is SHapley Additive exPlanations (SHAP).

The SHAP method is used to calculate influences of variables on the particular observation. This method is based on Shapley values, a technique borrowed from the game theory. SHAP was introduced by Scott M. Lundberg and Su-In Lee in A Unified Approach to Interpreting Model Predictions NIPS paper. Originally it was implemented in the Python library shap.

The R package shapper is a port of the Python library shap. In this post we show the functionalities of shapper. The examples are provided on titanic_train data set for classification.

While shapper is a port for Python library shap, there are also pure R implementations of the SHAP method, e.g. iml or shapleyR.

Installation

The shapper wraps up the Python library, therefore installation requires a bit more effort than installation of an ordinary R package.

Install the R package shapper

First of all we need to install shapper, this may be the stable release from CRAN

or the developer version form GitHub.

Install the Python library shap

Before you run shapper, make sure that you have installed Python.

Python library shap is required to use shapper. The shap can be installed both by Python or R. To install it through R, you an use function install_shap() from the shapper package.

If you experience any problems related to the installation of Python libraries or evaluation of Python code, see the reticulate documentation. The shapper access Python within reticulate, therefore the solution to the problem is likely to be in there ;-).

Would you survive sinking of the RMS Titanic?

The example usage is presented on the titanic_train dataset from the R package titanic. We will predict the Survived status. The other variables used by the model are: Pclass, Sex, Age, SibSp, Parch, Fare and Embarked.

Let’s build a model

Let’s see what are our chances assessed by the random forest model.

Prediction to be explained

Let’s assume that we want to explain the prediction of a particular observation (male, 8 years old, traveling 1-st class embarked at C, without parents and siblings.

Model prediction for this observation is .558 for survival.

Here shapper starts

To use the function shap() function (alias for individual_variable_effect()) we need four elements

  • a model,
  • a data set,
  • a function that calculated scores (predict function),
  • an instance (or instances) to be explained.

The shap() function can be used directly with these four arguments, but for the simplicity here we are using the DALEX package with preimplemented predict functions.

The explainer is an object that wraps up a model and meta-data. Meta data consists of, at least, the data set used to fit model and observations to explain.

And now it’s enough to generate SHAP attributions with explainer for RF model.

Plotting results

To generate a plot of Shapley values you can simply pass an object of class importance_variable_effect to a plot() function. Since we are interested in the class Survived = 1 we may add additional parameter that filter only selected classes.

Labels on y-axis show values of variables for this particular observation. Black arrows show predictions of model, in this case, probabilities of each status. Other arrows show effect of each variable on this prediction. Effects may be positive or negative and they sum up to the value of prediction.

On this plot we can see that model predicts that the passenger will survive. Changes are higher due to young age and 1st class, only the Sex = male decreases chances of survival for this observation.

More models

It is useful to contrast prediction of two models. Here we will show how to use shapper for such contrastive explanations.

We will compare randomForest with svm implemented in the e1071.

This model predict 32.5% chances of survival.

Shapley values plot may be modified. To show more than one model you can pass more individual_variable_plot objects.

To see only attributions use option show_predcited = FALSE.

More

Documentation and more examples are available at https://modeloriented.github.io/shapper/. The stable version of the package is on CRAN, the development version is on GitHub (https://github.com/ModelOriented/shapper). Shapper is a part of the DALEX universe.

Co nowego w ,,Przewodnik po pakiecie R – wydanie 4.0”?

Przewodnik1234okladka
Czwarte wydanie ,,Przewodnika po pakiecie R” trafiło do księgarń w połowie grudnia. Pierwszy nakład był mały i szybko się skończył, ale od połowy stycznia Przewodnik jest ponownie dostępny.
A ja mam trochę czasu by napisać co nowego można znaleźć w czwartym wydaniu.

Zmian jest wiele. Kierunki zmian są dwa. Po pierwsze, obniżyć próg wejścia dla osób, które dopiero zaczynają przygodę z analizą danych. Łagodnym wprowadzeniem w temat są pierwsze dwa rozdziały. W upraszczaniu tej części przydały się doświadczenia z Pogromców Danych (2000+ osób) i z różnych szkoleń dla nie-programistów.
Drugi kierunek zmian to szersze wprowadzenie do pakietów z grupy tidyverse oraz ułatwień, które oferuje RStudio. Weterani R mają różne ulubione edytory i rozwiązania codziennych problemów, ale dla osób rozpoczynających przygodę z pewnością najefektywniejszą drogą wejścia jest połączenie RStudio i pakietów z tidyverse. Również osoby pracujące z R od lat mogą z zaskoczeniem odkryć, że praca z datami jest bardzo prosta dzięki pakietowi lubridate (ok, lubridate ma już kilka lat na karku) lub że praca z czynnikami jest prosta dzięki pakietowi forcats.

Wzorem poprzednich wydań, pierwsze 3 rozdziały (150 stron) są dostępne bezpłatnie jako pdf online tutaj.

Rozdział 1 – Wprowadzenie
W pierwszym rozdziale znajduje się krótki opis narzędzia jakim jest język R i edytor RStudio.
Zaczynam od pytania ,,Dlaczego warto poznać R?”,
Czytelnik może zapoznać się z przykładowymi fragmentami kodu R do pobierania danych z internetu (z nadzieją na reakcję czytelnika ,,WoW, to się da zrobić w 5 linijkach! Ja też tak chcę!!!”), wizualizacji pobranych danych oraz prostego modelowania statystycznego. Wszystko w zaledwie kilku linijkach kodu, możliwe dzięki dużej ekspresji języka.
Jeżeli ktoś jeszcze nie wie, czy praca z R jest dla niego, ten rozdział pomoże podjąć decyzję.
Jest tutaj też krótka historia rozwoju R, od początków S po lawinowy rozwój R w ostatnich latach/miesiącach.

Czytaj dalej Co nowego w ,,Przewodnik po pakiecie R – wydanie 4.0”?

R + Kraków = eRka

erka11

We wtorek w Krakowie eRka. Poniżej zaproszenie od organizatorów. Za dwa tygodnie SER. Liczymy na sporą frekwencję, na ostatnie spotkanie przyszło ~100 osób (~66% z zapisanych).


Bartosz Sękiewicz

Zapraszamy na 11 spotkanie entuzjastów R w Krakowie, które odbędzie się w najbliższy wtorek (6 grudnia) o godz. 18:30 na Politechnice Krakowskiej (ul. Podchorążych 1, sala 101).

Rozpoczynamy wystąpieniem o bardzo istotnej tematyce, zwłaszcza z punktu widzenia Krakowa. Bartosz Czernecki opowie nam jak w analityczny sposób opisać zjawisko smogu, jakie metody możemy wykorzystywać oraz jakie są czynniki umożliwiające dokładne prognozowanie stężenia szkodliwych pyłów w powietrzu. Będzie to bardzo ciekawy przykład połączenia specjalistycznej wiedzy badacza z metodami machine learning, a wszystko to oczywiście w R.

W trakcie drugiego wystąpienia Zygmunt Zawadzki opowie o wadach i zaletach R w porównaniu do C++, przede wszystkim w kontekście szybkości działania. Programując w R nie zawsze zdajemy sobie sprawę z istnienia funkcjonalności, które z jednej strony pozwalają skupić się tylko na tym co z punktu widzenia analityka najważniejsze, ale za to z drugiej nakładają wiele ograniczeń. Zygmunt przedstawi bardzo ciekawe rozwiązanie tego dylematu i opowie jak możemy łączyć ze sobą te dwa języki programowania.

Zmieniamy delikatnie formułę spotkań i teraz przerwa na pizzę będzie pomiędzy pierwszym i drugim wystąpieniem.

Szczegółowy plan spotkania, bio prelegentów oraz abstrakty prezentacji można znaleźć na naszej stronie www.erkakrakow.pl.

Rejestracja na wydarzenie poprzez meetup oraz facebook.

Zapraszamy serdecznie!

Który telefon się wcześniej zepsuje?

Uczestniczyłem ostatnio w ciekawym (od strony metodologicznej) projekcie (patrz też komentarz na dole). A ponieważ przy okazji używałem kilku ciekawych narzędzi, więc poniżej krótko o nich opowiem.

Problem:

Przypuśćmy, że mamy taki problem. Mamy telefony dwóch marek, np A i B. Załóżmy też, że w tych telefonach może zepsuć się bateria lub ekran, innych uszkodzeń nie rozważamy. Jeżeli się cokolwiek zepsuje, telefon jest wymieniany na nowy, nie ma sensu go naprawiać.
Chcemy teraz sprawdzić czy tempo psucia się ekranu w telefonach marki A jest inne niż w telefonach marki B. Podobnie chcemy porównać działanie baterii.
Przypuśćmy, że w badaniu obserwujemy dużo (tysiące) telefonów każdej z marek. Badanie prowadzimy przez 3 lata, ale niektóre telefony uczestniczą w badaniu od roku, inne od 2 lub 3 lat.
Pytanie, jaką metoda porównać te marki?

Podejście brutalne 1:

Najbardziej brutalne podejście to policzenie jaki procent ekranów zepsuł się w marce A a jaki w marce B i porównanie tych procentów np. dokładnym testem Fishera (w R funkcja fisher.test {stats}, więcej o teście tutaj) lub regresją logistyczną.
Ale…
To nie jest dobre podejście, ponieważ telefony nie były obserwowane przez tak samo długi czas. Mogło się zdarzyć, że w jednej grupie pozyskano więcej telefonów na początku, 3 lata temu, więc ich czas obserwacji był dłuższy. A może w innej grupie pozyskano telefony niedawno, np. rok temu, czas obserwacji jest krótszy i dlatego mniej ich się zepsuło.

Podejście brutalne 2:

Skoro czas obserwacji jest różny, to porównajmy liczbę dni do czasu zepsucia się ekranu lub baterii. Możemy to zrobić np. testem Wilcoxona (w R funkcja wilcox.test {stats} więcej o teście np. tutaj).
Ale…
To nie jest dobre podejście, ponieważ nie wszystkie telefony się zepsuły, nie dla wszystkich znam czas do zepsucia. Być może w jednej grupie średni czas do uszkodzenia jest niższy, ale też mniej telefonów się zepsuło. Nie chcemy też odrzucać informacji o tych wszystkich telefonach, które się nie zepsuły.

Mniej brutalne podejście 3:

Skoro mamy różne czasy obserwacji i nie dla wszystkich obserwacji wiemy jaki był czas uszkodzenia ekranu lub baterii, to zastosujmy metody analizy przeżycia. Przedstawmy czas działania ekranu/baterii za pomocą krzywych przeżycia (np. Kaplan-Meier) i sprawdźmy czy marki różnią się używając np. testu log rank (w R funkcja survdiff {survival}, więcej o teście).
Ale…
Obserwujemy dwa zdarzenia, uszkodzenie ekranu lub uszkodzenie baterii. Jeżeli wystąpi jedno uszkodzenie to nie widzimy drugiego (telefon oddajemy i nie wiemy kiedy nastąpiłaby druga awaria). Ale analizując uszkodzenia baterii, nie możemy traktować uszkodzeń ekranu jako obserwacje cenzorowane ponieważ nie musi być to cenzorowanie nieinformatywne (brak spełnionych założeń). Być może uszkodzony ekran wpływa na sposób pracy baterii, a być może wadliwe ekrany i baterie produkowane są w jednej fabryce, itp.

Podejście najmniej brutalne 4:

Skoro różne telefony obserwujemy przez różny czas, nie dla wszystkich znamy całkowity czas działania, czyli mamy cenzorowanie, ale obserwujemy dwa (lub więcej) elementów, które mogą ulec uszkodzeniu, to możemy modelować je łącznie używając modelu konkurujących ryzyk.

W R można to zrobić używając pakietu mstate (artykuł z JSS) czy też od jakiegoś czasu można to zrobić pakietem survival (teraz funkcja Surv{survival} może przyjąć argument type=mstate).
Do porównania krzywych można wykorzystać test Gray’s log-rank test. A w pakiecie ggplotit jest nawet przeciążona funkcja ggplotit() rysujące skumulowane funkcje zdarzeń dla obiektów klasy survfitms lub Cuminc.

btw:
Projekt na którym pracowałem dotyczył biologii molekularnej i analizy przeżycia pacjentów, ale dla przykładu i wnioskowania to akurat bez znaczenia. Podobna metodologię można zastosować do analizy usterek w samochodach czy kliknięć internautów na stronach www.

european R users meeting 12-14.10.2016

6_edited

Po wakacjach, 12-14 października, w Poznaniu odbędzie się pierwsza europejska konferencja użytkowników R.

Na 12 października zaplanowane są warsztaty, a czwartek – piątek 13-14 października to regularna konferencja. Listę zaproszonych prelegentów i wiele innych informacji można znaleźć na stronie konferencji http://erum.ue.poznan.pl/.

Do zobaczenia w Poznaniu w październiku!

5_edited

Hack the Proton. A data-crunching game from the Beta and Bit series

logo_eng
I’ve prepared a short console-based data-driven R game named ,,The Proton Game’’. The goal of a player is to infiltrate Slawomir Pietraszko’s account on a Proton server. To do this, you have to solve four data-based puzzles.

The game can be played by beginners as well as heavy users of R. Survey completed by people who completed the beta version of this game shows that the game gives around 15 minutes of fun to people experienced in R and up to around 60 minutes to people that just start programming and using R. More details about the results from beta-version are presented on the plot on the bottom.

PieczaraPieraszki

Czytaj dalej Hack the Proton. A data-crunching game from the Beta and Bit series

The PISA2009lite package is released

This post introduces a new R package named PISA2009lite. I will show how to install this package, what is inside and how to use it.

Introduction

PISA (Programme for International Student Assessment) is a worldwide study focused on measuring performance of 15-year-old school pupils. More precisely, scholastic performance on mathematics, science and reading is measured in more than 500 000 pupils from 65 countries.

First PISA study was performed in 2000, second in 2003, and then in 2006, 2009 and the last one in 2012. Data from the last study will be made public on December 2013. Data from previous studies are accessible through the PISA website http://www.oecd.org/pisa/.
Note that this data set is quite large. This is why the PISA2009lite package will use more than 220MB of your disk [data sets are compressed on disk] and much more RAM [data sets will be decompressed after [lazy]loading to R].

Let's see some numbers. In PISA 2009 study, number of examined pupils: 515 958 (437 variables for each pupil), number of examined parents: 106 287 (no, their questions are not related to scholastic performance), number of schools from which pupils were sampled: 18 641. Pretty large, complex and interesting dataset!

On the official PISA webpage there are instructions how to read data from 2000-2009 studies into SAS and SPSS statistical packages. I am transforming these dataset to R packags, to make them easier to use for R users.
Right now, PISA2009lite is mature enough to share it. There are still many things to correct/improve/add. Fell free to point them [or fix them].

This is not the first attempt to get the PISA data available for R users. On the github you can find 'pisa' package maintained by Jason Bryer (https://github.com/jbryer/pisa) with data from PISA 2009 study.
But since I need data from all PISA editions, namely 2000, 2003, 2006, 2009 and 2012 I've decided to create few new packages, that will support consistent way to access data from different PISA studies.

Open the R session

The package is on github, so in order to install it you need just

now PISA2009lite is ready to be loaded

You will find five data sets in this package [actually ten, I will explain this later]. These are: data from student questionnaire, school questionnaire, parent questionnaire, cognitive items and scored cognitive items.

You can do a lot of things with these data sets. And I am going to show some examples in next posts.

Country ranking in just few lines of code

But as a warmer let's use it to calculate average performance in mathematics for each country.

Note that student2009$W_FSTUWT stands for sampling weights, student2009$PV1MATH stands for first plausible value from MATH scale while student2009$CNT stands for country

And plot it.

plot of chunk unnamed-chunk-6

Ponieważ możemy! czyli o mapie na której widać 38 511 800 Polaków

Czy chcielibyście zobaczyć mapę Polski, na której zaznaczony jest każdy Polak? Wizualizację ponad 38 milionów osób rozrzuconych mniej lub bardziej losowo na obszarze ponad 300 tysięcy km^2? Jeżeli w tej chwili w waszych głowach zapala się pytanie ‚po co?’, nie czytajcie dalej. Jeżeli zaś już widzicie taką mapę oczami wyobraźni, poniższy wpis bardzo Wam się spodoba.

Celem tego projektu było pokazanie w interesujący sposób informacji z Narodowego Spisu Powszechnego 2011. Jedną z inspiracji była interaktywna mapa przedstawiająca spis powszechny w Stanach Zjednoczonych o której przeczytać można tutaj. W ramach dzisiaj opisywanego projektu opracowano mapę Polski składającą się z ponad 38 milionów punktów, każdy punkt odpowiadający jednej osobie, rozmieszczenie punktów odpowiadające rzeczywistej gęstości zaludnienia.

Ten wpis jak i cały projekt został wykonany przez Pawła Wiechuckiego w ramach wolontariatu dla naszej fundacji.
Kiedyś podczas jakiejś prezentacji usłyszałem, że organizacje non-profit to ciekawe miejsca, ponieważ przyciągają osoby, które za darmo tworzą rzeczy bezcenne, dlatego że chcą by te rzeczy powstały. Projekt Pawła Wiechuckiego jest świetnym przykładem takiej aktywności. Nie było łatwo, trzeba było pokonać wiele trudności, których nie widać przy bardziej standardowych zastosowaniach [np. R nie potrafi wygenerować w formacie wektorowym kwadratu o boku mniejszym niż pół punktu tzn. 1/144 cala, nie można też wygenerować dowolnie dużego rysunku, a przynajmniej 100GB RAM to za mało, sama mapa to połowa sukcesu, bufor drukarki może być niewystarczający by taką mapę wydrukować itp], ale się udało.

Przy okazji dziękuję też firmie Iqor Polska za udostępnienie plotera drukującego na 42 calowej rolce, dzięki temu udało się wydrukować mapę na połączonych arkuszach o łącznej powierzchni ponad 6 metrów kwadratowych, powiesić na ścianie i zrobić zdjęcie młodej wskazującej kropkę przedstawiającą babcię.

Nieszablonowa mapa Polski

Paweł Wiechucki

Czytaj dalej Ponieważ możemy! czyli o mapie na której widać 38 511 800 Polaków

Pięć raportów od wielkanocnego zająca

Tak naprawdę to nie od zająca, ale od Krzyśka Trajkowskiego.
Ale zacznijmy od początku.

Jakiś czas temu [to już cztery lata?] pracowałem nad zbiorem ,,luźnych” notatek w okolicy eksploracji danych, czego wynikiem był dokument ,,Na przełaj przez Data Mining”. Do prac nad tym dokumentem dołączył Krzysiek i w aktualnej [jeszcze nie ukończonej] wersji, połowa rozdziałów jest jego autorstwa.

Poza tą pozycją wspomniany już autor przygotował pięć interesujących dokumentów, które być może zaciekawią czytelników tego bloga.
Poniżej lista dokumentów wraz z odnośnikami do plików pdf.
Miłej lektury.

  1. ,,Przegląd pakietów do optymalizacji liniowej”, Krzysztof Trajkowski (2011).
    Można przeczytać o funkcjach solveLP(linprog), lp(lpSolve), lp.transport(lpSolve), lp.assign(lpSolve).
  2. ,,Przegląd pakietów do optymalizacji nieliniowej”, Krzysztof Trajkowski (2012).
    Można przeczytać o funkcjach optim(stats) i fminsearch(pracma) oraz solve.QP(quadprog), constrOptim(stats), constrOptim.nl(alabama), solnp(Rsolnp) i nloptr(nloptr).
  3. ,,Przegląd pakietów do analizy zmiennych niezależnych”, Krzysztof Trajkowski (2012).
    Można przeczytać o testach parametrycznych i nieparametrycznych dla dwóch lub większej liczby grup.
  4. ,,Przegląd pakietów do analizy zmiennych zależnych”, Krzysztof Trajkowski (2012).
    Można przeczytać o testach parametrycznych i nieparametrycznych dla danych sparowanych/związanych.
  5. ,,Analiza i wizualizacja danych naukowych”, Krzysztof Trajkowski (2013).
    Można przeczytać o pakietach ggplot2 i shape.

Słowa uznania dla autora tych dokumentów można zostawiać w komentarzach.