思路
利用主成分分析，从变量中提取主成分，进而利用主成分做回归
关于主成分分析的基本步骤，可参考> https://www.cnblogs.com/hznudmh/p/16324990.html>
步骤
Step1. 原始数据标准化
Step2. 计算自变量之间的相关矩阵
Step3. 检验是否适合做主成分分析
Step4. 主成分提取
Step5. 计算主成分得分
Step6. 做主成分回归
实例
首先我们导入数据
ce <- read.table(file.choose(), header = T)
ce
## X1 X2 X3 X4 Y
## 1 7 26 6 60 78.5
## 2 1 29 15 52 74.3
## 3 11 56 8 20 104.3
## 4 11 31 8 47 87.6
## 5 7 52 6 33 95.9
## 6 11 55 9 22 109.2
## 7 3 71 17 6 102.7
## 8 1 31 22 44 72.5
## 9 2 54 18 22 93.1
## 10 21 47 4 26 115.9
## 11 1 40 23 34 83.8
## 12 11 66 9 12 113.3
## 13 10 68 8 12 109.4
Step1. 原始数据标准化
scale(ce) -> ce.s
ce.s
## X1 X2 X3 X4 Y
## 1 -0.07846099 -1.42368840 -0.9007209 1.7923096 -1.12492615
## 2 -1.09845379 -1.23089726 0.5044037 1.3143603 -1.40411237
## 3 0.60153422 0.50422298 -0.5884710 -0.5974365 0.59007490
## 4 0.60153422 -1.10236984 -0.5884710 1.0156421 -0.52002268
## 5 -0.07846099 0.24716813 -0.9007209 0.1792310 0.03170246
## 6 0.60153422 0.43995926 -0.4323460 -0.4779492 0.91579215
## 7 -0.75845619 1.46817866 0.8166536 -1.4338477 0.48371824
## 8 -1.09845379 -1.10236984 1.5972783 0.8364111 -1.52376360
## 9 -0.92845499 0.37569555 0.9727785 -0.4779492 -0.15442168
## 10 2.30152223 -0.07415044 -1.2129708 -0.2389746 1.36116064
## 11 -1.09845379 -0.52399643 1.7534033 0.2389746 -0.77261973
## 12 0.60153422 1.14686010 -0.4323460 -1.0753857 1.18833108
## 13 0.43153542 1.27538753 -0.5884710 -1.0753857 0.92908673
## attr(,"scaled:center")
## X1 X2 X3 X4 Y
## 7.461538 48.153846 11.769231 30.000000 95.423077
## attr(,"scaled:scale")
## X1 X2 X3 X4 Y
## 5.882394 15.560881 6.405126 16.738180 15.043723
Step2. 计算自变量之间的相关矩阵
cor(ce[,1:4]) -> ce.cor # 使自变量之间的相关矩阵，因此选前四列
Step3. 检验是否适合做主成分分析
- KMO法
KMO(ce.cor)
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = ce.cor)
## Overall MSA = 0.24
## MSA for each item =
## X1 X2 X3 X4
## 0.22 0.26 0.20 0.26
发现 \(MSA = 0.24 < 0.7\)
- bartlett法
cortest.bartlett(ce.cor, n = 13) # n为样本个数，即有几行
## $chisq
## [1] 67.28248
##
## $p.value
## [1] 1.473373e-12
##
## $df
## [1] 6
发现 \(p << 0.05\).
虽然 KMO 检验未通过，但 bartlett 检验通过了，我们主要看bartlett检验，因此认为样本适合做主成分.
Step4. 主成分提取
princomp(ce[,1:4], cor = T) -> ce.pca
summary(ce.pca, loadings = T)
## Importance of components:
## Comp.1 Comp.2 Comp.3
## Standard deviation 1.495227 1.2554147 0.43197934
## Proportion of Variance 0.558926 0.3940165 0.04665154
## Cumulative Proportion 0.558926 0.9529425 0.99959406
## Comp.4
## Standard deviation 0.0402957285
## Proportion of Variance 0.0004059364
## Cumulative Proportion 1.0000000000
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4
## X1 0.476 0.509 0.676 0.241
## X2 0.564 -0.414 -0.314 0.642
## X3 -0.394 -0.605 0.638 0.268
## X4 -0.548 0.451 -0.195 0.677
因此可以选两个或三个主成分，这里我们选三个. 那么保留了 99.96% 的样本信息.
Step5. 计算主成分得分
as.matrix(ce.s[, 1:4]) %*% ce.pca$loadings[, 1:3] -> Z1
Z1
## Comp.1 Comp.2 Comp.3
## 1 -1.4672378 1.9030357 -0.53000021
## 2 -2.1358287 0.2383537 -0.29018631
## 3 1.1298705 0.1838772 -0.01071260
## 4 -0.6598955 1.5767742 0.17920364
## 5 0.3587646 0.4835379 -0.74012228
## 6 0.9666396 0.1699440 0.08570239
## 7 0.9307051 -2.1348165 -0.17298608
## 8 -2.2321380 -0.6916707 0.45971977
## 9 -0.3515156 -1.4322451 -0.03156438
## 10 1.6625430 1.8280966 0.85119311
## 11 -1.6401800 -1.2951128 0.49417845
## 12 1.6925941 -0.3922488 -0.01981007
## 13 1.7456787 -0.4375255 -0.27461543
Step6. 做主成分回归
ZZ1 <- data.frame(Z1) # 化为数据框
lm(ce.s[, 5] ~ z1 + z2 + z3, data = ZZ1) -> ce.fit1 # data 必须要是数据框
summary(ce.fit1)
##
## Call:
## lm(formula = ce.s[, 5] ~ z1 + z2 + z3, data = ZZ1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.20228 -0.12766 0.02411 0.07911 0.25621
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.831e-16 4.281e-02 0.000 1.0000
## z1 6.570e-01 2.980e-02 22.045 3.84e-09 ***
## z2 -8.309e-03 3.549e-02 -0.234 0.8202
## z3 3.028e-01 1.031e-01 2.935 0.0166 *
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1544 on 9 degrees of freedom
## Multiple R-squared: 0.9821, Adjusted R-squared: 0.9762
## F-statistic: 164.9 on 3 and 9 DF, p-value: 3.5e-08
常数项明显未通过，剔除
lm(ce.s[, 5] ~ 0 + z1 + z2 + z3, data = ZZ1) -> ce.fit2 #data 必须要是数据框
summary(ce.fit2)
##
## Call:
## lm(formula = ce.s[, 5] ~ 0 + z1 + z2 + z3, data = ZZ1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.20228 -0.12766 0.02411 0.07911 0.25621
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z1 0.656958 0.028271 23.238 4.93e-10 ***
## z2 -0.008309 0.033671 -0.247 0.8101
## z3 0.302770 0.097856 3.094 0.0114 *
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1464 on 10 degrees of freedom
## Multiple R-squared: 0.9821, Adjusted R-squared: 0.9768
## F-statistic: 183.2 on 3 and 10 DF, p-value: 4.895e-09
z2 检验未通过，但是主成分回归不太在乎显著性检验，因此可以剔除也可以不剔除，这里我们不剔除.
于是
\[y = 0.656958z_1 - 0.008309z_2 + 0.302770z_3 \]下面得到 \(y\) 关于 \(x\) 的回归方程
rep(ce.fit2$coefficients, each = 4)*ce.pca$loadings[, 1:3] -> tempp
apply(tempp, 1, sum) -> ce.fit3
ce.fit3
## X1 X2 X3 X4
## 0.51297502 0.27868114 -0.06078483 -0.42288461
于是
\[y = 0.51297502x_1 + 0.27868114x_2 − 0.06078483x_3 − 0.42288461x_4 \]标签：##,回归,ce,成分,X2,X3,X1 来源： https://www.cnblogs.com/hznudmh/p/16330423.html
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。