Negatywno-dwumianowy GLM vs. transformacja logów dla danych zliczania: zwiększony poziom błędu typu I.

Niektórzy z was mogli przeczytać ten miły artykuł:

O'Hara RB, Kotze DJ (2010) Nie log-transform danych zliczania. Metody w ekologii i ewolucji 1: 118–122. klick .

W mojej dziedzinie badań (ekotoksykologia) mamy do czynienia ze słabo powielonymi eksperymentami, a GLM nie są szeroko stosowane. Zrobiłem więc podobną symulację jak O'Hara i Kotze (2010), ale naśladowałem dane ekotoksykologiczne.

Symulacje mocy :

Symulowałem dane z modelu czynnikowego z jedną grupą kontrolną ( ) i 5 grupami leczenia ( ). Obfitość w leczeniu 1 była identyczna z kontrolą ( ), obfitość w leczeniu 2-5 stanowiła połowę obfitości w kontroli ( ). Do symulacji zróżnicowałem wielkość próby (3,6,9,12) i liczebność w grupie kontrolnej (2, 4, 8, ..., 1024). Obficie pobrano z ujemnych rozkładów dwumianowych ze stałym parametrem dyspersji ( ). Wygenerowano i przeanalizowano 100 zestawów danych, stosując ujemny dwumianowy GLM i gaussowskie GLM + dane przekształcone logarytmicznie.μ 1 - 5 μ 1 = μ C μ 2 - 5 = 0,5 μ $\mu_c$$\mu_{1-5}$$\mu_1 = \mu_c$$\mu_{2-5} = 0.5 \mu_c$$\theta = 3.91$

Rezultaty są zgodne z oczekiwaniami: GLM ma większą moc, szczególnie gdy nie pobrano wielu zwierząt. Kod jest tutaj.

Błąd typu I :

Następnie spojrzałem na błąd pierwszego typu. Symulacje przeprowadzono jak wyżej, jednak wszystkie grupy miały taką samą liczebność ( ).$\mu_c = \mu_{1-5}$

Jednak wyniki nie są zgodne z oczekiwaniami: ujemny dwumianowy GLM wykazał większy błąd typu I w porównaniu do transformacji LM +. Zgodnie z oczekiwaniami różnica zniknęła wraz ze wzrostem wielkości próby. Kod jest tutaj.

Pytanie:

Dlaczego występuje zwiększony błąd typu I w porównaniu do transformacji lm +?

Jeśli mamy słabe dane (mała wielkość próby, mała liczebność (wiele zer)), czy powinniśmy zastosować transformację lm +? Małe wielkości próbek (2-4 na zabieg) są typowe dla takich eksperymentów i nie można ich łatwo zwiększyć.

Chociaż neg. kosz. GLM można uzasadnić jako odpowiednie dla tych danych, transformacja lm + może uchronić nas przed błędami typu 1.

To niezwykle interesujący problem. Po sprawdzeniu kodu nie mogę znaleźć oczywistej literówki.

Chciałbym, abyś powtórzył tę symulację, ale skorzystaj z testu maksymalnego prawdopodobieństwa, aby wyciągnąć wnioski na temat niejednorodności między grupami. Wymagałoby to dopasowania modelu zerowego, aby można było uzyskać oszacowania pod hipotezą zerową jednorodności stawek między grupami. Myślę, że jest to konieczne, ponieważ ujemny model dwumianowy nie jest modelem liniowym (szybkość jest parametryzowana liniowo, ale nie są). Dlatego nie jestem przekonany, że argumentacja zapewnia prawidłowe wnioskowanie.$\theta$$\theta$drop1

Większość testów dla modeli liniowych nie wymaga ponownego obliczania modelu w oparciu o hipotezę zerową. Wynika to z tego, że można obliczyć nachylenie geometryczne (test punktowy) i przybliżać szerokość (test Walda) przy użyciu oszacowań parametrów i szacowanej kowariancji przy samej hipotezie alternatywnej.

Ponieważ ujemny dwumian nie jest liniowy, myślę, że będziesz musiał dopasować model zerowy.

EDYTOWAĆ:

Zredagowałem kod i otrzymałem:

Edytowany kod tutaj: https://github.com/aomidpanah/simulations/blob/master/negativeBinomialML.r

