Jakich metod statystycznej analizy danych uczyć się?


The New England Journal of Medicine (dla bibliometroholików: 50 pkt MNiSW, impact factor ponad 50, pierwsza liga) opublikował dwa tygodnie temu krótki, ale bardzo ciekawy komentarz Statistical Methods. Przebadano 238 prac z roku 2015 pod kątem używanych technik statystycznych oraz zestawiono te wyniki z artykułami z poprzednich lat (od 1978, łącznie ~1000 artykułów).

Całość streszcza wykres po prawej (to link do oryginalnego artykułu). W dobrych czasopismach medycznych wykorzystywane są coraz bardziej zaawansowane metody statystyczne. Metody nieparametryczne i regresja wielokrotna w ponad 1/3 prac. Analiza przeżycia i analiza mocy testów, każde w ponad połowie artykułów (bardzo duży wzrost w ostatnich latach). Wielokrotne testowanie w co piątym artykule. Tylko co dwudziesty artykuł nie ma ,,statystyki”.

Przecinające się krzywe przeżycia

Spotkałem się ostatnio z ciekawym problemem.
Mamy dwie grupy pacjentów na dwóch różnych schematach leczenia i chcemy porównać ich dalsze losy, a konkretnie krzywe niepowodzenia leczenia (prawdopodobieństwo zgonu/wznowy).
Na pierwszy rzut oka klasyczna analiza przeżycia, test log-rank i po sprawie.

Ale szybko okazuje się, że krzywe przeżycia się przecinają, co więcej oczekiwać można tego po wcześniejszej rozmowie z lekarzem. Jeden schemat leczenia jest bardziej agresywny więc może prowadzić do gorszych rokowań krótkookresowych, ale lepszych w dalszej perspektywie.

Klasyczny test dla krzywych przeżycia oparty jest o odległość pomiędzy krzywymi, mierzoną jest jako ważona suma kwadratów odległości w poszczególnych punktach. Ale gdy krzywe się przecinają to taki test ma niską moc i nie ma sensu go stosować.

A więc co robić?
Ciekawe studium symulacyjne porównujące różne podejścia do testowania przecinających się krzywych zostało opublikowane dwa lata temu w Plos One (!).
Okazuje się, że dobrze sprawdza się rodzina testów Renyi, która jest oparta o supremum ważonych odległości pomiędzy krzywymi przeżycia.
W R te testy są zaimplementowane w pakiecie survMisc w funkcji comp. Jest to znacznie mocniejszy test dla przecinających się krzywych.

A przy okazji, okazuje się, że zmianę w hazardach w rozpatrywanym problemie dobrze ilustrują reszty Schonefelda. Poniższy wykres pokazuje, że hazard w jednej grupie jest znacznie wyższy do 12 miesiąca, a później gorsze losy czekają pacjentów drugiej grupy.

Oba wykresy wykonane pakietem survminer.

Opisy osi usunąłem ponieważ wyniki tych analiz jeszcze nie są opublikowane, ale też nazwy nie mają większego znaczenia.

1024 bity (i Bety)

c20

Miesiąc temu pisałem o akcji prowadzenia warsztatów ,,Jak zważyć psa linijką”.
Dzięki grantowi mPotęga planowaliśmy dotrzeć z warsztatami statystycznymi do 5 klas z różnych szkół.
Zainteresowanie okazało się jednak znacznie większe!
Do dziś udało się dotrzeć z warsztatami do ponad 30 nauczycieli i ponad 1000 (słownie: tysiąca!) dzieciaków.

Statystyka i programowanie w podstawówce? Za moich czasów tego nie było!

Czytaj dalej 1024 bity (i Bety)

Hakaton 'Puls miasta’ @ WhyR 2017

WhyR to Ogólnopolska Konferencja Użytkowników R, która odbędzie się 27-29 września 2017 na Politechnice Warszawskiej (więcej o WhyR). Dzień przed konferencją (26 września) planujemy przeprowadzić bardzo ciekawy hakaton z wykorzystaniem naprawdę dużych miejskich danych.

