lm.hw4: Linear Regression

test_that("multiplication works", {

#test case 1: simulation data - small sample size
  n = 10; p = 3; q = 2;
  x = matrix(rnorm(n * p), n, p) # no intercept
  y1 = rnorm(n)
  y2 = matrix(rnorm(n * q), n, q)

  # invalid cases
  # input response variable as a vector or something other than a matrix
  x_error1 = c(1, 2, 3, 4, 5)
  expect_error(lm_fit(x_error1, y1), "'x' must be a matrix")
  # input variables have different dimensions in rows.
  y_error1 = rnorm(9)
  y_error2 = matrix(rnorm(9 * q), 9, q)
  expect_error(lm_fit(x, y_error1), "incompatible dimensions")
  expect_error(lm_fit(x, y_error2), "incompatible dimensions")
  # matrix x doesn't have full column rank (X^TX is a singular matrix).
  x_error2 = matrix(rep(1,20), 10, 2)
  expect_error(lm_fit(x_error2, y1), "singular fit encountered")
  # use methods other than 'qr' and 'inv'
  expect_error(lm_fit(x,y1, method = "error"), "only support 'qr' and 'inv' method")


  # input response variable as a vector or a n*1 matrix
  z1 = lm_fit(x = x, y = y1) # no intercept, using method "qr"
  z2 = lm_fit(x = x, y = y1, method ="inv") # no intercept, using method "inv"
  z3 = lm_fit(x = x, y = y1, add.intercept = TRUE) # add intercept, using method "qr"
  z4 = lm_fit(x = x, y = y1, method = "inv", add.intercept = TRUE) # add intercept, using method "inv"

  test1 = stats::lm(y1 ~ 0 + x)
  test2 = stats::lm(y1 ~ x)

  expect_lt(sum(abs(z1$coefficients - test1$coefficients)), 1e-10)
  expect_lt(sum(abs(z1$residuals - test1$residuals)), 1e-10)
  expect_lt(sum(abs(z1$fitted.values - test1$fitted.values)), 1e-10)

  expect_lt(sum(abs(z2$coefficients - test1$coefficients)), 1e-10)
  expect_lt(sum(abs(z2$residuals - test1$residuals)), 1e-10)
  expect_lt(sum(abs(z2$fitted.values - test1$fitted.values)), 1e-10)

  expect_lt(sum(abs(z3$coefficients - test2$coefficients)), 1e-10)
  expect_lt(sum(abs(z3$residuals - test2$residuals)), 1e-10)
  expect_lt(sum(abs(z3$fitted.values - test2$fitted.values)), 1e-10)

  expect_lt(sum(abs(z4$coefficients - test2$coefficients)), 1e-10)
  expect_lt(sum(abs(z4$residuals - test2$residuals)), 1e-10)
  expect_lt(sum(abs(z4$fitted.values - test2$fitted.values)), 1e-10)


  # input response variable as a n* p matrix (p > 1)
  z5 = lm_fit(x = x, y = y2) # no intercept, using method "qr"
  z6 = lm_fit(x = x, y = y2, method ="inv") # no intercept, using method "inv"
  z7 = lm_fit(x = x, y = y2, add.intercept = TRUE) # add intercept, using method "qr"
  z8 = lm_fit(x = x, y = y2, method = "inv", add.intercept = TRUE) # add intercept, using method "inv"

  test3 = stats::lm(y2 ~ 0 + x)
  test4 = stats::lm(y2 ~ x)

  expect_lt(sum(abs(z5$coefficients - test3$coefficients)), 1e-10)
  expect_lt(sum(abs(z5$residuals - test3$residuals)), 1e-10)
  expect_lt(sum(abs(z5$fitted.values - test3$fitted.values)), 1e-10)

  expect_lt(sum(abs(z6$coefficients - test3$coefficients)), 1e-10)
  expect_lt(sum(abs(z6$residuals - test3$residuals)), 1e-10)
  expect_lt(sum(abs(z6$fitted.values - test3$fitted.values)), 1e-10)

  expect_lt(sum(abs(z7$coefficients - test4$coefficients)), 1e-10)
  expect_lt(sum(abs(z7$residuals - test4$residuals)), 1e-10)
  expect_lt(sum(abs(z7$fitted.values - test4$fitted.values)), 1e-10)

  expect_lt(sum(abs(z8$coefficients - test4$coefficients)), 1e-10)
  expect_lt(sum(abs(z8$residuals - test4$residuals)), 1e-10)
  expect_lt(sum(abs(z8$fitted.values - test4$fitted.values)), 1e-10)


#test case 2: simulation data - big sample size
  n = 10000; p = 100; q = 2;
  x = matrix(rnorm(n * p), n, p) # no intercept
  y1 = rnorm(n)
  y2 = matrix(rnorm(n * q), n, q)

  # invalid cases
  # input response variable as a vector or something other than a matrix
  x_error1 = c(1, 2, 3, 4, 5)
  expect_error(lm_fit(x_error1, y1), "'x' must be a matrix")
  # input variables have different dimensions in rows.
  y_error1 = rnorm(9)
  y_error2 = matrix(rnorm(9 * q), 9, q)
  expect_error(lm_fit(x, y_error1), "incompatible dimensions")
  expect_error(lm_fit(x, y_error2), "incompatible dimensions")
  # matrix x doesn't have full column rank (X^TX is a singular matrix).
  x_error2 = matrix(rep(1,20000), 10000, 2)
  expect_error(lm_fit(x_error2, y1), "singular fit encountered")
  # use methods other than 'qr' and 'inv'
  expect_error(lm_fit(x,y1, method = "error"), "only support 'qr' and 'inv' method")


  # input response variable as a vector or a n*1 matrix
  z1 = lm_fit(x = x, y = y1) # no intercept, using method "qr"
  z2 = lm_fit(x = x, y = y1, method ="inv") # no intercept, using method "inv"
  z3 = lm_fit(x = x, y = y1, add.intercept = TRUE) # add intercept, using method "qr"
  z4 = lm_fit(x = x, y = y1, method = "inv", add.intercept = TRUE) # add intercept, using method "inv"

  test1 = stats::lm(y1 ~ 0 + x)
  test2 = stats::lm(y1 ~ x)

  expect_lt(sum(abs(z1$coefficients - test1$coefficients)), 1e-10)
  expect_lt(sum(abs(z1$residuals - test1$residuals)), 1e-10)
  expect_lt(sum(abs(z1$fitted.values - test1$fitted.values)), 1e-10)

  expect_lt(sum(abs(z2$coefficients - test1$coefficients)), 1e-10)
  expect_lt(sum(abs(z2$residuals - test1$residuals)), 1e-10)
  expect_lt(sum(abs(z2$fitted.values - test1$fitted.values)), 1e-10)

  expect_lt(sum(abs(z3$coefficients - test2$coefficients)), 1e-10)
  expect_lt(sum(abs(z3$residuals - test2$residuals)), 1e-10)
  expect_lt(sum(abs(z3$fitted.values - test2$fitted.values)), 1e-10)

  expect_lt(sum(abs(z4$coefficients - test2$coefficients)), 1e-10)
  expect_lt(sum(abs(z4$residuals - test2$residuals)), 1e-10)
  expect_lt(sum(abs(z4$fitted.values - test2$fitted.values)), 1e-10)


  # input response variable as a n* p matrix (p > 1)
  z5 = lm_fit(x = x, y = y2) # no intercept, using method "qr"
  z6 = lm_fit(x = x, y = y2, method ="inv") # no intercept, using method "inv"
  z7 = lm_fit(x = x, y = y2, add.intercept = TRUE) # add intercept, using method "qr"
  z8 = lm_fit(x = x, y = y2, method = "inv", add.intercept = TRUE) # add intercept, using method "inv"

  test3 = stats::lm(y2 ~ 0 + x)
  test4 = stats::lm(y2 ~ x)

  expect_lt(sum(abs(z5$coefficients - test3$coefficients)), 1e-10)
  expect_lt(sum(abs(z5$residuals - test3$residuals)), 1e-10)
  expect_lt(sum(abs(z5$fitted.values - test3$fitted.values)), 1e-10)

  expect_lt(sum(abs(z6$coefficients - test3$coefficients)), 1e-10)
  expect_lt(sum(abs(z6$residuals - test3$residuals)), 1e-10)
  expect_lt(sum(abs(z6$fitted.values - test3$fitted.values)), 1e-10)

  expect_lt(sum(abs(z7$coefficients - test4$coefficients)), 1e-10)
  expect_lt(sum(abs(z7$residuals - test4$residuals)), 1e-10)
  expect_lt(sum(abs(z7$fitted.values - test4$fitted.values)), 1e-10)

  expect_lt(sum(abs(z8$coefficients - test4$coefficients)), 1e-10)
  expect_lt(sum(abs(z8$residuals - test4$residuals)), 1e-10)
  expect_lt(sum(abs(z8$fitted.values - test4$fitted.values)), 1e-10)

#test case 3: real data
  x = cbind(mtcars$hp, mtcars$wt, mtcars$hp * mtcars$wt)
  y = mtcars$mpg
  z1 = lm_fit(x, y, add.intercept = TRUE)
  z2 = lm_fit(x, y, add.intercept = TRUE, method = "inv")
  test = lm(mpg ~ hp + wt + hp:wt, data = mtcars)

  expect_lt(sum(abs(z1$coefficients - test$coefficients)), 1e-10)
  expect_lt(sum(abs(z1$residuals - test$residuals)), 1e-10)
  expect_lt(sum(abs(z1$fitted.values - test$fitted.values)), 1e-10)

  expect_lt(sum(abs(z2$coefficients - test$coefficients)), 1e-10)
  expect_lt(sum(abs(z2$residuals - test$residuals)), 1e-10)
  expect_lt(sum(abs(z2$fitted.values - test$fitted.values)), 1e-10)

})