É mais fácil usar o comando predict() e sobrepor os ajustados sobre o gráfico dos pontos.

da <- expand.grid(x=1:10, tr=gl(3,1))
da$y <- as.numeric(da$tr)+0.5*da$x+rnorm(da$x,0,0.1)
plot(y~x, da, col=da$tr)

m0 <- lm(y~tr+x, da)

pred <- expand.grid(tr=levels(da$tr), x=seq(1,10,length=30))
pred$y <- predict(m0, newdata=pred)

plot(y~x, da, col=da$tr)
by(pred, INDICES=pred$tr, function(i) lines(i$x, i$y, col=i$tr))

Um exemplo com gráfico melhor em

https://ridiculas.wordpress.com/2012/04/23/ajuste-simultaneo-de-regressaoanova-para-varias-respostas/

À disposição.
Walmes.

==========================================================================
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: walmes@ufpr.br
skype: walmeszeviani
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==========================================================================