¿Cómo probar / probar que los datos están inflados a cero?

9

Tengo un problema que creo que debería ser simple pero no puedo resolverlo. Estoy mirando la polinización de semillas, tengo plantas (n = 36) que florecen en racimos, tomo muestras de 3 racimos de flores de cada planta y 6 vainas de semillas de cada grupo (18 vainas de semillas en total de cada planta). Una vaina puede tener entre 0 y como máximo 4 semillas polinizadas. Entonces, los datos se cuentan, con un límite superior. Estoy encontrando que un promedio de ~ 10% de las semillas son polinizadas, pero entre 1 y 30% en una planta determinada, por lo que sobre datos dispersos, y por supuesto, faltan 4 réplicas de racimos en 3 plantas, por lo que no son perfectamente simétricas .

La pregunta que hago es si estos datos respaldan la idea de que esta planta requiere polinizadores para el conjunto de semillas.

Estoy descubriendo que la distribución de la cantidad de semillas en una vaina parece que hay más 0 vainas de semillas polinizadas (6-9 vainas de 16) y más 3 y 4 vainas de semillas polinizadas (2-4 para cada una) de lo que lo haría se esperaría si las semillas de la población se polinizaran al azar. Básicamente, creo que este es un ejemplo clásico para datos inflados cero, primero un insecto visita o no la flor (un generador cero) y si lo hace, luego poliniza 0-4 de las semillas en otra distribución. La hipótesis alternativa es que la planta se está autoajustando parcialmente, y entonces se esperaría que cada semilla tuviera la misma probabilidad de ser polinizada (estos datos sugieren una probabilidad de aproximadamente 0.1, lo que significa una probabilidad de 0.01 para dos semillas en la misma vaina, etc.) .

Pero simplemente quiero demostrar que los datos se ajustan mejor a una u otra distribución, no HACER realmente un ZIP o ZINB en los datos. Creo que cualquier método que use debe tener en cuenta la cantidad real de semillas polinizadas y la cantidad de vainas muestreadas en cada planta. Lo mejor que se me ocurrió es hacer algún tipo de correa de arranque donde simplemente asigno al azar el número de semillas polinizadas para una planta determinada en el número de vainas de semillas que probé, lo hago 10,000 veces y veo qué tan probable es los datos experimentales para la planta dada surgieron de esa distribución aleatoria.

Siento que hay algo en esto que debería ser mucho más fácil que el arranque de fuerza bruta, pero después de días de pensar y buscar, me doy por vencido. No puedo compararme con una distribución de Poisson porque es el límite superior, no es binomial porque necesito generar la distribución esperada de alguna manera primero. ¿Alguna idea? Y estoy usando R, así que los consejos allí (especialmente cómo generar de manera más elegante 10,000 distribuciones aleatorias de n bolas en 16 cajas que pueden contener cada una como máximo 4 bolas) serían bienvenidas.

AGREGADO 07/09/2012 Primero, gracias a todos por todo el interés y la ayuda. Leer las respuestas me ha hecho pensar en reformular un poco mi pregunta. Lo que digo es que tengo una hipótesis (que por ahora estoy pensando como nula) de que las semillas se polinizan aleatoriamente entre las vainas, y mi hipótesis alternativa es que una vaina de semillas con al menos 1 semilla polinizada tiene más probabilidades de tener múltiples semillas polinizadas de lo que se esperaría de un proceso aleatorio. He proporcionado datos reales de tres plantas como ejemplos para ilustrar de lo que estoy hablando. La primera columna es el número de semillas polinizadas en una vaina, la segunda columna es la frecuencia de las vainas con ese recuento de semillas.

planta 1 (total de 3 semillas: 4% de polinización)

num.seeds :: pod.freq

0 :: 16

1 :: 1

2 :: 1

3 :: 0

4 :: 0

planta 2 (total de 19 semillas: 26% de polinización)

num.seeds :: pod.freq

0 :: 12

1 :: 1

2 :: 1

3 :: 0

4 :: 4

planta 3 (total de 16 semillas: 22% de polinización)

num.seeds :: pod.freq

0 :: 9

1 :: 4

2 :: 3

3 :: 2

4 :: 0

