Quiero trazar la superficie de regresión de un modelo con un término de interacción usandorglEl sistema de trazado interactivo. Es fácil trazar un plano de regresión para un modelo sin un término de interacción usando:

plot3d(x=x1, y=x2, z=y1, type="s", col="yellow", size=1)
planes3d(a=coef(mod1)[2], b=coef(mod1)[3], c=-1, d=coef(mod1)[1], alpha=.5)

Sin embargo, cuando el avión gira, esto parece ser más difícil. Siguiendo con esta pregunta:3D equivalente de la función de curva en r, Estoy intentando:

f2 <- function(x, y) as.vector(coef(mod2)%*%c(1, x, y, x*y))

curve_3d <- function(f2, x_range=c(0, 40), y_range=c(0, 40)){ 
  if (!require(rgl) ) {stop("load rgl")}

  xvec <- seq(x_range[1], x_range[2], by=1)
  yvec <- seq(y_range[1], y_range[2], by=1)
  fz   <- outer(xvec, yvec, FUN=f2)
  persp3d(xvec, yvec, fz, alpha=.5)
}
open3d()
plot3d(x=x1, y=x2, z=y2, type="s", col="yellow", size=1)
curve_3d(f2)

Pero, no está funcionando. (He intentado algunas otras cosas también, pero me estoy quedando corto). Mi problema principal hasta ahora parece estar relacionado conf2; Sin embargo, también quiero que esto se vea comoplanes3d, y no estoy seguro de si esto me va a dar un alámbrico.

Aquí hay un ejemplo:

set.seed(897)
x1 = rep(c(0, 10, 20, 30, 40), times=25)
x2 = rep(c(0, 10, 20, 30, 40), each=25)
y2 = 37 + 0.7*x1 + 1.2*x2 - 0.05*x1*x2 + rnorm(125, mean=0, sd=5)
mod2 = lm(y2~x1*x2)
open3d()
plot3d(x=x1, y=x2, z=y2, type="s", col="yellow", size=1)
curve_3d(f2)