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 n*p matrix (p > 1)
  z_error1 = lm_fit(x = x, y = y2)
  expect_error(lm_anova(z_error1))
  # use methods other than "qr"
  z_error2 = lm_fit(x = x, y = y1, method = "inv")
  expect_error(lm_anova(z_error2), "only support 'qr' 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, add.intercept = TRUE) # add intercept, using method "qr"

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

  result1 = lm_anova(z1)
  result2 = lm_anova(z2)
  test.a1 = anova(test1)
  test.a2 = anova(test2)

  expect_equal(result1$Df, test.a1$Df)
  expect_lt(sum(abs(result1$`Sum Sq` - test.a1$`Sum Sq`)), 1e-10)
  expect_lt(sum(abs(result1$`Mean Sq` - test.a1$`Mean Sq`)), 1e-10)
  expect_lt(abs(result1$`F value` - test.a1$`F value`[1]), 1e-10)
  expect_lt(abs(result1$`Pr(>F)` - test.a1$`Pr(>F)`[1]), 1e-10)

  expect_equal(result2$Df, test.a2$Df)
  expect_lt(sum(abs(result2$`Sum Sq` - test.a2$`Sum Sq`)), 1e-10)
  expect_lt(sum(abs(result2$`Mean Sq` - test.a2$`Mean Sq`)), 1e-10)
  expect_lt(abs(result2$`F value` - test.a2$`F value`[1]), 1e-10)
  expect_lt(abs(result2$`Pr(>F)` - test.a2$`Pr(>F)`[1]), 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 n*p matrix (p > 1)
  z_error1 = lm_fit(x = x, y = y2)
  expect_error(lm_anova(z_error1))
  # use methods other than "qr"
  z_error2 = lm_fit(x = x, y = y1, method = "inv")
  expect_error(lm_anova(z_error2), "only support 'qr' 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, add.intercept = TRUE) # add intercept, using method "qr"

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

  result1 = lm_anova(z1)
  result2 = lm_anova(z2)
  test.a1 = anova(test1)
  test.a2 = anova(test2)

  expect_equal(result1$Df, test.a1$Df)
  expect_lt(sum(abs(result1$`Sum Sq` - test.a1$`Sum Sq`)), 1e-10)
  expect_lt(sum(abs(result1$`Mean Sq` - test.a1$`Mean Sq`)), 1e-10)
  expect_lt(abs(result1$`F value` - test.a1$`F value`[1]), 1e-10)
  expect_lt(abs(result1$`Pr(>F)` - test.a1$`Pr(>F)`[1]), 1e-10)

  expect_equal(result2$Df, test.a2$Df)
  expect_lt(sum(abs(result2$`Sum Sq` - test.a2$`Sum Sq`)), 1e-10)
  expect_lt(sum(abs(result2$`Mean Sq` - test.a2$`Mean Sq`)), 1e-10)
  expect_lt(abs(result2$`F value` - test.a2$`F value`[1]), 1e-10)
  expect_lt(abs(result2$`Pr(>F)` - test.a2$`Pr(>F)`[1]), 1e-10)

#test case 3: real data
  x = cbind(mtcars$hp, mtcars$wt, mtcars$hp * mtcars$wt)
  y = mtcars$mpg
  z = lm_fit(x, y, add.intercept = TRUE)
  test = lm(mpg ~ x, data = mtcars)

  result = lm_anova(z)
  test.a = anova(test)

  expect_equal(result$Df, test.a$Df)
  expect_lt(sum(abs(result$`Sum Sq` - test.a$`Sum Sq`)), 1e-10)
  expect_lt(sum(abs(result$`Mean Sq` - test.a$`Mean Sq`)), 1e-10)
  expect_lt(abs(result$`F value` - test.a$`F value`[1]), 1e-10)
  expect_lt(abs(result$`Pr(>F)` - test.a$`Pr(>F)`[1]), 1e-10)
})