# Simulate data for a regression model with two explanatory variables

################################################################################
# Simulate data

  set.seed(7129)
  n <- 10
  x1 <- round(runif(n = n, min = 1, max = 4), digits = 1)
  x2 <- round(rnorm(n = n, mean = 3, sd = 0.5), digits = 1)
  beta0 <- -5
  beta1 <- 2
  beta2 <- 2
  sigma <- 1
  epsilon <- rnorm(n = n, mean = 0, sd = sigma)

  # Form y = beta0 + beta1*x1 + beta2*x2 + epsilon
  #   While rounding should not be done, it is here so that
  #   future displays of calculations are easier.
  y <- round(beta0 + beta1*x1 + beta2*x2 + epsilon, digits = 2)

  set1 <- data.frame(x1, x2, y)
  set1

  # Estimate model
  mod.fit <- lm(formula = y ~ x1 + x2, data = set1)
  summary(mod.fit)



################################################################################
# 3D plot

  # Hard-coded minimum value of y to be 0 in the plots below.
  #   This should be changed if a minimum value is less than 0.
  min(y)

  library(rgl)

  # Scatter plot
  plot3d(x = x1, y = x2, z = y, xlab = expression(x[1]), ylab = expression(x[2]),
    zlab = "y", col = "red", size = 16, zlim = c(0, max(y)))
  grid3d(side = c("x", "y", "z"), col = "lightgray")

  # Scatter plot with needles
  plot3d(x = x1, y = x2, z = y, xlab = expression(x[1]), ylab = expression(x[2]),
    zlab = "y", col = "red", size = 16, zlim = c(0, max(y)))
  plot3d(x = x1, y = x2, z = y, xlab = expression(x[1]), ylab = expression(x[2]),
    zlab = "y", col = "red", type = "h", add = TRUE)
  grid3d(side = c("x", "y", "z"), col = "lightgray")

  # Scatter plot with regression model
  #  NOTE: plot3d() is a generic function! There is a method function plot3.lm()
  #        that automatically does the scatter plot with the estimated regression
  #        plane on it.
  methods(generic.function = "plot3d") # Notice plot3d.lm()
  plot3d(x = mod.fit, xlab = expression(x[1]), ylab = expression(x[2]),
    zlab = "y", col = "red", plane.col = "blue", size = 16, zlim = c(0, max(y)))
  grid3d(side = c("x", "y", "z"), col = "lightgray")

  # Ideas from https://www.geeksforgeeks.org/r-language/3d-multiple-regression-graph-with-rgl-package-in-r/
  #   were used here



################################################################################
# Matrix calculations

  options(width = 60)

  X <- as.matrix(data.frame(Ones = 1, set1[,1:2]))
  Y <- as.matrix(set1[,3])

  beta.hat <- solve(t(X)%*%X)%*%t(X)%*%Y
  beta.hat

  # Hat matrix
  H <- X%*%solve(t(X)%*%X)%*%t(X)
  H[1:2,1:2] # Part of H

  H%*%Y
  predict(object = mod.fit, newdata = set1) # R finds the x1 and x2 in set1 and ignores y
  # predict(object = mod.fit) # this works too, R defaults to data frame used in lm()