Jakich danych?
Hakaton realizowany jest w ramach projektu VaVeL (więcej o VaVeL) w którym partnerem jest Ośrodek Badań dla Biznesu PW (więcej o OBB), dzięki czemu będziemy mieli dostęp do danych z najróżniejszych sensorów w Warszawie. Jakich sensorów? Przykładowo dane o położeniu każdego tramwaju i autobusu w praktycznie każdej chwili (live + spora historia), dane o natężeniu ruchu pieszego w różnych punktach miasta, z publicznych kanałów informacyjnych i z wielu innych źródeł (rysunek po prawej to ślad z jednego dnia po sensorach z tramwaju 22). Masa danych. Mikołaj w tym roku przychodzi we wrześniu.

Jak to ogarnąć?
W ramach warsztatów poprowadzimy bezpłatne mini-wykłady z technologii BigData-owych, takich jak Hadoop czy Hive, dzięki czemu uczestnicy będą mogli i będą wiedzieć jak dostać się do tych gigantycznych zasobów. Ale nawet jeżeli ktoś nie przepada za żółtymi słoniami będzie mógł pracować na przetworzonych skrawkach danych lub też będzie mógł wesprzeć zespół od strony wizualizacji, burzy mózgów, tworzenia aplikacji mobilnych czy innych aplikacji.

Co będziemy robić?
Zbieramy różne pomysły na hackaton, ale liczymy też na burzę mózgów podczas samego wydarzenia. Analiza danych oceniających zatłoczenie przystanków na Mordorze? Aplikacja informująca ile się średnio spóźnia linia 10 w okolicach godziny 16? Wizualizacja transferu mieszkańców w różnych godzinach. Zobaczymy co z tego wyjdzie.

Jak się zarejestrować?
Więcej informacji o rejestracji pojawi się po feriach zimowych. Z pewnością warto śledzić stronę konferencji WhyR.

Is it a job offer for a Data Scientist?

screen-shot-2017-01-09-at-21-41-24
TL;DR
Konrad Więcko and Krzysztof Słomczyński (with tiny help from my side) have created a system that is tracing what skills are currently in demand among job offers for data scientists in Poland. What skills, how frequent and how the demand is changing over time.

The full description how this was done. static, +shiny.

Here: The shiny application for browsing skill sets.

Here: The R package that allows to access the live data.

Full version
A data science track (MSc level) will be very soon offered at MiNI department/Warsaw University of Technology and we (i.e. program committee) are spending a lot of time setting up the program. How cool it would be taking into account the current (and forecasted) demand on data science skills on the market? The project executed by Konrad Więcko and Krzysztof Słomczyński is dealing exactly with this problem.

So, the data is scrapped from pracuj.pl, one of the most popular (in Poland) websites with job offers.
Only offers interested to data scientists were used for further analyses.
How these offers were identified? In Poland the job title 'Data Scientist’ is (still) not that common (except linkedin). And there are different translations and also there is a lot of different job titles that may be interested for a data scientist.
So, Konrad and Krzysztof have developed a machine learning algorithm that scores how much a given offer matches a 'data scientist profile’ (whatever that means) based on the content of job offer description. And then, based on these predictions, various statistics may be calculated, like: number of offers, locations, demand on skills etc.

Here: Description of the dataset used for the model training.

Here: Description of the machine learning part of the project.

Trends for selected skill sets.

StatTuba w Twojej szkole (lub szkole Twoich dzieci)


StatTuba to projekt dotarcia do uczniów szkół podstawowych i średnich z wnioskowaniem opartym o dane. Pokażmy dzieciakom, że fajne rzeczy (czasem bardzo nieoczywiste) można robić, o ile tylko mamy odpowiednie dane.

Z akcją na PolakPotrafi.pl się nie udało, ale dzięki wsparciu z programu mPotęga (mBank), przygotowaliśmy materiały, pozwalające praktycznie każdemu zainteresowanemu nauczycielowi lub rodzicowi poprowadzenie warsztatów statystycznych w szkole. Materiały o których piszę niżej, są przeznaczone głównie dla klas 3-5.

Jak to działa?

