我有以下可重现的示例,它在混合回归模型上计算简单的邻域交叉验证。如图所示,如果我将这组命令转换为一个函数,它们的行为将不正确,拟合的 lmer 模型不会在每个 for 循环实例上更新,而是始终采用最终 for 循环实例的值。
这似乎是 R 编程的一个怪癖,我很想知道发生了什么,因为它将来很容易被忽视。
可重现的示例:
library(lme4)
set.seed(2002)
dataf <- data.frame( x = rnorm(15), p = rep(1:5, each = 3), rp = rep(rnorm(5),each = 3))
dataf$y <- rnorm(15, dataf$x + dataf$rp)
mod <- lmer( y ~ x + (1|p), data = dataf)
# Inner commands for function:
dataf <- getData(mod)
pvec <- unique(dataf$p)
r2vec <- numeric(0)
for ( pnow in pvec ) {
selp <- dataf$p == pnow
newmod <- update(mod, data = dataf[!selp,])
print(paste(pnow, fixef(newmod),collapse = ", "))
r2vec <- c( r2vec,
(dataf[selp,]$y - predict(newmod,newdata = dataf[selp,], allow.new.levels = T))^2
)
}
# [1] "1 0.285049274214266, 1 0.95048020397442"
# [1] "2 0.185721150305479, 2 0.901716963753631"
# [1] "3 0.810224187608013, 3 0.615735520886265"
# [1] "4 0.869392151461034, 4 0.968302207873863"
# [1] "5 0.454297693044713, 5 0.914047431881036"
mean(r2vec) # 2.241611
# Function with inner commands as above:
lopofun <- function(mod) {
dataf <- getData(mod)
pvec <- unique(dataf$p)
r2vec <- numeric(0)
for ( pnow in pvec ) {
selp <- dataf$p == pnow
newmod <- update(mod, data = dataf[!selp,])
print(paste(pnow, fixef(newmod),collapse = ", "))
r2vec <- c( r2vec,
(dataf[selp,]$y - predict(newmod,newdata = dataf[selp,], allow.new.levels = T))^2
)
}
mean(r2vec)
}
lopofun(mod) # 0.3388175
# [1] "1 0.454297693044713, 1 0.914047431881036"
# [1] "2 0.454297693044713, 2 0.914047431881036"
# [1] "3 0.454297693044713, 3 0.914047431881036"
# [1] "4 0.454297693044713, 4 0.914047431881036"
# [1] "5 0.454297693044713, 5 0.914047431881036"
正如插入的 print 函数所示,newmod 并未更新,而是采用 for 循环的每个实例中应为最终值的值。
我想继续使用 update(),因为尝试手动提取模型特征并重新拟合会很笨拙。毫无疑问,有一种更好的方法可以在不使用 for 循环的情况下进行编码,尽管如此,for 循环是一个简单的结构,我不确定为什么它会在这里失败。
编辑:我很欣赏审稿人的建议,认为这是一个惰性评估问题,使用 sapply() 而不是 for 循环可能会解决问题。这与我对惰性求值的(基本)理解不符。 print() 需要评估循环变量 pnow 并使用 newmod:它返回正确的 pnow 值,但不返回 newmod 值。
下面是使用 sapply() 而不是 for 循环来实现同样事情的丑陋尝试。但是,这里原始命令和函数都失败了。
# Inner commands using sapply
dataf <- getData(mod)
pvec <- unique(dataf$p)
onepfun <- function(pnow) {
selp <- dataf$p == pnow
newmod <- update(mod, data = dataf[!selp,])
print(paste(pnow, paste(fixef(newmod),collapse = ", ")))
(dataf[selp,]$y - predict(newmod,newdata = dataf[selp,],
allow.new.levels = T))^2
}
bob <- sapply(pvec,onepfun)
dim(bob) <- NULL
mean(bob) # 0.3388175
# [1] "1 0.454297693044713, 0.914047431881036"
# [1] "2 0.454297693044713, 0.914047431881036"
# [1] "3 0.454297693044713, 0.914047431881036"
# [1] "4 0.454297693044713, 0.914047431881036"
# [1] "5 0.454297693044713, 0.914047431881036"
# Function attempt with sapply
lopofun2 <- function(mod) {
dataf <- getData(mod)
pvec <- unique(dataf$p)
onepfun <- function(pnow) {
selp <- dataf$p == pnow
newmod <- update(mod, data = dataf[!selp,])
print(paste(pnow, paste(fixef(newmod),collapse = ", ")))
(dataf[selp,]$y - predict(newmod,newdata = dataf[selp,],
allow.new.levels = T))^2
}
bob <- sapply(pvec,onepfun)
dim(bob) <- NULL
mean(bob)
}
lopofun2(mod) # 0.3388175
# [1] "1 0.454297693044713, 0.914047431881036"
# [1] "2 0.454297693044713, 0.914047431881036"
# [1] "3 0.454297693044713, 0.914047431881036"
# [1] "4 0.454297693044713, 0.914047431881036"
# [1] "5 0.454297693044713, 0.914047431881036"
在我看来,在调用函数之前添加
selp<- NULL
可以解决问题
selp<- NULL
lopofun(mod)
[1] "1 0.285049274214266, 1 0.95048020397442"
[1] "2 0.185721150305479, 2 0.901716963753631"
[1] "3 0.810224187608013, 3 0.615735520886265"
[1] "4 0.869392151461034, 4 0.968302207873863"
[1] "5 0.454297693044713, 5 0.914047431881036"
[1] 2.241611
我可能是错的,但我认为这可能与函数的本地环境如何与函数外部同名的变量交互有关。