如何按组/子集进行 Leave One Out 交叉验证？

Question

这个问题是上一个问题 (Linear Regression prediction in R using Leave One out Approach) 的第二部分。

我正在尝试为每个国家/地区建立模型并使用留一法生成线性回归预测。换句话说，在下面的代码中，在构建模型 1 和模型 2 时，使用的“数据”不应该是整个数据集。相反，它应该是数据集（国家）的一个子集。每个国家/地区的数据都应使用使用该国家/地区特定数据构建的模型进行评估。

下面的代码返回错误。我该如何修改/修复下面的代码来做到这一点？或者有更好的方法吗？

library(modelr)
install.packages("gapminder")
library(gapminder)                           
data(gapminder) 

#CASE 1
model1 <- lm(lifeExp ~ pop, data = gapminder, subset = country)
model2 <- lm(lifeExp ~ pop + gdpPercap, data = gapminder, subset = country)

models <- list(fit_model1 = model1,fit_model2 = model2)

gapminder %>% nest_by(continent, country) %>%
  bind_cols(
    map(1:nrow(gapminder), function(i) {
      map_dfc(models, function(model) {
        training <- data[-i, ] 
        fit <- lm(model, data = training)
        
        validation <- data[i, ]
        predict(fit, newdata = validation)
        
      })
    }) %>%
      bind_rows()
  )
 

Answer 1

最简洁直接的解决方案是嵌套的for 循环方法，其中外循环是两个模型公式，内循环是我们想要省略的统一体。这也可以使用outer 来完成，我稍后也会展示。

为了清楚起见，我首先展示如何在每次迭代中省略一个观察值（即一行）（第一部分）。稍后我将展示如何省略一个集群（例如国家/地区）（第二部分）。我还使用了内置的iris 数据集，它更小，因此更易于处理。它包含一个"Species" 列，旨在对应于您数据中的"countries"。

第一部分

首先，我们将这两个公式放入一个列表中，并按照我们希望它们稍后出现在结果列中的方式命名它们。

FOAE <- list(fit1=Petal.Length ~ Sepal.Length, 
             fit2=Petal.Length ~ Sepal.Length + Petal.Width)

对于循环，我们要初始化一个矩阵im，其行对应于我们要省略的行数，列对应于模型公式的数量。

im <- matrix(NA, nrow=nrow(iris), ncol=length(FOAE), 
             dimnames=list(NULL, names(FOAE)))

看起来像这样：

head(im, n=3)
#      fit1 fit2
# [1,]   NA   NA
# [2,]   NA   NA
# [3,]   NA   NA

现在我们如上所述循环遍历公式和行。

for (i in seq(FOAE)) {
  for(j in seq(nrow(iris))) {
    train <- iris[-j,]  
    test <- iris[j,]    
    fit <- lm(FOAE[[i]], data=train)
    im[j, i] <- predict(fit, newdata=test)
  }
}

im现在已经被填满了，我们可以将cbind它到原来的iris数据集中得到我们的结果res1。

res1 <- cbind(iris, im)
head(res1)
#   Sepal.Length Sepal.Width Petal.Length Petal.Width Species     fit1     fit2
# 1          5.1         3.5          1.4         0.2  setosa 2.388501 1.611976
# 2          4.9         3.0          1.4         0.2  setosa 2.014324 1.501389
# 3          4.7         3.2          1.3         0.2  setosa 1.639805 1.392955
# 4          4.6         3.1          1.5         0.2  setosa 1.446175 1.333199
# 5          5.0         3.6          1.4         0.2  setosa 2.201646 1.556620
# 6          5.4         3.9          1.7         0.4  setosa 2.944788 2.127184

为了替代地遵循outer方法，我们将for循环内的代码放入formulaVectorize，以便它可以处理矩阵列（即向量） .

FUN1 <- Vectorize(function(x, y) {
  train <- iris[-x,]
  test <- iris[x,]
  fit <- lm(y, data=train)
  predict(fit, newdata=test)
})

现在我们将FOAE 和1:nrow(iris) 的行随后省略，与FUN1 一起放入outer()。这已经给了我们结果，我们可以用与上面相同的方式将cbind 到iris 得到我们的结果res2。

o1 <- outer(FOAE, 1:nrow(iris), FUN1)
res2 <- cbind(iris, o1)

head(res2)
#   Sepal.Length Sepal.Width Petal.Length Petal.Width Species     fit1     fit2
# 1          5.1         3.5          1.4         0.2  setosa 2.388501 1.611976
# 2          4.9         3.0          1.4         0.2  setosa 2.014324 1.501389
# 3          4.7         3.2          1.3         0.2  setosa 1.639805 1.392955
# 4          4.6         3.1          1.5         0.2  setosa 1.446175 1.333199
# 5          5.0         3.6          1.4         0.2  setosa 2.201646 1.556620
# 6          5.4         3.9          1.7         0.4  setosa 2.944788 2.127184

## test if results are different is negative 
stopifnot(all.equal(res1, res2))

第二部分

在省略集群（即物种或国家）时，我们可能会采用类似的方法。我在这里展示outer 方法。我们想要改变的是，我们现在想要忽略属于特定集群的观察，这里是"Species"（在你的情况下是"countries"），我们将unique 的值放入向量Species.u 。由于这些值采用"character" 或"factor" 格式，我们使用data[!data$cluster %in% x, ] 而不是data[-x, ] 对数据进行子集化。因为predict 会在集群中产生多个值，但我们希望在各个集群中具有相同的值，所以我们可能想要使用统计信息，例如每个集群的mean 预测。我们根据集群使用rownames。

