GLM binomial negativo vs. transformação de log para dados de contagem: maior taxa de erro do tipo I

18

Alguns de vocês podem ter lido este belo artigo:

O'Hara RB, Kotze DJ (2010) Não transforme dados de contagem de transformações. Métodos em Ecologia e Evolução 1: 118–122. Klick .

No meu campo de pesquisa (ecotoxicologia), estamos lidando com experimentos mal replicados e os GLMs não são amplamente utilizados. Então, fiz uma simulação semelhante à de O'Hara e Kotze (2010), mas imitei dados ecotoxicológicos.

Simulações de potência :

dados de um planejamento fatorial com um grupo controle ( ) e 5 grupos de tratamento ( ). A abundância no tratamento 1 foi idêntica ao controle ( ), as abundâncias nos tratamentos 2-5 foram metade da abundância no controle ( ). Para as simulações, variei o tamanho da amostra (3,6,9,12) e a abundância no grupo controle (2, 4, 8, ..., 1024). Abundâncias foram obtidas de distribuições binomiais negativas com parâmetro de dispersão fixo ( ). 100 conjuntos de dados foram gerados e analisados ​​usando um GLM binomial negativo e um dado gaussiano transformado em log GLM +.μ 1 - 5 μ 1 = μ c μ 2 - 5 = 0,5 μ c θ = 3,91μcμ15μ1=μcμ25=0.5μcθ=3.91

Os resultados são os esperados: o GLM tem maior poder, principalmente quando não foram amostrados muitos animais. insira a descrição da imagem aqui O código está aqui.

Erro tipo I :

Em seguida, observei o erro do tipo um. As simulações foram feitas como acima, porém todos os grupos tiveram a mesma abundância ( ).μc=μ15

No entanto, os resultados não são os esperados: insira a descrição da imagem aqui GLM binomial negativo mostrou um erro maior do Tipo I em comparação com a transformação LM +. Como esperado, a diferença desapareceu com o aumento do tamanho da amostra. O código está aqui.

Questão:

Por que há um aumento de erro tipo I comparado à transformação lm +?

Se tivermos dados insuficientes (tamanho pequeno da amostra, baixa abundância (muitos zeros)), devemos usar a transformação lm +? Amostras pequenas (2-4 por tratamento) são típicas para essas experiências e não podem ser aumentadas facilmente.

Embora o negativo. bin. O GLM pode ser justificado como apropriado para esses dados; a transformação lm + pode nos impedir de erros do tipo 1.

EDi
fonte
1
Não é uma resposta para sua pergunta principal, mas algo para os leitores observarem: a menos que você cometa o erro real do tipo I equivalente para os dois procedimentos, comparar poder não faz sentido; Eu sempre posso aumentar a potência da menor (neste caso, pegar registros e ajustar o normal) levantando seu erro do tipo I. Por outro lado, se você especificar a situação específica (tamanho da amostra, abundância), poderá obter a taxa de erro do tipo I (por exemplo, por simulação) e, então, determinar em que taxa nominal testar para atingir a taxa de erro do tipo I desejada , para que seu poder se torne comparável.
Glen_b -Reinstala Monica 9/09
Os valores do eixo y em seus gráficos estão em média entre os 100 conjuntos de dados?
shadowtalker
Devo esclarecer meu comentário: no caso em que as estatísticas são inerentemente discretas, você não tem o controle perfeito da taxa de erro tipo I, mas geralmente é possível aproximar bastante as taxas de erro tipo I. Em situações em que você não pode aproximá-los o suficiente para serem comparáveis, a única maneira de torná-los comparáveis ​​é com testes aleatórios.
Glen_b -Reinstate Monica
@ssdecontrol: Não, é apenas a proporção de conjuntos de dados (dos 100) em que p <α
EDi
1
Existem duas questões: (i) é que as aproximações são assintóticas, mas não é infinito; portanto, a aproximação é apenas isso, uma aproximação - isso seria uma questão de haver ou não discrição e levaria a um nível de significância diferente do nominal (mas se for contínuo, é algo que você pode ajustar); (ii) existe a questão da discrição, que impede que você obtenha um nível exato de significância se você o ajustar. n
Glen_b -Reinstate Monica

Respostas:

17

Este é um problema extremamente interessante. Analisei seu código e não consigo encontrar erros de digitação imediatamente óbvios.

Eu gostaria que você refizesse esta simulação, mas use o teste de máxima verossimilhança para inferir sobre a heterogeneidade entre os grupos. Isso envolveria a remontagem de um modelo nulo para que você possa obter estimativas dos s sob a hipótese nula de homogeneidade nas taxas entre os grupos. Eu acho que isso é necessário porque o modelo binomial negativo não é um modelo linear (a taxa é parametrizada linearmente, mas os s não são). Portanto, não estou convencido de que o argumento fornece inferência correta.θθθdrop1

A maioria dos testes para modelos lineares não requer que você recompute o modelo sob a hipótese nula. Isso ocorre porque você pode calcular a inclinação geométrica (teste de pontuação) e aproximar a largura (teste de Wald) usando estimativas de parâmetros e covariância estimada apenas sob a hipótese alternativa.

Como o binômio negativo não é linear, acho que você precisará ajustar um modelo nulo.

EDITAR:

Editei o código e obtive o seguinte: insira a descrição da imagem aqui

Código editado aqui: https://github.com/aomidpanah/simulations/blob/master/negativeBinomialML.r

AdamO
fonte
Mas eu acho que drop1() faz internamente re-ajustar o modelo nulo ...
Ben Bolker
4
glm.nbθdrop1logLikgetS3method('logLik', 'negbin'
gostaria de marcar com +1 novamente, mas não posso. Agradável.
Ben Bolker
Obrigado! Eu apenas olhei o código de ambos drop1()e lrtest(). Você está certo, drop1.glmusos glm.fitque dão o desvio errado. Não estava ciente de que não podemos usar drop1()com glm.nb()!
Edi # 10/14
Portanto, a pontuação típica e os testes de Wald são inválidos no modelo binomial negativo?
shadowtalker
8

O artigo de O'Hara e Kotze (Methods in Ecology and Evolution 1: 118–122) não é um bom ponto de partida para discussão. Minha preocupação mais séria é a afirmação no ponto 4 do resumo:

Descobrimos que as transformações tiveram um desempenho ruim, exceto. . .. Os modelos quase-Poisson e binomial negativo ... [mostraram] pouco viés.

λθλ

λ

O código R a seguir ilustra o ponto:

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 tente

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

A escala na qual os parâmetros são estimados é muito importante!

λ

Observe que o diagnóstico padrão funciona melhor em uma escala de log (x + c). A escolha de c pode não importar muito; frequentemente 0,5 ou 1,0 fazem sentido. Também é um melhor ponto de partida para investigar as transformações de Box-Cox, ou a variante Yeo-Johnson de Box-Cox. [Yeo, I. e Johnson, R. (2000)]. Consulte mais a página de ajuda do powerTransform () na embalagem do carro de R. O pacote gamlss de R permite ajustar binomiais negativos tipos I (a variedade comum) ou II ou outras distribuições que modelam a dispersão, bem como a média, com links de transformação de potência de 0 (= log, ou seja, link de log) ou mais . Os ajustes nem sempre podem convergir.

Exemplo: Dados de Mortes x Dano Base são para furacões no Atlântico que chegaram ao continente americano. Os dados estão disponíveis (nome hurricNamed ) em uma versão recente do pacote DAAG para R. A página de ajuda para os dados possui detalhes.

Ajuste robusto linear versus negativo binomial

O gráfico compara uma linha ajustada obtida usando um ajuste de modelo linear robusto, com a curva obtida pela transformação de um ajuste binomial negativo com link de log na escala de log (count + 1) usada para o eixo y no gráfico. (Observe que é necessário usar algo semelhante a uma escala de log (count + c), com c positivo, para mostrar os pontos e a "linha" ajustada do ajuste binomial negativo no mesmo gráfico.) Observe o grande viés que é evidente para o ajuste binomial negativo na escala logarítmica. O ajuste robusto do modelo linear é muito menos tendencioso nessa escala, se alguém assumir uma distribuição binomial negativa para as contagens. Um ajuste de modelo linear seria imparcial sob as suposições clássicas da teoria normal. Achei o viés surpreendente quando criei o que era essencialmente o gráfico acima! Uma curva ajustaria melhor os dados, mas a diferença está dentro dos limites usuais de variabilidade estatística. O ajuste robusto do modelo linear faz um trabalho ruim para contagens na extremidade inferior da balança.

Nota --- Estudos com dados de RNA-Seq: A comparação dos dois estilos de modelo tem sido interessante para a análise de dados de contagem de experimentos de expressão gênica. O artigo a seguir compara o uso de um modelo linear robusto, trabalhando com log (contagem + 1), com o uso de ajustes binomiais negativos (como no pacote Bioconductor edgeR ). A maioria das contagens, no aplicativo RNA-Seq, principalmente em mente, é grande o suficiente para que o modelo log-linear adequadamente ponderado se encaixe e trabalhe extremamente bem.

Law, CW, Chen, Y, Shi, W, Smyth, GK (2014). Voom: pesos de precisão desbloqueiam ferramentas de análise de modelos lineares para contagens de leitura de RNA-seq. Genome Biology 15, R29. http://genomebiology.com/2014/15/2/R29

NB também o artigo recente:

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). Quantas réplicas biológicas são necessárias em um experimento de RNA-seq e qual ferramenta de expressão diferencial você deve usar? RNA http://www.rnajournal.org/cgi/doi/10.1261/rna.053959.115

