Como testar / provar que os dados são inflados com zero?

9

Eu tenho um problema que acho que deveria ser simples, mas não consigo entender direito. Eu estou olhando para a polinização das sementes, tenho plantas (n = 36) que florescem em cachos, provo 3 cachos de flores de cada planta e 6 vagens de cada cluster (18 vagens de sementes no total de cada planta). Uma vagem pode ter entre 0 e no máximo 4 sementes polinizadas. Portanto, os dados são contados, com um limite superior. Estou descobrindo que uma média de ~ 10% das sementes são polinizadas, mas entre 1 e 30% em uma determinada planta, portanto, com dados dispersos e, é claro, existem 4 réplicas de cluster ausentes em 3 plantas, portanto não são perfeitamente simétricas .

A pergunta que faço é se esses dados corroboram a ideia de que esta planta requer polinizadores para o conjunto de sementes.

Estou descobrindo que a distribuição do número de sementes em uma vagem parece que existem mais 0 vagens de sementes polinizadas (6-9 vagens de 16) e mais 3 e 4 vagens de sementes polinizadas (2-4 para cada) do que seria seria esperado se as sementes da população fossem polinizadas aleatoriamente. Basicamente, acho que este é um exemplo clássico de dados inflacionados zero, primeiro um inseto visita ou não a flor (um gerador zero) e, se o fizer, poliniza 0-4 das sementes em outra distribuição. A hipótese alternativa é que a planta é parcialmente egoísta, e seria de esperar que cada semente tivesse a mesma probabilidade de ser polinizada (esses dados sugerem uma chance de aproximadamente 0,1, o que significa 0,01 de chance para duas sementes na mesma vagem, etc.) .

Mas eu simplesmente quero demonstrar que os dados se encaixam melhor em uma ou outra distribuição, na verdade, NÃO FAZ UM ZIP ou ZINB nos dados. Eu acho que qualquer método que eu use deve levar em conta o número real de sementes polinizadas e o número de vagens amostradas em cada planta. A melhor coisa que inventei é fazer algum tipo de tira de bota, onde apenas aleatoriamente atribuo o número de sementes polinizadas para uma determinada planta ao número de vagens de amostras que eu amostramos, faça isso 10.000 vezes e veja qual é a probabilidade os dados experimentais para a planta em questão saíram dessa distribuição aleatória.

Eu apenas sinto que há algo sobre isso que deve ser muito mais fácil do que o bootstrap de força bruta, mas depois de dias de pensamento e pesquisa, estou desistindo. Não posso apenas comparar com uma distribuição Poisson porque é o limite superior, não é binomial porque preciso gerar a distribuição esperada de alguma maneira primeiro. Alguma ideia? E eu estou usando R, então conselhos (especialmente como gerar mais de 10.000 distribuições aleatórias de n bolas em 16 caixas que podem conter no máximo 4 bolas) seriam bem-vindos.

ADICIONADO 9/07/2012 Primeiro, obrigado a todos por todo o interesse e ajuda. Ler as respostas me fez pensar em reformular um pouco minha pergunta. O que estou dizendo é que tenho uma hipótese (que por enquanto estou pensando em nula) de que as sementes são polinizadas aleatoriamente entre as vagens, e minha hipótese alternativa é que uma vagem de sementes com pelo menos 1 semente polinizada tem mais probabilidade de possuem múltiplas sementes polinizadas do que seria esperado por um processo aleatório. Forneci dados reais de três plantas como exemplos para ilustrar o que estou falando. A primeira coluna é o número de sementes polinizadas em uma vagem, a segunda coluna é a frequência de vagens com essa contagem de sementes.

planta 1 (total de 3 sementes: 4% de polinização)

num.seeds :: pod.freq

0 :: 16

1 :: 1

2 :: 1

3 :: 0

4 :: 0

planta 2 (total de 19 sementes: 26% de polinização)

num.seeds :: pod.freq

0 :: 12

1 :: 1

2 :: 1

3 :: 0

4 :: 4

planta 3 (total de 16 sementes: 22% de polinização)

num.seeds :: pod.freq

0 :: 9

1 :: 4

2 :: 3

3 :: 2

4 :: 0

