GLM binomial négatif ou transformation logarithmique pour les données de comptage: augmentation du taux d'erreur de type I

18

Certains d'entre vous ont peut-être lu ce bel article:

O'Hara RB, Kotze DJ (2010) Ne pas enregistrer les données de comptage par transformation. Méthodes en écologie et évolution 1: 118-122. Klick .

Dans mon domaine de recherche (écotoxicologie), nous avons affaire à des expériences mal reproduites et les GLM ne sont pas largement utilisés. J'ai donc fait une simulation similaire à O'Hara & Kotze (2010), mais j'ai imité les données écotoxicologiques.

Simulations de puissance :

J'ai simulé des données d'un plan factoriel avec un groupe témoin ( ) et 5 groupes de traitement ( ). L'abondance dans le traitement 1 était identique au témoin ( ), les abondances dans les traitements 2-5 étaient la moitié de l'abondance dans le contrôle ( ). Pour les simulations, j'ai fait varier la taille de l'échantillon (3,6,9,12) et l'abondance dans le groupe témoin (2, 4, 8, ..., 1024). Les abondances ont été tirées d'une distribution binomiale négative avec un paramètre de dispersion fixe ( ). 100 jeux de données ont été générés et analysés à l'aide d'un GLM binomial négatif et d'un GLM gaussien + données transformées en logarithme.μ 1 - 5 μ 1 = μ c μ 2 - 5 = 0,5 μ c θ = 3,91μcμ15μ1=μcμ25=0.5μcθ=3.91

Les résultats sont comme prévu: le GLM a une plus grande puissance, surtout lorsque peu d'animaux ont été échantillonnés. entrez la description de l'image ici Le code est ici.

Erreur de type I :

Ensuite, j'ai examiné l'erreur de type un. Les simulations ont été effectuées comme ci-dessus, mais tous les groupes avaient la même abondance ( ).μc=μ15

Cependant, les résultats ne sont pas ceux attendus: le entrez la description de l'image ici GLM binomial négatif a montré une plus grande erreur de type I par rapport à la transformation LM +. Comme prévu, la différence a disparu avec l'augmentation de la taille de l'échantillon. Le code est ici.

Question:

Pourquoi y a-t-il une erreur de type I accrue par rapport à la transformation lm +?

Si nous avons des données médiocres (petite taille d'échantillon, faible abondance (plusieurs zéros)), devrions-nous alors utiliser la transformation lm +? De petites tailles d'échantillon (2 à 4 par traitement) sont typiques de telles expériences et ne peuvent pas être augmentées facilement.

Bien que le neg. poubelle. GLM peut être justifié comme étant approprié pour ces données, la transformation lm + peut nous empêcher des erreurs de type 1.

EDi
la source
1
Pas une réponse à votre question principale, mais quelque chose que les lecteurs doivent noter: à moins que vous ne rendiez l'erreur de type I équivalente pour les deux procédures, comparer la puissance n'a pas de sens; Je peux toujours augmenter la puissance de la plus basse (dans ce cas, les journaux de prise et la taille normale) en levant son erreur de type I. D'un autre côté, si vous spécifiez la situation particulière (taille de l'échantillon, abondance), vous pouvez obtenir le taux d'erreur de type I (par exemple par simulation), et ainsi déterminer le taux nominal à tester pour atteindre le taux d'erreur de type I souhaité , de sorte que leur pouvoir devient comparable.
Glen_b -Reinstate Monica
Les valeurs de l'axe des y dans vos tracés sont-elles moyennes sur les 100 ensembles de données?
shadowtalker
Je devrais clarifier mon commentaire: dans le cas où les statistiques sont intrinsèquement discrètes, vous n'avez pas un contrôle parfait du taux d'erreur de type I, mais vous pouvez généralement rendre les taux d'erreur de type I assez proches. Dans les situations où vous ne pouvez pas les rapprocher suffisamment pour être comparables, la seule façon de les rendre comparables est d'utiliser des tests randomisés.
Glen_b -Reinstate Monica
@ssdecontrol: Non, c'est juste la proportion d'ensembles de données (sur 100) où p <α
EDi
1
Il y a deux problèmes: (i) est que les approximations sont asymptotiques mais n'est pas l'infini, donc l'approximation n'est que cela, une approximation - ce serait un problème s'il y avait de la discrétion ou non, et conduirait à un niveau de signification autre que la valeur nominale (mais si elle est continue, c'est quelque chose que vous pouvez ajuster); (ii) il y a le problème de la discrétion, qui vous empêche d'obtenir un niveau de signification exact si vous vous ajustez pour cela. n
Glen_b -Reinstate Monica

