Jak przetestować / udowodnić, że dane są zerowane?

9

Mam problem, który moim zdaniem powinien być prosty, ale nie mogę go pojąć. Patrzę na zapylanie nasion, mam rośliny (n = 36), które kwitną w klastrach, próbuję 3 kępy kwiatów z każdej rośliny i 6 strąków nasion z każdej grupy (łącznie 18 strąków nasion z każdej rośliny). Strąki mogą mieć zapylane od 0 do maksymalnie 4 nasion. Dane są więc liczone z górną granicą. Zauważyłem, że średnio ~ 10% nasion jest zapylanych, ale gdziekolwiek od 1 do 30% na danej roślinie, więc dane są zbyt rozproszone, i oczywiście brakuje 4 replik skupień na 3 roślinach, więc nie idealnie symetrycznych .

Pytanie, które zadaję, brzmi: czy dane te wspierają pomysł, że roślina wymaga zapylaczy do zbioru nasion.

Zauważyłem, że rozkład liczby nasion w strąku wygląda na to, że jest więcej 0 zapylonych strąków nasiennych (6-9 strąków na 16) i więcej 3 i 4 zapylonych strąków nasiennych (2-4 na każdy) niż należy się spodziewać, jeśli nasiona w populacji zostaną po prostu losowo zapylone. Zasadniczo uważam, że jest to klasyczny przykład danych o zerowym napełnieniu, najpierw owad albo wcale nie odwiedza kwiatu (jeden generator zerowy), a jeśli tak, zapyla 0-4 nasion w innym rozkładzie. Alternatywna hipoteza jest taka, że ​​roślina jest częściowo samozapylona, ​​a zatem można oczekiwać, że każde ziarno będzie miało takie samo prawdopodobieństwo zapylenia (dane te sugerują około 0,1 szansy, co oznacza 0,01 szansy na dwa nasiona w tym samym strąku itp.) .

Ale po prostu chcę wykazać, że dane najlepiej pasują do jednej lub drugiej dystrybucji, a nie DO ZAKRESU ZIP lub ZINB na danych. Myślę, że jakakolwiek metoda zastosuję powinna brać pod uwagę faktyczną liczbę zapylanych nasion i liczbę strąków pobranych z każdej rośliny. Najlepszą rzeczą, jaką wymyśliłem, jest zrobienie czegoś w rodzaju paska do butów, w którym losowo przypisuję liczbę zapylonych nasion dla danej rośliny do liczby strąków nasion, z których próbowałem, zrób to 10.000 razy i sprawdź, jakie jest prawdopodobieństwo dane eksperymentalne dla danej rośliny pochodzą z tego losowego rozkładu.

Po prostu czuję, że jest w tym coś o wiele łatwiejszego niż ładowanie brutalnej siły, ale po wielu dniach myślenia i poszukiwań rezygnuję. Nie mogę po prostu porównać z rozkładem Poissona, ponieważ jest on górną granicą, nie jest dwumianowy, ponieważ muszę jakoś wygenerować oczekiwany rozkład. jakieś pomysły? I używam R, więc radzę tam (szczególnie jak najbardziej elegancko wygenerować 10 000 losowych rozkładów n kulek w 16 pudełkach, z których każda może zawierać maksymalnie 4 piłki) byłaby mile widziana.

DODANO 9/07/2012 Po pierwsze, dziękuję wszystkim za zainteresowanie i pomoc. Przeczytanie odpowiedzi skłoniło mnie do przeredagowania mojego pytania. Mówię o tym, że mam jedną hipotezę (którą na razie uważam za zerową), że nasiona są zapylane losowo między strąkami, a moją alternatywną hipotezą jest to, że strąki nasienne z co najmniej 1 zapylonym ziarnem są bardziej prawdopodobne mieć wiele zapylanych nasion, niż można by oczekiwać w wyniku losowego procesu. Podałem prawdziwe dane z trzech zakładów jako przykłady ilustrujące to, o czym mówię. Pierwsza kolumna to liczba zapylonych nasion w strąku, druga kolumna to częstotliwość strąków o tej liczbie nasion.

roślina 1 (łącznie 3 nasiona: 4% zapylenie)

liczba nasion :: pod.freq

0 :: 16

1 :: 1

2 :: 1

3 :: 0

