Appropriate way to deal with a 3-level contingency table

Question

I have a three level contingency table, with count data for several species, the host plant from which they were collected and whether that collection happened on a rainy day (this actually matters!). Using R, fake data might be something like this:

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

Now, I want to know two things: Are species associated with host plants? Does "rain or not" effect this association. I used loglm() from MASS for this:

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

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

This is a little outside my comfort level, and wanted to check that I'd set models up right and that this was the best way to approach these questions.

Answer 1

There are two ways to interpret your first question, which are reflected in the two ways you asked it: “Are species associated with host plants?” and, “Are species independent to host plants, given the effect of rain?”

The first interpretation corresponds to a model of joint independence, which states that species and hosts are dependent, but jointly independent of whether it rained:

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

where $p_{shr}$ is the probability that an observation falls into the $(s,h,r)$ cell where $s$ indexes species, $h$ host type, and $r$ rain value, $p_{sh}$ is the marginal probability of the $(s,h,\cdot)$ cell where we collapse over the rain variable, and $p_r$ is the marginal probability of rain.

The second interpretation corresponds to a model of conditional independence, which states that species and hosts are independent given whether it rained:

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

where $p_{sh|r}$ is the conditional probability of the $(s,h,r)$ cell, given a value of $r$.

You can test these models in R (loglin would work fine too but I’m more familiar with 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

Above, mod1 corresponds to joint independence and mod2 corresponds to conditional independence, whereas mod0 corresponds to a mutual independence model $p_{shr} = p_s p_h p_r$. You can see the parameter estimates using summary(mod2), etc. As usual, you should check to see if model assumptions are met. In the data you provided, the null model actually fits adequately.

A different way of approaching your first question would be to perform Fischer’s exact test (fisher.test(xtabs(count ~ host + species))) on the collapsed 2x2 table (first interpretation) or the Mantel-Haenszel test (mantelhaen.test(xtabs(count ~ host + species + rain))) for 2-stratified 2x2 tables or to write a permutation test that respects the stratification (second interpretation).

To paraphrase your second question, Does the relationship between species and host depend on whether it rained?

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)

The full model mod4 is saturated, but you can test the effect in question by looking at the deviance of mod3 as I've done above.

Answer 2

Logistic regression seems appropriate for your problem. The variable you are trying to predict is the probability that an observation (which is either species A or species B) is species A. The covariates are $host$, $rain$ and optionally $host*rain$.

The R command would be:

glm(formula = species ~ host + rain, family = binomial(logit), weights=counts)

and you will be interested in the $p$-values of the slopes. Keep in mind that you are testing multiple hypotheses, though.

Answer 3

Initially I suggested to try one of the constrained ordination techniques from the vegan package, but on a second thought I doubt that this would be useful, as you actually have 2 contingency tables. I hope that the second part of this example [PDF: R Demonstration – Categorical Analysis] will be helpful.