r-如何使用循环更新公式进行正向选择



下午好。

我正在为正向变量选择构建循环。在每一步中,我都会在dynlm回归中依次添加一个项,看看解释变量中的哪个变量在最大化平方方面效果最好。由于我有100个解释变量,我最多有100个模型,我想存储每个模型的残差,这样我就可以比较100个模型中哪一个的AIC、BIC(手动(最低

到目前为止,我所完成的是

### M0 (Null Model)
fm0 <- tr.y.nm[,1]~L(tr.y.nm[,1],1:4) #AR(4) model
r0 <- summary(dynlm(fm))$r.squared
### loop to choose M1 ~ M100 with forward selection
fml <- vector("list", 100)
r <- rep(0, 100) 
ssr <- rep(0,100)
rr <- matrix(c(1:100, 1:100), nrow=100,ncol=100) # 100 by 100
x <- rep(0,100) # index to indicate which regressor should be added when selecting best model
### M1
for (i in 1:ncol(tr.x.nm)){
fml[[1]] <- update(fm0, .~. + L(tr.x.nm[,i], 0:1), evaluate=F)
model <- dynlm(fml[[1]])
rr[i,1] <- summary(model)$r.squared
}
x[1] <- which.max(rr[,1])
ssr[1] <- sum(model$resid^2, na.rm=T)
### M2
for (i in 1:ncol(tr.x.nm)){
fm <- fml[[1]]
fml[[2]] <- update(fm, .~. + L(tr.x.nm[,x[[1]]], 0:1), evaluate=F)
model <- dynlm(fml[[2]])
rr[i,2] <- summary(model)$r.squared
}
x[2] <- which.max(rr[,2])
ssr[2] <- sum(model$resid^2, na.rm=T)

### M3 
for (i in 1:ncol(tr.x.nm)){
fm <- fml[[2]]
fml[[3]] <- update(fm, .~. + L(tr.x.nm[,x[[2]]], 0:1), evaluate=F)
model <- dynlm(fml[[3]])
rr[i,3] <- summary(model)$r.squared
}
x[3] <- which.max(rr[,3])
ssr[3] <- sum(model$resid^2, na.rm=T)
# M4
for (i in 1:ncol(tr.x.nm)){
fm <- fml[[3]]
fml[[4]] <- update(fm, .~. + L(tr.x.nm[,x[[3]]], 0:1), evaluate=F)
model <- dynlm(fml[[4]])
rr[i,4] <- summary(model)$r.squared
}
x[4] <- which.max(rr[,4])
ssr[4] <- sum(model$resid^2, na.rm=T)

fml是一个100的列表,我想存储每个模型中使用的公式。

例如>磁头

[[1]]
tr.y.nm[, 1] ~ L(tr.y.nm[, 1], 1:4) + L(tr.x.nm[, i], 0:1)
[[2]]
tr.y.nm[, 1] ~ L(tr.y.nm[, 1], 1:4) + L(tr.x.nm[, i], 0:1) + 
L(tr.x.nm[, x[[1]]], 0:1)
[[3]]
tr.y.nm[, 1] ~ L(tr.y.nm[, 1], 1:4) + L(tr.x.nm[, i], 0:1) + 
L(tr.x.nm[, x[[1]]], 0:1) + L(tr.x.nm[, x[[2]]], 0:1)
[[4]]
tr.y.nm[, 1] ~ L(tr.y.nm[, 1], 1:4) + L(tr.x.nm[, i], 0:1) + 
L(tr.x.nm[, x[[1]]], 0:1) + L(tr.x.nm[, x[[2]]], 0:1) + L(tr.x.nm[, 
x[[3]]], 0:1)

考虑到这一点,由于M2~M100模型分别从M1~M99中增加了一个项我想使用类似于跟随的循环

for (j in 2:100) {
for (i in 1:ncol(tr.x.nm)){
fm <- fml[[j-1]]
fml[[j]] <- update(fm, .~. + L(tr.x.nm[,x[j-1]], 0:1), evaluate=F)
model <- dynlm(fml[[j]])
rr[i,j] <- summary(model)$r.squared
}
x[j] <- which.max(rr[,j])
}

然而,j在上面的循环中没有被识别,并且由于第4个索引,输出与前一个相同,不起作用。

> fml
[[1]]
tr.y.nm[, 1] ~ L(tr.y.nm[, 1], 1:4) + L(tr.x.nm[, i], 0:1)
[[2]]
tr.y.nm[, 1] ~ L(tr.y.nm[, 1], 1:4) + L(tr.x.nm[, i], 0:1) + 
L(tr.x.nm[, x[j - 1]], 0:1)
[[3]]
tr.y.nm[, 1] ~ L(tr.y.nm[, 1], 1:4) + L(tr.x.nm[, i], 0:1) + 
L(tr.x.nm[, x[j - 1]], 0:1)
[[4]]
tr.y.nm[, 1] ~ L(tr.y.nm[, 1], 1:4) + L(tr.x.nm[, i], 0:1) + 
L(tr.x.nm[, x[j - 1]], 0:1)

我想用j索引是不起作用的,也不知道为什么。如何用前面定义过值的字母j索引列表?

实际上我成功地运行了循环,但它需要很长时间才能完成。我相信sapply或lappy会更有效率,但不知道如何修改我的代码。

for (j in 2:100){
for (i in 1:ncol(tr.x.nm)){
fm <- update(fm0,.~.+L(tr.x.nm[,i], 0:1))
scope <- paste0("L(tr.x.nm[,x[", 1:(j-1), "]],0:1)")
fit<-fml[[j]]<- update.formula(fm, paste(".~. +", paste(scope, collapse="+")))
model <- dynlm(fit)
rr[i,j] <- summary(model)$r.squared
}
x[j] <- which.max(rr[,j])
ssr[j] <-sum(model$resid^2, na.rm=T)
}

最新更新