Interaction between variables

Y = β_{0}^{} + β_{1}^{} X_{1}^{} + β_{2}^{} X_{2}^{} + β_{3}^{} X_{1}^{2} + β_{4}^{} X_{2}^{2} + β_{5}^{} X_{1}^{} X_{2}^{} + ε

library(ggplot2)
  library(plot3D)

# Uses the function lm() to compute a linear regression model.
    # Uses the argument surf in scatter3D() function to add a regression surface.
      
    # X1, X2, Y variables
    x <- rep(-2:2, 5)      
    y <- rep(-2:2, each=5)
      
    z <- beta0 + beta1*x + beta2*y + beta3*x*x + beta4*y*y + beta5*x*y + rnorm(25, 0, sigma)
    fit <- lm(z ~ x + y + I(x^2) + I(y^2) + I(x*y)) 
      
    # predict values on regular xy grid
    grid.lines = 20
    
    x.pred <- seq(min(x), max(x), length.out = grid.lines)
    y.pred <- seq(min(y), max(y), length.out = grid.lines)
    
    xy <- expand.grid(x = x.pred, y = y.pred)
    
    z.pred <- matrix(predict(fit, newdata = xy), nrow = grid.lines, ncol = grid.lines)
    
    # fitted points for droplines to surface
    fitpoints <- predict(fit)
    
    # scatter plot with regression plane
    scatter3D(x, y, z, pch = 19, colvar = FALSE, col = 2,
              theta = 30, phi = 20, ticktype = "detailed",
              xlab = "X1", ylab = "X2", zlab = "Y",  
              surf = list(x = x.pred, y = y.pred, z = z.pred, col = "lightgrey", alpha = 0.4, facets = NA, fit = fitpoints), 
              main = "Red: measurements,   Grid: estimated function")

# ------------------------------------------ Y against the X1-variable (as a line-graph)

model <- lm( z ~ x )                                   # Simple analysis    
   
mat_coeff <- model$coefficients                        # The two estimated coefficients.

b0 <- mat_coeff[1]   # (Intercept)                     # Extracts coefficients.
b1 <- mat_coeff[2]   # x
 
myData <- data.frame(x, z)

myDiagram <- ggplot(myData, aes(x = x, y = z)) + labs(x = "X1-values", y = "Y-values") +
     annotate("text", x = 0.0, y = Inf, label = "Estimated model (X1)", vjust = 2, fontface = "bold", color = "blue", size = 4) + 
     geom_point(colour = "red") + geom_abline(intercept = b0, slope = b1, color = "lightgrey")
myDiagram

# ------------------------------------------ Y against the X2-variable (as a line-graph)

model <- lm( z ~ y )                                   # Simple analysis    
   
mat_coeff <- model$coefficients                        # The two estimated coefficients.

b0 <- mat_coeff[1]   # (Intercept)                     # Extracts coefficients.
b1 <- mat_coeff[2]   # y
 
myData <- data.frame(y, z)

myDiagram <- ggplot(myData, aes(x = y, y = z)) + labs(x = "X2-values", y = "Y-values") +
     annotate("text", x = 0.0, y = Inf, label = "Estimated model (X2)", vjust = 2, fontface = "bold", color = "blue", size = 4) + 
     geom_point(colour = "red") + geom_abline(intercept = b0, slope = b1, color = "lightgrey")
myDiagram

library(ggplot2)
  library(plot3D)

# Uses the function lm() to compute a linear regression model.
    # Uses the argument surf in scatter3D() function to add a regression surface.
      
    # X1, X2, Y variables
    
    x <- c(2, 3, 4, 1, 5, 6, 7, 8 );
    y <- c(1, 2, 5, 2, 6, 4, 3, 4 );
          
    z <- beta0 + beta1*x + beta2*y + beta3*x*x + beta4*y*y + beta5*x*y + rnorm(8, 0, sigma)
    fit <- lm(z ~ x + y + I(x^2) + I(y^2) + I(x*y)) 
      
    # predict values on regular xy grid
    grid.lines = 20
    
    x.pred <- seq(min(x), max(x), length.out = grid.lines)
    y.pred <- seq(min(y), max(y), length.out = grid.lines)
    
    xy <- expand.grid(x = x.pred, y = y.pred)
    
    z.pred <- matrix(predict(fit, newdata = xy), nrow = grid.lines, ncol = grid.lines)
    
    # fitted points for droplines to surface
    fitpoints <- predict(fit)
    
    # scatter plot with regression plane
    scatter3D(x, y, z, pch = 19, colvar = FALSE, col = 2,
              theta = 30, phi = 20, ticktype = "detailed",
              xlab = "X1", ylab = "X2", zlab = "Y",  
              surf = list(x = x.pred, y = y.pred, z = z.pred, col = "lightgrey", alpha = 0.4, facets = NA, fit = fitpoints), 
              main = "Red: measurements,   Grid: estimated function")

# ------------------------------------------ Y against the X1-variable (as a line-graph)

model <- lm( z ~ x )                                   # Simple analysis    
   
mat_coeff <- model$coefficients                        # The two estimated coefficients.

b0 <- mat_coeff[1]   # (Intercept)                     # Extracts coefficients.
b1 <- mat_coeff[2]   # x
 
myData <- data.frame(x, z)

myDiagram <- ggplot(myData, aes(x = x, y = z)) + labs(x = "X1-values", y = "Y-values") +
     annotate("text", x = 2.0, y = Inf, label = "Estimated model (X1)", vjust = 2, fontface = "bold", color = "blue", size = 4) + 
     geom_point(colour = "red") + geom_abline(intercept = b0, slope = b1, color = "lightgrey")
myDiagram

# ------------------------------------------ Y against the X2-variable (as a line-graph)

model <- lm( z ~ y )                                   # Simple analysis    
   
mat_coeff <- model$coefficients                        # The two estimated coefficients.

b0 <- mat_coeff[1]   # (Intercept)                     # Extracts coefficients.
b1 <- mat_coeff[2]   # y
 
myData <- data.frame(y, z)

myDiagram <- ggplot(myData, aes(x = y, y = z)) + labs(x = "X2-values", y = "Y-values") +
     annotate("text", x = 2.0, y = Inf, label = "Estimated model (X2)", vjust = 2, fontface = "bold", color = "blue", size = 4) + 
     geom_point(colour = "red") + geom_abline(intercept = b0, slope = b1, color = "lightgrey")
myDiagram