1. Zainteresowany nauczyciel, rodzic lub inna osoba, która chciałaby w wybranej szkole poprowadzić warsztaty, wysyła na adres stattuba@gmail.com informację gdzie i kiedy chce te warsztaty przeprowadzić,
2. Do prowadzącego, na wskazany adres, wysyłamy całkowicie bezpłatnie materiały i instrukcje jak poprowadzić warsztaty (na jedną grupę średnio 30 papierowych kopii tego komiksu),
3. Nauczyciel lub rodzic prowadzi warsztaty, przesyła nam krótką (2-3 zdania) informację jak poszły warsztaty (najlepiej z jednym-dwoma zdjęciami) a my odsyłamy dyplom z podziękowaniami dla nauczyciela i dyrektora za udział w programie StatTuba,
4. Przewidujemy, że na warsztaty prowadzone w styczniu wyślemy materiały do 30 nauczycieli/prowadzących (10 już jest, szukamy kolejnych).

Na stronie http://betabit.wiki/warsztaty/ umieściliśmy szczegółowy wideo-opis tego, jak można poprowadzić warsztaty, jakie ćwiczenia przeprowadzić i jak wykorzystać przesłane materiały. W razie jakichkolwiek pytań chętnie udzielę szczegółowych informacji. Proszę o email na adres stattuba@gmail.com.

Materiały do warsztatów zostały przygotowane razem z Agnieszką Tomczyk i Martyną Śpiewak.

Bardzo zachęcam do poprowadzenia takich warsztatów w szkole swoich dzieci. To dla dzieci wielka frajda. Warto przesłać to ogłoszenie nauczycielowi matematyki (lub nauczania początkowego) w swojej szkole.

Czasem diabeł tkwi w szczegółach.

Ostatnio trafił do mnie taki ciekawy problem.

Problem:
Mamy dwa zbiory monet (oznaczmy typ 1/typ 2) z końca XIII w., chcemy sprawdzić się czy zostały wyprodukowane tą samą techniką = czy rozkłady ciężarów monet są takie same.
Wiemy, że ich średnie wielkości wagowe są takie same ale mogą różnić się np. skośnością. Jeżeli uznamy, skośności są różne to można przypuszczać, że monety w różnych zbiorach wyprodukowano różnymi technikami.
Średnia wielkość monet jest taka sama, bo (o ile nie dochodziło do oszustwa) z określonej ilości surowca (jednej marki czy jednej grzywny) wiadomo było ile powinno powstać monet. Ale nawet jeżeli średnia jest taka sama to techniki produkcji monet mogą się różnić udziałem lżejszych i cięższych monet.
Monet jednego typu jest 40, drugiego 32 egz.

A więc jak porównać te dwie grupy?

Gdyby obudzić w nocy studenta statystyki i zapytać go o to jak porównywać rozkłady w dwóch grupach, to ten pewnie by tylko wymamrotał przez sen 'test Kołmogorowa Smirnowa’ i przewrócił się na drugi bok.
W programie R ten test dla dwóch prób można wykonać funkcją ks.test{stats} (wiki)

Wikipedia mówi ’The two-sample K–S test is one of the most useful and general nonparametric methods for comparing two samples, as it is sensitive to differences in both location and shape of the empirical cumulative distribution functions of the two samples.’, znaczy że ten test wykrywa różnice i w średnich i skośności i innych nierównościach pomiędzy próbami. Można więc mieć nadzieję, że nawet gdy nie ma różnicy w średnich to różnica w skośnościach będzie wykryta.
Ale…
Test K-S wykrywa różne odstępstwa od równości z różną czułością. Najbardziej czuły jest na przesunięcia ale już na różnice w skośności jest mniej czuły.
Gdyby szukać analogii, można by porównać ten test ze strategią poszukiwania okularów rano. Strategią w której przez 60% czasu szuka się okularów na parapecie, przez 30% na biurku a przez 10% pod łóżkiem. Ta strategia pozwala na znalezienie okularów w każdym z tych trzech miejsc o ile wystarczająco długo szukamy. Ale zgodzicie się, że jeżeli mamy dodatkową wiedzę, że okularów nie ma na parapecie to lepiej wykorzystać inną, bardziej specyficzną strategię.

Podobnie jest z testem K-S. Odpowiednikiem czasu szukania okularów jest tutaj liczba monet potrzebnych do identyfikacji różnicy. Jeżeli nie wiemy co może różnić porównywane grupy to test K-S jest być może nienajgorszym wyborem. Ale jeżeli wiemy gdzie szukać różnic (np. w skośności) to lepiej użyć innych narzędzi.

