用矩阵的方法计算回归模型参数

用矩阵的方法计算回归分析参数

1.1 数据来源：来源R语言默认的数据集women

这是一个描述女性身高和体重的数据，我们以height为X变量（自变量），以weight为Y变量（因变量），进行模型的计算。

计算方法参考：https://stats.idre.ucla.edu/r/library/r-library-matrices-and-matrix-computations-in-r/

1.2 查看数据

data(women)

head(women)

height	weight
58	115
59	117
60	120
61	123
62	126
63	129

1.3 理论模型

$ y = X\beta + \epsilon,\ E(\epsilon) = 0 ,\ Cov(\epsilon) = \sigma^2I_n \
回归系数估计： \widehat{\beta} = (X’X)^{-1}X’y \
拟合值：\widehat{y} = X\beta \
残差估计: \widehat{\epsilon}= y - \widehat{y} \
残差的平方：\sigma^2 = \sum{\widehat{\epsilon}^2}$

2.1 构建X矩阵

X <- as.matrix(cbind(1, women$height))  

n <- dim(X)[1]

p <- dim(X)[2]

head(X)

1	58
1	59
1	60
1	61
1	62
1	63

2.2 构建y矩阵

y <- matrix(women$weight,,1)

head(y)

115

117

120

123

126

129

2.3 计算 $\beta$ 参数

第一种方法，是直接根据公式计算：

beta.hat <- solve(t(X) %*% X) %*% t(X) %*% y

beta.hat

-87.51667

3.45000

第二种方法，是用crossprod函数，在计算大数据时有优势

beta.hat1 <- solve(crossprod(X), crossprod(X,y)) # solve(A,B) == solve(A)%*%B

beta.hat1

-87.51667

3.45000

2.4 计算拟合值Fitted_value

y.hat <- X %*% beta.hat

round(y.hat[1:5, 1],3) # 拟合值

y[1:5, 1] #原始值

1. 112.583
2. 116.033
3. 119.483
4. 122.933
5. 126.383

1. 115
2. 117
3. 120
4. 123
5. 126

plot(y.hat,y)

2.5 计算残差值（residual）

residual <- y - y.hat

head(residual)

2.41666667

0.96666667

0.51666667

0.06666667

-0.38333333

-0.83333333

2.6 计算残差方差组分（残差的方差）

sigma2 <- sum((y - y.hat)^2)/(n - p)

sigma2

2.32564102564103

2.7 计算方差协方差矩阵（var.beta）

v <- solve(t(X) %*% X) * sigma2
v

35.247305	-0.539880952
-0.539881	0.008305861

2.8 计算参数的标准误

#standard errors of the parameter estimates
sqrt(diag(v))

1. 5.93694404890295
2. 0.0911364954662

2.9 对参数进行T检验，计算T值

t.values <- beta.hat/sqrt(diag(v))

t.values

-14.74103

37.85531

2.10 根据T值，计算p值

2 * (1 - pt(abs(t.values), n - p))

1.711082e-09

1.088019e-14

3. 用lm函数和矩阵得到的结果对比

3.1 对比参数估计

mod <- lm(weight ~ height,data=women)
summary(mod)$coef

	Estimate	Std. Error	t value	Pr(>|t|)
(Intercept)	-87.51667	5.9369440	-14.74103	1.711082e-09
height	3.45000	0.0911365	37.85531	1.090973e-14

beta.hat

-87.51667

3.45000

3.2 拟合值

head(fitted(mod))

1

112.583333333333

2

116.033333333333

3

119.483333333333

4

122.933333333333

5

126.383333333333

6

129.833333333333

head(y.hat)

112.5833

116.0333

119.4833

122.9333

126.3833

129.8333

3.3 残差值

head(residuals(mod))

1

2.41666666666676

2

0.966666666666627

3

0.516666666666641

4

0.0666666666666444

5

-0.383333333333353

6

-0.833333333333348

head(residual)

2.41666667

0.96666667

0.51666667

0.06666667

-0.38333333

-0.83333333

summary(mod)

Call:
lm(formula = weight ~ height, data = women)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.7333 -1.1333 -0.3833  0.7417  3.1167 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -87.51667    5.93694  -14.74 1.71e-09 ***
height        3.45000    0.09114   37.85 1.09e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.525 on 13 degrees of freedom
Multiple R-squared:  0.991,	Adjusted R-squared:  0.9903 
F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14

sigma2

2.32564102564103

sqrt(sigma2)

1.52500525429948