Réponses:

17

C'est un problème extrêmement intéressant. J'ai examiné votre code et je ne trouve aucune faute de frappe immédiatement évidente.

Je voudrais vous voir refaire cette simulation mais utiliser le test du maximum de vraisemblance pour faire l'inférence sur l'hétérogénéité entre les groupes. Cela impliquerait de réaménager un modèle nul afin que vous puissiez obtenir des estimations des sous l'hypothèse nulle d'homogénéité des taux entre les groupes. Je pense que cela est nécessaire car le modèle binomial négatif n'est pas un modèle linéaire (le taux est paramétré linéairement, mais les ne le sont pas). Par conséquent, je ne suis pas convaincu que l' argument fournit une inférence correcte.θθθdrop1

La plupart des tests pour les modèles linéaires ne nécessitent pas de recalculer le modèle sous l'hypothèse nulle. En effet, vous pouvez calculer la pente géométrique (test de score) et approximer la largeur (test de Wald) en utilisant des estimations de paramètres et une covariance estimée sous l'hypothèse alternative seule.

Comme le binôme négatif n'est pas linéaire, je pense que vous devrez ajuster un modèle nul.

ÉDITER:

J'ai édité le code et obtenu ce qui suit: entrez la description de l'image ici

Code édité ici: https://github.com/aomidpanah/simulations/blob/master/negativeBinomialML.r

AdamO
la source
Mais je pense que cela drop1() correspond en interne au modèle nul ...
Ben Bolker
4
Nan! Regardez à quel point l'ajusteur negbinom est compliqué dans le glm.nbcode, les sont estimés à l'aide de l'algorithme EM. n'a pas de méthode par défaut pour le binôme négatif. Cependant, la fonction le fait (voir ). J'ai corrigé le code de l'OP et montré que cela donne une inférence presque identique au modèle de Poisson sous le zéro (qui a toujours un problème d'étalonnage important en raison d'autres problèmes). θdrop1logLikgetS3method('logLik', 'negbin'
AdamO
aimerait +1 à nouveau mais je ne peux pas. Agréable.
Ben Bolker
Merci! J'ai juste regardé le code des deux drop1()et lrtest(). Vous avez raison, des drop1.glmutilisations glm.fitqui donnent la mauvaise déviance. Je ne savais pas que nous ne pouvions pas l'utiliser drop1()avec glm.nb()!
EDi
Donc, le score typique et les tests de Wald ne sont pas valides dans le modèle binomial négatif?
shadowtalker
8

L'article O'Hara et Kotze (Methods in Ecology and Evolution 1: 118–122) n'est pas un bon point de départ pour la discussion. Ma préoccupation la plus sérieuse est l'allégation au point 4 du résumé:

Nous avons constaté que les transformations se sont mal déroulées, sauf. . .. Les modèles binomiaux quasi-Poisson et négatifs ... [montraient] peu de biais.

La moyenne pour une distribution de Poisson ou une distribution binomiale négative est pour une distribution qui, pour les valeurs de <= 2 et pour la plage de valeurs de la moyenne qui a été étudiée, est très positivement asymétrique. Les moyennes des distributions normales ajustées sont sur une échelle de log (y + c) (c est le décalage) et d'estimation E (log (y + c)]. Cette distribution est beaucoup plus proche de symétrique que la distribution de y .θ λλθλ

Les simulations d'O'Hara et de Kotze comparent E (log (y + c)], estimé par la moyenne (log (y + c)), avec log (E [y + c]). Ils peuvent être, et dans les cas notés sont très différents. Leurs graphiques ne comparent pas un binôme négatif avec un ajustement log (y + c), mais comparent plutôt la moyenne (log (y + c)] avec log (E [y + c]). ) sur leurs graphiques, ce sont en fait les ajustements binomiaux négatifs qui sont les plus biaisés! λ

Le code R suivant illustre le point:

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

Ou essayez

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

L'échelle sur laquelle les paramètres sont estimés est très importante!

