Maneira apropriada de lidar com uma tabela de contingência de três níveis

12

Eu tenho uma tabela de contingência de três níveis, com dados de contagem de várias espécies, a planta hospedeira da qual elas foram coletadas e se essa coleta ocorreu em um dia chuvoso (isso realmente importa!). Usando R, dados falsos podem ser algo assim:

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

Agora, quero saber duas coisas: as espécies estão associadas às plantas hospedeiras? "Chover ou não" afeta essa associação? Eu usei loglm()a partir MASSpara isso:

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

Isso está um pouco fora do meu nível de conforto e eu queria verificar se havia configurado os modelos corretamente e se essa era a melhor maneira de abordar essas questões.

david w
fonte

Respostas:

10

Há duas maneiras de interpretar sua primeira pergunta, que se refletem nas duas maneiras que você fez: "As espécies estão associadas às plantas hospedeiras?" e "As espécies são independentes das plantas hospedeiras, dado o efeito da chuva?"

A primeira interpretação corresponde a um modelo de independência conjunta , que afirma que espécies e hospedeiros são dependentes, mas conjuntamente independentes de chover:

pshr=pshpr

pshr(s,h,r)shrpsh(s,h,)pr

A segunda interpretação corresponde a um modelo de independência condicional , que afirma que espécies e hospedeiros são independentes, independentemente de chover:

psh|r=ps|rph|rpshr=psrphr/pr

psh|r(s,h,r) célula, dado um valor de r.

Você pode testar esses modelos em R ( loglintambém funcionaria bem, mas estou mais familiarizado 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

Acima, mod1corresponde à independência conjunta e mod2à independência condicional, enquanto mod0corresponde a um modelo de independência mútuapshr=psphpr. Você pode ver as estimativas de parâmetros usando summary(mod2), etc. Como sempre, verifique se as suposições do modelo são atendidas. Nos dados que você forneceu, o modelo nulo realmente se encaixa adequadamente.

Uma maneira diferente de abordar sua primeira pergunta seria executar o teste exato de Fischer ( fisher.test(xtabs(count ~ host + species))) na tabela 2x2 recolhida (primeira interpretação) ou o teste de Mantel-Haenszel ( mantelhaen.test(xtabs(count ~ host + species + rain))) para tabelas 2x2 estratificadas ou escrever um teste de permutação que respeite a estratificação (segunda interpretação).

Parafraseando sua segunda pergunta, a relação entre espécies e hospedeiro depende se choveu?

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)

O modelo completo mod4está saturado, mas você pode testar o efeito em questão observando o desvio de mod3como fiz acima.

bloqueado
fonte
Graças Lockedoff, especialmente para me ajudar a resolver o meu próprio pensamento wrt a diferença entre o modelo de independência condicional e conjunta
david w
1

A regressão logística parece apropriada para o seu problema. A variável que você está tentando prever é a probabilidade de uma observação (que é da espécie A ou da espécie B) ser da espécie A. As covariáveis ​​sãohost, rumaEun e opcionalmente hostrumaEun.

O comando R seria:

glm (fórmula = espécie ~ hospedeiro + chuva, família = binomial (logit), pesos = contagens)

e você estará interessado no p-valores das encostas. Lembre-se de que você está testando várias hipóteses.

charles.y.zheng
fonte
1
A regressão logística parece boa, mas tem a restrição adicional do total de linhas e colunas sendo corrigida. Pode não ser o caso dos dados de Poisson. Eu acredito que as respostas não serão muito diferentes.
suncoolsu
1

Inicialmente, sugeri tentar uma das técnicas de ordenação restrita do veganpacote, mas, pensando bem, duvido que isso seja útil, pois você realmente tem duas tabelas de contingência. Espero que a segunda parte deste exemplo [PDF: R Demonstration - Categorical Analysis] seja útil.

ils
fonte
Pense que o link está quebrado, você quis dizer este categórico aqui ? Isso foi útil, obrigado!
David w
Sim, parece que o espaço no URL o interrompe.
21711