En la planta # 1, solo se polinizaron 3 semillas en 18 vainas, una vaina tenía una semilla y una vaina tenía dos semillas. Pensando en un proceso de agregar una semilla a las vainas al azar, las dos primeras semillas van a su propia vaina, pero para la 3ra semilla, hay 6 lugares disponibles en las vainas que ya tienen una semilla pero 64 en las 16 vainas. sin semillas, entonces la probabilidad más alta de una vaina con 2 semillas aquí es 6/64 = 0.094. Eso es un poco bajo, pero no realmente extremo, por lo que diría que esta planta se ajusta a la hipótesis de la polinización aleatoria en todas las semillas con una probabilidad de ~ 4% de polinización. Pero la planta 2 me parece mucho más extrema, con 4 vainas completamente polinizadas, pero 12 vainas sin nada. No estoy muy seguro de cómo calcular las probabilidades de esta distribución directamente (de ahí mi idea inicial), pero supongo que las probabilidades de que esta distribución ocurra al azar si cada semilla tiene un ~ 25% de posibilidades de polinización son bastante bajas. Planta # 3 Realmente no tengo idea, creo que hay más 0 y 3 de lo que uno esperaría para una distribución aleatoria, pero mi intuición es que esta distribución para este número de semillas es mucho más probable que la distribución para la planta # 2, y puede que no sea tan improbable. Pero obviamente quiero saber con certeza, y en todas las plantas. Creo que hay más 0 y 3 de lo que uno esperaría para una distribución aleatoria, pero mi intuición es que esta distribución para este número de semillas es mucho más probable que la distribución para la planta # 2, y puede no ser tan improbable. Pero obviamente quiero saber con certeza, y en todas las plantas. Creo que hay más 0 y 3 de lo que uno esperaría para una distribución aleatoria, pero mi intuición es que esta distribución para este número de semillas es mucho más probable que la distribución para la planta # 2, y puede no ser tan improbable. Pero obviamente quiero saber con certeza, y en todas las plantas.

Al final, estoy buscando escribir una declaración como “La distribución de semillas polinizadas en vainas de semillas se ajusta (o no se ajusta) a la hipótesis de que las plantas no son simplemente parcialmente autocompatibles, sino que requieren la visita de un polinizador para efectuar el conjunto de semillas. (resultados de la prueba estadística) ". Esto es realmente solo una parte de mi sección prospectiva, donde estoy hablando sobre qué experimentos realizar a continuación, así que no estoy desesperado por que esto sea una cosa u otra, pero quiero saber por mí mismo, si es posible. Si no puedo hacer lo que intento hacer con estos datos, ¡también me gustaría saberlo!

