Eu uso esse tipo de função (do conjunto de dados da íris)

model_test <- lm(Sepal.Length ~( Petal.Length + Sepal.Width + Petal.Width +Species)^2,
data=iris)
gg<-expand.grid(Petal.Length=0:6,Species=levels(Species))
vv<-expand.grid(Sepal.Width=0:4,Petal.Width=1:4)

pd<-do.call(rbind,Map(function(Petal.Length,Species,Sepal.Width,Petal.Width){
              nd <- cbind(vv, Petal.Length=Petal.Length,Species=Species,
                               Sepal.Width=Sepal.Width, Petal.Width=Petal.Width)
              cbind(nd, pred=predict(model_test, nd, type="response"))},
              Petal.Length=iris$Petal.Length,Species=iris$Species,
               Sepal.Width=iris$Sepal.Width,Petal.Width=iris$Petal.Width))

wireframe(pred~Sepal.Width+Petal.Width|Species*Petal.Length,
           pd, drape=FALSE,scale=list(arrows=FALSE),subset=(Species=="setosa"),
           layout = c(3, 3))

Eu recebo esse enredo

Minhas perguntas): Se eu incluir mais fatores no meu modelo, as tiras (coloridas) tornariam os gráficos muito pequenos. Como remover a faixa vermelha verde e clara e adicionar um título de 1 linha comoSpecies==setosa & Petal.Length==0, ....,Species==setosa & Petal.Length==6.