摘要: 本文由digging4發表於：http://www.cnblogs.com/digging4/p/5079829.html
統計建模與R軟件-第六章 回歸分析
6.1為了估計山上積雪融化后對下游灌溉的影響，在山上建立一個觀測站，測量最大積雪深度X與當年灌溉面積Y，測得連續10年的數據如表6.17所示
表6.17 10年中最大積雪深度與當年灌溉面積的數據
序號 X(米) Y(公頃)
1 5.1 1907
2 3.5 1287
3 7.1 2700
4 6.2 2373
5 8.8 3260
6 7.8 3000
7 4.5 1947
8 5.6 2273
9 8.0 3113
10 6.4 2493
（1）試畫出相應的散點圖，判斷Y與X是否有線性關系；
（2）求出Y關於X的一元線性回歸方程；
（3）對方程作顯著性檢驗；
（4）現測得今年的數據是X=7米，給出今年灌溉面積的預測值和相應的區間估計（\(\alpha=0.05\)）
x <- c(5.1, 3.5, 7.1, 6.2, 8.8, 7.8, 4.5, 5.6, 8, 6.4)
y <- c(1907, 1287, 2700, 2373, 3260, 3000, 1947, 2273, 3113, 2493)
# 從散點圖可以看出，這些點基本上在一條直線上，x和y有線性關系
plot(x, y)
# y=141.0 + 364.2*x ，但是截距項不顯著
lm.sol <- lm(y ~ x)
summary(lm.sol)
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -128.59 -70.98 -3.73 49.26 167.23
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 141.0 125.1 1.13 0.29
## x 364.2 19.3 18.91 6.3e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 96.4 on 8 degrees of freedom
## Multiple R-squared: 0.978, Adjusted R-squared: 0.975
## F-statistic: 358 on 1 and 8 DF, p-value: 6.33e-08
# 采用無截距項回歸，回歸方程通過了顯著性檢驗，且R平方值提高到0.999
# 回歸方程 y=385.23*x
lm.sol <- lm(y ~ x - 1)
summary(lm.sol)
##
## Call:
## lm(formula = y ~ x - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -130.0 -52.0 -10.1 30.3 213.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## x 385.23 4.76 80.9 3.4e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 97.8 on 9 degrees of freedom
## Multiple R-squared: 0.999, Adjusted R-squared: 0.998
## F-statistic: 6.54e+03 on 1 and 9 DF, p-value: 3.42e-14
abline(lm.sol)
#
# 預測，就算只有一個值，也要采用數據框的形式，並且數據框的變量名和回歸方程的自變量名要一致
new <- data.frame(x = 7)
predict(lm.sol, new, interval = "prediction", level = 0.95)
## fit lwr upr
## 1 2697 2463 2930
# 得到預測值為2697，其95%的可信區間為 [2463, 2930]
6.2研究同一地區土壤所含可給態磷的情況，得到18組數據如表6.18所示，表中\(X_1\)為土壤內所含無機磷濃度，\(X_2\)為土壤內溶於\(K_2CO_3\)溶液並受溴化物水解的有機磷，\(X_3\)為土壤內溶於\(K_2CO_3\)溶液但不溶於溴化物水解的有機磷。
（1）求出Y關於X的多元線性回歸方程
（2）對方程做顯著性檢驗
（3）對變量做逐步回歸分析
序號 X1 X2 X3 Y
1 0.4 52 158 64
2 0.4 23 163 60
3 3.1 19 37 71
4 0.6 34 157 61
5 4.7 24 59 54
6 1.7 65 123 77
7 9.4 44 46 81
8 10.1 31 117 93
9 11.6 29 173 93
10 12.6 58 112 51
11 10.9 37 111 76
12 23.1 46 114 96
13 23.1 50 134 77
14 21.6 44 73 93
15 23.1 56 168 95
16 1.9 36 143 54
17 26.8 58 202 168
18 29.9 51 124 99
soil <- data.frame(X1 = c(0.4, 0.4, 3.1, 0.6, 4.7, 1.7, 9.4, 10.1, 11.6, 12.6,
10.9, 23.1, 23.1, 21.6, 23.1, 1.9, 26.8, 29.9), X2 = c(52, 23, 19, 34, 24,
65, 44, 31, 29, 58, 37, 46, 50, 44, 56, 36, 58, 51), X3 = c(158, 163, 37,
157, 59, 123, 46, 117, 173, 112, 111, 114, 134, 73, 168, 143, 202, 124),
Y = c(64, 60, 71, 61, 54, 77, 81, 93, 93, 51, 76, 96, 77, 93, 95, 54, 168,
99))
lm.sol <- lm(Y ~ X1 + X2 + X3, data = soil)
summary(lm.sol)
##
## Call:
## lm(formula = Y ~ X1 + X2 + X3, data = soil)
##
## Residuals:
## Min 1Q Median 3Q Max
## -28.35 -11.38 -2.66 12.10 48.81
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 43.6501 18.0544 2.42 0.0298 *
## X1 1.7853 0.5398 3.31 0.0052 **
## X2 -0.0833 0.4204 -0.20 0.8458
## X3 0.1610 0.1116 1.44 0.1710
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20 on 14 degrees of freedom
## Multiple R-squared: 0.549, Adjusted R-squared: 0.453
## F-statistic: 5.69 on 3 and 14 DF, p-value: 0.00923
# 可以看出X2和X3的系數不顯著，可以考慮從畫出X2～Y，X3·Y的圖進行分析
plot(residuals(lm.sol))
# 發現第17個點的殘差比其他點都大，去掉17點后重新計算
lm.sol2 <- lm(Y ~ X1 + X2 + X3, data = soil[-17, ])
summary(lm.sol2)
##
## Call:
## lm(formula = Y ~ X1 + X2 + X3, data = soil[-17, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.79 -5.68 2.27 4.89 16.78
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 64.4929 12.2929 5.25 0.00016 ***
## X1 1.3034 0.3576 3.64 0.00297 **
## X2 -0.1322 0.2669 -0.50 0.62863
## X3 0.0227 0.0767 0.30 0.77177
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.7 on 13 degrees of freedom
## Multiple R-squared: 0.529, Adjusted R-squared: 0.42
## F-statistic: 4.86 on 3 and 13 DF, p-value: 0.0176
# 和之前一樣，X2和X3的系數都沒有通過顯著性檢驗。
# 使用step函數，進行逐步回歸
lm.step <- step(lm.sol2)
## Start: AIC=89.77
## Y ~ X1 + X2 + X3
##
## Df Sum of Sq RSS AIC
## - X3 1 14 2101 87.9
## - X2 1 39 2126 88.1
## <none> 2087 89.8
## - X1 1 2132 4218 99.7
##
## Step: AIC=87.89
## Y ~ X1 + X2
##
## Df Sum of Sq RSS AIC
## - X2 1 31 2131 86.1
## <none> 2101 87.9
## - X1 1 2119 4220 97.7
##
## Step: AIC=86.13
## Y ~ X1
##
## Df Sum of Sq RSS AIC
## <none> 2131 86.1
## - X1 1 2295 4426 96.6
# 可以看到，去掉自變量X2和X3以后，AIC=86.13 達到最小
lm.fin <- lm(Y ~ X1, data = soil[-17, ])
summary(lm.step)
##
## Call:
## lm(formula = Y ~ X1, data = soil[-17, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.056 -3.061 0.939 5.038 18.016
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 62.569 4.452 14.05 4.8e-10 ***
## X1 1.229 0.306 4.02 0.0011 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.9 on 15 degrees of freedom
## Multiple R-squared: 0.519, Adjusted R-squared: 0.486
## F-statistic: 16.2 on 1 and 15 DF, p-value: 0.00111
截距項和X1的系數都通過了顯著性檢驗，因此回歸方程為 \(Y=62.569+1.229X_1\)
6.3已知如下數據，由表6.19所示
序號 X Y
1 1 0.6
2 1 1.6
3 1 0.5
4 1 1.2
5 2 2.0
6 2 1.3
7 2 2.5
8 3 2.2
9 3 2.4
10 3 1.2
11 4 3.5
12 4 4.1
13 4 5.1
14 5 5.7
15 6 3.4
16 6 9.7
17 6 8.6
18 7 4.0
19 7 5.5
20 7 10.5
21 8 17.5
22 8 13.4
23 8 4.5
24 9 30.4
25 11 12.4
26 12 13.4
27 12 26.2
28 12 7.4
（1）畫出數據的散點圖，求回歸直線\(y= \widehat{\beta_0}+\widehat{\beta_1}x\)，同時將回歸直線也畫在散點圖上。
（2）分析T檢驗和F檢驗是否通過
（3）畫出殘差（普通殘差和標准化殘差）與預測值的殘差圖，分析誤差是否是等方差的
（4）修正模型，對響應變量Y作開方，再完成（1）-（3）的工作
data <- data.frame(x = c(1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6,
7, 7, 7, 8, 8, 8, 9, 11, 12, 12, 12), y = c(0.6, 1.6, 0.5, 1.2, 2, 1.3,
2.5, 2.2, 2.4, 1.2, 3.5, 4.1, 5.1, 5.7, 3.4, 9.7, 8.6, 4, 5.5, 10.5, 17.5,
13.4, 4.5, 30.4, 12.4, 13.4, 26.2, 7.4))
# 求回歸方程
lm.sol <- lm(y ~ x + 1, data)
# 畫出散點圖
plot(data$x, data$y)
abline(lm.sol)
summary(lm.sol)
##
## Call:
## lm(formula = y ~ x + 1, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.841 -2.337 -0.021 1.059 17.832
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.452 1.835 -0.79 0.44
## x 1.558 0.281 5.55 7.9e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.17 on 26 degrees of freedom
## Multiple R-squared: 0.542, Adjusted R-squared: 0.525
## F-statistic: 30.8 on 1 and 26 DF, p-value: 7.93e-06
# 截距項的t value = -0.79， Pr(>|t|)= 0.44>0.05, 此項不顯著，沒有通過t檢驗
# F-statistic: 30.8 on 1 and 26 DF, p-value: 7.93e-06 <0.05，表示通過F檢驗
# 殘差屬於回歸診斷問題 預測值
y.fit <- predict(lm.sol)
# 普通殘差，殘差在不斷增加，因此誤差為異方差
y.res.nor <- residuals(lm.sol)
plot(y.res.nor ~ y.fit)
abline(1.5, 0.6)
abline(-0.1, -0.6)
# 標准化殘差（學生化殘差），殘差在不斷增加，因此誤差為異方差
y.res.nor <- rstandard(lm.sol)
plot(y.res.nor ~ y.fit)
abline(0.3, 0.1)
abline(-0.1, -0.1)
##################################################################################
################################################################################## 修正模型，對響應變量Y作開方，再完成（1）-（3）的工作
################################################################################## 求回歸方程
data$y <- sqrt(data$y)
lm.sol <- lm(y ~ x + 1, data)
# 畫出散點圖
plot(data$x, data$y)
abline(lm.sol)
summary(lm.sol)
##
## Call:
## lm(formula = y ~ x + 1, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5426 -0.4528 -0.0118 0.3493 2.1249
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7665 0.2559 3.00 0.006 **
## x 0.2914 0.0391 7.44 6.6e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.721 on 26 degrees of freedom
## Multiple R-squared: 0.681, Adjusted R-squared: 0.668
## F-statistic: 55.4 on 1 and 26 DF, p-value: 6.64e-08
# 兩個系數的 Pr(>|t|)都小於0.05, 通過t檢驗 F-statistic: 55.4 on 1 and 26
# DF, p-value: 6.64e-08，表示通過F檢驗
# 殘差屬於回歸診斷問題 預測值
y.fit <- predict(lm.sol)
# 普通殘差，殘差在一定范圍內變化，因此誤差為同方差
y.res.nor <- residuals(lm.sol)
plot(y.res.nor ~ y.fit)
abline(1, 0)
abline(-1, 0)
# 標准化殘差（學生化殘差），殘差在不斷增加，因此誤差為異方差
y.res.nor <- rstandard(lm.sol)
plot(y.res.nor ~ y.fit)
abline(1.3, 0)
abline(-1.3, 0)
得到最終的回歸模型為\(\sqrt{y}=0.7665+0.2914x\)
6.4對牙膏銷售數據（數據表見例6.9）得到的線性模型做回歸診斷，分析哪些樣本點需要作進一步的研究？哪些樣本點需要在回歸計算中刪去，如果有，刪去再做線性回歸模型的計算。
toothpaste <- data.frame(X1 = c(-0.05, 0.25, 0.6, 0, 0.25, 0.2, 0.15, 0.05,
-0.15, 0.15, 0.2, 0.1, 0.4, 0.45, 0.35, 0.3, 0.5, 0.5, 0.4, -0.05, -0.05,
-0.1, 0.2, 0.1, 0.5, 0.6, -0.05, 0, 0.05, 0.55), X2 = c(5.5, 6.75, 7.25,
5.5, 7, 6.5, 6.75, 5.25, 5.25, 6, 6.5, 6.25, 7, 6.9, 6.8, 6.8, 7.1, 7, 6.8,
6.5, 6.25, 6, 6.5, 7, 6.8, 6.8, 6.5, 5.75, 5.8, 6.8), Y = c(7.38, 8.51,
9.52, 7.5, 9.33, 8.28, 8.75, 7.87, 7.1, 8, 7.89, 8.15, 9.1, 8.86, 8.9, 8.87,
9.26, 9, 8.75, 7.95, 7.65, 7.27, 8, 8.5, 8.75, 9.21, 8.27, 7.67, 7.93, 9.26))
# 例子中得到的模型
lm.sol <- lm(Y ~ X1 + X2 + I(X2^2) + X1:X2, data = toothpaste)
summary(lm.sol)
##
## Call:
## lm(formula = Y ~ X1 + X2 + I(X2^2) + X1:X2, data = toothpaste)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.4372 -0.1175 0.0049 0.1226 0.3841
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.113 7.483 3.89 0.00066 ***
## X1 11.134 4.446 2.50 0.01915 *
## X2 -7.608 2.469 -3.08 0.00496 **
## I(X2^2) 0.671 0.203 3.31 0.00282 **
## X1:X2 -1.478 0.667 -2.21 0.03610 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.206 on 25 degrees of freedom
## Multiple R-squared: 0.921, Adjusted R-squared: 0.908
## F-statistic: 72.8 on 4 and 25 DF, p-value: 2.11e-13
# 計算殘差
y.res <- rstandard(lm.sol)
# 驗證殘差的正態性，p-value = 0.931>0.05，能通過正態性檢驗
shapiro.test(y.res)
##
## Shapiro-Wilk normality test
##
## data: y.res
## W = 0.9822, p-value = 0.8817
plot(y.res)
# 通過標准化殘差圖可以看出，點5和，點11的殘差最大，去掉這兩個點后重新計算
lm.sol.n <- lm(Y ~ X1 + X2 + I(X2^2) + X1:X2, data = toothpaste[c(-5, -11),
])
summary(lm.sol.n)
##
## Call:
## lm(formula = Y ~ X1 + X2 + I(X2^2) + X1:X2, data = toothpaste[c(-5,
## -11), ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3560 -0.1126 -0.0063 0.0965 0.3348
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 23.378 6.538 3.58 0.0016 **
## X1 10.292 3.790 2.72 0.0123 *
## X2 -5.667 2.163 -2.62 0.0153 *
## I(X2^2) 0.508 0.178 2.85 0.0090 **
## X1:X2 -1.324 0.570 -2.32 0.0295 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.173 on 23 degrees of freedom
## Multiple R-squared: 0.944, Adjusted R-squared: 0.935
## F-statistic: 97.6 on 4 and 23 DF, p-value: 4.42e-14
# 新的模型可以看到，Residual standard error: 從0.206下降到0.173 Multiple
# R-squared: 從0.921上升到 0.944，模型變好
6.5診斷水泥數據（數據見例6.10）是否存在多重共線性，分析例6.10中step()函數去掉的變量是否合理。
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))
# kappa(z, exact = FALSE, ...) 計算矩陣的條件數 κ < 100 不存在多重共線性
# 100 ≤ κ ≤ 1000, 存在一定程度多重共線性 κ > 1000, 存在嚴重多重共線性
XX <- cor(cement[, 1:4])
kappa(XX, exact = TRUE)
## [1] 1377
# κ =1377 > 1000, 存在嚴重多重共線性
# 下一步要找出哪些變量有多重共線性，計算矩陣的特征值和特征向量
eigen(XX)
## $values
## [1] 2.235704 1.576066 0.186606 0.001624
##
## $vectors
## [,1] [,2] [,3] [,4]
## [1,] 0.4760 0.5090 0.6755 0.2411
## [2,] 0.5639 -0.4139 -0.3144 0.6418
## [3,] -0.3941 -0.6050 0.6377 0.2685
## [4,] -0.5479 0.4512 -0.1954 0.6767
得到 $$\lambda _{min}=0.001624, \varphi =(0.2411,0.6418,0.2685,0.6767)^{T}$$
因此有 $$0.2411X1+0.6418X2+0.2685X3+0.6767X4\approx 0$$
#
# 分析例6.10中step()函數去掉的變量是否合理，例子中去掉了變量X3和X4，從多重共線性角度分析去掉是否合理
XX2 <- cor(cement[, 1:2])
kappa(XX2, exact = TRUE)
## [1] 1.593
# κ =1.593，已經消除了共線性問題
6.6為研究一些因素（如用抗生素、有無危險因子和事先是否有計划）對“破腹產后是否有感染”的影響，表6.20給出的是某醫院破腹產后的數據，試用logistic回歸模型對這些數據進行研究，分析感染與這些因素的關系。
表6.20 某醫院進行破腹產后的數據
---------------------------事先有計划--------------------臨時決定
-----------------------有感染----無感染-----------有感染-----無感染
用抗生素--有危險因子----1----------17----------------------11--87
用抗生素--無危險因子----0----------2-----------------------0---0
不用抗生素--有危險因子--28---------30---------------------23--3
不用抗生素--無危險因子--8----------32---------------------0---9
修改下表格的形式
是否用抗生素 有無危險因子 事先有無計划 有感染 無感染
用 有 有 1 17
用 有 無 11 87
用 無 有 0 2
用 無 無 0 0
不 有 有 28 30
不 有 無 23 3
不 無 有 8 32
不 無 無 0 9
# x1 是否使用抗生素 1用 0 不用 x2 有無危險因子 1有 0 無 x3 事先有無計划
# 1有 0 無
patient <- data.frame(x1 = c(1, 1, 1, 1, 0, 0, 0, 0), x2 = c(1, 1, 0, 0, 1,
1, 0, 0), x3 = c(1, 0, 1, 0, 1, 0, 1, 0), y.yes = c(1, 11, 0, 0, 28, 23,
8, 0), y.no = c(17, 87, 2, 0, 30, 3, 32, 9))
patient$Ymat <- cbind(patient$y.yes, patient$y.no)
glm.sol <- glm(Ymat ~ x1 + x2 + x3, family = binomial, data = patient)
summary(glm.sol)
##
## Call:
## glm(formula = Ymat ~ x1 + x2 + x3, family = binomial, data = patient)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## 0.2647 -0.0716 -0.1523 0.0000 -0.7852 1.4962 1.2156 -2.5623
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.821 0.495 -1.66 0.097 .
## x1 -3.254 0.481 -6.76 1.4e-11 ***
## x2 2.030 0.455 4.46 8.2e-06 ***
## x3 -1.072 0.425 -2.52 0.012 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 83.491 on 6 degrees of freedom
## Residual deviance: 10.997 on 3 degrees of freedom
## AIC: 36.18
##
## Number of Fisher Scoring iterations: 5
# 可以看到截距項不能通過驗證，去掉截距項重新回歸，回歸方程的都通過了檢驗
glm.sol <- glm(Ymat ~ x1 + x2 + x3 - 1, family = binomial, data = patient)
summary(glm.sol)
##
## Call:
## glm(formula = Ymat ~ x1 + x2 + x3 - 1, family = binomial, data = patient)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## 0.640 -0.156 -0.153 0.000 -0.384 0.794 0.396 -3.532
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## x1 -3.596 0.458 -7.85 4.0e-15 ***
## x2 1.577 0.352 4.49 7.2e-06 ***
## x3 -1.545 0.329 -4.70 2.7e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 132.439 on 7 degrees of freedom
## Residual deviance: 13.868 on 4 degrees of freedom
## AIC: 37.05
##
## Number of Fisher Scoring iterations: 4
得到破腹產后是否有感染的概率模型：$$P=\frac{exp(-3.596x_1+1.577x_2-1.545x_3)}{1+exp(-3.596x_1+1.577x_2-1.545x_3)}$$
6.7一位飲食公司的分析人員想調查自助餐館中的自動咖啡售貨機數量與咖啡銷售量之間的關系，她選擇了14家餐館來進行實驗，這14家餐館在營業額、顧客類型和地理位置方面都是相近的，放在試驗餐館的自動售貨機數量從0（這里由咖啡服務員端來）到6不等，並且是隨機分配到每個餐館的，表6.21所示的是關於試驗結果的數據。
（1）作線性回歸模型
（2）作多項式回歸模型
（3）畫出數據的散點圖和擬合曲線
表6.21: 自動咖啡售貨機數量與咖啡銷售量數據
餐館 售貨機數量 咖啡銷售量
1 0 508.1
2 0 498.4
3 1 568.2
4 1 577.3
5 2 651.7
6 2 657.0
7 3 713.4
8 3 697.5
9 4 755.3
10 4 758.9
11 5 787.6
12 5 792.1
13 6 841.4
14 6 831.8
coffee <- data.frame(x = c(0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6), y = c(508.1,
498.4, 568.2, 577.3, 651.7, 657, 713.4, 697.5, 755.3, 758.9, 787.6, 792.1,
841.4, 831.8))
# 做線性回歸模型
lm.sol <- lm(y ~ x, coffee)
summary(lm.sol)
##
## Call:
## lm(formula = y ~ x, data = coffee)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.40 -11.48 -3.78 14.63 24.92
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 523.80 8.47 61.8 < 2e-16 ***
## x 54.89 2.35 23.4 2.3e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.6 on 12 degrees of freedom
## Multiple R-squared: 0.978, Adjusted R-squared: 0.977
## F-statistic: 545 on 1 and 12 DF, p-value: 2.26e-11
plot(coffee$x, coffee$y)
abline(lm.sol)
# 雖然線性回歸模型通過了檢驗，但是從圖上看可以繼續做多項式回歸
lm.sol2 <- lm(y ~ 1 + x + I(x^2), coffee)
summary(lm.sol2)
##
## Call:
## lm(formula = y ~ 1 + x + I(x^2), data = coffee)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.664 -5.662 -0.465 5.500 10.668
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 502.556 4.850 103.62 < 2e-16 ***
## x 80.386 3.786 21.23 2.8e-10 ***
## I(x^2) -4.249 0.606 -7.01 2.2e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.86 on 11 degrees of freedom
## Multiple R-squared: 0.996, Adjusted R-squared: 0.995
## F-statistic: 1.39e+03 on 2 and 11 DF, p-value: 5.95e-14
# 多項式回歸的系數都通過了檢驗 標准話殘差由
# 17.6下降到7.86，R平方由0.978上升到0.996，可以認為此多項式回歸優於線性回歸
# 擬合曲線為
xfit <- seq(0, 6.5, len = 200)
yfit <- predict(lm.sol2, data.frame(x = xfit))
plot(coffee$x, coffee$y)
lines(xfit, yfit)
得到的回歸方程為 $$y=502.556+80.386x-4.249x^2$$
6.8 表6.22是40名肺癌病人的生存資料，其中\(X_1\)表示生活行動能力評分（1-100），\(X_2\)表示病人的年齡， \(X_3\)表示由診斷到直入研究時間（月）， \(X_4\)表示腫瘤的類型（“0”是磷癌，“1”是小型細胞癌，“2”是腺癌，“3”是大型細胞癌）；\(X_5\)表示兩種化療方法（“1”常規，“2”試驗新方法）；\(Y\)表示病人的生存時間（“0”是生存時間短，即生存時間小於200天，“1”表示生存時間長，即生存時間大於或等於200天）。
序號 X1 X2 X3 X4 X5 Y
1 70 64 5 1 1 1
2 60 63 9 1 1 0
3 70 65 11 1 1 0
4 40 69 10 1 1 0
5 40 63 58 1 1 0
6 70 48 9 1 1 0
7 70 48 11 1 1 0
8 80 63 4 2 1 0
9 60 63 14 2 1 0
10 30 53 4 2 1 0
11 80 43 12 2 1 0
12 40 55 2 2 1 0
13 60 66 25 2 1 1
14 40 67 23 2 1 0
15 20 61 19 3 1 0
16 50 63 4 3 1 0
17 50 66 16 0 1 0
18 40 68 12 0 1 0
19 80 41 12 0 1 1
20 70 53 8 0 1 1
21 60 37 13 1 1 0
22 90 54 12 1 0 1
23 50 52 8 1 0 1
24 70 50 7 1 0 1
25 20 65 21 1 0 0
26 80 52 28 1 0 1
27 60 70 13 1 0 0
28 50 40 13 1 0 0
29 70 36 22 2 0 0
30 40 44 36 2 0 0
31 30 54 9 2 0 0
32 30 59 87 2 0 0
33 40 69 5 3 0 0
34 60 50 22 3 0 0
35 80 62 4 3 0 0
36 70 68 15 0 0 0
37 30 39 4 0 0 0
38 60 49 11 0 0 0
39 80 64 10 0 0 1
40 70 67 18 0 0 1
（1）建立\(P(Y=1)\) 對 \(X_1~X_5\) 的logistic回歸模型，\(X_1~X_5\)對\(P(Y=1)\)綜合影響是否顯著？哪些變量是主要的影響因素，顯著水平如何？計算各病人生存時間大於等於200天的概率估計值。
（2——用逐步回歸法選取自變量，結果如何？在所選模型下，計算病人生存時間大於等於200天的概率估計值，並將計算結果與（1）中的模型做比較，差異如何？哪一個模型更合理。
patient <- data.frame(X1 = c(70, 60, 70, 40, 40, 70, 70, 80, 60, 30, 80, 40,
60, 40, 20, 50, 50, 40, 80, 70, 60, 90, 50, 70, 20, 80, 60, 50, 70, 40,
30, 30, 40, 60, 80, 70, 30, 60, 80, 70), X2 = c(64, 63, 65, 69, 63, 48,
48, 63, 63, 53, 43, 55, 66, 67, 61, 63, 66, 68, 41, 53, 37, 54, 52, 50,
65, 52, 70, 40, 36, 44, 54, 59, 69, 50, 62, 68, 39, 49, 64, 67), X3 = c(5,
9, 11, 10, 58, 9, 11, 4, 14, 4, 12, 2, 25, 23, 19, 4, 16, 12, 12, 8, 13,
12, 8, 7, 21, 28, 13, 13, 22, 36, 9, 87, 5, 22, 4, 15, 4, 11, 10, 18), X4 = c(1,
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 0, 0, 0, 0, 1, 1, 1, 1, 1,
1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 0, 0, 0, 0, 0), X5 = c(1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0), Y = c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 1))
# 建立$P(Y=1)$ 對 $X_1~X_5$ 的logistic回歸模型
glm.sol <- glm(Y ~ X1 + X2 + X3 + X4 + X5, data = patient)
summary(glm.sol)
##
## Call:
## glm(formula = Y ~ X1 + X2 + X3 + X4 + X5, data = patient)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.6427 -0.3012 -0.0658 0.2848 0.8126
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.322177 0.457570 -0.70 0.4862
## X1 0.010246 0.003549 2.89 0.0067 **
## X2 0.003468 0.006250 0.55 0.5826
## X3 0.000788 0.004217 0.19 0.8528
## X4 -0.118238 0.066063 -1.79 0.0824 .
## X5 -0.117326 0.125824 -0.93 0.3577
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.1495)
##
## Null deviance: 7.5000 on 39 degrees of freedom
## Residual deviance: 5.0838 on 34 degrees of freedom
## AIC: 45
##
## Number of Fisher Scoring iterations: 2
# 可見只有X1的系數滿足了顯著性要求，其他自變量對Y的影響不顯著。
yfit <- predict(glm.sol, patient[, 1:5])
p <- exp(yfit)/(1 + exp(yfit))
# 得到各病人生存時間大於等於200天的概率估計值
cbind(patient, p)
## X1 X2 X3 X4 X5 Y p
## 1 70 64 5 1 1 1 0.5952
## 2 60 63 9 1 1 0 0.5702
## 3 70 65 11 1 1 0 0.5972
## 4 40 69 10 1 1 0 0.5248
## 5 40 63 58 1 1 0 0.5291
## 6 70 48 9 1 1 0 0.5825
## 7 70 48 11 1 1 0 0.5829
## 8 80 63 4 2 1 0 0.5903
## 9 60 63 14 2 1 0 0.5420
## 10 30 53 4 2 1 0 0.4547
## 11 80 43 12 2 1 0 0.5750
## 12 40 55 2 2 1 0 0.4816
## 13 60 66 25 2 1 1 0.5467
## 14 40 67 23 2 1 0 0.4961
## 15 20 61 19 3 1 0 0.4103
## 16 50 63 4 3 1 0 0.4849
## 17 50 66 16 0 1 0 0.5779
## 18 40 68 12 0 1 0 0.5537
## 19 80 41 12 0 1 1 0.6299
## 20 70 53 8 0 1 1 0.6149
## 21 60 37 13 1 1 0 0.5488
## 22 90 54 12 1 0 1 0.6634
## 23 50 52 8 1 0 1 0.5643
## 24 70 50 7 1 0 1 0.6120
## 25 20 65 21 1 0 0 0.5016
## 26 80 52 28 1 0 1 0.6415
## 27 60 70 13 1 0 0 0.6053
## 28 50 40 13 1 0 0 0.5550
## 29 70 36 22 2 0 0 0.5746
## 30 40 44 36 2 0 0 0.5080
## 31 30 54 9 2 0 0 0.4858
## 32 30 59 87 2 0 0 0.5055
## 33 40 69 5 3 0 0 0.4941
## 34 60 50 22 3 0 0 0.5321
## 35 80 62 4 3 0 0 0.5893
## 36 70 68 15 0 0 0 0.6554
## 37 30 39 4 0 0 0 0.5309
## 38 60 49 11 0 0 0 0.6157
## 39 80 64 10 0 0 1 0.6742
## 40 70 67 18 0 0 1 0.6551
# 用逐步回歸法選取自變量
glm.new <- step(glm.sol)
## Start: AIC=45
## Y ~ X1 + X2 + X3 + X4 + X5
##
## Df Deviance AIC
## - X3 1 5.09 43.0
## - X2 1 5.13 43.4
## - X5 1 5.21 44.0
## <none> 5.08 45.0
## - X4 1 5.56 46.6
## - X1 1 6.33 51.8
##
## Step: AIC=43.04
## Y ~ X1 + X2 + X4 + X5
##
## Df Deviance AIC
## - X2 1 5.14 41.4
## - X5 1 5.23 42.2
## <none> 5.09 43.0
## - X4 1 5.57 44.6
## - X1 1 6.38 50.1
##
## Step: AIC=41.41
## Y ~ X1 + X4 + X5
##
## Df Deviance AIC
## - X5 1 5.26 40.3
## <none> 5.14 41.4
## - X4 1 5.61 42.9
## - X1 1 6.38 48.1
##
## Step: AIC=40.35
## Y ~ X1 + X4
##
## Df Deviance AIC
## <none> 5.26 40.3
## - X4 1 5.75 41.9
## - X1 1 6.51 46.9
# 最后只選擇了變量X1和X4 ，如果取顯著性水平為0.1，則所有系數都通過了檢驗
summary(glm.new)
##
## Call:
## glm(formula = Y ~ X1 + X4, data = patient)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5381 -0.3203 -0.0427 0.3250 0.7991
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.15086 0.23036 -0.65 0.5166
## X1 0.00984 0.00331 2.97 0.0052 **
## X4 -0.11941 0.06433 -1.86 0.0714 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.1421)
##
## Null deviance: 7.500 on 39 degrees of freedom
## Residual deviance: 5.258 on 37 degrees of freedom
## AIC: 40.35
##
## Number of Fisher Scoring iterations: 2
yfit.new <- predict(glm.sol, patient[, 1:5])
p.new <- exp(yfit.new)/(1 + exp(yfit.new))
# 得到各病人生存時間大於等於200天的概率估計值
cbind(patient, p.new)
## X1 X2 X3 X4 X5 Y p.new
## 1 70 64 5 1 1 1 0.5952
## 2 60 63 9 1 1 0 0.5702
## 3 70 65 11 1 1 0 0.5972
## 4 40 69 10 1 1 0 0.5248
## 5 40 63 58 1 1 0 0.5291
## 6 70 48 9 1 1 0 0.5825
## 7 70 48 11 1 1 0 0.5829
## 8 80 63 4 2 1 0 0.5903
## 9 60 63 14 2 1 0 0.5420
## 10 30 53 4 2 1 0 0.4547
## 11 80 43 12 2 1 0 0.5750
## 12 40 55 2 2 1 0 0.4816
## 13 60 66 25 2 1 1 0.5467
## 14 40 67 23 2 1 0 0.4961
## 15 20 61 19 3 1 0 0.4103
## 16 50 63 4 3 1 0 0.4849
## 17 50 66 16 0 1 0 0.5779
## 18 40 68 12 0 1 0 0.5537
## 19 80 41 12 0 1 1 0.6299
## 20 70 53 8 0 1 1 0.6149
## 21 60 37 13 1 1 0 0.5488
## 22 90 54 12 1 0 1 0.6634
## 23 50 52 8 1 0 1 0.5643
## 24 70 50 7 1 0 1 0.6120
## 25 20 65 21 1 0 0 0.5016
## 26 80 52 28 1 0 1 0.6415
## 27 60 70 13 1 0 0 0.6053
## 28 50 40 13 1 0 0 0.5550
## 29 70 36 22 2 0 0 0.5746
## 30 40 44 36 2 0 0 0.5080
## 31 30 54 9 2 0 0 0.4858
## 32 30 59 87 2 0 0 0.5055
## 33 40 69 5 3 0 0 0.4941
## 34 60 50 22 3 0 0 0.5321
## 35 80 62 4 3 0 0 0.5893
## 36 70 68 15 0 0 0 0.6554
## 37 30 39 4 0 0 0 0.5309
## 38 60 49 11 0 0 0 0.6157
## 39 80 64 10 0 0 1 0.6742
## 40 70 67 18 0 0 1 0.6551
# 可見新的模型效果更好
得到最終的模型為$$P=\frac{exp(-0.15086+0.00984x_1-0.11941x_4)}{1+exp(-0.15086+0.00984x_1-0.11941x_4)}$$
6.9一位醫院管理人員想建立一個回歸模型，對重傷病人出院后的長期恢復情況進行預測，自變量是病人住院的天數（X），應變量是病人出院后長期恢復的預后指數（Y），指數的數值越大表示預后結局越好，為此，研究了15個別病人的數據，這些數據列在表6.23中，根據經驗表明，病人住院的天數（X）和預后指數（Y）服從非線性模型
\[Y_i=\theta_0exp(\theta_1X_i)+\varepsilon_i,i=1,\cdots, 15 \]
表6.23: 關於重傷病人的數據
病號 住院天數(X) 預后指數(Y )
1 2 54
2 5 50
3 7 45
4 10 37
5 14 35
6 19 25
7 26 20
8 31 16
9 34 18
10 38 13
11 45 8
12 52 11
13 53 8
14 60 4
15 65 6
（1）用內在線性模型方法計算其各種參的估計值；
（2）用非線性方法（nls()函數和nlm()函數）計算其各種參的估計值
patient <- data.frame(x = c(2, 5, 7, 10, 14, 19, 26, 31, 34, 38, 45, 52, 53,
60, 65), y = c(54, 50, 45, 37, 35, 25, 20, 16, 18, 13, 8, 11, 8, 4, 6))
# 用內在線性模型方法計算其各種參的估計值 對自變量做2次正交
poly(patient$x, degree = 2)
## 1 2
## [1,] -0.365891 0.385604
## [2,] -0.327689 0.253961
## [3,] -0.302221 0.173806
## [4,] -0.264019 0.064985
## [5,] -0.213083 -0.058810
## [6,] -0.149412 -0.179320
## [7,] -0.060274 -0.284134
## [8,] 0.003396 -0.313357
## [9,] 0.041598 -0.312633
## [10,] 0.092534 -0.290367
## [11,] 0.181672 -0.192826
## [12,] 0.270810 -0.020734
## [13,] 0.283544 0.009937
## [14,] 0.372682 0.267231
## [15,] 0.436352 0.496657
## attr(,"degree")
## [1] 1 2
## attr(,"coefs")
## attr(,"coefs")$alpha
## [1] 30.73 33.95
##
## attr(,"coefs")$norm2
## [1] 1 15 6167 1727980
##
## attr(,"class")
## [1] "poly" "matrix"
lm.pol <- lm(y ~ 1 + poly(x, 2), data = patient)
summary(lm.pol)
##
## Call:
## lm(formula = y ~ 1 + poly(x, 2), data = patient)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.672 -1.290 0.219 1.384 4.074
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 23.333 0.634 36.81 1.0e-13 ***
## poly(x, 2)1 -59.094 2.455 -24.07 1.6e-11 ***
## poly(x, 2)2 19.464 2.455 7.93 4.1e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.46 on 12 degrees of freedom
## Multiple R-squared: 0.982, Adjusted R-squared: 0.979
## F-statistic: 321 on 2 and 12 DF, p-value: 3.81e-11
# 畫出散點圖和預測曲線
xfit <- seq(2, 65, len = 200)
yfit <- predict(lm.pol, data.frame(x = xfit))
plot(patient$x, patient$y)
lines(xfit, yfit)
使用內在線性模型方法，得到y關於x的二次回歸方程：\(y=23.333-59.094*x+19.464*x^2\)
# 用非線性方法（nls()函數和nlm()函數）計算其各種參的估計值
nls.sol <- nls(y ~ a * exp(b * x), data = patient, start = list(a = 2, b = 0.01))
summary(nls.sol)
##
## Formula: y ~ a * exp(b * x)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 58.60655 1.47216 39.8 5.7e-15 ***
## b -0.03959 0.00171 -23.1 6.0e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.95 on 13 degrees of freedom
##
## Number of iterations to convergence: 7
## Achieved convergence tolerance: 3.48e-06
# 畫出散點圖和預測曲線
xfit <- seq(2, 65, len = 200)
yfit <- predict(nls.sol, data.frame(x = xfit))
plot(patient$x, patient$y)
lines(xfit, yfit)
使用非線性方法方法，得到y關於x的回歸方程：\(y=58.60655*exp(-0.03959*x)\)