4 :: 0

roślina 2 (łącznie 19 nasion: 26% zapylania)

liczba nasion :: pod.freq

0 :: 12

1 :: 1

2 :: 1

3 :: 0

4 :: 4

roślina 3 (łącznie 16 nasion: 22% zapylenie)

liczba nasion :: pod.freq

0 :: 9

1 :: 4

2 :: 3

3 :: 2

4 :: 0

W roślinie nr 1 zapylono tylko 3 nasiona w 18 strąkach, jeden strąk miał jedno nasiona, a jeden strąk miał dwa nasiona. Myśląc o procesie losowego dodawania jednego ziarna do strąków, pierwsze dwa nasiona trafiają do własnego strąka, ale dla trzeciego nasienia dostępnych jest 6 miejsc w strąkach, które mają już jedno ziarno, ale 64 miejsca w 16 strąkach bez nasion, więc najwyższe prawdopodobieństwo, że strąk z 2 nasionami jest tutaj, wynosi 6/64 = 0,094. To trochę niskie, ale nie bardzo ekstremalne, więc powiedziałbym, że ta roślina jest zgodna z hipotezą losowego zapylania we wszystkich nasionach z ~ 4% szansą na zapylenie. Ale roślina 2 wygląda dla mnie o wiele bardziej ekstremalnie, z 4 strąkami całkowicie zapylonymi, a jednak 12 strąków bez niczego. Nie jestem do końca pewien, jak bezpośrednio obliczyć szanse na ten rozkład (stąd mój pomysł na bootstrap), ale sądzę, że szanse na to, że rozkład ten wystąpi losowo, jeśli każde ziarno ma ~ 25% szansy na zapylenie, jest dość niskie. Roślina nr 3 Naprawdę nie mam pojęcia, myślę, że dla losowej dystrybucji jest więcej zer i 3, niż należy się spodziewać, ale mam przeczucie, że ten rozkład dla tej liczby nasion jest znacznie bardziej prawdopodobny niż rozkład dla rośliny nr 2, i może nie być tak mało prawdopodobne. Ale oczywiście chcę wiedzieć na pewno i we wszystkich roślinach. Wydaje mi się, że dla losowego rozmieszczenia jest więcej zer i 3 niż się spodziewałem, ale mam przeczucie, że rozkład dla tej liczby nasion jest znacznie bardziej prawdopodobny niż rozkład dla rośliny nr 2 i może nie być tak mało prawdopodobny. Ale oczywiście chcę wiedzieć na pewno i we wszystkich roślinach. Wydaje mi się, że dla losowego rozmieszczenia jest więcej zer i 3 niż się spodziewałem, ale mam przeczucie, że rozkład dla tej liczby nasion jest znacznie bardziej prawdopodobny niż rozkład dla rośliny nr 2 i może nie być tak mało prawdopodobny. Ale oczywiście chcę wiedzieć na pewno i we wszystkich roślinach.

Na koniec chcę napisać stwierdzenie, takie jak: „Rozkład zapylonych nasion w strąkach nasiennych odpowiada (lub nie pasuje) hipotezie, że rośliny nie są po prostu częściowo samozgodne, ale wymagają wizyty zapylacza, aby wykonać zbiór nasion. (wyniki testu statystycznego). ” To tak naprawdę tylko część mojej wybiegającej w przyszłość sekcji, w której mówię o tym, jakie eksperymenty przeprowadzę w następnej kolejności, więc nie jestem zdesperowany, aby było to jedno lub drugie, ale chcę to wiedzieć, jeśli to możliwe. Jeśli nie mogę zrobić tego, co próbuję zrobić z tymi danymi, chciałbym to również wiedzieć!

