7
$\begingroup$

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.

Effect of treatment on the relationship between elephant number and plant density

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – biolSas Jan 19 at 23:09
7
$\begingroup$

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")
$\endgroup$
  • $\begingroup$ Ah yes I see the error now! Revised the code with your suggestions and everything is now looking much better. Thanks again Robert, really appreciate the help. $\endgroup$ – biolSas Jan 19 at 14:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.