Множественный регрессионный анализ в R с использованием QR-разложения

Question

Я пытаюсь написать функцию для решения множественной регрессии с помощью QR-разложения. Вход: вектор y и матрица X; выход: б, е, р ^2. До сих пор я получил это и ужасно застрял; Я думаю, я сделал все слишком сложным:

QR.regression <- function(y, X) {
X <- as.matrix(X)
y <- as.vector(y)
p <- as.integer(ncol(X))
if (is.na(p)) stop("ncol(X) is invalid")
n <- as.integer(nrow(X))
if (is.na(n)) stop("nrow(X) is invalid")
nr <- length(y)
nc <- NCOL(X)

# Householder 
for (j in seq_len(nc)) {
id <- seq.int(j, nr)
sigma <- sum(X[id, j]^2)
s <- sqrt(sigma)
diag_ej <- X[j, j]
gamma <- 1.0 / (sigma + abs(s * diag_ej))
kappa <- if (diag_ej < 0) s else -s
X[j,j] <- X[j, j] - kappa
if (j < nc)
for (k in seq.int(j+1, nc)) {
yPrime <- sum(X[id,j] * X[id,k]) * gamma
X[id,k] <- X[id,k] - X[id,j] * yPrime
}
yPrime <- sum(X[id,j] * y[id]) * gamma
y[id] <- y[id] - X[id,j] * yPrime
X[j,j] <- kappa
} # end of Householder transformation

rss <- sum(y[seq.int(nc+1, nr)]^2)  # residuals sum of squares
e <- rss/nr
e <- mean(residuals(QR.regression)^2)
beta <- solve(t(X) %*% X, t(X) %*% y)
for (i in seq_len(ncol(X))) # set zeros in the lower triangular side of X
X[seq.int(i+1, nr),i] <- 0
Rsq <- (X[1:nc,1:nc])^2
return(list(Rsq=Rsq, y = y, beta = beta, e = e))
}


UPDATE:
my.QR <- function(y, X) {
X <- as.matrix(X)
y <- as.vector(y)
p <- as.integer(ncol(X))
if (is.na(p)) stop("ncol(X) is invalid")
n <- as.integer(nrow(X))
if (is.na(n)) stop("nrow(X) is invalid")
qr.X <- qr(X)
b <- solve(t(X) %*% X, t(X) %*% y)
e <- as.vector(y - X %*% beta) #e
R2 <- (X[1:p, 1:p])^2
return(list(b = b, e= e, R2 = R2 ))
}

X <- matrix(c(1,2,3,4,5,6), nrow = 2, ncol = 3)
y <- c(1,2,3,4)
my.QR(X, y)

Answer 1

Все зависит от того, сколько встроенных средств R вам разрешено использовать для решения этой проблемы. я уже знаю, что lm не допускается, так что вот остальная часть истории.

---

Если вам разрешено использовать любые другие процедуры, кромеlm

Тогда вы можете просто использовать lm.fit, .lm.fit или же lsfit для решения обычных наименьших квадратов на основе QR.

lm.fit(X, y)
.lm.fit(X, y)
lsfit(X, y, intercept = FALSE)

Среди тех, .lm.fit является самым легким, в то время как lm.fit а также lsfit очень похожи Вот что мы можем сделать через .lm.fit:

f1 <- function (X, y) {
  z <- .lm.fit(X, y)
  RSS <- crossprod(z$residuals)[1]
  TSS <- crossprod(y - mean(y))[1]
  R2 <- 1 - RSS / TSS
  list(coefficients = z$coefficients, residuals = z$residuals, R2 = R2)
  }

В вопросе вашего одноклассника: функция Toy R для решения обычных наименьших квадратов путем разложения по сингулярным числам, я уже использовал это для проверки правильности подхода SVD.

---

Если вам не разрешено использовать встроенную процедуру QR-факторизации Rqr.default

Если .lm.fit не разрешено, но qr.default есть, то это тоже не так сложно.

f2 <- function (X, y) {
  ## QR factorization `X = QR`
  QR <- qr.default(X)
  ## After rotation of `X` and `y`, solve upper triangular system `Rb = Q'y` 
  b <- backsolve(QR$qr, qr.qty(QR, y))
  ## residuals
  e <- as.numeric(y - X %*% b)
  ## R-squared
  RSS <- crossprod(e)[1]
  TSS <- crossprod(y - mean(y))[1]
  R2 <- 1 - RSS / TSS
  ## multiple return
  list(coefficients = b, residuals = e, R2 = R2)
  }

Если вы хотите получить дисперсию-ковариацию оценочных коэффициентов, следуйте инструкциям Как рассчитать дисперсию оценки наименьших квадратов, используя QR-разложение в R?,

---

Если вам даже не разрешено использоватьqr.default

Затем мы должны написать QR-разложение сами. Это дает функция написания QR-кода домохозяйства в коде R.

Используя функцию myqr там мы можем написать

f3 <- function (X, y) {
  ## our own QR factorization
  ## complete Q factor is not required
  QR <- myqr(X, complete = FALSE)
  Q <- QR$Q
  R <- QR$R
  ## rotation of `y`
  Qty <- as.numeric(crossprod(Q, y))
  ## solving upper triangular system
  b <- backsolve(R, Qty)
  ## residuals
  e <- as.numeric(y - X %*% b)
  ## R-squared
  RSS <- crossprod(e)[1]
  TSS <- crossprod(y - mean(y))[1]
  R2 <- 1 - RSS / TSS
  ## multiple return
  list(coefficients = b, residuals = e, R2 = R2)
  }

f3 не очень эффективно, так как мы сформировали Q явно, хотя это тонкийQ фактор. В принципе, мы должны вращаться y наряду с QR-факторизацией Xтаким образом Q не должно быть сформировано.

---

Если вы хотите исправить свой существующий код

Это требует определенных усилий по отладке, поэтому может занять некоторое время. Я сделаю еще один ответ по этому поводу позже.