Na planta 1, apenas 3 sementes foram polinizadas em 18 vagens, uma vagem tinha uma semente e uma vagem tinha duas sementes. Pensando em um processo de adicionar uma semente às vagens aleatoriamente, as duas primeiras sementes entram em sua própria vagem, mas para a terceira semente, existem 6 manchas disponíveis nas vagens que já possuem uma semente, mas 64 manchas nas 16 vagens sem sementes, portanto, a maior probabilidade de uma vagem com 2 sementes aqui é 6/64 = 0,094. Isso é um pouco baixo, mas não muito extremo, então eu diria que esta planta se encaixa na hipótese de polinização aleatória em todas as sementes com uma chance de ~ 4% de ocorrer polinização. Mas a planta 2 parece muito mais extrema para mim, com 4 vagens completamente polinizadas, mas 12 vagens sem nada. Não tenho muita certeza de como calcular diretamente as chances dessa distribuição (daí a minha ideia de inicialização), mas acho que as chances dessa distribuição ocorrerem aleatoriamente se cada semente tiver uma chance de ~ 25% de polinização ser bastante baixa. Planta # 3 Eu realmente não tenho idéia, acho que existem mais 0 e 3 do que se poderia esperar para uma distribuição aleatória, mas meu pressentimento é que essa distribuição para esse número de sementes é muito mais provável que a distribuição para a planta # 2, e pode não ser tão improvável. Mas, obviamente, quero ter certeza, e em todas as plantas. Eu acho que existem mais 0 e 3 do que se poderia esperar para uma distribuição aleatória, mas meu pressentimento é que essa distribuição para esse número de sementes é muito mais provável que a distribuição para a planta nº 2, e pode não ser tão improvável. Mas, obviamente, quero ter certeza, e em todas as plantas. Eu acho que existem mais 0 e 3 do que se poderia esperar para uma distribuição aleatória, mas meu pressentimento é que essa distribuição para esse número de sementes é muito mais provável que a distribuição para a planta nº 2, e pode não ser tão improvável. Mas, obviamente, quero ter certeza, e em todas as plantas.

No final, estou procurando escrever uma declaração como “A distribuição de sementes polinizadas em vagens se encaixa (ou não) com a hipótese de que as plantas não são apenas parcialmente autocompatíveis, mas requerem a visita de um polinizador para efetuar o conjunto de sementes. (resultados do teste estatístico). ” Isso é realmente apenas parte da minha seção prospectiva, onde estou falando sobre quais experimentos realizar a seguir, então não estou desesperado para que isso seja uma coisa ou outra, mas quero saber por mim mesmo, se possível. Se não posso fazer o que estou tentando fazer com esses dados, também gostaria de saber!

