#案例7.1 回归诊断Anscombe数据
rm(list=ls())
Anscombe<-data.frame(
  X=c(10.0, 8.0, 13.0, 9.0, 11.0, 14.0, 6.0, 4.0, 12.0, 7.0, 5.0),
  Y1=c(8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68),
  Y2=c(9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74),
  Y3=c(7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.44, 5.73),
  X4=c(rep(8,7), 19, rep(8,3)),
  Y4=c(6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89)
)
summary(lm(Y1~X, data=Anscombe))
summary(lm(Y2~X, data=Anscombe))
summary(lm(Y3~X, data=Anscombe))
summary(lm(Y4~X4,data=Anscombe))
par(mar = c(3, 2, 0.5, 0.5), mfrow = c(2, 2))
plot(c(3,20), c(3,13), type="n", xlab = "X", ylab = "Y1")
points(Anscombe$X,Anscombe$Y1)
abline(lm(Anscombe$Y1~Anscombe$X))
text(11,12.5,expression((a)),cex=1)
plot(c(3,20), c(3,13), type="n", xlab = "X", ylab = "Y2")
points(Anscombe$X,Anscombe$Y2)
abline(lm(Anscombe$Y2~Anscombe$X))
text(11,12.5,expression((b)),cex=1)
plot(c(3,20), c(3,13), type="n", xlab = "X", ylab = "Y3")
points(Anscombe$X,Anscombe$Y3)
abline(lm(Anscombe$Y3~Anscombe$X))
text(11,12.5,expression((c)),cex=1)
plot(c(3,20), c(3,13), type="n", xlab = "X", ylab = "Y4")
points(Anscombe$X4,Anscombe$Y4)
abline(lm(Anscombe$Y4~Anscombe$X4))
text(11,12.5,expression((d)),cex=1)

#图7.2 二次函数拟合图
X2<-Anscombe$X^2
lm2.sol<-lm(Anscombe$Y2~Anscombe$X+X2)
summary(lm2.sol)
x<-seq(min(Anscombe$X), max(Anscombe$X), by=0.1)
b<-coef(lm2.sol)
y<-b[1]+b[2]*x+b[3]*x^2
plot(c(3,20), c(3,13), type="n", xlab = "X", ylab = "Y2")
points(Anscombe$X,Anscombe$Y2)
lines(x,y)

#图7.3 剔除异常点后的拟合图
i<-1:11; Y31<-Anscombe$Y3[i!=3]; X3<-Anscombe$X[i!=3]
lm3.sol<-lm(Y31~X3)
summary(lm3.sol)
plot(c(3,20), c(3,13), type="n", xlab = "X", ylab = "Y3")
points(Anscombe$X,Anscombe$Y3)
abline(lm3.sol)

#案例7.2 残差的正态性检验
rm(list=ls())
DA=read.csv("szzs.csv",header=T)
DA=DA[2:28]
lm.sz=lm(DA$cf~DA$jr+DA$dc)
y.res<-residuals(lm.sz)#普通残差的计算
shapiro.test(y.res)

#案例7.4 
rm(list=ls())
DA=read.csv("szzs.csv",header=T)
lm.sz=lm(cf~0+.,data=DA[,2:28])
summary(lm.sz)

#案例7.3 残差
par(mar = c(4, 4, 0.5, 0.5), mfrow = c(1, 2))
y.rst=rstandard(lm.sz)
y.fit=predict(lm.sz)
plot(y.rst~y.fit)#以拟合值为横坐标的残差图（标准化残差散点图）
y.rst<-rstudent(lm.sz); 
plot(y.rst~y.fit)#以拟合值为横坐标的残差图（外学生化残差散点图）

#图7.6 修正模型后的残差图
lm.new<-update(lm.sz, sqrt(sqrt(.))~.)
coef(lm.new)
par(mar = c(4, 4, 0.5, 0.5), mfrow = c(1, 2))
y.rst<-rstandard(lm.new); y.fit<-predict(lm.sz)
plot(y.rst~y.fit)#以拟合值为横坐标的残差图（标准化残差散点图）
y.rst<-rstudent(lm.new)
plot(y.rst~y.fit)#以拟合值为横坐标的残差图（外学生化残差散点图）

#图7.7
plot(lm.sz,2)#对数正态QQ残差图

#条件数的计算
XX<-cor(DA[3:28])#变量相关系数矩阵
kappa(XX,exact=TRUE)#求矩阵的条件数
eigen(XX)#求矩阵的特征值

#模型（7.15）均方误差计算
sum(eigen(solve(XX))$values)

#案例7.5 
p=1
n=nrow(DA)
d1=dffits(lm.sz)
cf=1:n
cf[abs(d1)>2*sqrt((p+1)/n)]

#Cook距离计算
infl=lm.influence(lm.sz,do.coef=FALSE)
D=cooks.distance(lm.sz,infl=lm.influence(lm.sz,do.coef=FALSE),
rs=weighted.residuals(lm.sz),
sd=sqrt(deviance(lm.sz)/df.residual(lm.sz)),
hat=infl$hat)
sort(D)#按从小到大的顺序排列

#协方差比诊断
D=abs(1-covratio(lm.sz, infl=lm.influence(lm.sz, do.coef = FALSE),res = 
weighted.residuals(lm.sz)))
sort(D)

#回归诊断函数Reg_Diag
Reg_Diag<-function(fm){
   n<-nrow(fm$model); df<-fm$df.residual
   p<-n-df-1; s<-rep(" ", n);  
   res<-residuals(fm); s1<-s; s1[abs(res)==max(abs(res))]<-"*"
   sta<-rstandard(fm); s2<-s; s2[abs(sta)>2]<-"*"
   stu<-rstudent(fm); s3<-s; s3[abs(sta)>2]<-"*"
   h<-hatvalues(fm); s4<-s; s4[h>2*(p+1)/n]<-"*"
   d<-dffits(fm); s5<-s; s5[abs(d)>2*sqrt((p+1)/n)]<-"*"
   c<-cooks.distance(fm); s6<-s; s6[c==max(c)]<-"*"
   co<-covratio(fm);   abs_co<-abs(co-1)
   s7<-s; s7[abs_co==max(abs_co)]<-"*"
   data.frame(residual=res, s1, standard=sta, s2, 
              student=stu, s3,  hat_matrix=h, s4, 
              DFFITS=d, s5,cooks_distance=c, s6, 
              COVRATIO=co, s7)}
Reg_Diag(lm.sz)