Si l'on échantillonne à partir d'un Poisson, bien sûr, on s'attend à ce que le Poisson fasse mieux, si l'on en juge par les critères utilisés pour ajuster le Poisson. Idem pour un binôme négatif. La différence n'est peut-être pas si grande, si la comparaison est juste. Les données réelles (par exemple, peut-être, dans certains contextes génétiques) peuvent parfois être très proches de Poisson. Quand ils s'écartent de Poisson, le binôme négatif peut ou peut ne pas bien fonctionner. De même, surtout si est de l'ordre de peut-être 10 ou plus, pour modéliser log (y + 1) en utilisant la théorie normale standard.λ

Notez que les diagnostics standard fonctionnent mieux sur une échelle de journal (x + c). Le choix de c n'a peut-être pas trop d'importance; souvent 0,5 ou 1,0 ont du sens. C'est également un meilleur point de départ pour étudier les transformations de Box-Cox, ou la variante Yeo-Johnson de Box-Cox. [Yeo, I. et Johnson, R. (2000)]. Voir plus loin la page d'aide de powerTransform () dans le package de voiture de R. Le package gamlss de R permet d'adapter les types binomiaux négatifs I (la variété commune) ou II, ou d'autres distributions qui modélisent la dispersion ainsi que la moyenne, avec des liens de transformation de puissance de 0 (= log, c.-à-d. Lien log) ou plus . Les ajustements peuvent ne pas toujours converger.

Exemple: Les décès par rapport aux dommages de base concernent les ouragans de l'Atlantique nommés qui ont atteint le continent américain. Les données sont disponibles (nom hurricNamed ) à partir d'une version récente du package DAAG pour R. La page d'aide pour les données contient des détails.

Ajustement log-linéaire robuste vs ajustement binomial négatif

Le graphique compare une ligne ajustée obtenue à l'aide d'un ajustement de modèle linéaire robuste, avec la courbe obtenue en transformant un ajustement binomial négatif avec lien logarithmique sur l'échelle logarithmique (comptage + 1) utilisée pour l'axe des y sur le graphique. (Notez qu'il faut utiliser quelque chose qui ressemble à une échelle logarithmique (nombre + c), avec c positif, pour montrer les points et la "ligne" ajustée de l'ajustement binomial négatif sur le même graphique.) Notez le biais important qui est évident pour l'ajustement binomial négatif sur l'échelle logarithmique. L'ajustement robuste du modèle linéaire est beaucoup moins biaisé sur cette échelle, si l'on suppose une distribution binomiale négative pour les dénombrements. Un ajustement de modèle linéaire serait non biaisé selon les hypothèses de la théorie normale classique. J'ai trouvé le biais étonnant lorsque j'ai créé ce qui était essentiellement le graphique ci-dessus! Une courbe correspondrait mieux aux données, mais la différence se situe dans les limites des normes habituelles de variabilité statistique. L'ajustement robuste du modèle linéaire fait un mauvais travail pour les dénombrements à l'extrémité inférieure de l'échelle.

Remarque --- Études avec les données RNA-Seq: La comparaison des deux styles de modèle a été intéressante pour l'analyse des données de comptage des expériences d'expression génique. L'article suivant compare l'utilisation d'un modèle linéaire robuste, fonctionnant avec log (count + 1), avec l'utilisation d'ajustements binomiaux négatifs (comme dans le package R du boîtier du bioconducteur ). La plupart des dénombrements, dans l'application RNA-Seq qui est principalement à l'esprit, sont suffisamment grands pour que les ajustements log-linéaires correctement pesés fonctionnent extrêmement bien.

Law, CW, Chen, Y, Shi, W, Smyth, GK (2014). Voom: les poids de précision débloquent des outils d'analyse de modèle linéaire pour les comptages de lecture ARN-seq. Genome Biology 15, R29. http://genomebiology.com/2014/15/2/R29

NB également le récent article:

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). Combien de répétitions biologiques sont nécessaires dans une expérience d'ARN-seq et quel outil d'expression différentielle devez-vous utiliser? RNA http://www.rnajournal.org/cgi/doi/10.1261/rna.053959.115

Il est intéressant de noter que le modèle linéaire s'adapte en utilisant le package limma (comme edgeR , du groupe WEHI) résiste extrêmement bien (dans le sens de montrer peu de preuves de biais), par rapport aux résultats avec de nombreuses répétitions, car le nombre de répétitions est réduit.

