• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    公众号
登陆 注册

R语言与概率统计(三) 多元统计分析（中）

原作者: [db:作者] 来自: [db:来源] 收藏 邀请

模型修正

#但是，回归分析通常很难一步到位，需要不断修正模型
###############################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 ‘ ’ 1

Residual 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 ‘ ’ 1

Residual 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 ‘ ’ 1

Residual 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 ‘ ’ 1

Residual 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 ‘ ’ 1

Residual 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 + X4

       Df 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.576

Step:  AIC=24.97
Y ~ X1 + X2 + X4

       Df 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 ‘ ’ 1

Residual 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 deletions

Model:
Y ~ X1 + X2 + X4
       Df 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 ‘ ’ 1

Residual 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 ‘ ’ 1

Residual 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 test

data:  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

 

鲜花
握手
雷人
路过
鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
More+
上一篇：
R语言与概率统计(三) 多元统计分析（下）广义线性回归发布时间：2022-07-18
下一篇：
R语言多元分析系列之五：聚类分析（完） | R blogs / lang发布时间：2022-07-18
热门推荐
More+
热门话题
More+
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服（服务时间 9:00～18:00）

在线QQ客服
地址：深圳市南山区西丽大学城创智工业园
电邮：jeky_zhao#qq.com
移动电话：139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap