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
Os resultados são os esperados: o GLM tem maior poder, principalmente quando não foram amostrados muitos animais. 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 ( ).
No entanto, os resultados não são os esperados: 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.
Respostas:
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:
Código editado aqui: https://github.com/aomidpanah/simulations/blob/master/negativeBinomialML.r
fonte
drop1()
faz internamente re-ajustar o modelo nulo ...glm.nb
drop1
logLik
getS3method('logLik', 'negbin'
drop1()
elrtest()
. Você está certo,drop1.glm
usosglm.fit
que dão o desvio errado. Não estava ciente de que não podemos usardrop1()
comglm.nb()
!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:
O código R a seguir ilustra o ponto:
Ou tente
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.
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.
NB também o artigo recente:
É 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:
fonte
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
fonte