1.3线性回归之线性回归实例

2021年6月26日 5点热度 0条评论 来源: qq_16365849

线性回归实例

junjun

2016年2月7日

实例一、一元线性回归

测试沸点与气压的关系:Forbes数据共有四列,分别表示 沸点(F)、气压(英寸汞柱)、log(气压)、100*log(气压)

#1、读取数据
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")))

#2、数据预处理
#1)将数据转化为数据框格式
forbes <- as.data.frame(X)

#2)画图分析因变量与自变量之间的关系
plot(forbes$F, forbes$log100)

#3、建模:从散点图上可以看出,这些点基本上落在一条一直上,所以做回归分析
model_lm <- lm(log100~F, data=forbes)

#4、模型评估
#1)分析summary函数结果
summary(model_lm)
## 
## 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
#从上面可以看出:对应于两个系数的P值均<2.17*10-9,是非常显著的;关于方程的检验,残差的标准差为0.3789,相关系数(R2)为0.995,关于F分布的P值<2.2*10-16也是分成显著的。所以,该模型能通过t检验和F检验。因此,回归方程为 Y-42.13087+0.89546X 

#2)分析残差:函数residuals函数计算回归方程的残差。计算残差并画出关于残差的散点图
y.residuals <- residuals(model_lm)
plot(y.residuals)

#从图中可以看出第12个样本可能会有问题,比其他的样本点的残差大得多。因此,这个点可能不正确,或者模型的差假设不正确,或者是方差不是常数等。需要对这个问题进行分析。

#5、模型优化
#1)把每个点都标记出来 
a <- as.data.frame(y.residuals)
a$id <- 1:nrow(a)
names(a) <- list("value", "id")
b <- a[, c("id", "value")]
plot(value~id, b, xlab=b$id)
for(i in 1:nrow(b))
  text(b[i, 1], b[i, 2], labels = i, adj=1.2)

#2)去掉第12号样本点
i <- 1:17
forbes <- as.data.frame(X[i!=12, ])
model_lm12 <- lm(log100~F, data=forbes)
summary(model_lm12)
## 
## Call:
## lm(formula = log100 ~ F, data = forbes)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.21175 -0.06194  0.01590  0.09077  0.13042 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -41.30180    1.00038  -41.29 5.01e-16 ***
## F             0.89096    0.00493  180.73  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1133 on 14 degrees of freedom
## Multiple R-squared:  0.9996, Adjusted R-squared:  0.9995 
## F-statistic: 3.266e+04 on 1 and 14 DF,  p-value: < 2.2e-16
#去掉第12号样本后,回归方程的系数没有太大的变化,但系数的标准差和残差的标准差有很大的变化,小了3倍左右,相关系数R2也有提高。

#6、预测
x <- data.frame(F=195)
pred <- predict(model_lm12, x)
pred
##        1 
## 132.4347

实例二、多元线性回归

试根据提供的数据建立一个数学模型,分析牙膏销售量与其他因素的关系,为定制加个人策略提供数量依据。 牙膏销量为Y,价格差为X1,公司的广告费为X2,假设基本模型为线性模型 Y=b0+b1X1+b2X2+e 输入数据,调用R软件中的lm()函数求解,并用summary()显示计算结果

#1、加载数据
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)
)

#2、建立模型
model_lm <- lm(Y~X1+X2, data=toothpaste)

#3、模型评估
#1)查看summary函数的结果
summary(model_lm)
## 
## 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
#由此得到销售量与价格差与广告费之间的关系为 Y=4.4075+1.5883X1+0.5635X2

#2)分别做y与x1和y与x2的散点图,可以看出y与x1的散点图拟合为一条直线,y与x2的散点图拟合为一条曲线。所以,修改模型为:
attach(toothpaste)
plot(Y~X1)

plot(Y~X2)

model_lm2 <- update(model_lm, .~.+I(X2^2))
summary(model_lm2)
## 
## 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
#此时,模型残差的标准差有所下降,相关系数的平方有所上升,这说明修正是合理的。但是,对应于b2的P值>0.05。进一步分析,作b的区间估计

