R语言与概率统计(三) 多元统计分析（中）
2024-04-18 16:51:20
模型修正
#但是，回归分析通常很难一步到位，需要不断修正模型
###############################6.9通过牙膏销量模型学习模型修正
toothpaste<-data.frame(X1=c(-0.05, 0.25,0.60,0, 0.25,0.20, 0.15,0.05,-0.15, 0.15,0.20, 0.10,0.40,0.45,0.35,0.30, 0.50,0.50, 0.40,-0.05,-0.05,-0.10,0.20,0.10,0.50,0.60,-0.05,0, 0.05, 0.55),X2=c( 5.50,6.75,7.25,5.50,7.00,6.50,6.75,5.25,5.25,6.00,6.50,6.25,7.00,6.90,6.80,6.80,7.10,7.00,6.80,6.50,6.25,6.00,6.50,7.00,6.80,6.80,6.50,5.75,5.80,6.80),Y =c( 7.38,8.51,9.52,7.50,9.33,8.28,8.75,7.87,7.10,8.00,7.89,8.15,9.10,8.86,8.90,8.87,9.26,9.00,8.75,7.95,7.65,7.27,8.00,8.50,8.75,9.21,8.27,7.67,7.93,9.26)
)plot(toothpaste)#先画出三点图观察变量之间的关系
cor(toothpaste)
> lm.sol<-lm(Y~X1+X2,data=toothpaste)
> summary(lm.sol)Call:
lm(formula = Y ~ X1 + X2, data = toothpaste)Residuals:Min 1Q Median 3Q Max
-0.49779 -0.12031 -0.00867 0.11084 0.58106 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.4075 0.7223 6.102 1.62e-06 ***
X1 1.5883 0.2994 5.304 1.35e-05 ***
X2 0.5635 0.1191 4.733 6.25e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.2383 on 27 degrees of freedom
Multiple R-squared: 0.886, Adjusted R-squared: 0.8776
F-statistic: 105 on 2 and 27 DF, p-value: 1.845e-13
> lm.new<-update(lm.sol, .~.+I(X2^2)) #.~.表示原来的Y~X1+X2，+I(X2^2)多加了一个平方项
> summary(lm.new)Call:
lm(formula = Y ~ X1 + X2 + I(X2^2), data = toothpaste)Residuals:Min 1Q Median 3Q Max
-0.40330 -0.14509 -0.03035 0.15488 0.46602 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.3244 5.6415 3.071 0.00495 **
X1 1.3070 0.3036 4.305 0.00021 ***
X2 -3.6956 1.8503 -1.997 0.05635 .
I(X2^2) 0.3486 0.1512 2.306 0.02934 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.2213 on 26 degrees of freedom
Multiple R-squared: 0.9054, Adjusted R-squared: 0.8945
F-statistic: 82.94 on 3 and 26 DF, p-value: 1.944e-13
> confint(lm.new)2.5 % 97.5 %
(Intercept) 5.72818421 28.9205529
X1 0.68290927 1.9310682
X2 -7.49886317 0.1076898
I(X2^2) 0.03786354 0.6593598
> #x2可能为0，考虑去掉
> lm2.new<-update(lm.new, .~.-X2)
> summary(lm2.new)Call:
lm(formula = Y ~ X1 + I(X2^2), data = toothpaste)Residuals:Min 1Q Median 3Q Max
-0.4859 -0.1141 -0.0046 0.1053 0.5592 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.07667 0.35531 17.102 5.17e-16 ***
X1 1.52498 0.29859 5.107 2.28e-05 ***
I(X2^2) 0.04720 0.00952 4.958 3.41e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.2332 on 27 degrees of freedom
Multiple R-squared: 0.8909, Adjusted R-squared: 0.8828
F-statistic: 110.2 on 2 and 27 DF, p-value: 1.028e-13
> #考虑两个变量之间不独立
> lm3.new<-update(lm2.new, .~.+X1*X2-X2)
> summary(lm3.new)Call:
lm(formula = Y ~ X1 + I(X2^2) + X1:X2, data = toothpaste)Residuals:Min 1Q Median 3Q Max
-0.48652 -0.11434 -0.00502 0.10524 0.55927 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.076101 0.364641 16.663 2.15e-15 ***
X1 1.573061 3.667273 0.429 0.671
I(X2^2) 0.047219 0.009786 4.825 5.33e-05 ***
X1:X2 -0.007064 0.536935 -0.013 0.990
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.2377 on 26 degrees of freedom
Multiple R-squared: 0.8909, Adjusted R-squared: 0.8783
F-statistic: 70.76 on 3 and 26 DF, p-value: 1.234e-12
逐步回归
#########################6.10逐步回归
cement<-data.frame(X1=c( 7, 1, 11, 11, 7, 11, 3, 1, 2, 21, 1, 11, 10),X2=c(26, 29, 56, 31, 52, 55, 71, 31, 54, 47, 40, 66, 68),X3=c( 6, 15, 8, 8, 6, 9, 17, 22, 18, 4, 23, 9, 8),X4=c(60, 52, 20, 47, 33, 22, 6, 44, 22, 26, 34, 12, 12),Y =c(78.5, 74.3, 104.3, 87.6, 95.9, 109.2, 102.7, 72.5, 93.1,115.9, 83.8, 113.3, 109.4)
)
plot(cement)
> lm.sol<-lm(Y ~ X1+X2+X3+X4, data=cement)
> summary(lm.sol)Call:
lm(formula = Y ~ X1 + X2 + X3 + X4, data = cement)Residuals:Min 1Q Median 3Q Max
-3.1750 -1.6709 0.2508 1.3783 3.9254 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) 62.4054 70.0710 0.891 0.3991
X1 1.5511 0.7448 2.083 0.0708 .
X2 0.5102 0.7238 0.705 0.5009
X3 0.1019 0.7547 0.135 0.8959
X4 -0.1441 0.7091 -0.203 0.8441
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 2.446 on 8 degrees of freedom
Multiple R-squared: 0.9824, Adjusted R-squared: 0.9736
F-statistic: 111.5 on 4 and 8 DF, p-value: 4.756e-07
> lm.step<-step(lm.sol)
Start: AIC=26.94
Y ~ X1 + X2 + X3 + X4Df Sum of Sq RSS AIC
- X3 1 0.1091 47.973 24.974
- X4 1 0.2470 48.111 25.011
- X2 1 2.9725 50.836 25.728
<none> 47.864 26.944
- X1 1 25.9509 73.815 30.576Step: AIC=24.97
Y ~ X1 + X2 + X4Df Sum of Sq RSS AIC
<none> 47.97 24.974
- X4 1 9.93 57.90 25.420
- X2 1 26.79 74.76 28.742
- X1 1 820.91 868.88 60.629
>
> summary(lm.step)Call:
lm(formula = Y ~ X1 + X2 + X4, data = cement)Residuals:Min 1Q Median 3Q Max
-3.0919 -1.8016 0.2562 1.2818 3.8982 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) 71.6483 14.1424 5.066 0.000675 ***
X1 1.4519 0.1170 12.410 5.78e-07 ***
X2 0.4161 0.1856 2.242 0.051687 .
X4 -0.2365 0.1733 -1.365 0.205395
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 2.309 on 9 degrees of freedom
Multiple R-squared: 0.9823, Adjusted R-squared: 0.9764
F-statistic: 166.8 on 3 and 9 DF, p-value: 3.323e-08
> drop1(lm.step)
Single term deletionsModel:
Y ~ X1 + X2 + X4Df Sum of Sq RSS AIC
<none> 47.97 24.974
X1 1 820.91 868.88 60.629
X2 1 26.79 74.76 28.742
X4 1 9.93 57.90 25.420
> lm.opt<-lm(Y ~ X1+X2, data=cement); summary(lm.opt)Call:
lm(formula = Y ~ X1 + X2, data = cement)Residuals:Min 1Q Median 3Q Max
-2.893 -1.574 -1.302 1.363 4.048 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) 52.57735 2.28617 23.00 5.46e-10 ***
X1 1.46831 0.12130 12.11 2.69e-07 ***
X2 0.66225 0.04585 14.44 5.03e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 2.406 on 10 degrees of freedom
Multiple R-squared: 0.9787, Adjusted R-squared: 0.9744
F-statistic: 229.5 on 2 and 10 DF, p-value: 4.407e-09
> library(car)#多重共线性，vif方差膨胀因子
> vif(lm.opt)X1 X2
1.055129 1.055129
> sqrt(vif(lm.opt))>2#problem?若存在共线性怎么办？X1 X2
FALSE FALSE
> AIC(lm.sol,lm.opt)df AIC
lm.sol 6 65.83669
lm.opt 4 64.31239
> BIC(lm.sol,lm.opt)df BIC
lm.sol 6 69.22639
lm.opt 4 66.57219做线性回归要做：t检验，F检验，残差分析，模型解释（灵敏度分析）
残差分析（回归诊断）
> ######6.5forbes数据分析的残差分析
> X<-matrix(c(
+ 194.5, 20.79, 1.3179, 131.79,
+ 194.3, 20.79, 1.3179, 131.79,
+ 197.9, 22.40, 1.3502, 135.02,
+ 198.4, 22.67, 1.3555, 135.55,
+ 199.4, 23.15, 1.3646, 136.46,
+ 199.9, 23.35, 1.3683, 136.83,
+ 200.9, 23.89, 1.3782, 137.82,
+ 201.1, 23.99, 1.3800, 138.00,
+ 201.4, 24.02, 1.3806, 138.06,
+ 201.3, 24.01, 1.3805, 138.05,
+ 203.6, 25.14, 1.4004, 140.04,
+ 204.6, 26.57, 1.4244, 142.44,
+ 209.5, 28.49, 1.4547, 145.47,
+ 208.6, 27.76, 1.4434, 144.34,
+ 210.7, 29.04, 1.4630, 146.30,
+ 211.9, 29.88, 1.4754, 147.54,
+ 212.2, 30.06, 1.4780, 147.80),
+ ncol=4, byrow=T,
+ dimnames = list(1:17, c("F", "h", "log", "log100")))
> forbes<-data.frame(X)
> plot(forbes$F, forbes$log100)#画出两个变量之间的散点图，观察是否存在线性趋势
> lm.sol<-lm(log100~F, data=forbes)
> summary(lm.sol)Call:
lm(formula = log100 ~ F, data = forbes)Residuals:Min 1Q Median 3Q Max
-0.32261 -0.14530 -0.06750 0.02111 1.35924 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) -42.13087 3.33895 -12.62 2.17e-09 ***
F 0.89546 0.01645 54.45 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.3789 on 15 degrees of freedom
Multiple R-squared: 0.995, Adjusted R-squared: 0.9946
F-statistic: 2965 on 1 and 15 DF, p-value: < 2.2e-16> abline(lm.sol)#在散点图上添加直线
> #残差正态性检验
> y.res=residuals(lm.sol)
> shapiro.test(y.res)Shapiro-Wilk normality testdata: y.res
W = 0.54654, p-value = 3.302e-06
所有代码：
#但是，回归分析通常很难一步到位，需要不断修正模型
###############################6.9通过牙膏销量模型学习模型修正
toothpaste<-data.frame(X1=c(-0.05, 0.25,0.60,0, 0.25,0.20, 0.15,0.05,-0.15, 0.15,0.20, 0.10,0.40,0.45,0.35,0.30, 0.50,0.50, 0.40,-0.05,-0.05,-0.10,0.20,0.10,0.50,0.60,-0.05,0, 0.05, 0.55),X2=c( 5.50,6.75,7.25,5.50,7.00,6.50,6.75,5.25,5.25,6.00,6.50,6.25,7.00,6.90,6.80,6.80,7.10,7.00,6.80,6.50,6.25,6.00,6.50,7.00,6.80,6.80,6.50,5.75,5.80,6.80),Y =c( 7.38,8.51,9.52,7.50,9.33,8.28,8.75,7.87,7.10,8.00,7.89,8.15,9.10,8.86,8.90,8.87,9.26,9.00,8.75,7.95,7.65,7.27,8.00,8.50,8.75,9.21,8.27,7.67,7.93,9.26)
)plot(toothpaste)#先画出三点图观察变量之间的关系
cor(toothpaste)lm.sol<-lm(Y~X1+X2,data=toothpaste)
summary(lm.sol)
#根据散点图做如下修正
lm.new<-update(lm.sol, .~.+I(X2^2)) #.~.表示原来的Y~X1+X2，+I(X2^2)多加了一个平方项
summary(lm.new)confint(lm.new)
#x2可能为0，考虑去掉
lm2.new<-update(lm.new, .~.-X2)
summary(lm2.new)
#考虑两个变量之间不独立
lm3.new<-update(lm2.new, .~.+X1*X2-X2)
summary(lm3.new)#########################6.10逐步回归
cement<-data.frame(X1=c( 7, 1, 11, 11, 7, 11, 3, 1, 2, 21, 1, 11, 10),X2=c(26, 29, 56, 31, 52, 55, 71, 31, 54, 47, 40, 66, 68),X3=c( 6, 15, 8, 8, 6, 9, 17, 22, 18, 4, 23, 9, 8),X4=c(60, 52, 20, 47, 33, 22, 6, 44, 22, 26, 34, 12, 12),Y =c(78.5, 74.3, 104.3, 87.6, 95.9, 109.2, 102.7, 72.5, 93.1,115.9, 83.8, 113.3, 109.4)
)
plot(cement)
lm.sol<-lm(Y ~ X1+X2+X3+X4, data=cement)
summary(lm.sol)lm.step<-step(lm.sol)summary(lm.step)drop1(lm.step)lm.opt<-lm(Y ~ X1+X2, data=cement); summary(lm.opt)library(car)#多重共线性，vif方差膨胀因子
vif(lm.opt)
sqrt(vif(lm.opt))>2#problem?若存在共线性怎么办？AIC(lm.sol,lm.opt)
BIC(lm.sol,lm.opt)######6.5forbes数据分析的残差分析
X<-matrix(c(194.5, 20.79, 1.3179, 131.79,194.3, 20.79, 1.3179, 131.79,197.9, 22.40, 1.3502, 135.02,198.4, 22.67, 1.3555, 135.55,199.4, 23.15, 1.3646, 136.46,199.9, 23.35, 1.3683, 136.83,200.9, 23.89, 1.3782, 137.82,201.1, 23.99, 1.3800, 138.00,201.4, 24.02, 1.3806, 138.06,201.3, 24.01, 1.3805, 138.05,203.6, 25.14, 1.4004, 140.04,204.6, 26.57, 1.4244, 142.44,209.5, 28.49, 1.4547, 145.47,208.6, 27.76, 1.4434, 144.34,210.7, 29.04, 1.4630, 146.30,211.9, 29.88, 1.4754, 147.54,212.2, 30.06, 1.4780, 147.80),ncol=4, byrow=T,dimnames = list(1:17, c("F", "h", "log", "log100")))forbes<-data.frame(X)
plot(forbes$F, forbes$log100)#画出两个变量之间的散点图，观察是否存在线性趋势
lm.sol<-lm(log100~F, data=forbes)
summary(lm.sol)
abline(lm.sol)#在散点图上添加直线#残差正态性检验
y.res=residuals(lm.sol)
shapiro.test(y.res)y.res<-residuals(lm.sol);plot(y.res)#画出残差图
text(12,y.res[12], labels=12,adj=1.2)#或者
plot(lm.sol)plot(lm.sol,2)#残差的qq图#异常值的判断
library(car)
outlierTest(lm.sol)#去除异常值
i<-1:17; forbes12<-data.frame(X[i!=12, ])
lm12<-lm(log100~F, data=forbes12)
summary(lm12)#残差正态性检验
y.res=residuals(lm12)
shapiro.test(y.res)View Code
R语言与概率统计(三) 多元统计分析（中）相关推荐
- R语言与概率统计(六) 主成分分析 因子分析
超高维度分析,N*P的矩阵,N为样本个数,P为指标,N<<P PCA:抓住对y对重要的影响因素 主要有三种:PCA,因子分析,回归方程+惩罚函数(如LASSO) 为了降维,用更少的变量解决 ...
- R语言与概率统计(四) 判别分析(分类)
Fisher就是找一个线L使得组内方差小,组间距离大.即找一个直线使得d最大. ####################################1.判别分析,线性判别:2.分层抽样#insta ...
- R语言：优雅、卓越的统计分析及绘图环境
文 / 刘思喆 历史 R语言由新西兰奥克兰大学的Ross Ihaka和Robert Gentleman两人共同发明,其词法和语法分别源自Scheme和S语言,一般认为R语言是S语言[注:John Ch ...
- R语言的各种统计分布函数
转载自品略图书馆 http://www.pinlue.com/article/2018/09/1613/487222559948.html R语言的各种统计分布函数 1.二项分布Binomial di ...
- R语言循环函数编写三境界
1. 三境界 R语言写循环有三境界: 手动for循环 apply系列 purrr泛函式编程 其中,手动for循环我最常用,apply系列半吊子,purrr函数越用越趁手,这里介绍一下purrr编写循环 ...
- R语言使用lm函数拟合多元线性回归模型、假定预测变量没有交互作用（Multiple linear regression）
R语言使用lm函数拟合多元线性回归模型.假定预测变量没有交互作用(Multiple linear regression) 目录
- R语言数据描述性统计（Descriptive statistics）实战：数据全局描述信息、数值数据的描述性统计（Numerical data）、离散型数据的描述性统计（Categorical）
R语言数据描述性统计(Descriptive statistics)实战:数据全局描述信息.数值数据的描述性统计(Numerical data).离散型数据的描述性统计(Categorical) 目录
- R语言Wilcoxon Signed-rank统计分布函数（dsignrank, psignrank, qsignrank rsignrank ）实战
R语言Wilcoxon Signed-rank统计分布函数(dsignrank, psignrank, qsignrank & rsignrank )实战 目录 R语言Wilcoxon Sig ...
- R语言nchar函数统计字符串中字符个数实战
R语言nchar函数统计字符串中字符个数实战 目录 R语言nchar函数统计字符串中字符个数实战 #基础语法
最新文章
- HDU 2034 人见人爱A-B【STL/set】
- linux修改进程优先级
- python爬虫反爬-python爬虫--爬虫与反爬
- dispay的flex属性
- 【Java】关于Java的一些基础知识点
- Delphi的MessageBox对话框使用
- 异常处理第一讲(SEH),筛选器异常,以及__asm的扩展,寄存器注入简介
- golang中apend_golang的append()为什么不会影响slice的地址？
- docker没有下载完全_一个时代的结束：苹果正式关闭iTunes，但歌曲下载并没有完全消失...
- kettle使用数据库来生成序列_kettle 生成 ktr
- Rust或C#,Python 等如何封装C++的接口 (比如CTP)？
- 安装svn + vs code配置svn
- qqkey获取原理_QQKEY获取多功能软件+【附源码】
- ionic-vue 开发app移动端
- 模拟高清和数字高清摄像机的区别，全局快门CMOS 图像传感器，Interline Transfer CCD图像传感器
- mysql dump hbase_mysqldump 导出部分数据库
- 增长量计算n+1原则_资料分析几大常用公式，增速、A/B型公式.....
- 有才的人全败给“傲”，平庸的人皆输在“懒”！
- 免费SSL证书申请和部署
- GUI使用2——总结NGUI、tookit2D、GUI比较