FUN2 <- Vectorize(function(x, y) {
  train <- iris[!iris$Species %in% x,]
  test <- iris[iris$Species %in% x,]
  fit <- lm(y, data=train)
  mean(predict(fit, newdata=test))
})
Species.u <- unique(iris$Species)

o2 <- `rownames<-`(outer(Species.u, FOAE, FUN2), Species.u)

这给了我们一个比我们的数据集更小的矩阵。感谢rownames，我们可以match 预测它们所属的集群。

o2
#                fit1     fit2
# setosa     3.609943 2.662609
# versicolor 3.785760 3.909919
# virginica  4.911009 5.976922

res3 <- cbind(iris, o2[match(iris$Species, rownames(o2)), ])

head(res3)
#          Sepal.Length Sepal.Width Petal.Length Petal.Width Species     fit1     fit2
# setosa            5.1         3.5          1.4         0.2  setosa 3.609943 2.662609
# setosa.1          4.9         3.0          1.4         0.2  setosa 3.609943 2.662609
# setosa.2          4.7         3.2          1.3         0.2  setosa 3.609943 2.662609
# setosa.3          4.6         3.1          1.5         0.2  setosa 3.609943 2.662609
# setosa.4          5.0         3.6          1.4         0.2  setosa 3.609943 2.662609
# setosa.5          5.4         3.9          1.7         0.4  setosa 3.609943 2.662609

tail(res3)
#              Sepal.Length Sepal.Width Petal.Length Petal.Width   Species     fit1     fit2
# virginica.44          6.7         3.3          5.7         2.5 virginica 4.911009 5.976922
# virginica.45          6.7         3.0          5.2         2.3 virginica 4.911009 5.976922
# virginica.46          6.3         2.5          5.0         1.9 virginica 4.911009 5.976922
# virginica.47          6.5         3.0          5.2         2.0 virginica 4.911009 5.976922
# virginica.48          6.2         3.4          5.4         2.3 virginica 4.911009 5.976922
# virginica.49          5.9         3.0          5.1         1.8 virginica 4.911009 5.976922

编辑

在这个版本的FUN2、FUN3中，每个集群的模型的输出都是rbinded（当然是两列，因为有两个模型）。

FUN3 <- Vectorize(function(x, y) {
  train <- iris[!iris$Species %in% x,]
  test <- iris[iris$Species %in% x,]
  fit <- lm(y, data=train)
  (predict(fit, newdata=test))
}, SIMPLIFY=F)
Species.u <- unique(iris$Species)

o3 <- `rownames<-`(outer(Species.u, FOAE, FUN3), Species.u)

res32 <- cbind(iris, apply(o3, 2, unlist))
head(res32)
#          Sepal.Length Sepal.Width Petal.Length Petal.Width Species     fit1     fit2
# setosa.1          5.1         3.5          1.4         0.2  setosa 3.706940 2.678255
# setosa.2          4.9         3.0          1.4         0.2  setosa 3.500562 2.547587
# setosa.3          4.7         3.2          1.3         0.2  setosa 3.294183 2.416919
# setosa.4          4.6         3.1          1.5         0.2  setosa 3.190994 2.351586
# setosa.5          5.0         3.6          1.4         0.2  setosa 3.603751 2.612921
# setosa.6          5.4         3.9          1.7         0.4  setosa 4.016508 3.073249

编辑 2

正如我在您的评论中了解到的那样，您想要 1. 集群中的数据子集。这将是下面的FUN4 中的ss。然后ss 也通过在子集ss 的行上省略一行z 来进行子集化。

FUN4 <- Vectorize(function(x, y) {
  ## subsets first by cluster then by row
  ss <- iris[iris$Species %in% x,]  ## cluster subset
  sapply(1:nrow(ss), function(z) {  ## subset rows using `sapply`
    train <- ss[-z,]  ## train data w/o row z
    test <- ss[z,]    ## test data for `predict`, just row z
    fit <- lm(y, data=train)
    predict(fit, newdata=test)
  })
}, SIMPLIFY=F)

## the two models
FOAE <- list(fit1=Petal.Length ~ Sepal.Length, 
             fit2=Petal.Length ~ Sepal.Length + Petal.Width)

## unique cluster names
Species.u <- unique(iris$Species)

## with the `outer` we iterate over all the permutations of clusters and models `FOAE`.
o4 <- `rownames<-`(outer(Species.u, FOAE, FUN4), Species.u)

## `unlist`ed result is directly `cbind`able to original data
res4 <- cbind(iris, apply(o4, 2, unlist))

## result
head(res4)
# Sepal.Length Sepal.Width Petal.Length Petal.Width Species     fit1     fit2
# setosa.1          5.1         3.5          1.4         0.2  setosa 1.476004 1.451029
# setosa.2          4.9         3.0          1.4         0.2  setosa 1.449120 1.431737
# setosa.3          4.7         3.2          1.3         0.2  setosa 1.426185 1.416492
# setosa.4          4.6         3.1          1.5         0.2  setosa 1.404040 1.398103
# setosa.5          5.0         3.6          1.4         0.2  setosa 1.462460 1.441295
# setosa.6          5.4         3.9          1.7         0.4  setosa 1.504990 1.559045