#3)编写函数:求拟合函数fm的区间估计(a=0.05)
beta.int <- function(fm, alpha=0.05) { 
    A <- summary(fm)$coefficients      
    df <- fm$df.residual
    left <- A[,1]-A[,2]*qt(1-alpha/2, df)
    right <- A[,1]+A[,2]*qt(1-alpha/2, df)
    rowname <- dimnames(A)[[1]]
    colname <- c("Estimate", "Left", "Right")
    matrix(c(A[,1], left, right), ncol=3, dimnames=list(rowname, colname))
}

beta.int(model_lm2)
##               Estimate        Left      Right
## (Intercept) 17.3243685  5.72818421 28.9205529
## X1           1.3069887  0.68290927  1.9310682
## X2          -3.6955867 -7.49886317  0.1076898
## I(X2^2)      0.3486117  0.03786354  0.6593598
#由此可知,b2的区间估计是[-7.49886317, 0.1076898],包含了0,也就是说b2的值可能会使0。所以,去掉X2的一次项,再进行分析
model_lm3 <- update(model_lm2, .~.-X2)
summary(model_lm3)
## 
## 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
#由此可知,此模型虽然通过了F检验和T检验,但是与上一模型对比来看,残差的标准差上升了,相关系数下降了,所以模型有不足之处,需要进一步修正,考虑x1与x2的交互作用

#4)添加X1与X2的相互作用
model_lm4 <- update(model_lm3, .~.+X1*X2)
summary(model_lm4)
## 
## Call:
## lm(formula = Y ~ X1 + I(X2^2) + X2 + X1:X2, data = toothpaste)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.43725 -0.11754  0.00489  0.12263  0.38410 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  29.1133     7.4832   3.890 0.000656 ***
## X1           11.1342     4.4459   2.504 0.019153 *  
## I(X2^2)       0.6712     0.2027   3.312 0.002824 ** 
## X2           -7.6080     2.4691  -3.081 0.004963 ** 
## X1:X2        -1.4777     0.6672  -2.215 0.036105 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2063 on 25 degrees of freedom
## Multiple R-squared:  0.9209, Adjusted R-squared:  0.9083 
## F-statistic: 72.78 on 4 and 25 DF,  p-value: 2.107e-13
#由上可知:模型通过t检验和F检验,并且残差的标准差减少了,相关系数增加了。因此,模型最终为Y=29.1133+11.1342X2-7.6080X2+0.6712X2^2-1.4777X1X2+e

实例三、逐步回归

某种水泥在凝固时放出的热量Y(卡/克)与水泥中四种化学成分X1、X2、X3、X4有关,现测得13组数据,选出主要的变量,建立Y关于他们的线性回归方程

#1、加载数据
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)
)

#2、建模
model_lm <- lm(Y~., data=cement)

#3、模型评估
summary(model_lm)
## 
## Call:
## lm(formula = Y ~ ., 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
#从上可以看到,如果选择全部变量作回归方程,效果是不好的。因为,回归方程的系数全部都没有通过检验。

#4、用函数step()作逐步回归
model_step <- step(model_lm)
## 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
##这里去掉x3后,再去掉任意变量,都会使AIC值增大。

#分析计算结果
model_update <- update(model_lm, .~.-X3)
summary(model_update)
## 
## 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
#由上可知:回归系数检验的显著性水平有很大提高,但变量X2、X4系数检验的显著水平仍不理想。在R软件中,还有add1()和drop1()方法,尝试去掉一个变量。

