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

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

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β+ϵ, E(ϵ)=0, Cov(ϵ)=σ2In 回归系数估计：βˆ=(X′X)−1X′y 拟合值：yˆ=Xβ 残差估计:ϵˆ=y−yˆ 残差的
d3d1
平方：σ2=∑ϵˆ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.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] #原始值

112.583

116.033

119.483

122.933

126.383

115

117

120

123

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))

5.93694404890295

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