What statistical model can incorporate my categorical variable (2 levels) and 2 continuous variables? I have some data which has a categorical variable (treatment) with 2 levels, control and 'mice added'. I have 2 continuous variables of plant density and elephant number/density. I have examined the data visually (see graph) and the addition of mice seems to affect the relationship between plant density and elephant density. What statistical test can I use which 1. incorporates both categorical and continuous variables and 2. which will allow me to examine not only the relationship between elephant density and plant density, but also test if treatment affects this relationship? Would using an ANCOVA model be acceptable? And in this case based on the research question, would I need to treat the categorical variable as the covariate? Advice would be much appreciated!

 A: An ANCOVA model would be one way. It would look something like
ElephantDensity ~ Control * PlantDensity

This will fit fixed effects for Control and PlantDensity, along with an interaction term between them. The interaction will quantify the extent to which the association between PlantDensityand ElephantDensity varies at the different levels of treatment.
Alternatively it may be a better idea to model elephant numbers rather than density, if the denominator that you use to calculate the densities also varies. This is because the two densities will then be linked and bias due to mathematical coupling may be invoked. If the denominator (an area, I assume) is fixed, then this is not a concern.
Here is a very simple simulation that shows what can go wrong when dividing the response and the exposure by a third variable. In this case the response is the number of elephants observed, the exposure is the number of plants, and the third variable is the size / area of where the samples were taken:
> set.seed(15)

> N <- 100    # number of sites sampled

> x <- rpois(N, 5) # number of plants

> y <- round(5 - 0.2 * x + rnorm(N, 0, 2))   # number of elephants

> m0 <- lm(y ~ x)
> summary(m0)

Here I have fitted a linear model, where a Poisson would be better, for the sake of simplicity. We obtain:
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.61810    0.56270   8.207 9.16e-13 ***
x           -0.17232    0.09384  -1.836   0.0693 .  

So we obtain results that are quite close to the "truth" of 5 and -0.2. All is rosey !
Now, let's introduce an area variable and create density variables for both elephants and plants:
> area <- runif(N,2,5)
> m1 <- lm( I(y/area ~ I(x/area)))
> summary(m1)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.8681     0.1738   4.994 2.57e-06 ***
I(x/area)     0.1603     0.0913   1.756   0.0822 . 

Whoops ! We now have an estimate of similar size, but opposite sign!! So be careful when dividing two variables by a third one, as severe bias due to mathematical coupling may be invoked.
