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!

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.

• That's very helpful thank you Robert. Will this structure effect how I test/interpret the assumptions? e.g. linearity – biolSas Jan 18 '20 at 12:34
• You could test assumptions in a number of ways , for example by inspecting residuals, or adding nonlinear terms to see if the fit improves, but that is really another question so please post a new question about that,.... also please see the extra paragraph I just added. – Robert Long Jan 18 '20 at 12:43
• Just seen the extra paragraph, thank you. Just to clarify further with your point about the elephant values, the values for elephant 'density' are actually a measures of elephants per sq kilometer in each area. – biolSas Jan 18 '20 at 13:06
• Was each value calculated on the same size area, or different sizes ? If the latter then you will have mathematical coupling, although the extent to which it will bias results depends on the relative sizes and variability in your data. I could add a small example to my answer, of how this happens, if you like ? – Robert Long Jan 18 '20 at 13:13
• Thank you very much for the explanation and further help, I'm working with theoretical data and wasn't given details on the areas used unfortunately! I ran into a few issues later on after examining and visualizing the results. I'm unsure if it's personal error or model error. If you are interested here is my newer related post: stats.stackexchange.com/questions/445465/… – biolSas Jan 19 '20 at 9:58