如何從R中的方差協方差矩陣得到迴歸係數？

Question

我想通過使用矩陣代數來計算迴歸係數，找出一個多元迴歸的例子。

#create vectors -- these will be our columns 
y <- c(3,3,2,4,4,5,2,3,5,3) 
x1 <- c(2,2,4,3,4,4,5,3,3,5) 
x2 <- c(3,3,4,4,3,3,4,2,4,4) 

#create matrix from vectors 
M <- cbind(y,x1,x2) 
k <- ncol(M) #number of variables 
n <- nrow(M) #number of subjects 

#create means for each column 
M_mean <- matrix(data=1, nrow=n) %*% cbind(mean(y),mean(x1),mean(x2)); M_mean 

#creates a difference matrix which gives deviation scores 
D <- M - M_mean; D 

#creates the covariance matrix, the sum of squares are in the diagonal and the sum of cross products are in the off diagonals. 
C <- t(D) %*% D; C 

我可以看到最終值應該是什麼（-.19，-.01）以及這個計算之前的矩陣是什麼樣的。

E<-matrix(c(10.5,3,3,4.4),nrow=2,ncol=2) 
F<-matrix(c(-2,-.6),nrow=2,ncol=1) 

但我不知道如何從方差協方差矩陣創建這些矩陣來使用矩陣代數來得到係數。

希望你能幫上忙。

Answer 1

我可以看到你在做中心迴歸：

通過sandipan答案是不是你所需的東西，因爲它按照通常的正常方程式估算：

後者已經有一個線程：Solving normal equation gives different coefficients from using lm?這裏我把重點放在前者。

---

其實你已經相當接近。你已經得到的混合協方差C：

#  y x1 x2 
#y 10.4 -2.0 -0.6 
#x1 -2.0 10.5 3.0 
#x2 -0.6 3.0 4.4 

從你的E和F定義，你知道你需要的子矩陣進行。事實上，你可以做矩陣子集，而不是手動歸咎於：

E <- C[2:3, 2:3] 

#  x1 x2 
#x1 10.5 3.0 
#x2 3.0 4.4 

F <- C[2:3, 1, drop = FALSE] ## note the `drop = FALSE` 

#  y 
#x1 -2.0 
#x2 -0.6 

然後估計是剛剛，並在R中可以做（讀?solve）：

c(solve(E, F)) ## use `c` to collapse matrix into a vector 
# [1] -0.188172043 -0.008064516 

---

其他建議

• 你可以找到列的意思是colMeans，而不是矩陣乘法（讀取?colMeans）;
• 您可以使用sweep（閱讀?sweep）執行對中;
• 使用crossprod(D)而不是t(D) %*% D（讀取?crossprod）。

這裏是一個會議上，我會做：

y <- c(3,3,2,4,4,5,2,3,5,3) 
x1 <- c(2,2,4,3,4,4,5,3,3,5) 
x2 <- c(3,3,4,4,3,3,4,2,4,4) 

M <- cbind(y,x1,x2) 
M_mean <- colMeans(M) 
D <- sweep(M, 2, M_mean) 
C <- crossprod(D) 

E <- C[2:3, 2:3] 
F <- C[2:3, 1, drop = FALSE] 
c(solve(E, F)) 
# [1] -0.188172043 -0.008064516 

Answer 2

也許你想是這樣的：

X <- cbind(1, x1, x2) 
C <- t(X) %*% X # No need of centering the columns with means 
D <- t(X) %*% y 

coef <- t(solve(C) %*% D) 
coef 
#      x1   x2 
# [1,] 4.086022 -0.188172 -0.008064516 

lm(y~x1+x2)$coef # coefficients with R lm() 
# (Intercept)   x1   x2 
# 4.086021505 -0.188172043 -0.008064516