Artykuł O'Hara i Kotze (Methods in Ecology and Evolution 1: 118–122) nie jest dobrym punktem wyjścia do dyskusji. Moje największe zaniepokojenie budzi twierdzenie zawarte w punkcie 4 streszczenia:

Stwierdziliśmy, że transformacje przebiegały słabo, z wyjątkiem. . .. Modele quasi-Poissona i ujemne dwumianowe ... [pokazały] małe odchylenie.

$\lambda$$\theta$$\lambda$

$\lambda$

Poniższy kod R ilustruje ten punkt:

x <- rnbinom(10000, 0.5, mu=2)  
## NB: Above, this 'mu' was our lambda. Confusing, is'nt it?
log(mean(x+1))
[1] 1.09631
log(2+1)  ## Check that this is about right
[1] 1.098612

mean(log(x+1))
[1] 0.7317908

Lub spróbuj

log(mean(x+.5))
[1] 0.9135269
mean(log(x+.5))
[1] 0.3270837

Skala szacowania parametrów ma ogromne znaczenie!

$\lambda$

Zauważ, że standardowa diagnostyka działa lepiej na skali log (x + c). Wybór c może nie mieć większego znaczenia; często 0,5 lub 1,0 ma sens. Jest to także lepszy punkt wyjścia do badania transformacji Box-Cox lub wariant Box-Cox Yeo-Johnsona. [Yeo, I. and Johnson, R. (2000)]. Zobacz także stronę pomocy dla powerTransform () w pakiecie samochodowym R. Pakiet gamlss R umożliwia dopasowanie ujemnych dwumianowych typów I (wspólna odmiana) lub II lub innych rozkładów, które modelują dyspersję, a także średnią, z łączami transformacji mocy równymi 0 (= log, tj. Log log) lub więcej . Pasowania nie zawsze mogą się zbiegać.

Przykład: Dane dotyczące zgonów w porównaniu z podstawowymi uszkodzeniami dotyczą nazwanych huraganów atlantyckich, które dotarły na kontynent USA. Dane są dostępne (nazwa hurricNamed ) z najnowszej wersji pakietu DAAG dla R. Strona pomocy dla danych zawiera szczegółowe informacje.

Wykres porównuje dopasowaną linię uzyskaną przy użyciu solidnego dopasowania modelu liniowego, z krzywą uzyskaną przez przekształcenie ujemnego dopasowania dwumianowego z łącznikiem log na skalę log (liczenie + 1) zastosowaną dla osi y na wykresie. (Zauważ, że trzeba użyć czegoś podobnego do skali log (liczenie + c), z dodatnim c, aby pokazać punkty i dopasowaną „linię” z ujemnego dopasowania dwumianowego na tym samym wykresie.) Zwróć uwagę na duże odchylenie, które jest widoczne dla ujemnego dopasowania dwumianowego na skali logarytmicznej. Solidne dopasowanie modelu liniowego jest znacznie mniej tendencyjne na tej skali, jeśli przyjmie się ujemny rozkład dwumianowy dla zliczeń. Dopasowanie modelu liniowego byłoby bezstronne przy założeniach klasycznej teorii normalnej. Zdecydowanie zaskoczyło mnie to, kiedy po raz pierwszy stworzyłem coś, co było zasadniczo powyższym wykresem! Krzywa lepiej pasowałaby do danych, ale różnica mieści się w granicach zwykłych standardów zmienności statystycznej. Solidne dopasowanie do modelu liniowego źle sprawdza się w przypadku liczników na dolnym końcu skali.

Uwaga --- Badania z danymi RNA-Seq: Porównanie dwóch stylów modelu było interesujące dla analizy danych zliczeniowych z eksperymentów ekspresji genów. Poniższy artykuł porównuje użycie solidnego modelu liniowego, pracującego z log (liczba + 1), z wykorzystaniem ujemnych dopasowań dwumianowych (jak w EdgeR pakietu Bioconductor ). Większość zliczeń, w aplikacji RNA-Seq, o której przede wszystkim chodzi, jest wystarczająco duża, aby odpowiednio zważone pasujące logarytmiczne dopasowanie modelu działało wyjątkowo dobrze.

