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

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


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.

• Oh that was a genuine mistake, didn't notice that. I ran into some issues earlier and was halfway through adding a new graph to show the commenter when I realized I'd made an error in the code and it wasn't necessary. I must have deleted the original graph accidentally in the process. I will re-add it. Jan 19, 2020 at 23:09

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