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.
λθλ
λ
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!
λ
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"))
Kod jest tutaj.
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.