Fiz uma pergunta bastante ampla a princípio, pois estou curioso para saber se existem ou não bons testes para mostrar se os dados devem entrar em um modelo inflado zero em primeiro lugar. Todos os exemplos que eu vi parecem dizer - "veja, há muitos zeros aqui e há uma explicação razoável para isso, então vamos usar um modelo inflado zero". É o que estou fazendo agora neste fórum, mas tive uma experiência no meu último capítulo em que usei um Poisson glm para dados de contagem e um dos meus supervisores disse: “Não, os glms são muito complexos e desnecessários, esses dados devem entrar em uma tabela de contingência ”e depois me enviou um despejo de dados da enorme tabela de contingência gerada pelo seu caro pacote de estatísticas que forneceu os mesmos valores de p para todos os meus fatores + interações com três dígitos significativos !! Então, estou tentando manter as estatísticas claras e simples, e certifique-se de entendê-las bem o suficiente para defender minhas escolhas com firmeza, o que não acho que possa fazer por um modelo inflado com zero no momento. Usei um quase -IBMOMIAL (para plantas inteiras se livrar da pesudoreplicação) e um modelo misto para os dados acima para comparar tratamentos e responder às minhas principais perguntas experimentais, ou parece fazer o mesmo trabalho, mas também vou brinque com o ZINB hoje à noite, para ver como é o desempenho. Estou pensando que, se puder demonstrar explicitamente que esses dados estão fortemente agrupados (ou inflados com zero) a princípio, depois forneça uma boa razão biológica para isso, estarei muito melhor configurado para retirar um ZINB do que para basta comparar um com um modelo quasibinomial / misto e argumentar, uma vez que fornece melhores resultados, é o que devo usar. o que acho que não posso fazer por um modelo inflado zero agora. Usei um quase -IBMOMIAL (para plantas inteiras se livrar da pesudoreplicação) e um modelo misto para os dados acima para comparar tratamentos e responder às minhas principais perguntas experimentais, ou parece fazer o mesmo trabalho, mas também vou brinque com o ZINB hoje à noite, para ver como é o desempenho. Estou pensando que, se puder demonstrar explicitamente que esses dados estão fortemente agrupados (ou inflados com zero) a princípio, depois forneça uma boa razão biológica para isso, estarei muito melhor configurado para retirar um ZINB do que para basta comparar um com um modelo quasibinomial / misto e argumentar, uma vez que fornece melhores resultados, é o que devo usar. o que acho que não posso fazer por um modelo inflado zero agora. Usei um quase -IBMOMIAL (para plantas inteiras se livrar da pesudoreplicação) e um modelo misto para os dados acima para comparar tratamentos e responder às minhas principais perguntas experimentais, ou parece fazer o mesmo trabalho, mas também vou brinque com o ZINB hoje à noite, para ver como é o desempenho. Estou pensando que, se puder demonstrar explicitamente que esses dados estão fortemente agrupados (ou inflados com zero) a princípio, depois forneça uma boa razão biológica para isso, estarei muito melhor configurado para retirar um ZINB do que para basta comparar um com um modelo quasibinomial / misto e argumentar, uma vez que fornece melhores resultados, é o que devo usar. Usei um quase -IBMOMIAL (para plantas inteiras se livrar da pesudoreplicação) e um modelo misto para os dados acima para comparar tratamentos e responder às minhas principais perguntas experimentais, ou parece fazer o mesmo trabalho, mas também vou brinque com o ZINB hoje à noite, para ver como é o desempenho. Estou pensando que, se puder demonstrar explicitamente que esses dados estão fortemente agrupados (ou inflados com zero) a princípio, depois forneça uma boa razão biológica para isso, estarei muito melhor configurado para retirar um ZINB do que para basta comparar um com um modelo quasibinomial / misto e argumentar, uma vez que fornece melhores resultados, é o que devo usar. Usei um quase -IBMOMIAL (para plantas inteiras se livrar da pesudoreplicação) e um modelo misto para os dados acima para comparar tratamentos e responder às minhas principais perguntas experimentais, ou parece fazer o mesmo trabalho, mas também vou brinque com o ZINB hoje à noite, para ver como é o desempenho. Estou pensando que, se puder demonstrar explicitamente que esses dados estão fortemente agrupados (ou inflados com zero) a princípio, depois forneça uma boa razão biológica para isso, estarei muito melhor configurado para retirar um ZINB do que para basta comparar um com um modelo quasibinomial / misto e argumentar, uma vez que fornece melhores resultados, é o que devo usar.

Mas não quero me distrair muito da minha pergunta principal. Como posso determinar se meus dados realmente são mais inflados em zero do que o esperado em uma distribuição aleatória? No meu caso, a resposta para isso é o que é de real interesse para mim, com o possível benefício da justificação de modelo sendo um bônus.

Obrigado novamente por todo o seu tempo e ajuda!

Saúde, BWGIA

BWGIA
fonte
por que você não deseja ajustar o modelo binomial inflado a zero?
atiretoo - restabelece monica
a hipótese do "eu parcial" é exclusiva da hipótese do "polinizador"? Nesse caso, seu segundo modelo seria simplesmente um modelo binomial com probabilidade pe tamanho = 4.
atiretoo - restabelecer monica 07/07

Respostas:

5

Parece-me um modelo misto relativamente direto (não-linear) para mim. Você tem vagens de sementes aninhadas em clusters aninhados em plantas e pode ajustar um modelo binomial com efeitos aleatórios em cada estágio:

    library(lme4)
    binre <- lmer( pollinated ~ 1 + (1|plant) + (1|cluster), data = my.data, family = binomial)

ou com covariáveis, se você os tiver. Se as flores se polinizarem, você poderá observar alguns efeitos leves devido à variabilidade natural da viabilidade das plantas por si mesmas. No entanto, se a maior parte da variabilidade na resposta é direcionada pela variabilidade de cluster, por exemplo, você tem uma evidência mais forte de polinização por insetos que podem visitar apenas os clusters selecionados de uma planta. Idealmente, você desejaria uma distribuição não paramétrica dos efeitos aleatórios em vez de gaussiana: uma massa pontual em zero, sem visitas a insetos e uma massa pontual em um valor positivo - este é essencialmente o modelo de mistura que Michael Chernick pensou. Você pode ajustar isso com o pacote GLLAMM Stata. Eu ficaria surpreso se isso não fosse possível em R.

Provavelmente, para um experimento limpo, você gostaria de ter as plantas dentro, ou pelo menos em um local sem acesso a insetos, e ver quantas sementes seriam polinizadas. Isso provavelmente responderia a todas as suas perguntas de uma maneira mais metodologicamente rigorosa.

