如何从 R 中的方差协方差矩阵中获取回归系数?
How do I get regression coefficients from a variance covariance matrix in R?
我想通过使用矩阵代数计算回归系数来计算一个多元回归的例子。
#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)
但我不确定如何从方差-协方差矩阵创建这些以使用矩阵代数获得系数。
希望能帮到你。
可能你想要这样的东西:
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
我看到你在做中心回归:
sandipan的答案不是你想要的,因为它通过通常的正规方程来估计:
已经有关于后者的跟帖:这里重点说一下前者。
其实你已经很接近了。您已获得混合协方差 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
的定义,您知道需要 sub-matrices 才能继续。事实上,你可以做矩阵子集而不是手动插补:
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
我想通过使用矩阵代数计算回归系数来计算一个多元回归的例子。
#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)
但我不确定如何从方差-协方差矩阵创建这些以使用矩阵代数获得系数。
希望能帮到你。
可能你想要这样的东西:
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
我看到你在做中心回归:
sandipan的答案不是你想要的,因为它通过通常的正规方程来估计:
已经有关于后者的跟帖:
其实你已经很接近了。您已获得混合协方差 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
的定义,您知道需要 sub-matrices 才能继续。事实上,你可以做矩阵子集而不是手动插补:
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
那么估计就是?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