R语言beta1hat <- as.matrix(solve(t(x1) %*% x1) %*% t(x1) %*% y) beta2hat <- as.matrix(solve(t(x2) %*% x2) %*% t(x2) %*% y) sig<- as.numeric(sqrt(sum((y-x1%*%beta1hat-x2%*%beta2hat)^2)/(n-2-1))) sig <- as.matrix(sig) hlo1 <- c(beta1hat-sigma*solve(t(x1)%*%x1)*qnorm(0.025)) hup1 <- c(beta1hat+sigma*solve(t(x1)%*%x1)*qnorm(0.025))报错
时间: 2023-12-16 08:03:37 浏览: 86
这个错误通常是由于维度不匹配导致的。在计算hlo1和hup1时,使用了sigma作为矩阵或数组,并且与其他矩阵或数组进行乘法运算。但是,sigma应该是一个数值型变量,而不是矩阵或数组。建议将sigma转换为数值型变量,以避免这个错误。另外,由于beta1hat和beta2hat已经是矩阵或数组,不需要再使用as.matrix()函数进行转换。
下面是修改后的代码:
```
beta1hat <- solve(t(x1) %*% x1) %*% t(x1) %*% y
beta2hat <- solve(t(x2) %*% x2) %*% t(x2) %*% y
sig <- as.numeric(sqrt(sum((y - x1 %*% beta1hat - x2 %*% beta2hat)^2) / (n - 2 - 1)))
hlo1 <- c(beta1hat - sig * solve(t(x1) %*% x1) * qnorm(0.025))
hup1 <- c(beta1hat + sig * solve(t(x1) %*% x1) * qnorm(0.025))
```
这样就可以避免这个错误了。
相关问题
hup1 <- vector(beta1hat+sigma1*solve(t(x1)%*%x1)*qnorm(0.025)) Error in vector(beta1hat + sigma1 * solve(t(x1) %*% x1) * qnorm(0.025)) : invalid 'mode' argument In addition: Warning message: In beta1hat + sigma1 * solve(t(x1) %*% x1) * qnorm(0.025) : Recycling array of length 1 in vector-array arithmetic is deprecated. Use c() or as.vector() instead.
在你提供的代码中,`hup1 <- vector(beta1hat+sigma1*solve(t(x1)%*%x1)*qnorm(0.025))` 出现了两个问题。
第一个问题是 `vector()` 函数的错误。这个函数是用来创建由重复元素组成的向量的,而你需要用 `c()` 函数将计算结果连接成向量。因此,应该将 `vector()` 替换为 `c()`。
第二个问题是一个警告信息,它表明在向量间的运算中,向量长度不匹配会导致数组重复。警告信息建议使用 `c()` 或 `as.vector()` 函数来解决这个问题。因此,你需要使用 `c()` 函数来将计算结果连接成一个长度为 k 的向量。
综上所述,修改后的代码如下:
```
hup1 <- c(beta1hat+sigma1*solve(t(x1)%*%x1)*qnorm(0.025))
hup1 <- c(hup1, rep(0, k-length(hup1)))
```
同样的,对于 `hlo2` 和 `hup2` 的计算也需要做相应的修改。
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)解释这段程序
这段程序的作用是进行线性回归模型的诊断和优化,并绘制相关的图形。
首先,程序设置数据文件的路径为 "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 距离图。
阅读全文