StasK
fonte
Vou tentar isso, acho que ajudará a responder minhas próprias perguntas, mas não tenho tanta certeza de como convencerá os outros. Você está bem na segunda parte, estou tentando pensar em como esses dados informam um experimento futuro mais direcionado.
BWGIA
1

Parece-me que esta é uma distribuição de mistura para cada inseto individual. Com probabilidade p, o inseto pousa com probabilidade 1-p pousa e distribui 0 a 4 sementes. Mas se você não tem informações sobre se o inseto pousa na planta ou não, não é possível distinguir as duas maneiras de obter 0. Portanto, você pode deixar p ser a probabilidade de 0 e ter a distribuição multinomial (p1, p2, p3, p4) onde pi é a probabilidade de i sementes, dado o inseto polinizar sujeito à restrição p1 + p2 + p3 + p4 = 1. O modelo possui cinco incógnitas p, p1, p2, p3, p4 com a restrição 0 = 0 para cada i. Com dados suficientes, você poderá estimar esses parâmetros, talvez usando uma abordagem de máxima verossimilhança restrita.

Michael R. Chernick
fonte
Concordo, mas a questão não é se encaixar nesse modelo, mas gerar distribuições previstas sob duas hipóteses biológicas diferentes. Talvez a resposta seja encaixar um ZIB e "algum outro modelo" que corresponda à hipótese do ego e compará-los.
atiretoo - restabelece monica
@atiretoo, o modelo não fornece uma distribuição estimada para o número de sementes polinizadas que você poderia comparar com sua distribuição hipotética?
22912 Michael Michael Chernick
Concordo - se você tem os modelos certos para as 2 hipóteses.
atiretoo - restabelecer monica 07/07
1

Esta é uma resposta para a última parte da sua pergunta, como gerar rapidamente os dados que você deseja para a hipótese do polinizador:

n = 16
max = 4
p1 = 0.1
p2 = 0.9
Y1 = rbinom(10000*n,1,p1)
Y2 = matrix(Y1*rbinom(10000*n,4,p2),ncol=16)

Você também pode usar rzibinom()no pacote VGAM. Embora eu não tenha certeza do que você quer fazer com isso. Você tem 2 parâmetros livres, p1 e p2, que precisam ser estimados. Por que não usar um modelo binomial inflado com zero para estimar a partir dos dados?

Você deve procurar o pacote VGAM, que se encaixa nos modelos ZIB, entre outros. De fato, você pode obter a distribuição esperada para um ZIB a partir da função VGAM dzibinom(), com a qual você poderia comparar sua distribuição observada se conhecer os parâmetros de visitação e polinização. Novamente, você realmente deve se encaixar no modelo ZIB.

Se a sua hipótese de auto parcial é exclusiva da polinização por insetos, a distribuição esperada é simplesmente binomial, e você pode estimar os parâmetros com uma família binomial glm ou talvez um glmm com identificação da planta como um efeito aleatório. No entanto, se eles podem se auto parcial e receber polinização por insetos, você volta a precisar de uma mistura de duas distribuições binomiais. Nesse caso, eu investigaria o uso do OpenBUGS ou JAGS para ajustar o modelo usando o MCMC.

Depois de ajustar os dois modelos aos seus dados, você os compara para ver qual deles se encaixa melhor, usando AIC ou BIC ou alguma outra métrica de sua escolha.

atiretoo - restabelecer monica
fonte
Obrigado por isso, mas executar esse código parece gerar um número aleatório de sementes, bem como uma distribuição aleatória. Eu estava pensando que eu queria que o nubmer de sementes a ser fixado (digamos 19 sementes, veja abaixo) e depois ver como provável uma dada distribuição era para que nubmer exata
BWGIA
Opps, clique no post muito cedo e eu quis dizer "veja acima", pois adicionei algumas informações à minha pergunta. Estou intrigado com o seu comentário sobre o uso da AIC para comparar modelos. Posso fazer isso entre modelos (com a mesma variável de resposta) com diferentes distribuições? Eu pensei que a comparação da AIC só era válida quando você adiciona / descarta termos em um modelo, mas com a mesma família de distribuição especificada?
BWGIA
Não, essa é a principal vantagem da AIC sobre, por exemplo, a seleção para trás. Desde que os dados sejam os mesmos, é possível comparar o AIC entre modelos diferentes, mesmo que não estejam aninhados. Você precisa ter cuidado para que o software esteja calculando probabilidades sem deixar constantes de fora, mas em uma única função você pode comparar modelos não aninhados com facilidade.
Atiretoo - reinstate monica