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”.