Jakich?

Niestety [?] teoria oparta o rozkłady gaussowskie się nie przyda, bo te nie są skośne. Ale zwolennicy rozwiązań parametrycznych mogą wykorzystać rozkłady z parametrem położenia, skali i kształtu (GAMLSS) aby testować skośność. Więcej o tych modelach można przeczytać np. w pracy magisterskiej mojej dyplomantki Anny Lis.

Jednak, jeżeli lubicie dźwięk rozpędzającego się wiatraczka chłodzącego procesor, to interesującą alternatywą będzie test oparty o metodę bootstrap lub test permutacyjny.

Technika bootstrap (wiki) pozwala na szacowanie przedziałów ufności dla skośności w obu próbach. Szacując dokładność oszacowania skośności możemy już prosto skonstruować test dla dwóch prób. Schemat jest następujący: 1. z próby pierwszej losujemy ze zwracaniem 40 monet, 2. na otrzymanej próbie wyznaczamy skośność, 3. powtarzamy 1-2 wiele razy, np. 10000 razy, 4. robimy 1-3 dla drugiej próby, 5. porównujemy oszacowania skośności.

Test permutacyjny przemiesza obie grupy (wiki) i sprawdzi, czy po resamplingu połączonej próby i po ponownym podzieleniu wymieszanych obserwacji na dwie podpróby, skośności w nowo utworzonych podpróbach różnią się bardziej niż różniły się w oryginalnych populacjach. Tutaj schemat jest następujący: 1. Łączymy monety w jeden wektor 72 elementowy, 2. permutujemy wektor, 3. dzielimy wektor na dwie podpróby o wielkości 40 i 32, 4. wyznaczamy skośności na podpróbach, 5. wyznaczamy moduł różnicy w skośnościach z punktu 4, 6. powtarzamy 1-5 wiele razy, np. 10000 razy, 7. sprawdzamy czy ta różnica w skośnościach którą widzimy w oryginalnych danych jest większa niż różnice w takich przepermutowanych próbach.

Akurat dla tych konkretnych danych obie te techniki dają bardzo zbliżone wyniki, p-wartości w okolicach 0.04, co przy tak małych próbach można uznać nawet za wyraźny sygnał pochodzenia z różnych populacji.
K-S test nie widział żadnych różnic, raportował p-wartość powyżej 0.8.

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.

Statystyk na wakacjach

Miejsce: Park rozrywki w Szklarskiej Porębie, kolejka do kina 6D.
Aktorzy: [S]taystyk na wakacjach i [B]ileterka.
Czas: 14:55, na 5 minut przed seansem w ww. kinie. Seanse odbywają się co 30 minut. Przed wejściem ustawia się kolejka. 10 minut przed seansem osoby z kolejki zaczynają wchodzić do kina. Wchodzi pierwsze 25 osób.

– Na ten seans już nie ma miejsc, proszę przyjść na kolejny o 15:30 – informuje Bileterka.
– A ile minut przed seansem powinienem przyjść by były jeszcze miejsca? – grzecznie pyta Statystyk.
– 5 minut przed seansem, tak jak jest napisane w regulaminie – Bileterka wskazuje palcem regulamin.
– Ale teraz jestem 5 minut przed seansem i już nie ma miejsc – zauważa Statytyk. – Więc ile minut wcześniej powinienem przyjść aby były jeszcze miejsca? – docieka.
– To zależy od tego ile osób przyjdzie. Musi być Pan najpóźniej 5 minut przed seansem. – powtarza Bileterka zniecierpliwionym głosem.
– A ile minut przed seansem się zazwyczaj kończą bilety? – dopytuje Statystyk.

Mniej więcej tutaj dla mojej interlokutorki staje się jasne, że trafił się jej wyjątkowo dociekliwy/upierdliwy (strony mogą różnie określać tę cechę) osobnik. Jej odpowiedź jest już bardziej stanowcza.

– Różnie się kończą. To zależy ile osób przyjdzie na kolejny seans. A tego nikt nie wie – rozmówczyni niesłusznie zakłada, że odstraszy mnie brak precyzyjnych szacunków.

W tym miejscu przerwę relacjonowanie naszej rozmowy. Na kolejny seans przyszedłem 10 minut przed czasem i wszedłem mniej więcej w połowie kolejki.