#5、下面,用drop1()计算:
drop1(model_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
#去掉变量x4残差平方和增大9.93,AIC增长也是最小的,所以去掉x4

#6、建立最终的模型
model <- lm(Y~X1+X2, data=cement)
summary(model)
## 
## 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
#从上可知:所有的检验均是显著的,最后得到“最优”的回归方程为 Y=52.57735+1.46831X1+0.66225X2

实例四、回归诊断

在建立模型前后,需要做回归诊断:异常样品的存在往往会给回归模型带来不稳定,所以提出了回归诊断的问题,主要内容有

(1)关于误差项是否满足:①独立性;②等方差性;③正态性 (2)选择线性模型是否合适? (3)是否存在异常样本? (4)回归分析的结果是否对某些样本的依赖过重?即回归模型是否具备稳定性? (5)自变量之间是否存在高度相关?即是否有多重共线性问题存在?

#1、残差:判断模型的参数是否符合正态性分布
#residuals()或resid()函数提供了模型残差的计算,得到残差后,可以对残差进行检验,如正态性检验等。shapiro.test()函数为正态性检验函数。
#实例:对Forbes数据得到回归模型,得到的残差做W正态性检验
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 <- as.data.frame(X)
model_lm <- lm(log100~F, data=forbes)
y.residuals <- residuals(model_lm)
shapiro.test(y.residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  y.residuals
## W = 0.54654, p-value = 3.302e-06
#此时,P值小于0.05拒绝原假设,即该数据不符合正态分布。因此,在这个例子中我们也可发现较大的W统计量的情况下,我们依然可能拒绝其符合正态分布。

i <- 1:17
forbes <- as.data.frame(X[i!=12, ])
model_lm2 <- lm(log100~F, data=forbes)
y.residuals <- residuals(model_lm2)
shapiro.test(y.residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  y.residuals
## W = 0.92215, p-value = 0.1827
#此时,其W统计量接近1,p值显著大于0.05,所以我们没办法拒绝其符合正态分布.所以,通过正态性检验,去掉 第12号样本点是合理的。

#2、残差图:rstandard()计算回归模型的标准化残差
#利用上面的血压与身高的数据:
blood<-data.frame(
X1=c(76.0, 91.5, 85.5, 82.5, 79.0, 80.5, 74.5, 79.0, 85.0, 76.5, 82.0, 95.0, 92.5),
X2=c(50, 20, 20, 30, 30, 50, 60, 50, 40, 55, 40, 40, 20),
Y= c(120, 141, 124, 126, 117, 125, 123, 125, 132, 123, 132, 155, 147)
 )

#画残差图:
model_lm <- lm(Y~X1+X2, data=blood)
y.residuals <- resid(model_lm)
y.fit <- predict(model_lm)
plot(y.residuals~y.fit)

#画标准化残差图:
dev.new()
y.rst <- rstandard(model_lm)
plot(y.rst~y.fit)
#从上可知:当残差服从正态分布的假设成立时,标准化残差应近似的服从标准正态分布。对于标准化残差,应该有95%的样本点落在区间[-2, 2]中,通过标准化残差图,更容易诊断出回归模型是否出现问题。

#3、残差的Q-Q图:正态分布的检验方法——QQ图,可以检验残差的正态性。
a <- c(1:100)
qqnorm(a)
qqline(a)
#qqplot()函数提供了更为精确的正态假设检验方法,画出了n-p-1个自由度的t分布下的学生化残差图形,n为样本大小,p是回归参数的数目(包括截距项) 用法:qqplot(lm(),data=))

#注意:car包中的qqPlot(model, labels, id.method, simulate)函数:分位比较图。labels为文本字符串向量;id.method为点识别方法;simulate为T,计算置信区间。用于检测正态性。

#4、多重共线性:一个变量可以由其他变量求出,例如,学生的总成绩可以由各科成绩求出。度量多重共线性严重程度的一个重要指标是矩阵的条件数,可以由函数kappa()求出。在R中,函数kappa()计算矩阵的条件数,其使用方法为:kappa(z, exact=F, …) 其中,z是矩阵(相关系数的矩阵,可由cor函数求得),exact是逻辑变量,当exact=T时,精确计算条件数,否则近似计算条件数。
#注意:一般条件数K<100,则认为多重共线性的程度很小;若100<=K<=1000则认为存在中等程度或较强的多重共线性;若K>1000则认为存在严重的多重共线性。
collinear<-data.frame(
Y=c(10.006, 9.737, 15.087, 8.422, 8.625, 16.289, 5.958, 9.313, 12.960, 5.541, 8.756, 10.937),
X1=rep(c(8, 0, 2, 0), c(3, 3, 3, 3)),
X2=rep(c(1, 0, 7, 0), c(3, 3, 3, 3)),
X3=rep(c(1, 9, 0), c(3, 3, 6)),
X4=rep(c(1, 0, 1, 10), c(1, 2, 6, 3)),
X5=c(0.541, 0.130, 2.116, -2.397, -0.046, 0.365, 1.996, 0.228, 1.38, -0.798, 0.257, 0.440),
X6=c(-0.099, 0.070, 0.115, 0.252, 0.017, 1.504, -0.865, -0.055, 0.502, -0.399, 0.101, 0.432) )

XX <- cor(collinear[2:7])
kappa(XX, exact = T)
## [1] 2195.908
#由上可知:得到条件数 K=2195.908>1000,认为有严重的多重共线性。

#找出哪些变量是多重共线性的,计算矩阵的特征值和相应的特征向量:
eigen(XX)
## $values
## [1] 2.428787365 1.546152096 0.922077664 0.793984690 0.307892134 0.001106051
## 
## $vectors
##            [,1]        [,2]        [,3]        [,4]       [,5]
## [1,] -0.3907189  0.33968212  0.67980398 -0.07990398  0.2510370
## [2,] -0.4556030  0.05392140 -0.70012501 -0.05768633  0.3444655
## [3,]  0.4826405  0.45332584 -0.16077736 -0.19102517 -0.4536372
## [4,]  0.1876590 -0.73546592  0.13587323  0.27645223 -0.0152087
## [5,] -0.4977330  0.09713874 -0.03185053  0.56356440 -0.6512834
## [6,]  0.3519499  0.35476494 -0.04864335  0.74817535  0.4337463
##              [,6]
## [1,] -0.447679719
## [2,] -0.421140280
## [3,] -0.541689124
## [4,] -0.573371872
## [5,] -0.006052127
## [6,] -0.002166594

实例五、广义线性回归模型

Norell实验,高压电线对牲畜的影响

#1、加载数据
norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )

#2、数据预处理
norell$Ymat<-cbind(norell$success, norell$n-norell$success)

#3、建模
glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)

#4、模型评估
summary(glm.sol)
## 
## Call:
## glm(formula = Ymat ~ x, family = binomial, data = norell)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6  
## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
## x             1.2459     0.1119   11.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 250.4866  on 5  degrees of freedom
## Residual deviance:   9.3526  on 4  degrees of freedom
## AIC: 34.093
## 
## Number of Fisher Scoring iterations: 4
#5、预测:与线性回归模型相同,在得到回归模型后,可以作预测:电流强度为3.5毫安时,有响应的牛的概率
pre <- predict(glm.sol, data.frame(x=3.5))
p <- exp(pre)/(1+exp(pre))

#6、求有50%的牛响应时的电流强度:当P=0.5时,ln(P/(1-P))=0,所以X=-b0/b1
glm.sol$coefficients
## (Intercept)           x 
##   -3.301035    1.245937
X <- -glm.sol$coefficients[[1]]/glm.sol$coefficients[[2]]

#7、画出响应比例与logistic回归曲线
d <- seq(0, 5, length=100)
pre <- predict(glm.sol, data.frame(x=d))
p <- exp(pre)/(1+exp(pre))
norell$y <- norell$success/norell$n
plot(norell$x, norell$y)
lines(d, p)

#其中,d是给出曲线横坐标的点,pre是计算预测值,p是相应的预测概率。用plot函数和lines给出散点图和对应的预测曲线。
    原文作者:qq_16365849
    原文地址: https://blog.csdn.net/qq_16365849/article/details/50642970
    本文转自网络文章,转载此文章仅为分享知识,如有侵权,请联系管理员进行删除。