The lines on my scatterplot for ANCOVA results doesn't look right, personal error or model error?

Question

I ran an ANCOVA model, to test the effect of treatment on the relationship between two continuous variables (elephant number and plant density) - please see more details in my last question (What statistical model can incorporate my categorical variable (2 levels) and 2 continuous variables?) but when visualizing my results on a scatterplot something isn't quite right. Here are my ANCOVA summary results:

lm(formula = eledens ~ treat * plants, data = elemice)

Residuals:
   Min     1Q Median     3Q    Max 
-7.102 -2.715  0.264  1.814  9.235 

Coefficients:
                       Estimate Std. Error t value    Pr(>|t|)    
(Intercept)              8.3028     1.8952   4.381 0.000097820 ***
treatMice added         -0.7066     2.7105  -0.261    0.795810    
plants                   0.7368     0.1232   5.978 0.000000743 ***
treatMice added:plants  -0.6840     0.1613  -4.241    0.000148 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.856 on 36 degrees of freedom
Multiple R-squared:  0.7393,    Adjusted R-squared:  0.7176 
F-statistic: 34.04 on 3 and 36 DF,  p-value: 0.0000000001308

And here is my plotting code trying to visualize the result:

plot(eledens~plants, data=elemice, type="n", xlab="Plant Density", ylab="Elephant Density")
points(elemice$plants[elemice$treat=="Control"], elemice$eledens[elemice$treat=="Control"], col="skyblue3", pch=16)
points(elemice$plants[elemice$treat=="Mice added"], elemice$eledens[elemice$treat=="Mice added"], col="salmon", pch=16)
abline(fit.mice$coefficients[1:2], col="skyblue3")
abline(fit.mice$coefficients[1]+fit.mice$coefficients[3],fit.mice$coefficients[2], col="salmon")
new.x <- rep(seq(min(elemice$plants), max(elemice$plants), len=100),2)
new.s <- rep(c("Control","Mice added"), each=100)
pred <- predict(fit.mice, new=data.frame(plants=new.x, treat=new.s), interval="conf")
pred <- data.frame(pred, treat=new.s, plants=new.x)
head(pred)
lines(new.x[1:100],pred[1:100,"lwr"],lty=2, col="skyblue3")
lines(new.x[1:100],pred[1:100,"upr"],lty=2, col="skyblue3")
lines(new.x[101:200],pred[101:200,"lwr"],lty=2, col="salmon")
lines(new.x[101:200],pred[101:200,"upr"],lty=2, col="salmon")
legend("topleft", pch=16, col=c("skyblue3","salmon"), legend=c("Control","Mice added"))

As you can see below (see output graph) the fitted lines from the coefficients look a bit odd, so does this mean my model is incorrect? Or have I assigned the coefficients incorrectly in the coding? I am fairly new to R coding, so advice on interpreting this would be much appreciated.

Answer 1

You appear to be misunderstanding the output from your model. In your code, the line:

abline(fit.mice$coefficients[1:2], col="skyblue3")

plots a line with the correct intercept, but the wrong slope. fit.mice$coefficients[1] is the intercept, but fit.mice$coefficients[2] is the estimate for treat in the Mice added group, hence this is an offset to the intercept for the Mice added group. What you want for the slope in the Control group is simply the estimate for plants which is fit.mice$coefficients[3]. So:

abline(fit.mice$coefficients[1], fit.mice$coefficients[3] , col="skyblue3")

Then to plot the line for the Mice added group, the intercept will be the global intercept plus the estimate for treat and the slope will be the estimate for plants plus the interaction term. So:

abline(fit.mice$coefficients[1] + fit.mice$coefficients[2], fit.mice$coefficients[3] + fit.mice$coefficients[4], col="salmon")