處理3級列聯表的適當方法

我有一個三級列聯表，其中包含幾個物種的計數數據、收集它們的宿主植物以及該收集是否發生在下雨天（這實際上很重要！）。使用 R，假數據可能是這樣的：

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

現在，我想知道兩件事：物種是否與寄主植物有關？“下雨與否”會影響這種關聯嗎？我為此使用loglm()了MASS：

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

這有點超出我的舒適度，我想檢查我是否正確設置了模型，這是解決這些問題的最佳方法。

有兩種方式可以解釋您的第一個問題，這反映在您提出的兩種方式中：“物種是否與寄主植物有關？” 並且，“考慮到雨水的影響，物種是否獨立於寄主植物？”

第一種解釋對應於聯合獨立模型，該模型指出物種和宿主是相互依賴的，但共同獨立於是否下雨：

在哪裡是一個觀測值落入單元格在哪裡索引物種，主機類型，以及雨值，是邊際概率我們在雨變量上崩潰的單元格，以及是下雨的邊際概率。

第二種解釋對應於條件獨立模型，該模型表明物種和宿主是獨立的，無論是否下雨：

或者

在哪裡是條件概率單元格，給定一個值.

您可以在 R 中測試這些模型（loglin也可以正常工作，但我更熟悉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

上，mod1對應聯合獨立，mod2對應條件獨立，而mod0對應相互獨立模型. 您可以使用 等查看參數估計值summary(mod2)。像往常一樣，您應該檢查是否滿足模型假設。在您提供的數據中，空模型實際上非常適合。

處理第一個問題的另一種方法是對折疊的 2x2 表（第一種解釋）執行 Fischer 精確檢驗fisher.test(xtabs(count ~ host + species))( mantelhaen.test(xtabs(count ~ host + species + rain))（第二種解釋）。

套用你的第二個問題，物種和宿主之間的關係是否取決於是否下雨？

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)

完整的模型已經飽和，但是您可以通過查看我在上面所做mod4的偏差來測試有問題的效果。mod3

引用自：https://stats.stackexchange.com/questions/8342