Na początku zadałem dość szerokie pytanie, ponieważ jestem ciekawy, czy istnieją jakieś dobre testy, które mogłyby pokazać, czy dane powinny w pierwszej kolejności przejść do modelu o zerowym napełnieniu. Wydaje się, że wszystkie przykłady, które widziałem, mówią - „spójrz, jest tu wiele zer i jest na to rozsądne wytłumaczenie, więc zastosujmy model z napompowaniem zerowym”. Właśnie to robię teraz na tym forum, ale miałem doświadczenie w ostatnim rozdziale, w którym użyłem glm Poissona do danych zliczania, a mój jeden z moich przełożonych powiedział: „Nie, glm są zbyt złożone i niepotrzebne, dane powinny przejdź do tabeli zdarzeń ”, a następnie przesłałem mi zrzut danych z ogromnej tabeli zdarzeń wygenerowanej przez ich drogi pakiet statystyk, który dał te same wartości p dla wszystkich moich czynników + interakcje do trzech cyfr znaczących !! Staram się więc, aby statystyki były jasne i proste, i upewnij się, że rozumiem je wystarczająco dobrze, aby solidnie bronić moich wyborów, czego nie sądzę, żebym mógł zrobić dla modelu z napompowaniem zerowym. Użyłem zarówno quasibinomial (dla całych roślin, aby pozbyć się pesudoreplicaitonu), jak i mieszanego modelu dla powyższych danych, aby porównać zabiegi i odpowiedzieć na moje główne pytania eksperymentalne, wydaje się, że wykonują tę samą pracę, ale zamierzam również baw się dzisiaj z ZINB, aby zobaczyć, jak dobrze to działa. Myślę, że jeśli mogę jednoznacznie wykazać, że te dane są na początku silnie skupione w klastrze (lub są zawyżone do zera), a następnie podać dobry biologiczny powód takiego wystąpienia, będę znacznie lepiej przygotowany do późniejszego wyciągnięcia ZINB, niż do po prostu porównaj jeden z quasibinomial / mixed model i argumentuj, ponieważ daje lepsze wyniki, tego powinienem użyć. co nie wydaje mi się, że mogę teraz zrobić dla modelu z napompowaniem zerowym. Użyłem zarówno quasibinomial (dla całych roślin, aby pozbyć się pesudoreplicaitonu), jak i mieszanego modelu dla powyższych danych, aby porównać zabiegi i odpowiedzieć na moje główne pytania eksperymentalne, wydaje się, że wykonują tę samą pracę, ale zamierzam również baw się dzisiaj z ZINB, aby zobaczyć, jak dobrze to działa. Myślę, że jeśli mogę jednoznacznie wykazać, że te dane są na początku silnie skupione w klastrze (lub są zawyżone do zera), a następnie podać dobry biologiczny powód takiego wystąpienia, będę znacznie lepiej przygotowany do późniejszego wyciągnięcia ZINB, niż do po prostu porównaj jeden z quasibinomial / mixed model i argumentuj, ponieważ daje lepsze wyniki, tego powinienem użyć. co nie wydaje mi się, że mogę teraz zrobić dla modelu z napompowaniem zerowym. Użyłem zarówno quasibinomial (dla całych roślin, aby pozbyć się pesudoreplicaitonu), jak i mieszanego modelu dla powyższych danych, aby porównać zabiegi i odpowiedzieć na moje główne pytania eksperymentalne, wydaje się, że wykonują tę samą pracę, ale zamierzam również baw się dzisiaj z ZINB, aby zobaczyć, jak dobrze to działa. Myślę, że jeśli mogę jednoznacznie wykazać, że te dane są na początku silnie skupione w klastrze (lub są zawyżone do zera), a następnie podać dobry biologiczny powód takiego wystąpienia, będę znacznie lepiej przygotowany do późniejszego wyciągnięcia ZINB, niż do po prostu porównaj jeden z quasibinomial / mixed model i argumentuj, ponieważ daje lepsze wyniki, tego powinienem użyć. Użyłem zarówno quasibinomial (dla całych roślin, aby pozbyć się pesudoreplicaitonu), jak i mieszanego modelu dla powyższych danych, aby porównać zabiegi i odpowiedzieć na moje główne pytania eksperymentalne, wydaje się, że wykonują tę samą pracę, ale zamierzam również baw się dzisiaj z ZINB, aby zobaczyć, jak dobrze to działa. Myślę, że jeśli mogę jednoznacznie wykazać, że te dane są na początku silnie skupione w klastrze (lub są zawyżone do zera), a następnie podać dobry biologiczny powód takiego wystąpienia, będę znacznie lepiej przygotowany do późniejszego wyciągnięcia ZINB, niż do po prostu porównaj jeden z quasibinomial / mixed model i argumentuj, ponieważ daje lepsze wyniki, tego powinienem użyć. Użyłem zarówno quasibinomial (dla całych roślin, aby pozbyć się pesudoreplicaitonu), jak i mieszanego modelu dla powyższych danych, aby porównać zabiegi i odpowiedzieć na moje główne pytania eksperymentalne, wydaje się, że wykonują tę samą pracę, ale zamierzam również baw się dzisiaj z ZINB, aby zobaczyć, jak dobrze to działa. Myślę, że jeśli mogę jednoznacznie wykazać, że te dane są na początku silnie skupione w klastrze (lub są zawyżone do zera), a następnie podać dobry biologiczny powód takiego wystąpienia, będę znacznie lepiej przygotowany do późniejszego wyciągnięcia ZINB, niż do po prostu porównaj jeden z quasibinomial / mixed model i argumentuj, ponieważ daje lepsze wyniki, tego powinienem użyć.

