Comment obtenir des coefficients de régression à partir d'une matrice de covariance de variance dans R?

Question

Je veux élaborer un exemple de régression multiple tout au long de l'utilisation de l'algèbre matricielle pour calculer les coefficients de régression.

#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 

Je peux voir ce que les valeurs finales devraient être (-.19,, 01) et ce que les matrices avant ce regard de calcul comme.

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

Mais je ne sais pas comment créer ces de la matrice de variance-covariance pour obtenir les coefficients en utilisant l'algèbre matricielle.

J'espère que vous pouvez aider.

Answer 1

Je peux voir que vous faites une régression centrée:

La réponse par sandipan est pas tout à fait ce que vous voulez, comme passe par l'équation normale habituelle pour estimer:

Il y a déjà un fil sur ce dernier: Solving normal equation gives different coefficients from using lm? Ici, je me concentre sur le premier.

---

En fait, vous êtes déjà assez proche. Vous avez obtenu la covariance mixte 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 

À partir de votre définition de E et F, vous savez que vous avez besoin des sous-matrices pour continuer. En fait, vous pouvez le faire matrice subsetting plutôt que imputant manuellement:

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 

l'estimation est juste , et vous pouvez le faire en R (lire ?solve):

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

---

Autre suggestions

• vous pouvez trouver des moyens de colonne par colMeans, au lieu d'une multiplication matricielle (lire ?colMeans); Vous pouvez effectuer le centrage en utilisant sweep (lire ?sweep);
• utiliser crossprod(D) que t(D) %*% D (lire ?crossprod).

Voici une session que je ferais:

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

Probablement vous voulez quelque chose comme ceci:

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