Linear Regression Analysis
第二章
2.12
(1)拟合模型：
> library(openxlsx) #加载包 openxlsx
> data = read.xlsx("22_data.xlsx",sheet = 2) #read.xlsx 函数读入数据
> x = data[,1]
> y = data[,2]
> res = lm(y~x) #构造线性回归模型函数
> res #结果
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept) x #得出线性回归模型 y = -6.332 + 9.208 x
-6.332 9.208
> summary(res) #打印方差分析，系数的估计值及其检验。
Call:
lm(formula = y ~ x)
Residuals: #残差统计量，残差第一四分位数(1Q)和第三分位数(3Q)有大约相同的幅度，意味着有较对称的钟形分布
Min 1Q Median 3Q Max
-2.5629 -1.2581 -0.2550 0.8681 4.0581
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.33209 1.67005 -3.792 0.00353 ** #截距的点估计值及其检验
x 9.20847 0.03382 272.255 < 2e-16 *** #自变量系数的点估计值及其检验
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.946 on 10 degrees of freedom
Multiple R-squared: 0.9999, Adjusted R-squared: 0.9999 #相关系数与调整的相关系数
F-statistic: 7.412e+04 on 1 and 10 DF, p-value: < 2.2e-16 #模型的显著性检验（F检验）
(2)根据上面程序结果，自变量具有显著性，模型具有显著性。
(3) 不能支持管理员的观点，根据构造的线性回归模型，平均环境温度增加1°F，平均月水蒸气消耗量将增加 9208+lb ，达不到10000lb.
(4) 使用58°F的平均环境温度构造一个月中水蒸气消耗量的99%预测区间：
> library(openxlsx)
> data = read.xlsx("22_data.xlsx",sheet = 2)
> x = data[,1]
> y = data[,2]
> fun = function(x) #计算预测值函数
+ {
+ y = -6.332 + 9.208*x
+ }
> y_pred = fun(x) #计算所有数的预测值
> s_y0_pred = function(x0,x,y,n) #构造计算预测值标准差的函数
+ {
+ n = 12
+ y_pred = fun(x)
+ sse = sum((y_pred - y)*(y_pred - y))
+ se = sqrt(sse/(n-2))
+ se * sqrt(1/n + ((x0-mean(x))^2)/sum((x-rep(mean(x),n))^2))
+ }
> x0 = 58 ; n = 12
> y0_pred = fun(x0) #当环境温度为58°F,对应的因变量预测值
> s = s_y0_pred(x0,x,y,n)
> print(c(y0_pred-qt(0.995,n-2)*s,y0_pred+qt(0.995,n-2)*s)) #输出结果
[1] 525.5666 529.8974
2.13
a.做出散点图
> data = read.xlsx("22_data.xlsx",sheet = 1)
> x = data[,2]
> y = data[,1]
> plot(x,y,main = "散点图",xlab = "index",ylab = "days")
> abline(lm(y~x))
b.估计预测方程
> lm(y~x)
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept) x
-193.0 15.3
预测方程为：y = *-193.0 + 15.3 x
c.进行回归显著性检验
> summary(lm(y~x))
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-41.70 -21.54 2.12 18.56 36.42
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -192.984 163.503 -1.180 0.258 #p值大于0.05
x 15.296 9.421 1.624 0.127 #p值大于0.05 ， 回归变量与响应变量没有显著相关性
Residual standard error: 23.79 on 14 degrees of freedom
Multiple R-squared: 0.1585, Adjusted R-squared: 0.09835
F-statistic: 2.636 on 1 and 14 DF, p-value: 0.1267
根据上述结果，指数与天数并没有显著相关性。
d.计算并画出95%置信带与95%预测带
> sx = sort(x)
> #计算置信区间 > conf = predict(fm,data.frame(x = sx),interval = "confidence") > #计算预测区间 > pred = predict(fm,data.frame(x=sx),interval = "prediction") > plot(x,y,ylim = c(0,150),xlab = "index",ylab = "days",main = "95%预测带、置信带") > abline(fm) > lines(sx,conf[,2],col = "red") > lines(sx,conf[,3],col = "red") > lines(sx,pred[,2],col = "blue") > lines(sx,pred[,3],col = "blue")
2.14
a.散点图
b.估计预测方程
> fm = lm(y~x)
> fm
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept) x
0.6714 -0.2964
预测方程为：y = 0.6714 - 0.2964 x
c.数据分析
> summary(fm)
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-0.20464 -0.10634 0.02196 0.08527 0.20643
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.6714 0.1595 4.209 0.00563 **
x -0.2964 0.2314 -1.281 0.24754 #p值大于0.05 ，该自变量没有显著相关
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.15 on 6 degrees of freedom
Multiple R-squared: 0.2147, Adjusted R-squared: 0.08382 R^2 = 0.2147
F-statistic: 1.64 on 1 and 6 DF, p-value: 0.2475 #整个模型不具有显著性。
d.计算并画出95%置信带和95%预测带
> plot(x,y,main = "散点图",xlab = "比率",ylab = "黏度",ylim = c(-0.1,1)) > sx = sort(x) > conf = predict(fm,data.frame(sx),interval = "confidence") > pred = predict(fm,data.frame(sx),interval = "prediction") > abline(fm) > lines(sx,conf[,2],col = "red") #绘制置信下限 > lines(sx,conf[,3],col = "red") #绘制置信上限 > lines(sx,pred[,2],col = "blue") #绘制预测下限 > lines(sx,pred[,3],col = "blue") #绘制预测上限
2.15
a.估计预测方程
Call: lm(formula = y ~ x) Coefficients: (Intercept) x 1.281511 -0.008758
预测方程为： y = 1.281511 - 0.008758 x
b.全面分析此模型
> fm = lm(y~x)
> summary(fm)
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-0.043955 -0.035863 -0.009305 0.019900 0.069559
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.2815107 0.0468683 27.34 1.58e-07 ***
x -0.0087578 0.0007284 -12.02 2.01e-05 *** #根据 p 值，自变量温度极显著
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.04743 on 6 degrees of freedom
Multiple R-squared: 0.9602, Adjusted R-squared: 0.9535
F-statistic: 144.6 on 1 and 6 DF, p-value: 2.007e-05 #根据 p 值，整个回归模型是显著的
c.画95%置信带、预测带
2.16
先画出散点图：
从散点图可以看出容量与压力之间具有明显的线性关系，我们构造一元线性模型：
> fm = lm(y~x) > fm Call: lm(formula = y ~ x) Coefficients: (Intercept) x -290.707 2.346
估计预测模型为： y = -290.707 + 2.346x
再对模型进行检验：
> summary(fm)
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-4.3276 -0.9227 0.0773 1.2676 2.9577
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.907e+02 1.355e+00 -214.6 <2e-16 ***
x 2.346e+00 4.007e-04 5855.4 <2e-16 *** #该自变量极显著
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.741 on 31 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 3.429e+07 on 1 and 31 DF, p-value: < 2.2e-16 #整个回归模型极显著
2.17
> x = data[,2]
> y = data[,1]
> n = length(x)
> plot(x,y)
>
> fm = lm(y~x) #一元回归模型
> abline(fm)
> coef(fm)
(Intercept) x
163.930734 1.579647
> summary(fm)
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-1.41483 -0.91550 -0.05148 0.76941 2.72840
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 163.9307 2.6551 61.74 < 2e-16 ***
x 1.5796 0.1051 15.04 1.88e-10 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.291 on 15 degrees of freedom
Multiple R-squared: 0.9378, Adjusted R-squared: 0.9336
F-statistic: 226 on 1 and 15 DF, p-value: 1.879e-10
> anova(fm) #方差分析
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x 1 376.92 376.92 226.04 1.879e-10 ***
Residuals 15 25.01 1.67
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> rm(list = ls())
2.18
library(openxlsx)
data = read.xlsx("2.18.xlsx",sheet = 1)
x = data[,2]
y = data[,3]
n = length(x)
plot(x,y)
fm = lm(y~x) #一元回归
coef(fm) #输出回归系数
summary(fm)
anova(fm)
#构造此数据的95%置信带与预测带
sx = sort(x)
conf = predict(fm,data.frame(x = sx),interval = "confidence")
pred = predict(fm,data.frame(x = sx),interval = "prediction")
abline(fm)
lines(sx,conf[,2],col = 'red')
lines(sx,conf[,3],col = 'red')
lines(sx,pred[,2],col = 'blue')
lines(sx,pred[,3],col = 'blue')
plot(x,y,main = "%95置信带与95%预测带",xlab = "花费钱数",ylab = "每周挣回的印象",ylim=c(-100,200))
rm(list = ls())
第三章
3.8
#a.拟合co2产量y与总溶剂量x6和氢消耗量x7关系的多元回归模型
library(openxlsx) data = read.xlsx("3.8.xlsx",sheet = 1) data y = data[,1] #响应变量 x = data[,c(7,8)] #回归变量 fm = lm(y~.,x) #多元线性回归 summary(fm) anova(fm) #检验显著性 summary(fm) #d confint(fm) #e x6 = data[,7] fm1 = lm(y~x6) summary(fm1) anova(fm1) confint(fm1,level = 0.95) rm(list = ls())
3.9
library(openxlsx)
data = read.xlsx("3.9.xlsx",sheet = 1)
y = data[,1]
x = data[,c(2,5)]
#a.拟合多元回归模型
fm = lm(y~.,x)
coef(fm)
#b,c 检验回归显著性（）
anova(fm)
summary(fm)
#e
#检验是否有潜在的多重共线性
r2 = 0.6367
vif = 1/(1-r2)
rm(list = ls())
3.10
> #3.10
> library(openxlsx)
> data = read.xlsx("3.10.xlsx",sheet = 1)
> #a
> y = data[,6]
> x = data[,c(1:5)]
> fm = lm(y~.,x)
> coef(fm)
(Intercept) Clarity Aroma Body Flavor
3.9968648 2.3394535 0.4825505 0.2731612 1.1683238
Oakiness
-0.6840102
> #b,c
> summary(fm)
Call:
lm(formula = y ~ ., data = x)
Residuals:
Min 1Q Median 3Q Max
-2.85552 -0.57448 -0.07092 0.67275 1.68093
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.9969 2.2318 1.791 0.082775 .
Clarity 2.3395 1.7348 1.349 0.186958
Aroma 0.4826 0.2724 1.771 0.086058 .
Body 0.2732 0.3326 0.821 0.417503
Flavor 1.1683 0.3045 3.837 0.000552 ***
Oakiness -0.6840 0.2712 -2.522 0.016833 *
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.163 on 32 degrees of freedom
Multiple R-squared: 0.7206, Adjusted R-squared: 0.6769
F-statistic: 16.51 on 5 and 32 DF, p-value: 4.703e-08
> #d
> xx = data[,c(2,4)]
> fm1 = lm(y~.,xx)
> summary(fm1)
Call:
lm(formula = y ~ ., data = xx)
Residuals:
Min 1Q Median 3Q Max
-2.19048 -0.60300 -0.03203 0.66039 2.46287
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.3462 1.0091 4.307 0.000127 ***
Aroma 0.5180 0.2759 1.877 0.068849 .
Flavor 1.1702 0.2905 4.027 0.000288 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.229 on 35 degrees of freedom
Multiple R-squared: 0.6586, Adjusted R-squared: 0.639
F-statistic: 33.75 on 2 and 35 DF, p-value: 6.811e-09
> AIC(fm) #优先考虑的模型应是AIC值最小的那一个
[1] 126.7552
> AIC(fm1)
[1] 128.3761
> #e
> conf = confint(fm)
> conf1 = confint(fm1)
> conf = as.matrix(conf)
> conf1 = as.matrix(conf1)
>
> conf[5,2]-conf[5,1]
[1] 1.240414
> conf1[3,2]-conf[3,1]
[1] 1.83241
>
> rm(list = ls())
>
3.11
> #3.11
> library(openxlsx)
> data = read.xlsx("3.11.xlsx",sheet = 1)
> y = data[,6]
> x = data[,c(1:5)]
> #a
> fm = lm(y~.,x)
> coef(fm)
(Intercept) x1 x2 x3
5.207905e+01 5.555556e-02 2.821429e-01 1.250000e-01
x4 x5
1.776357e-16 -1.606498e+01
> #b,c
> summary(fm)
Call:
lm(formula = y ~ ., data = x)
Residuals:
Min 1Q Median 3Q Max
-12.250 -4.438 0.125 5.250 9.500
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.208e+01 1.889e+01 2.757 0.020218 *
x1 5.556e-02 2.987e-02 1.860 0.092544 .
x2 2.821e-01 5.761e-02 4.897 0.000625 ***
x3 1.250e-01 4.033e-01 0.310 0.762949
x4 1.776e-16 2.016e-01 0.000 1.000000
x5 -1.606e+01 1.456e+00 -11.035 6.4e-07 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 8.065 on 10 degrees of freedom
Multiple R-squared: 0.9372, Adjusted R-squared: 0.9058
F-statistic: 29.86 on 5 and 10 DF, p-value: 1.055e-05
> #d
> xx = data[,c(2,5)]
> fm1 = lm(y~.,xx)
> summary(fm1)
Call:
lm(formula = y ~ ., data = xx)
Residuals:
Min 1Q Median 3Q Max
-15.375 -4.188 -0.875 3.438 12.625
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 80.13461 5.69146 14.080 3.01e-09 ***
x2 0.28214 0.05883 4.796 0.000349 ***
x5 -16.06498 1.48659 -10.807 7.26e-08 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 8.236 on 13 degrees of freedom
Multiple R-squared: 0.9149, Adjusted R-squared: 0.9018
F-statistic: 69.89 on 2 and 13 DF, p-value: 1.107e-07
> AIC(fm)
[1] 118.6885
> AIC(fm1)
[1] 117.5552
>
> #e
> confint(fm)
2.5 % 97.5 %
(Intercept) 9.99688896 94.1612109
x1 -0.01100273 0.1221138
x2 0.15378045 0.4105053
x3 -0.77353688 1.0235369
x4 -0.44926844 0.4492684
x5 -19.30879739 -12.8211665
> confint(fm1) #温度:x2
2.5 % 97.5 %
(Intercept) 67.8389462 92.4302647
x2 0.1550559 0.4092298
x5 -19.2765650 -12.8533989
>
第四章
例4.1 根据例3.1数据，输出残差，标准化残差，学生化残差，press残差，外部学生化残差 表格
library(openxlsx)
#例4.1---------------------------------------------------------
#处理数据
data = read.xlsx("3.1.xlsx",sheet = 1)
data = data[,c(2,3,4)]
names(data)=c("time","cases","distance")
y = data$time
x1 = data$cases
x2 = data$distance
#线性回归
fm = lm(y~x1+x2)
#残差
ei = residuals(fm)
View(ei)
#标准化残差(1)
di = rstandard(fm)
View(di)
#计算mse的函数
mse = function(ei,p) #ei是残差向量，p是回归参数个数
{
n = length(ei)
sse = sum(ei**2)
mse = sse/(n-p)
return(mse)
}
di_ = ei/sqrt(mse(ei,3))#标准化残差(2)
View(di_)
#学生化残差(1)
ri = rstudent(fm)
View(ri)
#计算帽子矩阵，并提取对角线元素
H = function(X) #X是回归向量矩阵
{
h = X%*%solve(t(X)%*%X)%*%t(X)
hii = diag(h)
return(hii)
}
X = cbind(1,x1,x2)
hii = H(X) #计算hii
View(hii)
ri_ = ei/sqrt(mse(ei,3)*(1-hii)) #学生化残差(2)
View(ri_)
#计算PRESS统计量
e_i = ei/(1-hii) #计算e(i)
View(e_i)
#外部学生化残差
ti = function(ei,X) #输入残差回归变量矩阵
{
p = ncol(X) #回归参数个数
n = length(ei) #数据个数
hii = H(X) #帽子矩阵主对角线元素
s2_i = ((n-p)*mse(ei,p) -(ei**2)/(1-hii)) / (n-p-1) #计算S(i)^2
ans = ei / sqrt(s2_i*(1-hii))
return(ans)
}
ti = ti(ei,X)
View(ti)
#计算PRESS统计量
press = function(ei,X)
{
hii = H(X)
res = sum((ei/(1-hii))**2)
#View(res)
}
Press = (ei/(1-hii))**2
View(Press)
PRESS = press(ei,X) #输出PRESS统计量
#将所有残差数据写入表格
Num = seq(1,length(ei))
mydata = cbind(Num,ei,di_,ri_,hii,e_i,ti,Press)
class(mydata)
View(mydata)
write.xlsx(mydata,"C:\\Users\\86130\\Desktop\\mydata.xlsx")
#例4.2
ti #外部学生化残差
View(ti)
n = length(ti) #数据个数
order = rank(ti) #rank函数返回ti按升序排序之后的序号
Pi = (order-1/2)/n #累积概率
plot(ti,Pi,xlim=c(-3,5)) #画正态概率图
fm_tP = lm(Pi~ti) #线性回归模型
abline(fm_tP) #添加回归线
#例4.3
#画残差与拟合值y_i的残差图
plot(fitted(fm),ti) #fitted(fm)返回模型fm的预测值
abline(h = 0) #添加直线y=0
#例4.4
#画残差与回归变量的残差图
par(mfrow =c(1,2))
plot(x1,ti,xlab = "箱数",ylab = "外部学生化残差")
abline(h=0) #h:y轴 v:x轴
plot(x2,ti,xlab = "距离",ylab = "外部学生化残差")
abline(h=0)
#例4.5
#画偏回归图
#回归变量x1的偏回归图
lm.y_x2 = lm(y~x2)
lm.x1_x2 = lm(x1~x2)
plot(resid(lm.x1_x2),resid(lm.y_x2),xlab = "箱数",ylab = "时间")
#回归变量x2的偏回归图
lm.y_x1 = lm(y~x1)
lm.x2_x1 = lm(x2~x1)
plot(resid(lm.x2_x1),resid(lm.y_x1),xlab = "距离",ylab = "时间",pch = 10)
#例4.6
#计算PRESS的预测R^2
R_pred = function(X,y)
{
hii = H(X)
ei = resid(lm(y~X[,2]+X[,3]))
PRESS = sum((ei/(1-hii))**2)
sst = sum((y-mean(y))**2)
ans = 1-PRESS/sst
return(ans)
}
R_pred(X,y)
#例4.7
data = read.xlsx("2.1.xlsx",sheet = 1)
names(data) = c("order","y","x")
x = data$x
y = data$y
X = cbind(1,x)
fm = lm(y~x)
#绘制正态概率图
plot_ZP = function(ti) #输入外部学生化残差
{
n = length(ti)
order = rank(ti) #按升序排列，t(i)是第order个
Pi = (order-1/2)/n #累积概率
plot(ti,Pi,xlim=c(-3,3),xlab = "学生化残差",ylab = "百分比") #画正态概率图
}
ei = resid(fm)
ti = ti(ei,X) #计算外部学生化残差ti
plot_ZP(ti) #绘制正态概率图
plot(fitted(fm),ti) #绘制残差与所预测y_pred的残差图
abline(h = 0)
#绘制除去5,6两点的正态概率图
data = data[-c(5,6),]
x = data$x
y = data$y
X = cbind(1,x)
fm1 = lm(y~x) #线性模型
ei = resid(fm1)
ti = ti(ei,X) #计算外部学生化残差ti
plot(fitted(fm1),ti) #绘制残差与所预测y_pred的残差图
abline(h = 0)
#例4.8
data = read.xlsx("4.8.xlsx",sheet = 1)
x = data$x
y = data$y
fm = lm(y~x) #线性回归
a = anova(fm) #方差分析
sst = sum(a[2]) #总平方和
ssg = a[1,2] #回归平方和
sse = a[2,2] #残差平方和
level_x = c(table(x)>1) #查看哪些自变量重复
#进行失拟检验
library(rsm) #加载rsm包用于失拟检验
lm.rsm<-rsm(y~FO(x))
loftest(lm.rsm) #调用失拟检验函数loftest
rm(list = ls())
例4.10 通过近邻点估计纯误差
#例4.10
data = read.xlsx("3.1.xlsx",sheet = 1) #导入数据
names(data)=c("order","time","cases","distance")
y = data$time #准备数据
x1 = data$cases
x2 = data$distance
fm = lm(y~x1 + x2) #线性回归
coef(fm)
b_cases = coef(fm)[2] #beta1
b_distance = coef(fm)[3] #beta2
y_pred = predict(fm) #计算预测值
ei = resid(fm) #残差
new_data = cbind(data,y_pred,ei) #构建新数据
new_data = new_data[order(new_data$y_pred),] #按照预测值升序排序
a = anova(fm) #方差分析
sse = a[3,2] #残差平方和
mse = a[3,3] #残差均方和
#计算Dii'
Di_i = function(i,i_,mse,beta1,beta2,new_data) #i第i个点，i_第i_个点,data数据集
{
one = beta1*(new_data$cases[i]-new_data$cases[i_])/sqrt(mse)
two = beta2*(new_data$distance[i]-new_data$distance[i_])/sqrt(mse)
ans = one**2 + two**2
return(ans)
}
#定义一个数据框用来存储结果
σ_ans = data.frame(
i = numeric(0), #观测值i
i_ = numeric(0), #观测值i_
Dii = numeric(0), #Di_i
delta = numeric(0) #E|ei-ei_|
)
#计算相邻k个点的两点的 Di_i,i,i_,Delta残差
for (k in c(1:4))
{
for (i in c(1:24))
{
if (i+k>25)
break
D = Di_i(i,i+k,mse,b_cases,b_distance,new_data) #计算相邻k个点的两点的Di_i
E = abs(new_data$ei[i]-new_data$ei[i+k]) #计算相邻k个点的两点的Delta残差
another = data.frame(
i = new_data$order[i],
i_ = new_data$order[i+k],
Dii = D,
delta = E
)
σ_ans = rbind(σ_ans,another) #合并两个数据框
}
}
names(σ_ans) = c("i","i_","Dii^2","Delta") #重命名最后的数据框
σ_ans = σ_ans[order(σ_ans$Dii^2),] #按照Di_i对数据框进行排序
row.names(σ_ans) = c(1:90) #对每一行进行编号
#计算累计标准差
std = numeric(0) #存储累计标准差
sum_Delta = 0 #存储累计Delta残差
for (i in 1:90)
{
sum_Delta = sum_Delta + σ_ans$Delta[i] #0.886/m*Σ(Delta)
res = 0.886/i*sum_Delta
std[i] = res
}
σ_ans = cbind(std,σ_ans)
4.16
#######################自己编的函数，运行后直接调用#######################
#计算mse的函数
mse = function(ei,p) #ei是残差向量，p是回归参数个数
{
n = length(ei)
sse = sum(ei**2)
mse = sse/(n-p)
return(mse)
}
#计算帽子矩阵，并提取对角线元素
H = function(X) #X是回归向量矩阵
{
h = X%*%solve(t(X)%*%X)%*%t(X)
hii = diag(h)
return(hii)
}
#外部学生化残差
ti = function(ei,X) #输入残差回归变量矩阵
{
p = ncol(X) #回归参数个数
n = length(ei) #数据个数
hii = H(X) #帽子矩阵主对角线元素
s2_i = ((n-p)*mse(ei,p) -(ei**2)/(1-hii)) / (n-p-1) #计算S(i)^2
ans = ei / sqrt(s2_i*(1-hii))
return(ans)
}
#计算PRESS统计量
press = function(ei,X) #X是自变量的设计矩阵
{
hii = H(X)
res = sum((ei/(1-hii))**2)
#View(res)
}
#计算PRESS的预测R^2
R_pred = function(X,y)
{
hii = H(X)
ei = resid(lm(y~X[,2]+X[,3]))
PRESS = sum((ei/(1-hii))**2)
sst = sum((y-mean(y))**2)
ans = 1-PRESS/sst
return(ans)
}
#绘制正态概率图
plot_ZP = function(ti) #输入外部学生化残差
{
n = length(ti)
order = rank(ti) #按升序排列，t(i)是第order个
Pi = (order-1/2)/n #累积概率
plot(ti,Pi,xlim=c(-3,3),xlab = "学生化残差",ylab = "百分比") #画正态概率图
}
#进行失拟检验
library(rsm) #加载rsm包用于失拟检验
lm.rsm<-rsm(y~FO(x))
loftest(lm.rsm) #调用失拟检验函数loftest
#计算Dii'
Di_i = function(i,i_,mse,beta1,beta2,new_data) #i第i个点，i_第i_个点,data数据集
{
one = beta1*(new_data$cases[i]-new_data$cases[i_])/sqrt(mse)
two = beta2*(new_data$distance[i]-new_data$distance[i_])/sqrt(mse)
ans = one**2 + two**2
return(ans)
}
#4.16
#a.残差的正态概率图
data = read.xlsx('3.12.xlsx',sheet = 1) #导入数据
y = data$y
x1 = data$x1
x2 = data$x2
X = cbind(1,x1,x2) #处理数据
fm = lm(y~x1+x2) #线性回归模型
ei = resid(fm) #计算残差
ti = ti(ei,X) #ti()自制求外部学生化残差函数
plot_ZP(ti) #plot_zp()自制绘制正态概率图函数
#为什么要编写函数？
#1.这些题目都是重复的代码操作
#2.如果是想多次重复打代码来熟悉，大可不必，因为会忘的。
#正态概率图有一个异常点，order(ti) 返回第一小的是第28号点
#b.残差与响应变量预测值的残差图
plot(fitted(fm),ti)
#残差图表明残差包含在一条水平带中，模型不存在明显的缺点。
#c.
#模型fm的PRESS统计量
press_fm = press(ei,X)
#新模型fm1的PRESS统计量
fm1 = lm(y~x2)
ei = resid(fm1)
X = cbind(1,x2)
press_fm1 = press(ei,X) #press()自制求press统计量函数
#选择press统计量小的模型
