dir <- "G:/实用回归分析/" li5.1 <- read.csv(paste(dir, "li5.1.csv", sep=""),header = T) head(li5.1) li5.1 <- li5.1[,-(1:2)] y <- li5.1[,1] x <- cbind(1, as.matrix(li5.1[,-1])) # sample size n <- length(y) #dim of variable p <- dim(x)[2]-1 #original least square lm5.1 <- lm(y~., data=li5.1) #x <- scale(x)/sqrt(n-1) # anova(lm(y~1), lm5.1) # #residual vector e <- lm5.1$residuals #fitted value yhat <- lm5.1$fitted.values #残差图 plot(yhat, e) plot(e~yhat) #leverage vector # #HX=X # #H1=1 H <- x %*% solve(t(x) %*% x) %*% t(x) apply(H, 1, sum) apply(H, 2, sum) h <- diag(H) round(h, 4) hatvalues(lm5.1) #deleting residual ei <- as.vector(e/(1-h)) #the estimator of sigma sigma <- sqrt(sum(e^2)/(n-p-1)) #standard residual zre <- e/sigma #student residual sre <- e/(sigma*sqrt(1-h)) rstandard(lm5.1) library(MASS) stdres(lm5.1) #学生化残差图 plot(stdres(lm5.1)~yhat) #delete student residual srei <- sre*(1/sqrt((n-p-1)/(n-p-2)-sre^2/(n-p-2))) rstudent(lm5.1) studres(lm5.1) #MASS is.vector(studres(lm5.1)) which(abs(studres(lm5.1))>3) #center leverage vector centx <- scale(x[,-1],center=T,scale=F) ch <- diag(centx %*% solve(t(centx) %*% centx) %*% t(centx)) sum(ch) which(ch>3*mean(ch)) ch1 <- h-1/n #cook distance cookd <- cooks.distance(lm5.1) which(cookd>1) cook <- e^2/((p+1)*sigma^2)*h/(1-h)^2 plot(lm5.1)解释这段程序
时间: 2024-02-26 17:52:31 浏览: 68
这段程序的作用是进行线性回归模型的诊断和优化,并绘制相关的图形。
首先,程序设置数据文件的路径为 "G:/实用回归分析/",并读取名为 "li5.1.csv" 的 CSV 文件,将其转换为数据框 li5.1。然后删除数据框的前两列,将第一列赋值给 y,将其余列作为自变量构成的矩阵并添加截距项,赋值给 x。
接下来,程序计算样本量 n 和自变量个数 p。使用 lm() 函数拟合线性回归模型 lm5.1,计算残差 e 和拟合值 yhat。绘制残差图,使用 plot() 函数,分别绘制 yhat 和 e 的散点图以及 e 关于 yhat 的散点图。
然后,程序计算杠杆值和帽子值,使用 hatvalues() 函数计算帽子值,使用 apply() 函数计算杠杆值,并将其打印出来。计算标准化残差 zre 和学生化残差 sre,使用 rstandard() 函数和 stdres() 函数计算标准化残差。绘制学生化残差图,使用 plot() 函数,绘制 sre 关于 yhat 的散点图。
接下来,程序计算修正后的学生化残差 srei,使用 rstudent() 函数和 studres() 函数计算修正后的学生化残差,并使用 which() 函数找出绝对值大于3的学生化残差的位置。
然后,程序计算中心化杠杆值 ch,使用 diag() 函数计算中心化矩阵的对角元素,使用 apply() 函数计算中心化杠杆值,并打印出其总和以及绝对值大于 3 倍平均值的中心化杠杆值的位置。计算 Cook 距离 cookd 和 cook 值,使用 cooks.distance() 函数计算 Cook 距离,使用 e^2/((p+1)*sigma^2)*h/(1-h)^2 计算 cook 值。
最后,使用 plot() 函数对 lm5.1 进行绘图,绘制回归诊断图,包括标准化残差图、杠杆值图和 Cook 距离图。
阅读全文