Ale historia dopiero tutaj się zaczyna.

Przez kolejne dwie godziny moje szkraby szalały na dmuchańcach obok kina. Miałem trochę czasu by poobserwować kolejkę do kina, zebrać trochę danych i zastanowić się, jak sam bym odpowiedział na pytanie, które zadałem Bileterce.

Zagadnienie:

Oszacować ile minut przed seansem należy przyjść aby mieć 90% pewności, że wystarczy dla nas miejsc w kinie.

Dane:

Dla 4 seansów (dwie godziny obserwacji) mamy informację ile osób (najczęściej przychodzą całe rodziny) i ile minut przed seansem dołączyło do kolejki.

Model 1:

Rozwiązanie brutalne, praktycznie bez modelowania.
Dla każdego seansu liczymy ile minut przed seansem przyszła ostatnia osoba, która zmieściła się na sali. Dla naszych seansów było to odpowiednio 8,9,7,8 minut.

Rozwiązanie proste o uroku cepa. Bez modelu parametrycznego z czterech liczb trudno wyznaczyć 90% kwantyl. (Ok, można jeżeli jest się ultrasem bootstrapowcem).

Szukamy więc czegoś parametrycznego.

Model 2:

Zakładamy, że liczba osób dochodzących do kolejki opisana jest jednorodnym procesem Poissona.
Oznacza to, że zakładamy, że w pewnym okresie czasu, np. -15 do -5 minut przed seansem, chętni przychodzą pojedynczo ze stałą intensywnością (=nie w stałych odstępach czasu ale ze stałem prawdopodobieństwem pojawienia się).
Więcej o procesie Poissona np. tutaj.

I co dalej? Szacujemy intensywność przychodzenia osób (w tym modelu to średnia) i liczymy czas oczekiwania na przekroczenie przez proces Poissona bariery 22 osób (jeszcze my się musimy zmieścić).

Piękny parametryczny model.
Drażniące są tylko te nierealne założenia.
Może da się je osłabić.

Model 3:

Zakładamy, że liczba osób dochodzących do kolejki opisana jest złożonym procesem Poissona.
Złożony proces Poissona to połączenie zwykłego procesu Poissona (opisuje momenty, w których do kolejki dochodzi rodzina) oraz skoków o różnej wielkości (wielkość skoku to liczba osób w rodzinie, które dołączyły do kolejki, z obserwacji od 1 do 5, najczęściej 2-3). Jest to rozszerzenie modelu 2, w którym uwzględniamy to, że do kolejki na raz dołączyć może kilka osób.
Więcej o złożonym procesie Poissona np. tutaj.

I co dalej? Osobno szacujemy intensywność pojawiania się rodzin (podobnie model z jednym parametrem szacowanym średnią), osobno szacujemy rozkład wielkości rodziny. Mając te składowe, wyznaczamy (np. symulacyjnie) rozkład czasu przekroczenia bariery 22 osób.

Model coraz piękniejszy, wymaga estymacji parametrów dwóch rozkładów (czasu przyjścia i wielkości rodziny). Drażni jedynie to założenie o stałej intensywności pojawiania się rodzin na odcinku -15 min do -5 min przed seansem.

Model 4:

Wykorzystajmy złożony niejednorodny proces Poissona. Czyli to co powyżej, ale tym razem intensywność pojawiania się rodzin jest nieujemną funkcją na odcinku -30 min do -5 min. Na początku raczej bliska zera (kto ustawia się w kolejce na 20 minut przed seansem gdy nikt inny w kolejce nie stoi?) a później szybko skacząca w czasie -15 min do -5 min przed seansem (nauczeni doświadczeniem wiedzą, że warto zjawić się wcześniej).
To już jest TEN model. W zmiennej intensywności możemy nawet uwzględnić porę dnia, liczbę osób przebywających w parku rozrywki i kilka innych parametrów. Samą intensywność można szacować np. estymatorem jądrowym.

Jedynym problemem okazało się to, że o 18 zamykali park i nie dało się zebrać więcej danych.

Więcej o niejednorodnym procesie Poissona można przeczytać tutaj.

Inne pomysły na modele?

[*] Ilustracja pochodzi z opowiadania ,,Jak długo żyją Muffinki”.