Code R pour le graphique ci-dessus:

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"))    
John Maindonald
la source
2
Merci pour votre commentaire M. Maindonald. Au cours des deux dernières années, il y a eu aussi d'autres articles (se concentrant davantage sur les tests d'hypothèses, puis sur les biais): Ives 2015, Warton et al 2016, Szöcs 2015.
EDi
c'est peut-être un bon point de départ pour la discussion, même si ce point particulier est problématique? (Je dirais plus généralement que c'est une raison de ne pas trop se concentrer sur les préjugés, mais plutôt de considérer quelque chose comme RMSE ... [avertissement, je n'ai pas relu ces articles récemment, et je n'ai lu que le résumé de le journal Warton ...]
Ben Bolker
1
Le point de Warton et al (2016), selon lequel les propriétés des données devraient être les motifs de choix, est crucial. Les tracés quantile-quantile sont un bon moyen de comparer les détails des ajustements. En particulier, l'ajustement à l'un ou l'autre ou aux deux extrêmes peut être important pour certaines applications. Les modèles à gonflement nul ou à obstacle peuvent être un raffinement efficace pour obtenir le bon compte zéro. À l'extrême supérieur, n'importe lequel des modèles examinés peut être gravement compromis. Warton et al en ont, louable, un exemple. J'aimerais voir des comparaisons entre un large éventail d'ensembles de données écologiques.
John Maindonald
Mais dans les ensembles de données écologiques, les espèces de la partie inférieure (= espèces rares) ne sont-elles pas intéressantes? Il ne devrait pas être trop difficile de compiler des ensembles de données écologiques et de comparer ...
EDi
En fait, c'est pour le bas de la catégorie des dommages que le modèle binomial négatif semble, pour les données de décès par ouragan, être le moins satisfaisant. Le package gamlss de R a une fonction qui permet de comparer facilement les centiles de la distribution ajustée avec les centiles des données:
John Maindonald
6

Le message original reflète l'article de Tony Ives: Ives (2015) . Il est clair que les tests de signification donnent des résultats différents à l'estimation des paramètres.

John Maindonald explique pourquoi les estimations sont biaisées, mais son ignorance du contexte est ennuyeuse - il nous critique pour avoir montré qu'une méthode dont nous convenons tous qu'elle est défectueuse est défectueuse. De nombreux écologistes enregistrent aveuglément la transformation, et nous essayions de souligner les problèmes liés à cette opération.

Il y a une discussion plus nuancée ici: Warton (2016)

Ives, AR (2015), Pour tester la signification des coefficients de régression, allez-y et transformez les données de décompte. Méthodes Ecol Evol, 6: 828–835. doi: 10.1111 / 2041-210X.12386

Warton, DI, Lyon, M., Stoklosa, J. et Ives, AR (2016), Trois points à considérer lors du choix d'un test LM ou GLM pour les données de comptage. Méthodes Ecol Evol. doi: 10.1111 / 2041-210X.12552

Bob O'Hara
la source
Bienvenue sur CV. Bien qu'utile, cette réponse est principalement de type "lien uniquement". Les liens changent et se dissocient. Il serait plus utile pour le CV si vous deviez développer les points clés de chacun.
Mike Hunter
Merci pour la réponse. Je pense que l'article de Warton et al. inventorie l'état actuel de la discussion.
EDi
Merci, bienvenue! J'ai pris la liberté d'ajouter les références dans leur intégralité.
Scortchi - Réintégrer Monica
1
Veuillez décrire les principaux points soulevés dans les nouvelles références et, là où cela a du sens, reliez-les également à la question d'origine. Il s'agit d'une contribution précieuse, mais elle est plus proche actuellement d'un commentaire sur une autre réponse que d'une réponse à la question (qui devrait fournir un contexte pour les liens , par exemple). Quelques phrases de contexte supplémentaires aideraient considérablement le poste.
Glen_b -Reinstate Monica
3
Plus précisément, mes commentaires portent sur le point 4 de l'article d'O'Hara et de Kotze: "Nous avons constaté que les transformations se sont mal déroulées, sauf ... Les modèles binomiaux quasi-Poisson et négatifs ... [montraient] peu de biais." Les simulations sont un commentaire sur la comparaison entre la moyenne attendue sur une échelle de y (les dénombrements) et la moyenne attendue sur une échelle de log (y + c), pour une distribution asymétrique très positive, rien de plus. Le paramètre binomial négatif lambda est sans biais sur l'échelle de y, tandis que la moyenne log-normale est sans biais (sous normalité à cette échelle) sur une échelle de log (y + c).
John Maindonald