Law, CW, Chen, Y, Shi, W, Smyth, GK (2014). Voom: odważniki precyzyjne odblokowują narzędzia do analizy modeli liniowych dla zliczeń odczytu sekwencji RNA. Genome Biology 15, R29. http://genomebiology.com/2014/15/2/R29

NB także najnowszy artykuł:

Schurch NJ, Schofield P, Gierliński M, Cole C, Sherstnev A, Singh V, Wrobel N, Gharbi K, Simpson GG, Owen-Hughes T, Blaxter M, Barton GJ (2016). Ile kopii biologicznych jest potrzebnych w eksperymencie z sekwencją RNA i jakiego narzędzia ekspresji różnicowej należy użyć? RNA http://www.rnajournal.org/cgi/doi/10.1261/rna.053959.115

Interesujące jest to, że model liniowy ataki wykorzystujące limma pakiet (jak Edger z grupy WEHI) wstać bardzo dobrze (w sensie pokazując niewiele dowodów stronniczości) w porównaniu do wyników z wielu powtórzeń, jak liczba powtórzeń jest zredukowany.

Kod R dla powyższego wykresu:

library(latticeExtra, quietly=TRUE)
hurricNamed <- DAAG::hurricNamed
ytxt <- c(0, 1, 3, 10, 30, 100, 300, 1000)
xtxt <- c(1,10, 100, 1000, 10000, 100000, 1000000 )
funy <- function(y)log(y+1)
gph <- xyplot(funy(deaths) ~ log(BaseDam2014), groups= mf, data=hurricNamed,
   scales=list(y=list(at=funy(ytxt), labels=paste(ytxt)),
           x=list(at=log(xtxt), labels=paste(xtxt))),
   xlab = "Base Damage (millions of 2014 US$); log transformed scale",
   ylab="Deaths; log transformed; offset=1",
   auto.key=list(columns=2),
   par.settings=simpleTheme(col=c("red","blue"), pch=16))
gph2 <- gph + layer(panel.text(x[c(13,84)], y[c(13,84)],
           labels=hurricNamed[c(13,84), "Name"], pos=3,
           col="gray30", cex=0.8),
        panel.text(x[c(13,84)], y[c(13,84)],
           labels=hurricNamed[c(13,84), "Year"], pos=1, 
           col="gray30", cex=0.8))
ab <- coef(MASS::rlm(funy(deaths) ~ log(BaseDam2014), data=hurricNamed))

gph3 <- gph2+layer(panel.abline(ab[1], b=ab[2], col="gray30", alpha=0.4))
## 100 points that are evenly spread on a log(BaseDam2014) scale
x <- with(hurricNamed, pretty(log(BaseDam2014),100))
df <- data.frame(BaseDam2014=exp(x[x>0])) 
hurr.nb <- MASS::glm.nb(deaths~log(BaseDam2014), data=hurricNamed[-c(13,84),])
df[,'hatnb'] <- funy(predict(hurr.nb, newdata=df, type='response'))
gph3 + latticeExtra::layer(data=df,
       panel.lines(log(BaseDam2014), hatnb, lwd=2, lty=2, 
           alpha=0.5, col="gray30"))    

Oryginalny post odzwierciedla artykuł Tony'ego Ivesa: Ives (2015) . Oczywiste jest, że testowanie istotności daje różne wyniki w szacowaniu parametrów.

John Maindonald wyjaśnia, dlaczego szacunki są stronnicze, ale jego ignorancja w tle jest denerwująca - krytykuje nas za wykazanie, że metoda, którą wszyscy zgadzamy się, jest wadliwa. Wielu ekologów ślepo rejestruje transformację, a my próbowaliśmy wskazać problemy z tym związane.

Tutaj jest bardziej szczegółowa dyskusja: Warton (2016)

Ives, AR (2015), Aby przetestować istotność współczynników regresji, śmiało i zlicz dane transformacji logarytmicznej. Methods Ecol Evol, 6: 828–835. doi: 10.1111 / 2041-210X.12386

Warton, DI, Lyons, M., Stoklosa, J. and Ives, AR (2016), Trzy punkty do rozważenia przy wyborze testu LM lub GLM dla danych zliczania. Metody Ecol Evol. doi: 10.1111 / 2041-210X.12552