É interessante que o modelo linear se ajuste usando o pacote limma (como edgeR , do grupo WEHI) se mantenha extremamente bem (no sentido de mostrar pouca evidência de viés), em relação aos resultados com muitas repetições, pois o número de repetições é reduzido.

Código R para o gráfico acima:

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
fonte
2
Obrigado pelo seu comentário, Sr. Maindonald. Nos últimos dois anos, também houve mais alguns trabalhos (focando mais no teste de hipóteses, em seguida, no viés): Ives 2015, Warton et al 2016, Szöcs 2015.
EDi
talvez seja um bom ponto de partida para discussão, mesmo que esse ponto em particular seja problemático? (Eu argumentaria de maneira mais geral que essa é uma razão para não focar demais no viés, mas considerar algo como o RMSE ... the Warton paper ...]
Ben Bolker
1
O argumento de Warton et al (2016), de que as propriedades dos dados devem ser o fundamento da escolha, é crucial. Os gráficos quantil-quantil são uma boa maneira de comparar os detalhes dos ajustes. Em particular, o ajuste em um ou outro ou em ambos os extremos pode ser importante para algumas aplicações. Modelos inflados a zero ou com obstáculos podem ser um aprimoramento eficaz para acertar a contagem zero. No extremo superior, qualquer um dos modelos em discussão pode estar seriamente comprometido. Warton et al., Louvávelmente, têm um exemplo. Eu gostaria de ver comparações em uma ampla variedade de conjuntos de dados ecológicos.
John Maindonald
Mas, nos conjuntos de dados ecológicos, as espécies da parte inferior (= espécies raras) não são interessantes? Não deve ser muito difícil de compilar alguns conjuntos de dados ecológicos e comparar ...
EDi
Na verdade, é para o extremo inferior da categoria de dano que o modelo binomial negativo parece, para os dados de mortes por furacões, ser menos satisfatório. De R gamlss pacote tem uma função que torna mais fácil para comparar percentis da distribuição equipada com percentis dos dados:
John Maindonald
6

O post original reflete o artigo de Tony Ives: Ives (2015) . É claro que o teste de significância fornece resultados diferentes para a estimativa de parâmetros.

John Maindonald explica por que as estimativas são tendenciosas, mas sua ignorância sobre o cenário é irritante - ele nos critica por mostrar que um método que todos nós concordamos que é falho é falho. Muitos ecologistas registram transformações às cegas, e estávamos tentando apontar os problemas para fazer isso.

Há uma discussão mais sutil aqui: Warton (2016)

Ives, AR (2015), Para testar a significância dos coeficientes de regressão, vá em frente e registre os dados da contagem de transformação. Métodos Ecol Evol, 6: 828–835. doi: 10.1111 / 2041-210X.12386

Warton, DI, Lyons, M., Stoklosa, J. e Ives, AR (2016), três pontos a serem considerados ao escolher um teste LM ou GLM para dados de contagem. Métodos Ecol Evol. doi: 10.1111 / 2041-210X.12552

Bob O'Hara
fonte
Bem-vindo ao CV. Embora útil, essa resposta é principalmente da resposta do tipo "somente link". Os links mudam e desassociam. Seria mais útil para o CV se você elaborasse os pontos principais de cada um.
Mike Hunter
Obrigado pela resposta. Eu acho que o artigo de Warton et al. inventa o estado atual da discussão.
EDi
Obrigado e bem-vindo! Tomei a liberdade de adicionar as referências na íntegra.
Scortchi - Restabelece Monica
1
Descreva os principais pontos mencionados nas novas referências e para onde faz sentido, relacione-os novamente com a pergunta original. Essa é uma contribuição valiosa, mas atualmente está mais próxima de um comentário em outra resposta do que uma resposta à pergunta (que deve fornecer contexto para os links , por exemplo). Algumas frases de contexto adicionais ajudariam o post substancialmente.
Glen_b -Reinstate Monica
3
Especificamente, meus comentários abordam o ponto 4 no artigo de O'Hara e Kotze: "Descobrimos que as transformações tiveram um desempenho fraco, exceto ... Os modelos quase-Poisson e binomial negativo ... [mostraram] pouco viés". As simulações são um comentário sobre a comparação entre a média esperada em uma escala de y (as contagens) versus a média esperada em uma escala de log (y + c), para uma distribuição altamente inclinada positivamente, nada mais. O parâmetro binomial negativo lambda é imparcial na escala de y, enquanto a média log-normal é imparcial (abaixo da normalidade nessa escala) em uma escala de log (y + c).
John Maindonald