Al principio hice una pregunta bastante amplia, ya que tengo curiosidad por saber si hay buenas pruebas para mostrar si los datos deben entrar en un modelo inflado cero en primer lugar. Todos los ejemplos que he visto parecen decir: "mira, hay muchos ceros aquí, y hay una explicación razonable para eso, así que usemos un modelo inflado cero". Eso es lo que estoy haciendo ahora en este foro, pero tuve una experiencia en mi último capítulo en el que utilicé una película de Poisson para datos de conteo y uno de mis supervisores dijo: "No, las películas son demasiado complejas e innecesarias, estos datos deberían ir a una tabla de contingencia "y luego me envió un volcado de datos de la tabla de contingencia masiva generada por su costoso paquete de estadísticas que dio los mismos valores de p para todos mis factores + interacciones a tres dígitos significativos !! Entonces, trato de mantener las estadísticas claras y simples, y me aseguro de comprenderlos lo suficientemente bien como para defender con firmeza mis elecciones, lo que no creo que pueda hacer por un modelo inflado cero en este momento. He utilizado tanto un cuasibinomio (para plantas enteras para deshacerse de la pesudoreplicación) como un modelo mixto para los datos anteriores para comparar tratamientos y responder a mis preguntas experimentales principales, o parece que hago el mismo trabajo, pero también voy a juega con ZINB's esta noche, para ver qué tan bien funciona. Estoy pensando que si puedo demostrar explícitamente que estos datos están fuertemente agrupados (o cero inflados) al principio, luego proporcionar una buena razón biológica para que eso ocurra, estaré mucho mejor preparado para sacar posteriormente un ZINB, que simplemente compare uno con un modelo cuasibinomial / mixto y discuta ya que da mejores resultados, eso es lo que debería usar. lo que no creo que pueda hacer por un modelo inflado cero en este momento. He utilizado tanto un cuasibinomio (para plantas enteras para deshacerse de la pesudoreplicación) como un modelo mixto para los datos anteriores para comparar tratamientos y responder a mis preguntas experimentales principales, o parece que hago el mismo trabajo, pero también voy a juega con ZINB's esta noche, para ver qué tan bien funciona. Estoy pensando que si puedo demostrar explícitamente que estos datos están fuertemente agrupados (o cero inflados) al principio, luego proporcionar una buena razón biológica para que eso ocurra, estaré mucho mejor preparado para sacar posteriormente un ZINB, que simplemente compare uno con un modelo cuasibinomial / mixto y discuta ya que da mejores resultados, eso es lo que debería usar. lo que no creo que pueda hacer por un modelo inflado cero en este momento. He utilizado tanto un cuasibinomio (para plantas enteras para deshacerse de la pesudoreplicación) como un modelo mixto para los datos anteriores para comparar tratamientos y responder a mis preguntas experimentales principales, o parece que hago el mismo trabajo, pero también voy a juega con ZINB's esta noche, para ver qué tan bien funciona. Estoy pensando que si puedo demostrar explícitamente que estos datos están fuertemente agrupados (o cero inflados) al principio, luego proporcionar una buena razón biológica para que eso ocurra, estaré mucho mejor preparado para sacar posteriormente un ZINB, que simplemente compare uno con un modelo cuasibinomial / mixto y discuta ya que da mejores resultados, eso es lo que debería usar. He utilizado tanto un cuasibinomio (para plantas enteras para deshacerse de la pesudoreplicación) como un modelo mixto para los datos anteriores para comparar tratamientos y responder a mis preguntas experimentales principales, o parece que hago el mismo trabajo, pero también voy a juega con ZINB's esta noche, para ver qué tan bien funciona. Estoy pensando que si puedo demostrar explícitamente que estos datos están fuertemente agrupados (o cero inflados) al principio, luego proporcionar una buena razón biológica para que eso ocurra, estaré mucho mejor preparado para sacar posteriormente un ZINB, que simplemente compare uno con un modelo cuasibinomial / mixto y discuta ya que da mejores resultados, eso es lo que debería usar. He utilizado tanto un cuasibinomio (para plantas enteras para deshacerse de la pesudoreplicación) como un modelo mixto para los datos anteriores para comparar tratamientos y responder a mis preguntas experimentales principales, o parece que hago el mismo trabajo, pero también voy a juega con ZINB's esta noche, para ver qué tan bien funciona. Estoy pensando que si puedo demostrar explícitamente que estos datos están fuertemente agrupados (o cero inflados) al principio, luego proporcionar una buena razón biológica para que eso ocurra, estaré mucho mejor preparado para sacar posteriormente un ZINB, que simplemente compare uno con un modelo cuasibinomial / mixto y discuta ya que da mejores resultados, eso es lo que debería usar.

Pero no quiero distraer demasiado de mi pregunta principal, ¿cómo puedo determinar si mis datos están realmente más inflados de lo esperado de una distribución aleatoria? En mi caso, la respuesta a eso es lo que realmente me interesa, y el posible beneficio de la justificación del modelo es una bonificación.

Gracias de nuevo por todo su tiempo y ayuda!

Saludos, BWGIA

BWGIA
fuente
¿Por qué no quieres ajustar el modelo binomial inflado a cero?
atiretoo - reinstalar monica
¿Es la hipótesis del "selfing parcial" exclusiva de la hipótesis del "polinizador"? Si es así, entonces su segundo modelo sería simplemente un modelo binomial con probabilidad p y tamaño = 4.
atiretoo - reinstalar a monica

Respuestas:

5

Esto me parece un modelo mixto relativamente sencillo (no lineal). Tiene vainas de semillas anidadas en racimos anidados en plantas, y puede ajustar un modelo binomial con efectos aleatorios en cada etapa:

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

o con covariables si los tienes. Si las flores se autopolinizan, es posible que vea algunos efectos leves debido a la variabilidad natural en la viabilidad de las plantas por sí mismas. Sin embargo, si la mayor parte de la variabilidad en la respuesta es impulsada por la variabilidad del grupo, tendría una evidencia más fuerte de polinización por insectos que podrían visitar solo grupos seleccionados en una planta. Idealmente, desearía una distribución no paramétrica de los efectos aleatorios en lugar de gaussiana: una masa puntual a cero, para ninguna visita de insectos, y una masa puntual a un valor positivo; este es esencialmente el modelo de mezcla en el que pensó Michael Chernick. Puede encajar esto con el paquete GLLAMM Stata, me sorprendería si esto no fuera posible en R.

