Как получить коэффициенты регрессии из матрицы ковариации дисперсии в 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