Ale nie chcę zbytnio odwracać uwagi od mojego głównego pytania, jak mogę ustalić, czy moje dane są bardziej zera bardziej niż oczekiwane z losowego rozkładu? W moim przypadku odpowiedź na to pytanie jest tym, co naprawdę mnie interesuje, a ewentualna korzyść z uzasadnienia modelu jest premią.

Jeszcze raz dziękuję za poświęcony czas i pomoc!

Pozdrawiam, BWGIA

BWGIA
źródło
dlaczego nie chcesz dopasować dwumianowego modelu z napompowaniem zerowym?
atiretoo
czy hipoteza „częściowego samozapylenia” dotyczy wyłącznie hipotezy „zapylacza”? Jeśli tak, to twój drugi model byłby po prostu modelem dwumianowym z prawdopodobieństwem p i rozmiarem = 4.
atiretoo

Odpowiedzi:

5

Wydaje mi się to stosunkowo prostym (nieliniowym) modelem mieszanym. Masz strąki nasienne zagnieżdżone w klastrach zagnieżdżonych w roślinach i możesz dopasować model dwumianowy z losowymi efektami na każdym etapie:

    library(lme4)
    binre <- lmer( pollinated ~ 1 + (1|plant) + (1|cluster), data = my.data, family = binomial)

lub ze zmiennymi towarzyszącymi, jeśli je masz. Jeśli kwiaty zapylą się samoczynnie, możesz zauważyć łagodne skutki wynikające z naturalnej zmienności w tym, jak same rośliny są żywotne. Jeśli jednak większość zmienności odpowiedzi wynika z powiedzmy zmienności skupisk, miałbyś silniejszy dowód zapylania przez owady, które mogą odwiedzać tylko wybrane skupiska na roślinie. Idealnie byłoby chcieć nieparametrycznego rozkładu efektów losowych niż Gaussa: masa punktowa zero, bez wizyt owadów, i masa punktowa o wartości dodatniej - jest to zasadniczo model mieszany, o którym myślał Michael Chernick. Możesz dopasować to do pakietu GLLAMM Stata, byłbym zaskoczony, gdyby nie było to możliwe w R.

Prawdopodobnie dla czystego eksperymentu chciałbyś mieć rośliny w środku, a przynajmniej w miejscu bez dostępu owadów i zobaczyć, ile nasion zostanie zapylonych. To prawdopodobnie odpowiadałoby na wszystkie pytania w bardziej rygorystyczny metodologicznie sposób.

StasK
źródło
Spróbuję tego, myślę, że pomoże mi to odpowiedzieć na własne pytania, ale nie jestem pewien, jak przekonać innych. Zajmujesz się drugą częścią, staram się myśleć o tym, w jaki sposób te dane kształtują przyszły bardziej ukierunkowany eksperyment.
BWGIA,
1

Wydaje mi się, że jest to rozkład mieszanki dla każdego owada. Z prawdopodobieństwem p owad ląduje z prawdopodobieństwem 1-p ląduje i rozdaje 0 do 4 nasion. Ale jeśli nie masz informacji o tym, czy owad ląduje na roślinie, nie możesz rozróżnić dwóch sposobów uzyskania 0. Więc możesz pozwolić p być prawdopodobieństwem dla 0, a następnie masz rozkład wielomianowy (p1, p2, p3, p4) gdzie pi jest prawdopodobieństwem nasion i, biorąc pod uwagę, że owady zapylają z zastrzeżeniem ograniczenia p1 + p2 + p3 + p4 = 1. Model ma pięć niewiadomych p, p1, p2, p3, p4 z ograniczeniem 0 = 0 dla każdego i. Przy wystarczającej ilości danych powinieneś być w stanie oszacować te parametry, być może stosując metodę ograniczonego maksymalnego prawdopodobieństwa.