Probablemente para un experimento limpio, desearía tener las plantas adentro, o al menos en un lugar sin acceso para insectos, y ver cuántas semillas se polinizarían. Eso probablemente respondería a todas sus preguntas de una manera metodológicamente más rigurosa.

StasK
fuente
Voy a intentar esto, creo que me ayudará a responder mis propias preguntas, pero no estoy tan seguro de cómo convencerá a los demás. Estás en lo cierto con la segunda parte, estoy tratando de pensar en cómo estos datos informan un futuro experimento más dirigido.
BWGIA
1

Me parece que esta es una distribución de mezcla para cada insecto individual. Con probabilidad p, el insecto aterriza con probabilidad 1-p aterriza y distribuye de 0 a 4 semillas. Pero si no tiene información sobre si el insecto aterriza o no en la planta, no puede distinguir las dos formas de obtener 0. Entonces podría dejar que p sea la probabilidad de 0 y luego tener la distribución multinomial (p1, p2, p3, p4) donde pi es la probabilidad de i semillas dado que el insecto poliniza sujeto a la restricción p1 + p2 + p3 + p4 = 1. El modelo tiene cinco incógnitas p, p1, p2, p3, p4 con la restricción 0 = 0 para cada i. Con suficientes datos, debería poder estimar estos parámetros tal vez utilizando un enfoque restringido de máxima verosimilitud.

Michael R. Chernick
fuente
Estoy de acuerdo, pero la pregunta no es ajustar ese modelo, sino generar distribuciones predichas bajo dos hipótesis biológicas diferentes. Tal vez la respuesta sea ajustar un ZIB y "algún otro modelo" que coincida con la hipótesis del selfing, y compararlos.
atiretoo - restablecer monica
@atiretoo ¿no le proporciona el modelo una distribución estimada para el número de semillas polinizadas que podría comparar con su distribución hipotética?
Michael R. Chernick
De acuerdo: si tiene los modelos correctos para las 2 hipótesis.
atiretoo - restablecer monica
1

Esta es una respuesta a la última parte de su pregunta, cómo generar rápidamente los datos que desea para la hipótesis del 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)

También puede usar rzibinom()en el paquete VGAM. Aunque no estoy seguro de qué quieres hacer con él. Tiene 2 parámetros libres, p1 y p2, que deben estimarse. ¿Por qué no utilizar un modelo binomial inflado a cero para estimarlos a partir de los datos?

Debe mirar el paquete VGAM, que se ajusta a los modelos ZIB entre otros. De hecho, puede obtener la distribución esperada para un ZIB de la función VGAM dzibinom(), que podría usar para comparar su distribución observada si conoce los parámetros de visitas y polinización. Una vez más, realmente debería ajustarse al modelo ZIB.

Si su hipótesis de selfing parcial es exclusiva de la polinización de insectos, entonces la distribución esperada es simplemente binomial, y podría estimar los parámetros con una glm de la familia binomial o tal vez una glmm con identificación de la planta como efecto aleatorio. Sin embargo, si pueden autoinyectarse Y recibir polinización de insectos, entonces volverá a necesitar una mezcla de dos distribuciones binomiales. En ese caso, investigaría usando OpenBUGS o JAGS para ajustar el modelo usando MCMC.

Una vez que tenga los dos modelos ajustados a sus datos, compare los modelos para ver cuál se ajusta mejor, utilizando AIC o BIC o alguna otra métrica de su elección.

atiretoo - restablecer monica
fuente
Gracias por ese atiretoo, pero ejecutar ese código parece generar un número aleatorio de semillas, así como una distribución aleatoria. Estaba pensando que quería que se arreglara el nubmer de semillas (digamos 19 semillas, ver abajo) y luego ver cuán probable era una distribución dada para ese
nubmer
Opps, presioné la publicación demasiado pronto y quise decir "ver arriba" ya que agregué información a mi pregunta. Estoy intrigado por su comentario sobre el uso de AIC para comparar modelos, ¿puedo hacerlo a través de modelos (con la misma variable de respuesta) con diferentes distribuciones? ¿Pensé que la comparación de AIC solo era válida cuando agregaba / quitaba términos a un modelo pero con la misma familia de distribución especificada?
BWGIA
No, esa es la ventaja clave de AIC sobre, por ejemplo, la selección hacia atrás. Mientras los datos sean los mismos, puede comparar AIC entre diferentes modelos, incluso si no están anidados. Debe tener cuidado de que el software esté calculando las probabilidades sin omitir constantes, pero dentro de una sola función puede comparar modelos no anidados con facilidad.
atiretoo - reinstalar monica