Реализация функции в R для нахождения коэффициентов модели линейной регрессии с нуля

Question

# but cannot handle categorical variables
my_lm <- function(explanatory_matrix, response_vec) {
  exp_mat <- as.matrix(explanatory_matrix)
  intercept <- rep(1, nrow(exp_mat))
  exp_mat <- cbind(exp_mat, intercept)
  solve(t(exp_mat) %*% exp_mat) %*% (t(exp_mat) %*% response_vec)
}

Приведенный выше код не будет работать, если в файле объяснительной_матрицы есть категориальные переменные.Как я могу это реализовать?

Answer 1

Для категориальной обработки вы должны использовать однократное кодирование.Сделайте что-то вроде

  formula <- dep_var ~ indep_var
  exp_mat <- model.matrix(formula, explanatory_matrix)
  solve(t(exp_mat) %*% exp_mat) %*% (t(exp_mat) %*% response_vec)

Answer 2

Вот пример для набора данных с одной категориальной переменной:

set.seed(123)

x <- 1:10
a <- 2
b <- 3

y <- a*x + b + rnorm(10)
# categorical variable
x2 <- sample(c("A", "B"), 10, replace = T)
# one-hot encoding
x2 <- as.integer(d$x2 == "A")


xm <- matrix(c(x, x2, rep(1, length(x))), ncol = 3, nrow = 10)
ym <- matrix(y, ncol = 1, nrow = 10)

beta_hat <- MASS::ginv(t(xm) %*% xm) %*% t(xm) %*% ym
beta_hat

Это дает (обратите внимание на порядок коэффициентов - он соответствует порядку столбцов предиктора):

      [,1]
[1,]  1.9916754
[2,] -0.7594809
[3,]  3.2723071

, что соответствует выходу lm:

d <- data.frame(x = x,
                x2 = x2,
                y = y)
lm(y ~ ., data = d)

Выход

# Call:
#   lm(formula = y ~ ., data = d)
# 
# Coefficients:
#   (Intercept)            x           x2  
# 3.2723       1.9917      -0.7595