Właściwy sposób radzenia sobie z 3-poziomową tabelą awaryjną

Mam trzypoziomową tabelę zdarzeń awaryjnych, z danymi zliczeń dla kilku gatunków, rośliny żywicielskiej, z której zostały zebrane, oraz tego, czy zbiór ten miał miejsce w deszczowy dzień (to naprawdę ma znaczenie!). Przy użyciu R fałszywe dane mogą wyglądać mniej więcej tak:

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

Teraz chcę wiedzieć dwie rzeczy: Czy gatunki są związane z roślinami żywicielskimi? Czy „deszcz czy nie” wpływa na to skojarzenie? Użyłem loglm()z MASStego do:

 # 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)

Jest to nieco poza moim poziomem komfortu i chciałem sprawdzić, czy odpowiednio ustawiłem modele i że był to najlepszy sposób, aby odpowiedzieć na te pytania.

Istnieją dwa sposoby interpretacji pierwszego pytania, które znajdują odzwierciedlenie w dwóch zadanych przez Ciebie pytaniach: „Czy gatunki są związane z roślinami żywicielskimi?” oraz „Czy gatunki są niezależne od roślin żywicielskich, biorąc pod uwagę efekt deszczu?”

Pierwsza interpretacja odpowiada modelowi wspólnej niezależności , który stwierdza, że ​​gatunki i żywiciele są zależne, ale wspólnie niezależne od tego, czy padało:

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

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

Druga interpretacja odpowiada modelowi warunkowej niezależności , który stwierdza, że ​​gatunki i żywiciele są niezależni, biorąc pod uwagę, czy padało:

$\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$

Możesz przetestować te modele w R (też loglinby działało dobrze, ale jestem bardziej zaznajomiony 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)

Innym sposobem podejścia do pierwszego pytania byłoby wykonanie dokładnego testu Fischera ( fisher.test(xtabs(count ~ host + species))) na zwiniętej tabeli 2x2 (pierwsza interpretacja) lub testu Mantela-Haenszela ( mantelhaen.test(xtabs(count ~ host + species + rain))) dla 2-warstwowych tabel 2x2 lub napisanie testu permutacyjnego uwzględniającego stratyfikację (druga interpretacja).

Parafrazując swoje drugie pytanie: Czy związek między gatunkiem a żywicielem zależy od tego, czy padał deszcz?

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)

Pełny model mod4jest nasycony, ale można przetestować dany efekt, patrząc na dewiację, mod3tak jak to zrobiłem powyżej.

$host$$rain$$host*rain$

Poleceniem R byłoby:

glm (formuła = gatunek ~ gospodarz + deszcz, rodzina = dwumianowy (logit), wagi = liczy)

$p$

Początkowo zasugerowałem wypróbowanie jednej z ograniczonych technik święceń z veganpakietu, ale z drugiej strony wątpię, czy byłoby to przydatne, ponieważ faktycznie masz 2 tabele zdarzeń awaryjnych. Mam nadzieję, że druga część tego przykładu [PDF: R Demonstracja - Analiza kategoryczna] będzie pomocna.