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

来源：互联网发布：java获取本地磁盘文件编辑：程序博客网时间：2024/04/27 06:32

实例一、一元线性回归

测试沸点与气压的关系: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:17forbes <- 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:17forbes <- 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/b1glm.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$nplot(norell$x, norell$y)lines(d, p)

#其中，d是给出曲线横坐标的点，pre是计算预测值，p是相应的预测概率。用plot函数和lines给出散点图和对应的预测曲线。

0 0