¿Comparación simple de dos medias de Poisson en R?

Question

Tengo dos muestras y me gustaría determinar si la diferencia entre ellas es estadísticamente significativa o no:

Dado que se trata claramente de datos de recuentos pequeños que no pueden ser aproximados por una distribución normal, no creo que una prueba t sea apropiada. En su lugar, creo que puedo asumir que estos recuentos siguen (dos distintas) distribuciones de Poisson. Como se ha dicho aquí La mejor manera de resumir cada una de estas distribuciones de Poisson es mediante las respectivas medias muestrales, como se indica a continuación:

$$ \overline{y} = \frac{1}{N} \sum_{i = 1}^n{y_i} $$

Obviamente, esto da lugar a medias no enteras ( 0.4 y 2.6 respectivamente), por lo que no puedo utilizar funciones como poisson.test() o cualquiera de las funciones del exactci biblioteca directamente.

Soy consciente de esta respuesta con respecto a la prueba C y la prueba E - pero ¿hay una implementación directa en R que haga esto?

La única forma razonable de hacerlo utilizando poisson.test() Lo que se me ocurre sería sumar los índices de cada condición y utilizar los números de los recuentos como los respectivos T parámetros:

# Count data for each respective Condition
Cond1 <- c(0, 0, 0, 1, 1)
Cond2 <- c(1, 2, 3, 3, 4)

poisson.test(
  x = c(sum(Cond1), sum(Cond2)),
  T = c(length(Cond1), length(Cond2)),
  alternative = "two.sided"
)

lo que arroja el siguiente resultado:

    Comparison of Poisson rates

data:  c(sum(Cond1), sum(Cond2)) time base: c(length(Cond1), length(Cond2))
count1 = 2, expected count1 = 7.5, p-value = 0.007385
alternative hypothesis: true rate ratio is not equal to 1
95 percent confidence interval:
 0.01685531 0.67955077
sample estimates:
rate ratio 
 0.1538462

¿Es esto válido o hay una forma mejor de hacerlo en R?

Answer 1

Esto es lo que yo haría que creo que es lo que @Dave está insinuando en el comentario.

Ajuste el GLM con la distribución de la familia Poisson y vea el efecto de la "condición" en la media de cada grupo.

Prepara el conjunto de datos:

Cond1 <- c(0, 0, 0, 1, 1)
Cond2 <- c(1, 2, 3, 3, 4)

dat <- data.frame(
    cond= c(rep('C1', length(Cond1)), rep('C2', length(Cond2))),
    count= c(Cond1, Cond2)
)

Ajuste el modelo y evalúe la importancia de la "condición":

fit <- glm(count ~ cond, data= dat, family= 'poisson')
summary(fit)

Call:
glm(formula = count ~ cond, family = "poisson", data = dat)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.13533  -0.89443  -0.07296   0.65703   0.80391  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -0.9163     0.7071  -1.296   0.1950  
condC2        1.8718     0.7595   2.464   0.0137 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 14.8823  on 9  degrees of freedom
Residual deviance:  5.8682  on 8  degrees of freedom
AIC: 27.731

Number of Fisher Scoring iterations: 5

Así, la media de Cond 1 se estima como exp(-0.9163) = 0.4 la media de Cond 2 es exp(-0.9163 + 1.8718) = 2.6 y la diferencia tiene p = 0,0137 para la hipótesis nula de ser cero.