title: Multiple regression analysis in R using QR decomposition
body: <p>I am trying to write a function for solving multiple regression using QR decomposition. Input: y vector and X matrix; output: b, e, R^2. So far I`ve got this and am terribly stuck; I think I have made everything way too complicated:</p>
  
  <pre><code>QR.regression &lt;- function(y, X) {
  X &lt;- as.matrix(X)
  y &lt;- as.vector(y)
  p &lt;- as.integer(ncol(X))
  if (is.na(p)) stop("ncol(X) is invalid")
  n &lt;- as.integer(nrow(X))
  if (is.na(n)) stop("nrow(X) is invalid")
  nr &lt;- length(y)
  nc &lt;- NCOL(X)
  
  # Householder 
  for (j in seq_len(nc)) {
  id &lt;- seq.int(j, nr)
  sigma &lt;- sum(X[id, j]^2)
  s &lt;- sqrt(sigma)
  diag_ej &lt;- X[j, j]
  gamma &lt;- 1.0 / (sigma + abs(s * diag_ej))
  kappa &lt;- if (diag_ej &lt; 0) s else -s
  X[j,j] &lt;- X[j, j] - kappa
  if (j &lt; nc)
  for (k in seq.int(j+1, nc)) {
  yPrime &lt;- sum(X[id,j] * X[id,k]) * gamma
  X[id,k] &lt;- X[id,k] - X[id,j] * yPrime
  }
  yPrime &lt;- sum(X[id,j] * y[id]) * gamma
  y[id] &lt;- y[id] - X[id,j] * yPrime
  X[j,j] &lt;- kappa
  } # end of Householder transformation
  
  rss &lt;- sum(y[seq.int(nc+1, nr)]^2)  # residuals sum of squares
  e &lt;- rss/nr
  e &lt;- mean(residuals(QR.regression)^2)
  beta &lt;- 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] &lt;- 0
  Rsq &lt;- (X[1:nc,1:nc])^2
  return(list(Rsq=Rsq, y = y, beta = beta, e = e))
  }
  
  
  UPDATE:
  my.QR &lt;- function(y, X) {
  X &lt;- as.matrix(X)
  y &lt;- as.vector(y)
  p &lt;- as.integer(ncol(X))
  if (is.na(p)) stop("ncol(X) is invalid")
  n &lt;- as.integer(nrow(X))
  if (is.na(n)) stop("nrow(X) is invalid")
  qr.X &lt;- qr(X)
  b &lt;- solve(t(X) %*% X, t(X) %*% y)
  e &lt;- as.vector(y - X %*% beta) #e
  R2 &lt;- (X[1:p, 1:p])^2
  return(list(b = b, e= e, R2 = R2 ))
  }
  
  X &lt;- matrix(c(1,2,3,4,5,6), nrow = 2, ncol = 3)
  y &lt;- c(1,2,3,4)
  my.QR(X, y)
  </code></pre>