Michael R. Chernick
źródło
Zgadzam się, ale nie chodzi o dopasowanie tego modelu, ale o wygenerowanie przewidywanych rozkładów w oparciu o dwie różne hipotezy biologiczne. Być może odpowiedzią jest dopasowanie ZIB i „innego modelu”, który pasuje do hipotezy samozaparcia, i porównanie ich.
atiretoo
@atiretoo czy model nie zapewnia szacunkowego rozkładu liczby zapylonych nasion, które można porównać do hipotetycznego rozkładu?
Michael R. Chernick,
Uzgodnione - jeśli masz odpowiednie modele dla 2 hipotez.
atiretoo
1

To jest odpowiedź na ostatnią część twojego pytania, jak szybko wygenerować dane, które chcesz dla hipotezy zapylacza:

n = 16
max = 4
p1 = 0.1
p2 = 0.9
Y1 = rbinom(10000*n,1,p1)
Y2 = matrix(Y1*rbinom(10000*n,4,p2),ncol=16)

Możesz także użyć rzibinom()w pakiecie VGAM. Chociaż nie jestem pewien, co chcesz z tym zrobić. Masz 2 wolne parametry, p1 i p2, które należy oszacować. Dlaczego nie zastosować modelu dwumianowego z napompowaniem zerowym do oszacowania ich na podstawie danych?

Powinieneś spojrzeć na pakiet VGAM, który pasuje między innymi do modeli ZIB. W rzeczywistości można uzyskać oczekiwany rozkład dla ZIB z funkcji VGAM dzibinom(), którego można użyć do porównania zaobserwowanego rozkładu z, jeśli znasz parametry odwiedzin i zapylania. Ponownie, naprawdę powinieneś pasować do modelu ZIB.

Jeśli twoja hipoteza częściowego samozapylenia dotyczy wyłącznie zapylania owadów, to oczekiwany rozkład jest po prostu dwumianowy i możesz oszacować parametry za pomocą dwumianowej rodziny glm lub może glmm z identyfikatorem rośliny jako efektem losowym. Jeśli jednak potrafią częściowo zapanować nad sobą ORAZ otrzymać zapylenie owadów, to wrócisz do konieczności mieszania dwóch rozkładów dwumianowych. W takim przypadku zbadałbym użycie OpenBUGS lub JAGS w celu dopasowania modelu za pomocą MCMC.

Po dopasowaniu dwóch modeli do danych należy porównać modele, aby sprawdzić, który z nich pasuje lepiej, za pomocą AIC lub BIC lub innej wybranej przez Ciebie metryki.

atiretoo - przywróć monikę
źródło
Dzięki za ten tatuaż, ale uruchomienie tego kodu wydaje się generować losową liczbę nasion, a także losową dystrybucję. Myślałem, że chcę, aby nubmer nasion został naprawiony (powiedzmy 19 nasion, patrz poniżej), a potem przekonam się, jak prawdopodobne jest dane rozłożenie tego dokładnego
ziarna
Opps, trafiłem za wcześnie i miałem na myśli „patrz wyżej”, ponieważ dodałem trochę informacji do mojego pytania. Intryguje mnie twój komentarz na temat używania AIC do porównywania modeli. Czy mogę to zrobić w różnych modelach (z tą samą zmienną odpowiedzi) o różnych rozkładach? Myślałem, że porównanie AIC jest poprawne tylko wtedy, gdy dodajesz / upuszczasz warunki do modelu, ale z tą samą rodziną dystrybucji określoną?
BWGIA,
Nie, to kluczowa przewaga AIC nad np. Wstecznym wyborem. O ile dane są takie same, możesz porównać AIC między różnymi modelami, nawet jeśli nie są zagnieżdżone. Trzeba uważać, aby oprogramowanie obliczało prawdopodobieństwa bez pomijania stałych, ale w ramach jednej funkcji można łatwo porównywać modele zagnieżdżone.
atiretoo