Manera apropiada de lidiar con una tabla de contingencia de 3 niveles

Tengo una tabla de contingencia de tres niveles, con datos de conteo para varias especies, la planta huésped de la cual fueron recolectadas y si esa recolección ocurrió en un día lluvioso (¡esto realmente importa!). Usando R, los datos falsos podrían ser algo como esto:

count    <- rpois(8, 10)
species  <- rep(c("a", "b"), 4)
host     <- rep(c("c","c", "d", "d"), 2)
rain     <- c(rep(0,4), rep(1,4))
my.table <- xtabs(count ~ host + species + rain)


, , rain = 0

    species
host  a  b
   c 12 15
   d 10 13

, , rain = 1

    species
host  a  b
   c 11 12
   d 12  7

Ahora, quiero saber dos cosas: ¿Están las especies asociadas con las plantas hospederas? ¿La lluvia o no afecta esta asociación? Solía loglm()de MASSpara esto:

 # Are species independent to host plants, given the effect of rain?
loglm(~species + host + rain + species*rain + host*rain, data=my.table)

 # Given any relationship between host plants and species, does rain change it?
loglm(~species + host + rain + species*host)

Esto está un poco fuera de mi nivel de comodidad, y quería comprobar que había configurado correctamente los modelos y que esta era la mejor manera de abordar estas preguntas.

Hay dos formas de interpretar su primera pregunta, que se reflejan en las dos formas en que la hizo: "¿Las especies están asociadas con las plantas hospederas?" y, "¿Son las especies independientes de las plantas hospederas, dado el efecto de la lluvia?"

La primera interpretación corresponde a un modelo de independencia conjunta , que establece que las especies y los huéspedes son dependientes, pero conjuntamente independientes de si llovió:

$\quad p_{shr} = p_{sh} p_r$

$p_{shr}$$(s,h,r)$$s$$h$$r$$p_{sh}$$(s,h,\cdot)$$p_r$

La segunda interpretación corresponde a un modelo de independencia condicional , que establece que las especies y los huéspedes son independientes dado que llovió:

$\quad p_{sh|r} = p_{s|r}p_{h|r}$$p_{shr} = p_{sr}p_{hr} / p_r$

$p_{sh|r}$$(s,h,r)$$r$

Puede probar estos modelos en R ( loglintambién funcionaría bien, pero estoy más familiarizado con esto glm):

count <- c(12,15,10,13,11,12,12,7)
species <- rep(c("a", "b"), 4)
host <- rep(c("c","c", "d", "d"), 2)
rain <- c(rep(0,4), rep(1,4))
my.table <- xtabs(count ~ host + species + rain)
my.data <- as.data.frame.table(my.table)
mod0 <- glm(Freq ~ species + host + rain, data=my.data, family=poisson())
mod1 <- glm(Freq ~ species * host + rain, data=my.data, family=poisson())
mod2 <- glm(Freq ~ (species + host) * rain, data=my.data, family=poisson())
anova(mod0, mod1, test="Chi") #Test of joint independence
anova(mod0, mod2, test="Chi") #Test of conditional independence

mod1mod2mod0$p_{shr} = p_s p_h p_r$summary(mod2)

Una forma diferente de abordar su primera pregunta sería realizar la prueba exacta de Fischer ( fisher.test(xtabs(count ~ host + species))) en la tabla colapsada 2x2 (primera interpretación) o la prueba Mantel-Haenszel ( mantelhaen.test(xtabs(count ~ host + species + rain))) para tablas 2x2 estratificadas 2 o escribir una prueba de permutación que respete la estratificación (segunda interpretación).

Parafraseando su segunda pregunta, ¿la relación entre las especies y el huésped depende de si llovió?

mod3 <- glm(Freq ~ species*host*rain - species:host:rain, data=my.data, family=poisson())
mod4 <- glm(Freq ~ species*host*rain, data=my.data, family=poisson())
anova(mod3, mod4, test=”Chi”)
pchisq(deviance(mod3), df.residual(mod3), lower=F)

El modelo completo mod4está saturado, pero puede probar el efecto en cuestión observando la desviación de mod3lo que he hecho anteriormente.

$host$$rain$$host*rain$

El comando R sería:

glm (fórmula = especie ~ huésped + lluvia, familia = binomio (logit), pesos = recuentos)

$p$

Inicialmente sugerí probar una de las técnicas de ordenación restringida del veganpaquete, pero pensándolo bien, dudo que esto sea útil, ya que en realidad tienes 2 tablas de contingencia. Espero que la segunda parte de este ejemplo [PDF: Demostración R - Análisis categórico] sea útil.